From: Cezary Czaplewski Date: Fri, 17 Jul 2015 10:48:27 +0000 (+0200) Subject: the first version of READ2SIGMA - the same domains in all *sco files, X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=commitdiff_plain;ds=sidebyside;h=4764dff6d4a43d5ae03320fbd0d2dbef3bb37e3f;p=unres.git the first version of READ2SIGMA - the same domains in all *sco files, readpdb_template cannot read pdb files with gaps --- diff --git a/source/unres/src_MD/COMMON.CONTROL b/source/unres/src_MD/COMMON.CONTROL index e9b05fe..e512e53 100644 --- a/source/unres/src_MD/COMMON.CONTROL +++ b/source/unres/src_MD/COMMON.CONTROL @@ -6,14 +6,14 @@ & energy_dec, & sideadd,lsecondary,read_cart,unres_pdb, & vdisulf,searchsc,lmuca,dccart,extconf,out1file, - & gnorm_check,gradout,split_ene + & gnorm_check,gradout,split_ene,read2sigma common /cntrl/ modecalc,iscode,indpdb,indback,indphi,iranconf, & icheckgrad,minim,i2ndstr,refstr,pdbref,outpdb,outmol2,outx, & iprint, & overlapsc,energy_dec,sideadd,lsecondary,read_cart,unres_pdb & ,vdisulf,searchsc,lmuca,dccart,mucadyn,extconf,out1file, & constr_dist,gnorm_check,gradout,split_ene,constr_homology, - & homol_nset + & homol_nset,read2sigma common /homol/ waga_homology(maxprocs/20), & waga_dist, waga_angle, waga_theta, waga_d, dist_cut C... minim = .true. means DO minimization. diff --git a/source/unres/src_MD/COMMON.MD b/source/unres/src_MD/COMMON.MD index bc0014b..7f8a20f 100644 --- a/source/unres/src_MD/COMMON.MD +++ b/source/unres/src_MD/COMMON.MD @@ -26,7 +26,8 @@ c modified by FP (Nov.,2014) & sigma_d(max_template,maxres) c - integer ires_homo(maxdim),jres_homo(maxdim) + integer ires_homo(maxdim), + & jres_homo(maxdim),idomain(max_template,maxres) double precision v_ini,d_time,d_time0,t_bath,tau_bath, & EK,potE,potEcomp(0:n_ene+4),totE,totT,amax,kinetic_T,dvmax,damax, @@ -62,9 +63,11 @@ c common /homrestr/ odl,dih,sigma_dih,sigma_odl, & lim_odl,lim_dih,ires_homo,jres_homo,link_start_homo, & link_end_homo,idihconstr_start_homo,idihconstr_end_homo, + & idomain c c FP (30/10/2014,04/03/2015) -c +c + common /homrestr_double/ & xxtpl,yytpl,zztpl,thetatpl,sigma_theta,sigma_d,sigma_odlir c common /qmeas/ qfrag,qpair,qinfrag,qinpair,wfrag,wpair,eq_time, diff --git a/source/unres/src_MD/readpdb.F b/source/unres/src_MD/readpdb.F index eaf41c8..745b256 100644 --- a/source/unres/src_MD/readpdb.F +++ b/source/unres/src_MD/readpdb.F @@ -496,3 +496,243 @@ c print *,vbld(i+1),vbld(i+1+nres) enddo return end + subroutine readpdb_template(k) +C Read the PDB file with gaps for read_constr_homology with read2sigma +C and convert the peptide geometry into virtual-chain geometry. + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.LOCAL' + include 'COMMON.VAR' + include 'COMMON.CHAIN' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.GEO' + include 'COMMON.NAMES' + include 'COMMON.CONTROL' + include 'COMMON.DISTFIT' + include 'COMMON.SETUP' + include 'COMMON.MD' + integer i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity, + & ishift_pdb + logical lprn /.false./,fail + double precision e1(3),e2(3),e3(3) + double precision dcj,efree_temp + character*3 seq,res + character*5 atom + character*80 card + double precision sccor(3,20) + integer rescode + efree_temp=0.0d0 + ibeg=1 + ishift1=0 + ishift=0 +c write (2,*) "UNRES_PDB",unres_pdb + ires=0 + ires_old=0 + iii=0 + lsecondary=.false. + nhfrag=0 + nbfrag=0 + do i=1,10000 + read (ipdbin,'(a80)',end=10) card +c write (iout,'(a)') card + if (card(:3).eq.'END' .or. card(:3).eq.'TER') goto 10 +C Fish out the ATOM cards. + if (index(card(1:4),'ATOM').gt.0) then + read (card(12:16),*) atom +c write (iout,*) "! ",atom," !",ires +c if (atom.eq.'CA' .or. atom.eq.'CH3') then + read (card(23:26),*) ires + read (card(18:20),'(a3)') res +c write (iout,*) "ires",ires,ires-ishift+ishift1, +c & " ires_old",ires_old +c write (iout,*) "ishift",ishift," ishift1",ishift1 +c write (iout,*) "IRES",ires-ishift+ishift1,ires_old + if (ires-ishift+ishift1.ne.ires_old) then +C Calculate the CM of the preceding residue. +c if (ibeg.eq.0) call sccenter(ires,iii,sccor) + if (ibeg.eq.0) then +c write (iout,*) "Calculating sidechain center iii",iii + if (unres_pdb) then + do j=1,3 + dc(j,ires)=sccor(j,iii) + enddo + else + call sccenter(ires_old,iii,sccor) + endif + iii=0 + endif +C Start new residue. + if (res.eq.'Cl-' .or. res.eq.'Na+') then + ires=ires_old + cycle + else if (ibeg.eq.1) then +c write (iout,*) "BEG ires",ires + ishift=ires-1 + if (res.ne.'GLY' .and. res.ne. 'ACE') then + ishift=ishift-1 + itype(1)=21 + endif + ires=ires-ishift+ishift1 + ires_old=ires +c write (iout,*) "ishift",ishift," ires",ires, +c & " ires_old",ires_old + ibeg=0 + else + ishift=ishift-(ires-ishift+ishift1-ires_old-1) + ires=ires-ishift+ishift1 + ires_old=ires + endif + if (res.eq.'ACE' .or. res.eq.'NHE') then + itype(ires)=10 + else + itype(ires)=rescode(ires,res,0) + endif + else + ires=ires-ishift+ishift1 + endif +c write (iout,*) "ires_old",ires_old," ires",ires + if (card(27:27).eq."A" .or. card(27:27).eq."B") then +c ishift1=ishift1+1 + endif +c write (2,*) "ires",ires," res ",res," ity",ity + if (atom.eq.'CA' .or. atom.eq.'CH3' .or. + & res.eq.'NHE'.and.atom(:2).eq.'HN') then + read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3) +c write (iout,*) "backbone ",atom +#ifdef DEBUG + write (iout,'(2i3,2x,a,3f8.3)') + & ires,itype(ires),res,(c(j,ires),j=1,3) +#endif + iii=iii+1 + do j=1,3 + sccor(j,iii)=c(j,ires) + enddo + if (ishift.ne.0) then + ires_ca=ires+ishift-ishift1 + else + ires_ca=ires + endif +c write (*,*) card(23:27),ires,itype(ires) + else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. + & atom.ne.'N' .and. atom.ne.'C' .and. + & atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. + & atom.ne.'OXT' .and. atom(:2).ne.'3H') then +c write (iout,*) "sidechain ",atom + iii=iii+1 + read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3) + endif + endif + enddo + 10 continue +#ifdef DEBUG + write (iout,'(a,i5)') ' Number of residues found: ',ires +#endif + if (ires.eq.0) return +C Calculate the CM of the last side chain. + if (iii.gt.0) then + if (unres_pdb) then + do j=1,3 + dc(j,ires)=sccor(j,iii) + enddo + else + call sccenter(ires,iii,sccor) + endif + endif + nres=ires + nsup=nres + nstart_sup=1 + if (itype(nres).ne.10) then + nres=nres+1 + itype(nres)=21 + do j=1,3 + dcj=c(j,nres-2)-c(j,nres-3) + c(j,nres)=c(j,nres-1)+dcj + c(j,2*nres)=c(j,nres) + enddo + endif + do i=2,nres-1 + do j=1,3 + c(j,i+nres)=dc(j,i) + enddo + enddo + do j=1,3 + c(j,nres+1)=c(j,1) + c(j,2*nres)=c(j,nres) + enddo + if (itype(1).eq.21) then + nsup=nsup-1 + nstart_sup=2 + if (unres_pdb) then +C 2/15/2013 by Adam: corrected insertion of the first dummy residue + call refsys(2,3,4,e1,e2,e3,fail) + if (fail) then + e2(1)=0.0d0 + e2(2)=1.0d0 + e2(3)=0.0d0 + endif + do j=1,3 + c(j,1)=c(j,2)-3.8d0*e2(j) + enddo + else + do j=1,3 + dcj=c(j,4)-c(j,3) + c(j,1)=c(j,2)-dcj + c(j,nres+1)=c(j,1) + enddo + endif + endif +C Calculate internal coordinates. + if (lprn) then + write (iout,'(/a)') + & "Cartesian coordinates of the reference structure" + write (iout,'(a,3(3x,a5),5x,3(3x,a5))') + & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)" + do ires=1,nres + write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') + & restyp(itype(ires)),ires,(c(j,ires),j=1,3), + & (c(j,ires+nres),j=1,3) + enddo + endif +C Calculate internal coordinates. + if(me.eq.king.or..not.out1file)then + write (iout,'(a)') + & "Backbone and SC coordinates as read from the PDB" + do ires=1,nres + write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') + & ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3), + & (c(j,nres+ires),j=1,3) + enddo + endif + call int_from_cart(.true.,.false.) + call sc_loc_geom(.false.) + do i=1,nres + thetaref(i)=theta(i) + phiref(i)=phi(i) + enddo + do i=1,nres-1 + do j=1,3 + dc(j,i)=c(j,i+1)-c(j,i) + dc_norm(j,i)=dc(j,i)*vbld_inv(i+1) + enddo + enddo + do i=2,nres-1 + do j=1,3 + dc(j,i+nres)=c(j,i+nres)-c(j,i) + dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres) + enddo +c write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3), +c & vbld_inv(i+nres) + enddo +c call chainbuild +C Copy the coordinates to reference coordinates + do i=1,2*nres + do j=1,3 + cref(j,i)=c(j,i) + enddo + enddo + + + ishift_pdb=ishift + return + end diff --git a/source/unres/src_MD/readrtns.F b/source/unres/src_MD/readrtns.F index 918c761..3aa5fd6 100644 --- a/source/unres/src_MD/readrtns.F +++ b/source/unres/src_MD/readrtns.F @@ -2724,8 +2724,9 @@ c & sigma_odl_temp(maxres,maxres,max_template) c c FP - Nov. 2014 Temporary specifications for new vars c - double precision rescore_tmp,x12,y12,z12 + double precision rescore_tmp,x12,y12,z12,rescore2_tmp double precision, dimension (max_template,maxres) :: rescore + double precision, dimension (max_template,maxres) :: rescore2 character*24 pdbfile,tpl_k_rescore c ----------------------------------------------------------------- c Reading multiple PDB ref structures and calculation of retraints @@ -2743,6 +2744,7 @@ c Alternative: reading from input call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma call readi(controlcard,"HOMOL_NSET",homol_nset,1) + read2sigma=(index(controlcard,'READ2SIGMA').gt.0) if (homol_nset.gt.1)then call card_concat(controlcard) read(controlcard,*) (waga_homology(i),i=1,homol_nset) @@ -2770,30 +2772,8 @@ c lim_theta=0 lim_xx=0 c -c Reading HM global scores (prob not required) -c -c open (4,file="HMscore") -c do k=1,constr_homology -c read (4,*,end=521) hmscore_tmp -c hmscore(k)=hmscore_tmp ! Another transformation can be used -c write(*,*) "Model", k, ":", hmscore(k) -c enddo -c521 continue - -c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d do k=1,constr_homology - read(inp,'(a)') pdbfile -c Next stament causes error upon compilation (?) -c if(me.eq.king.or. .not. out1file) -c write (iout,'(2a)') 'PDB data will be read from file ', -c & pdbfile(:ilen(pdbfile)) - open(ipdbin,file=pdbfile,status='old',err=33) - goto 34 - 33 write (iout,'(a)') 'Error opening PDB file.' - stop - 34 continue -c print *,'Begin reading pdb data' c c Files containing res sim or local scores (former containing sigmas) c @@ -2801,24 +2781,7 @@ c write(kic2,'(bz,i2.2)') k tpl_k_rescore="template"//kic2//".sco" -c tpl_k_sigma_odl="template"//kic2//".sigma_odl" -c tpl_k_sigma_dih="template"//kic2//".sigma_dih" -c tpl_k_sigma_theta="template"//kic2//".sigma_theta" -c tpl_k_sigma_d="template"//kic2//".sigma_d" - unres_pdb=.false. - call readpdb -c -c Distance restraints -c -c ... --> odl(k,ii) -C Copy the coordinates from reference coordinates (?) - do i=1,2*nres - do j=1,3 - c(j,i)=cref(j,i) -c write (iout,*) "c(",j,i,") =",c(j,i) - enddo - enddo c c From read_dist_constr (commented out 25/11/2014 <-> res sim) c @@ -2829,23 +2792,53 @@ c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore rescore(k,irec)=0.0d0 goto 1301 endif - read (ientin,*,end=1401) rescore_tmp + if (read2sigma) then + read (ientin,*,end=1401) rescore2_tmp,rescore_tmp, + & idomain(k,irec) + rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) + rescore2(k,irec)=0.5d0*(rescore2_tmp+0.5d0) + else + idomain(k,irec)=1 + read (ientin,*,end=1401) rescore_tmp + c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values - rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores + rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec) + endif 1301 continue enddo 1401 continue close (ientin) -c open (ientin,file=tpl_k_sigma_odl,status='old') -c do irec=1,maxdim ! loop for reading sigma_odl -c read (ientin,*,end=1401) i, j, -c & sigma_odl_temp(i+nnt-1,j+nnt-1,k) ! new variable (?) -c sigma_odl_temp(j+nnt-1,i+nnt-1,k)= ! which purpose? -c & sigma_odl_temp(i+nnt-1,j+nnt-1,k) -c enddo -c 1401 continue -c close (ientin) + + read(inp,'(a)') pdbfile +c Next stament causes error upon compilation (?) +c if(me.eq.king.or. .not. out1file) +c write (iout,'(2a)') 'PDB data will be read from file ', +c & pdbfile(:ilen(pdbfile)) + open(ipdbin,file=pdbfile,status='old',err=33) + goto 34 + 33 write (iout,'(a)') 'Error opening PDB file.' + stop + 34 continue +c print *,'Begin reading pdb data' + + unres_pdb=.false. + if (read2sigma) then + call readpdb_template(k) + else + call readpdb + endif +c +C Copy the coordinates from reference coordinates (?) + do i=1,2*nres + do j=1,3 + c(j,i)=cref(j,i) +c write (iout,*) "c(",j,i,") =",c(j,i) + enddo + enddo + + + if (waga_dist.ne.0.0d0) then ii=0 do i = nnt,nct-2 ! right? without parallel. @@ -2854,6 +2847,9 @@ c do i = 1,nres ! alternative for bounds as used to set initial values i c do j=i+2,nres ! ibid c do i = nnt,nct-2 ! alternative for bounds as used to assign dist restraints in orig. read_constr_homology (s. above) c do j=i+2,nct ! ibid + + if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0) then + ii=ii+1 c write (iout,*) "k",k c write (iout,*) "i",i," j",j," constr_homology", @@ -2882,7 +2878,16 @@ c if (odl(k,ii).le.6.0d0) then c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j) c Other functional forms possible depending on odl(k,ii), eg. c - if (odl(k,ii).le.dist_cut) then + if (read2sigma) then + sigma_odl(k,ii)=0 + do ik=i,j + sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik) + enddo + sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1) + if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = + & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0) + else + if (odl(k,ii).le.dist_cut) then sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) ! other exprs possible c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j) else @@ -2893,25 +2898,24 @@ c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j) sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0) #endif - + endif + endif c Following expr replaced by a positive exp argument c sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* c & dexp(-0.5d0*(odl(k,ii)/dist_cut)**2) c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)* c & dexp(-0.5d0*(odl(k,ii)/dist_cut)**2) - endif + c sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) ! rescore ~ error c sigma_odl(k,ii)=sigma_odl(k,ii)*sigma_odl(k,ii) c c sigma_odl(k,ii)=sigma_odl_temp(i,j,k)* ! new var read from file (?) c & sigma_odl_temp(i,j,k) ! not inverse because of use of res. similarity + endif enddo -c read (ientin,*) sigma_odl(k,ii) ! 1st variant enddo -c lim_odl=ii -c if (constr_homology.gt.0) call homology_partition endif c c Theta, dihedral and SC retraints @@ -2930,6 +2934,10 @@ c close (ientin) c do i=1,nres ! alternative for bounds acc to readpdb? c do i=1,nres-3 ! alternative for bounds as used to set initial values in orig. read_constr_homology c do i=idihconstr_start_homo,idihconstr_end_homo ! with FG parallel. + if (idomain(k,i).eq.0) then + sigma_dih(k,i)=0.0 + cycle + endif dih(k,i)=phiref(i) ! right? c read (ientin,*) sigma_dih(k,i) ! original variant c write (iout,*) "dih(",k,i,") =",dih(k,i) @@ -2966,6 +2974,10 @@ c close (ientin) do i = nnt+2,nct ! right? without parallel. c do i = i=1,nres ! alternative for bounds acc to readpdb? c do i=ithet_start,ithet_end ! with FG parallel. + if (idomain(k,i).eq.0) then + sigma_theta(k,i)=0.0 + cycle + endif thetatpl(k,i)=thetaref(i) c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i) c write(iout,*) "rescore(",k,i,") =",rescore(k,i), @@ -2991,12 +3003,16 @@ c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use c & sigma_d(k,i+nnt-1) c enddo c1404 continue - close (ientin) + do i = nnt,nct ! right? without parallel. c do i=2,nres-1 ! alternative for bounds acc to readpdb? c do i=loc_start,loc_end ! with FG parallel. - if (itype(i).eq.10) goto 1 ! right? + if (itype(i).eq.10) cycle + if (idomain(k,i).eq.0 ) then + sigma_d(k,i)=0.0 + cycle + endif xxtpl(k,i)=xxref(i) yytpl(k,i)=yyref(i) zztpl(k,i)=zzref(i) @@ -3011,10 +3027,8 @@ c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ? c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i) c read (ientin,*) sigma_d(k,i) ! 1st variant if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right? - 1 continue enddo endif - close(ientin) enddo if (waga_dist.ne.0.0d0) lim_odl=ii if (constr_homology.gt.0) call homology_partition