From 6e5f3f4866eb982f17fb6d5c2bb30bd7bf0585be Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Fri, 3 Jul 2015 13:16:13 +0200 Subject: [PATCH] WeFold Lorentzian like potentials introduction --- source/unres/src_MD-M/COMMON.SBRIDGE | 3 +- source/unres/src_MD-M/energy_p_new_barrier.F | 58 ++++++++++++++++++++++++-- source/unres/src_MD-M/gnmr1.f | 30 +++++++++++++ source/unres/src_MD-M/readrtns_CSA.F | 23 ++++++++-- 4 files changed, 107 insertions(+), 7 deletions(-) diff --git a/source/unres/src_MD-M/COMMON.SBRIDGE b/source/unres/src_MD-M/COMMON.SBRIDGE index 1d69501..3279b28 100644 --- a/source/unres/src_MD-M/COMMON.SBRIDGE +++ b/source/unres/src_MD-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 common /links/ dhpb(maxdim),dhpb1(maxdim),forcon(maxdim), + & fordepth(maxdim), & ihpb(maxdim),jhpb(maxdim),ibecarb(maxdim),nhpb double precision weidis common /restraints/ weidis 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 700fb07..cf9ddcb 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -4084,6 +4084,7 @@ C include 'COMMON.VAR' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' + include 'COMMON.CONTROL' dimension ggg(3) ehpb=0.0D0 cd write(iout,*)'edis: nhpb=',nhpb,' fbr=',fbr @@ -4116,21 +4117,72 @@ C 15/02/13 CC dynamic SSbond - additional check ehpb=ehpb+2*eij endif cd write (iout,*) "eij",eij +cd & ' waga=',waga,' fac=',fac + 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 + ehpb=ehpb+fordepth(i)**4 + & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) + fac=fordepth(i)**4 + & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(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,*) "beta nmr", +c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + else + dd=dist(ii,jj) + 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,*) "beta 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 + do j=1,3 + ghpbx(j,iii)=ghpbx(j,iii)-ggg(j) + ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j) + enddo + do k=1,3 + ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k) + ghpbc(k,iii)=ghpbc(k,iii)-ggg(k) + enddo else C Calculate the distance between the two points and its difference from the C target distance. dd=dist(ii,jj) + if (constr_dist.eq.11) then + ehpb=ehpb+fordepth(i)**4*rlornmr1(dd,dhpb(i),dhpb1(i)) + fac=fordepth(i)**4*rlornmr1prim(dd,dhpb(i),dhpb1(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 -cd print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd, -cd & ' waga=',waga,' fac=',fac + endif + endif do j=1,3 ggg(j)=fac*(c(j,jj)-c(j,ii)) enddo @@ -4154,7 +4206,7 @@ cgrad enddo enddo endif enddo - ehpb=0.5D0*ehpb + if (constr_dist.ne.11) ehpb=0.5D0*ehpb return end C-------------------------------------------------------------------------- diff --git a/source/unres/src_MD-M/gnmr1.f b/source/unres/src_MD-M/gnmr1.f index 905e746..f44e404 100644 --- a/source/unres/src_MD-M/gnmr1.f +++ b/source/unres/src_MD-M/gnmr1.f @@ -41,3 +41,33 @@ c------------------------------------------------------------------------------- return end c--------------------------------------------------------------------------------- + double precision function rlornmr1(y,ymin,ymax,sigma) + implicit none + double precision y,ymin,ymax,sigma + double precision wykl /4.0d0/ + if (y.lt.ymin) then + rlornmr1=(ymin-y)**wykl/((ymin-y)**wykl+sigma**wykl) + else if (y.gt.ymax) then + rlornmr1=(y-ymax)**wykl/((y-ymax)**wykl+sigma**wykl) + else + rlornmr1=0.0d0 + endif + return + end +c------------------------------------------------------------------------------ + double precision function rlornmr1prim(y,ymin,ymax,sigma) + implicit none + double precision y,ymin,ymax,sigma + double precision wykl /4.0d0/ + if (y.lt.ymin) then + rlornmr1prim=-(ymin-y)**(wykl-1)*sigma**wykl*wykl/ + & ((ymin-y)**wykl+sigma**wykl) + else if (y.gt.ymax) then + rlornmr1prim=(y-ymax)**(wykl-1)*sigma**wykl*wykl/ + & ((y-ymax)**wykl+sigma**wykl) + else + rlornmr1prim=0.0d0 + endif + return + end + diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index 3987a4d..f53da12 100644 --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@ -2352,10 +2352,27 @@ c write (iout,*) "j",j," k",k endif enddo do i=1,ndist_ - read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1) + 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) + else + read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), + & ibecarb(i),forcon(nhpb+1) + endif if (forcon(nhpb+1).gt.0.0d0) then nhpb=nhpb+1 - dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) + 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)) + endif +C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1) +C if (forcon(nhpb+1).gt.0.0d0) then +C nhpb=nhpb+1 +C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) #ifdef MPI if (.not.out1file .or. me.eq.king) & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ", @@ -2364,7 +2381,7 @@ c write (iout,*) "j",j," k",k write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #endif - endif + enddo call flush(iout) return -- 1.7.9.5