Lorentzian restraints on distances in cluster_wham

diff --git a/source/cluster/wham/src/energy_p_new.F b/source/cluster/wham/src/energy_p_new.F
index 7ee0e64..f0d81cd 100644 (file)
--- 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 (file)
--- 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 (file)
--- 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 (file)
--- 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)