From c9dddd7cdc2cbf55f6fa8df9404eea7ecfe59605 Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Tue, 28 Jul 2015 12:21:39 +0200 Subject: [PATCH] Intoduction Lorentizan like restrains to cluster --- source/cluster/wham/src-M/COMMON.CONTROL | 9 +- source/cluster/wham/src-M/COMMON.SBRIDGE | 3 +- source/cluster/wham/src-M/energy_p_new.F | 87 ++++++++++++- source/cluster/wham/src-M/readrtns.F | 178 ++++++++++++++++++++++++++ source/unres/src_MD-M/energy_p_new_barrier.F | 9 +- source/wham/src-M/energy_p_new.F | 8 +- 6 files changed, 281 insertions(+), 13 deletions(-) diff --git a/source/cluster/wham/src-M/COMMON.CONTROL b/source/cluster/wham/src-M/COMMON.CONTROL index 9a2bd18..992affa 100644 --- a/source/cluster/wham/src-M/COMMON.CONTROL +++ b/source/cluster/wham/src-M/COMMON.CONTROL @@ -1,7 +1,10 @@ double precision betaT - integer iscode,indpdb,outpdb,outmol2,iopt,nstart,nend,symetr + integer iscode,indpdb,outpdb,outmol2,iopt,nstart,nend,symetr, + & constr_dist logical refstr,pdbref,punch_dist,print_dist,caonly,lside, - & lprint_cart,lprint_int,from_cart,efree,from_bx,from_cx + & lprint_cart,lprint_int,from_cart,efree,from_bx,from_cx, + & with_dihed_constr common /cntrl/ betaT,iscode,indpdb,refstr,pdbref,outpdb,outmol2, & punch_dist,print_dist,caonly,lside,lprint_cart,lprint_int, - & from_cart,from_bx,from_cx,efree,iopt,nstart,nend,symetr + & from_cart,from_bx,from_cx,efree,iopt,nstart,nend,symetr, + & constr_dist diff --git a/source/cluster/wham/src-M/COMMON.SBRIDGE b/source/cluster/wham/src-M/COMMON.SBRIDGE index bea57e4..028f9ae 100644 --- a/source/cluster/wham/src-M/COMMON.SBRIDGE +++ b/source/cluster/wham/src-M/COMMON.SBRIDGE @@ -2,9 +2,10 @@ integer ns,nss,nfree,iss common /sbridge/ ss_depth,ebr,d0cm,akcm,akth,akct,v1ss,v2ss,v3ss, & ns,nss,nfree,iss(maxss) - double precision dhpb,dhpb1,forcon + double precision dhpb,dhpb1,forcon,fordepth integer ihpb,jhpb,nhpb,idssb,jdssb,ibecarb common /links/ dhpb(maxdim),dhpb1(maxdim),forcon(maxdim), + & fordepth(maxdim), & ihpb(maxdim),jhpb(maxdim),nhpb double precision weidis common /restraints/ weidis diff --git a/source/cluster/wham/src-M/energy_p_new.F b/source/cluster/wham/src-M/energy_p_new.F index 642cd29..a65c91c 100644 --- a/source/cluster/wham/src-M/energy_p_new.F +++ b/source/cluster/wham/src-M/energy_p_new.F @@ -2929,11 +2929,41 @@ C iii and jjj point to the residues for which the distance is assigned. endif C 24/11/03 AL: SS bridges handled separately because of introducing a specific C distance and angle dependent SS bond potential. - if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. - & iabs(itype(jjj)).eq.1) then +C if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. +C & iabs(itype(jjj)).eq.1) then +C call ssbond_ene(iii,jjj,eij) +C ehpb=ehpb+2*eij +C else + if (.not.dyn_ss .and. i.le.nss) then + if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. + & iabs(itype(jjj)).eq.1) then call ssbond_ene(iii,jjj,eij) ehpb=ehpb+2*eij - else + endif !ii.gt.neres + else if (ii.gt.nres .and. jj.gt.nres) then +c Restraints from contact prediction + dd=dist(ii,jj) + if (constr_dist.eq.11) then +C ehpb=ehpb+fordepth(i)**4.0d0 +C & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) + ehpb=ehpb+fordepth(i)**4.0d0 + & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) + fac=fordepth(i)**4.0d0 + & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd +C write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj, +C & ehpb,fordepth(i),dd +C write(iout,*) ehpb,"atu?" +C ehpb,"tu?" +C fac=fordepth(i)**4.0d0 +C & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd + else !constr_dist.eq.11 + if (dhpb1(i).gt.0.0d0) then + ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd +c write (iout,*) "beta nmr", +c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + else !dhpb(i).gt.0.00 + C Calculate the distance between the two points and its difference from the C target distance. dd=dist(ii,jj) @@ -2946,6 +2976,8 @@ C C Evaluate gradient. C fac=waga*rdis/dd + endif !dhpb(i).gt.0 + endif cd print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd, cd & ' waga=',waga,' fac=',fac do j=1,3 @@ -2960,6 +2992,53 @@ C Cartesian gradient in the SC vectors (ghpbx). ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j) enddo endif + else !ii.gt.nres +C write(iout,*) "before" + dd=dist(ii,jj) +C write(iout,*) "after",dd + if (constr_dist.eq.11) then + ehpb=ehpb+fordepth(i)**4.0d0 + & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) + fac=fordepth(i)**4.0d0 + & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd +C ehpb=ehpb+fordepth(i)**4*rlornmr1(dd,dhpb(i),dhpb1(i)) +C fac=fordepth(i)**4*rlornmr1prim(dd,dhpb(i),dhpb1(i))/dd +C print *,ehpb,"tu?" +C write(iout,*) ehpb,"btu?", +C & dd,dhpb(i),dhpb1(i),fordepth(i),forcon(i) +C write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj, +C & ehpb,fordepth(i),dd + else + if (dhpb1(i).gt.0.0d0) then + ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd +c write (iout,*) "alph nmr", +c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + else + rdis=dd-dhpb(i) +C Get the force constant corresponding to this distance. + waga=forcon(i) +C Calculate the contribution to energy. + ehpb=ehpb+waga*rdis*rdis +c write (iout,*) "alpha reg",dd,waga*rdis*rdis +C +C Evaluate gradient. +C + fac=waga*rdis/dd + endif + endif + do j=1,3 + ggg(j)=fac*(c(j,jj)-c(j,ii)) + enddo +cd print '(i3,3(1pe14.5))',i,(ggg(j),j=1,3) +C If this is a SC-SC distance, we need to calculate the contributions to the +C Cartesian gradient in the SC vectors (ghpbx). + if (iii.lt.ii) then + do j=1,3 + ghpbx(j,iii)=ghpbx(j,iii)-ggg(j) + ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j) + enddo + endif do j=iii,jjj-1 do k=1,3 ghpbc(k,j)=ghpbc(k,j)+ggg(k) @@ -2967,7 +3046,7 @@ C Cartesian gradient in the SC vectors (ghpbx). enddo endif enddo - ehpb=0.5D0*ehpb + if (constr_dist.ne.11) ehpb=0.5D0*ehpb return end C-------------------------------------------------------------------------- diff --git a/source/cluster/wham/src-M/readrtns.F b/source/cluster/wham/src-M/readrtns.F index 6e88ef3..06777ec 100644 --- a/source/cluster/wham/src-M/readrtns.F +++ b/source/cluster/wham/src-M/readrtns.F @@ -35,6 +35,11 @@ C pdbref=(index(controlcard,'PDBREF').gt.0) iscode=index(controlcard,'ONE_LETTER') tree=(index(controlcard,'MAKE_TREE').gt.0) + with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0 + call readi(controlcard,'CONSTR_DIST',constr_dist,0) + write (iout,*) "with_dihed_constr ",with_dihed_constr, + & " CONSTR_DIST",constr_dist + call flush(iout) min_var=(index(controlcard,'MINVAR').gt.0) plot_tree=(index(controlcard,'PLOT_TREE').gt.0) punch_dist=(index(controlcard,'PUNCH_DIST').gt.0) @@ -87,6 +92,7 @@ C include 'COMMON.CONTROL' include 'COMMON.CONTACTS' include 'COMMON.TIME1' + include 'COMMON.TORCNSTR' #ifdef MPL include 'COMMON.INFO' #endif @@ -278,6 +284,26 @@ C Convert sequence to numeric code print *,'Call Read_Bridge.' call read_bridge +C this fragment reads diheadral constrains + if (with_dihed_constr) then + + read (inp,*) ndih_constr + if (ndih_constr.gt.0) then + read (inp,*) ftors + write (iout,*) 'FTORS',ftors +C ftors is the force constant for torsional quartic constrains + read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr) + write (iout,*) + & 'There are',ndih_constr,' constraints on phi angles.' + do i=1,ndih_constr + write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i) + enddo + do i=1,ndih_constr + phi0(i)=deg2rad*phi0(i) + drange(i)=deg2rad*drange(i) + enddo + endif ! endif ndif_constr.gt.0 + endif ! with_dihed_constr nnt=1 nct=nres print *,'NNT=',NNT,' NCT=',NCT @@ -391,6 +417,11 @@ C & 'The chain contains',ns,' disulfide-bridging cysteines.' dyn_ss_mask(iss(i))=.true. enddo endif +c Read distance restraints + if (constr_dist.gt.0) then + call read_dist_constr + call hpb_partition + endif return end c----------------------------------------------------------------------------- @@ -563,6 +594,25 @@ c---------------------------------------------------------------------------- read (rekord(iread:),*) wartosc return end +C---------------------------------------------------------------------- + subroutine multreadi(rekord,lancuch,tablica,dim,default) + implicit none + integer dim,i + integer tablica(dim),default + character*(*) rekord,lancuch + character*80 aux + integer ilen,iread + external ilen + do i=1,dim + tablica(i)=default + enddo + iread=index(rekord,lancuch(:ilen(lancuch))//"=") + if (iread.eq.0) return + iread=iread+ilen(lancuch)+1 + read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim) + 10 return + end + c---------------------------------------------------------------------------- subroutine card_concat(card) include 'DIMENSIONS' @@ -664,3 +714,131 @@ C #endif return end + subroutine read_dist_constr + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' +#ifdef MPI + include 'mpif.h' +#endif + include 'COMMON.CONTROL' + include 'COMMON.CHAIN' + include 'COMMON.IOUNITS' + include 'COMMON.SBRIDGE' + integer ifrag_(2,100),ipair_(2,100) + double precision wfrag_(100),wpair_(100) + character*500 controlcard + logical lprn /.true./ + write (iout,*) "Calling read_dist_constr" +C write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup +C call flush(iout) + write(iout,*) "TU sie wywalam?" + call card_concat(controlcard) + write (iout,*) controlcard + call flush(iout) + call readi(controlcard,"NFRAG",nfrag_,0) + call readi(controlcard,"NPAIR",npair_,0) + call readi(controlcard,"NDIST",ndist_,0) + call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) + call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0) + call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0) + call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0) + call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0) + write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_ + write (iout,*) "IFRAG" + do i=1,nfrag_ + write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i) + enddo + write (iout,*) "IPAIR" + do i=1,npair_ + write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i) + enddo + call flush(iout) + do i=1,nfrag_ + if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup + if (ifrag_(2,i).gt.nstart_sup+nsup-1) + & ifrag_(2,i)=nstart_sup+nsup-1 +c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i) + call flush(iout) + if (wfrag_(i).gt.0.0d0) then + do j=ifrag_(1,i),ifrag_(2,i)-1 + do k=j+1,ifrag_(2,i) + write (iout,*) "j",j," k",k + ddjk=dist(j,k) + if (constr_dist.eq.1) then + nhpb=nhpb+1 + ihpb(nhpb)=j + jhpb(nhpb)=k + dhpb(nhpb)=ddjk + forcon(nhpb)=wfrag_(i) + else if (constr_dist.eq.2) then + if (ddjk.le.dist_cut) then + nhpb=nhpb+1 + ihpb(nhpb)=j + jhpb(nhpb)=k + dhpb(nhpb)=ddjk + forcon(nhpb)=wfrag_(i) + endif + else + nhpb=nhpb+1 + ihpb(nhpb)=j + jhpb(nhpb)=k + dhpb(nhpb)=ddjk + forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2) + endif + if (lprn) + & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ", + & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) + enddo + enddo + endif + enddo + do i=1,npair_ + if (wpair_(i).gt.0.0d0) then + ii = ipair_(1,i) + jj = ipair_(2,i) + if (ii.gt.jj) then + itemp=ii + ii=jj + jj=itemp + endif + do j=ifrag_(1,ii),ifrag_(2,ii) + do k=ifrag_(1,jj),ifrag_(2,jj) + nhpb=nhpb+1 + ihpb(nhpb)=j + jhpb(nhpb)=k + forcon(nhpb)=wpair_(i) + dhpb(nhpb)=dist(j,k) + write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ", + & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) + enddo + enddo + endif + enddo + do i=1,ndist_ + if (constr_dist.eq.11) then + read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), + & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1) + fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1) +C write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ", +C & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) + else + read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1) + endif + if (forcon(nhpb+1).gt.0.0d0) then + nhpb=nhpb+1 + if (ibecarb(i).gt.0) then + ihpb(i)=ihpb(i)+nres + jhpb(i)=jhpb(i)+nres + endif + if (dhpb(nhpb).eq.0.0d0) + & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) +C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) + write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ", + & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) + endif +C endif + enddo + call hpb_partition + call flush(iout) + return + end diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F index 131e817..662a796 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -5776,9 +5776,12 @@ c do i=1,ndih_constr else difi=0.0 endif -cd write (iout,'(2i5,4f8.3,2e14.5)') i,itori,rad2deg*phii, -cd & rad2deg*phi0(i), rad2deg*drange(i), -cd & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg) + if (energy_dec) then + write (iout,'(a6,2i5,4f8.3,2e14.5)') "edihc", + & i,itori,rad2deg*phii, + & rad2deg*phi0(i), rad2deg*drange(i), + & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg) + endif enddo cd write (iout,*) 'edihcnstr',edihcnstr return diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F index 8cf29b2..26e2a9b 100644 --- a/source/wham/src-M/energy_p_new.F +++ b/source/wham/src-M/energy_p_new.F @@ -4485,8 +4485,9 @@ c write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg) edihcnstr=edihcnstr+0.25d0*ftors*difi**4 gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 endif -! write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii, -! & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg) + write (iout,'(a6,2i5,2f8.3,2e14.5)') "edih", + & i,itori,rad2deg*phii, + & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg) enddo ! write (iout,*) 'edihcnstr',edihcnstr return @@ -4584,6 +4585,9 @@ c write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg) else difi=0.0d0 endif + write (iout,'(a6,2i5,2f8.3,2e14.5)') "edih", + & i,itori,rad2deg*phii, + & rad2deg*difi,0.25d0*ftors*difi**4 c write (iout,'(2i5,4f10.5,e15.5)') i,itori,phii,phi0(i),difi, c & drange(i),edihi ! write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii, -- 1.7.9.5