c Main loop to go over possible covalent bond and ligand c orientations, should be called once per each ligand in c the db. c NL 01/2012 subroutine cov_sampler() implicit none include 'max.h' !always necessary include 'sphere.h' !read in receptor 'spheres' include 'molecule.h' !for ligand coordinates and colors include 'match.h' !for rot matrix include 'status.h' !for docking status include 'mscore.h' !for mol_write include 'molsav.h' !for t_ variables include 'dock.h' c local variables integer ligi1,ligi2,ligi3,ligi4 !colored ligand atoms indices common ligi1 real cov_loc(3) !tmp vector to hold the covalent location real a,fi,d !cov bond parameters see cov_bond.f real la,lfi,ld !loop variables to iterate these real lb,lfi1,lfi2 !loop variables to iterate the lig params real lfi1mag !magnitude of bins for fi1 integer i,j,k,m,n,s,o !reg. loop variables integer cli,clj integer sp ! geometry of covalent atom in ligand 1/2/3 real p1(3),p2(3),p3(3) !receptor spheres - see cov_bond.f real best_energy !best energy for any of the flexiable !conformations sampled when score_flexible is called integer dock_status !see status.h for enumerated types real telaps, time2 !time keeping real rigid_part(3,matm(1)) !array to test the rigid part position !before applying to the rest of the !cl array (which could be huge if there !are a lot of conformations) real tmp_cl(3,MAXEATM) real top100scores(100) real tmp_worse integer indexofworse real rvdwa, rvdwb, relec, rsp, rsa integer gbump integer lensteps, ang1steps, ang2steps real curr_energy,curr_a,curr_d,curr_fi real curr_cl(3,MAXEATM) c Declare local constant Pi REAL, PARAMETER :: Pi = 3.1415927 REAL, PARAMETER :: deg_to_rad = 0.0174532925 c Main code c ========= telaps = 0.0 ! time elapsed call get_time(time1) !get receptor atoms from spheres file p1(1)=spcorr(1,1); p1(2)=spcorr(2,1); p1(3)=spcorr(3,1) p2(1)=spcorr(1,2); p2(2)=spcorr(2,2); p2(3)=spcorr(3,2) p3(1)=spcorr(1,3); p3(2)=spcorr(2,3); p3(3)=spcorr(3,3) !iterate over ligand atoms to detect cov1 and cov2 !TODO to allow several covalent mod site on a ligand ! change this loop to run outside cov_sampler and ! call it over again with different atoms !detecting covalent attachment point and its neighbours !assumptions: !- cov attach point color is 8 (see inhier) !- cov n1 to n3 are 9,10,11 !- if 9,10,11 present - assume sp3 geometry !- if only 9,10 present - assume sp2 geometry ligi1=0; ligi2=0; ligi3=0; ligi4=0; do i=1,matm(1) !number of atoms in rigid component if (lcolor(i) .eq. 8) then !assumes covalent atom color is 8 ligi1=i endif if (lcolor(i) .eq. 9) then !assumes its neighbour is 9/10/11 ligi2=i endif if (lcolor(i) .eq. 10) then !assumes its neighbour is 9/10/11 ligi3=i endif if (lcolor(i) .eq. 11) then !assumes its neighbour is 9/10/11 ligi4=i endif enddo !determine geometry of covalent attachmed point by number of colored !neighbours if ((ligi1 .eq. 0) .OR. (ligi2 .eq. 0)) then nsaved = nsaved-1 write (6,*) 'ERROR: could not find two colored ligand & atoms in the rigid part, check db preparation' return else if (ligi3 .eq. 0) then sp=1 else if (ligi4 .eq. 0) then sp=2 else sp=3 endif !setup for scoring: t_variables are used as tmps to store a !current ligs scores and coordinates requiered by mol_sav t_escore = 1.0e9 t_molname = molname t_MFCD = refcod best_energy = 1e15 curr_energy = 1e15 dock_status = NOMATCH lfi1mag=20.0 !store original ligand coordinates read !restore ligand orig coordinates do cli=1,3 do clj=1,eatm tmp_cl(cli,clj)=cl(cli,clj) enddo enddo do cli=1,3 do clj=1,matm(1) rigid_part(cli,clj)=cl(cli,clj) enddo enddo j=1 !counter for created conformations in this call used !only for debug writing of the poses pdbs. d=bondlen-bondrange !set initial cov bond length and ang steps lensteps=int((bondrange/bondstep)*2+1) ang1steps=int((bondang1range/bondang1step)*2+1) ang2steps=int((bondang2range/bondang2step)*2+1) do s=1,lensteps !iterate over bond length d la=bondang1-bondang1range !initial cov bond angle do k=1,ang1steps !iterate over la do lfi=0,340,20 !105,105 !iterate over lfi !convert to radians TODO this can be iterated in RADs a =Pi-(la*deg_to_rad) fi=lfi*deg_to_rad call cov_bond(p1,p2,p3,a,fi,d,cov_loc) !!!!!!Inner loop - for each bond position !!!!!! sample ligand's orientations do lb=bondang2-bondang2range, & bondang2+bondang2range, & bondang2step !restore ligand orig coordinates do cli=1,3 do clj=1,matm(1) rigid_part(cli,clj)=tmp_cl(cli,clj) enddo enddo n=1; m=1 do n=1,18 !iterate over lfi1 lfi1=lfi1mag; lfi2=0 !this call just changes lfi1 call cov_lig(p3,cov_loc,ligi1,ligi2, & ligi3,ligi4,rigid_part,matm(1), & lb,lfi1,lfi2,sp, & ((n.eq.1) .and. (m.eq.1))) ! adjust lb only the first time ! 'time saver' !update xatm do clj=1,matm(1) xatm(1,clj)=rigid_part(1,clj) xatm(2,clj)=rigid_part(2,clj) xatm(3,clj)=rigid_part(3,clj) enddo !check if rigid part is bumping bumping gbump=0; call bump_check(gbump, 0, 1, nhvy(1),bumped) ! write (6,*) "this many bumping: ",gbump if (gbump .eq. 0) then !.eq. 0 !call score_group on rigid part call score_group(1, rvdwa, rvdwb, relec, rsp, rsa) if (rvdwa-rvdwb+relec-rsp-rsa < 200) then !200) then !update cl to current rigid position by applying !same cov_lig transformation, and score hierarchy do cli=1,3 do clj=1,eatm cl(cli,clj)=tmp_cl(cli,clj) enddo enddo call cov_lig(p3,cov_loc,ligi1,ligi2,ligi3,ligi4, & cl,eatm,lb,n*lfi1mag,m*0.0,sp,.true.) !scores the entire hierarchy and returns the !energy of the best conformation call score_flexible(best_energy,dock_status) !save this pose in the correct index call mol_sav(nsaved) endif !rigid score < 0 endif !bumpcheck enddo !n enddo !lb enddo !lfi la=la+bondang1step !increment la enddo !k d = d+bondstep !length update enddo !s end 'length' while !write out results call get_time(time2) telaps = time2 - time1 !if no pose matched - reduce nsaved. if (t_escore .eq. 1.0e9) then nsaved = nsaved-1 write (6,*) 'no poses found for ligand ',refcod else !endif !if (nsaved.eq.0) then !if no ligands were saved (all bump) write(6,'(i6,x,a9,i8,x,i10,i6,i10,1x,f8.1,1x,f8.2)') & nsaved,refcod, nsaved, nsaved, nhvy(0), & tconf, 0.9, telaps write(6,1515) refcod, nbr, t_eessum , t_vdw, & -1*t_pds, -1*t_ads, t_escore 1515 format (" E ",a9,2x,i6,2x,f9.2,f13.2,2f9.2,f13.2) endif return end