From 39b10050c05749106eb2c7ca50deea0b59505964 Mon Sep 17 00:00:00 2001 From: Cezary Czaplewski Date: Wed, 24 Aug 2016 22:20:32 +0200 Subject: [PATCH] Lorentzian restraints on distances in cluster_wham --- source/cluster/wham/src/energy_p_new.F | 27 ++++++++++++++++-- source/cluster/wham/src/gnmr1.f | 30 ++++++++++++++++++++ .../cluster/wham/src/include_unres/COMMON.SBRIDGE | 3 +- source/cluster/wham/src/readrtns.F | 12 ++++++++ 4 files changed, 69 insertions(+), 3 deletions(-) diff --git a/source/cluster/wham/src/energy_p_new.F b/source/cluster/wham/src/energy_p_new.F index 7ee0e64..f0d81cd 100644 --- a/source/cluster/wham/src/energy_p_new.F +++ b/source/cluster/wham/src/energy_p_new.F @@ -2876,6 +2876,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 @@ -2908,6 +2909,12 @@ cd write (iout,*) "eij",eij 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.0d0 + & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) + fac=fordepth(i)**4.0d0 + & *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 @@ -2925,7 +2932,8 @@ C C Evaluate gradient. C fac=waga*rdis/dd - endif + endif !end dhpb1(i).gt.0 + endif !end const_dist=11 do j=1,3 ggg(j)=fac*(c(j,jj)-c(j,ii)) enddo @@ -2941,6 +2949,20 @@ C C Calculate the distance between the two points and its difference from the C target distance. 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 @@ -2958,6 +2980,7 @@ C Evaluate gradient. C fac=waga*rdis/dd endif + endif cd print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd, cd & ' waga=',waga,' fac=',fac do j=1,3 @@ -2978,7 +3001,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/gnmr1.f b/source/cluster/wham/src/gnmr1.f index 905e746..8bfc43a 100644 --- a/source/cluster/wham/src/gnmr1.f +++ b/source/cluster/wham/src/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)**2 + else if (y.gt.ymax) then + rlornmr1prim=(y-ymax)**(wykl-1)*sigma**wykl*wykl/ + & ((y-ymax)**wykl+sigma**wykl)**2 + else + rlornmr1prim=0.0d0 + endif + return + end + diff --git a/source/cluster/wham/src/include_unres/COMMON.SBRIDGE b/source/cluster/wham/src/include_unres/COMMON.SBRIDGE index f866aa7..b68d0e3 100644 --- a/source/cluster/wham/src/include_unres/COMMON.SBRIDGE +++ b/source/cluster/wham/src/include_unres/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),ibecarb(maxdim),nhpb double precision weidis common /restraints/ weidis diff --git a/source/cluster/wham/src/readrtns.F b/source/cluster/wham/src/readrtns.F index 4051df0..fac2508 100644 --- a/source/cluster/wham/src/readrtns.F +++ b/source/cluster/wham/src/readrtns.F @@ -880,8 +880,14 @@ c call flush(iout) 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) + 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 if (ibecarb(i).gt.0) then @@ -893,8 +899,14 @@ c call flush(iout) endif enddo do i=1,nhpb + if (constr_dist.eq.11) then + write (iout,'(a,3i5,2f8.2,i2,2f10.1)') "+dist.constr11 ", + & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i), + & fordepth(i) + else write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ", & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i) + endif enddo #ifdef AIX call flush_(iout) -- 1.7.9.5