Lorenz restraints on distances

diff --git a/source/unres/src_MD/COMMON.SBRIDGE b/source/unres/src_MD/COMMON.SBRIDGE
index 91dd2cd..a146c54 100644 (file)
--- a/source/unres/src_MD/COMMON.SBRIDGE
+++ b/source/unres/src_MD/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/energy_p_new_barrier.F b/source/unres/src_MD/energy_p_new_barrier.F
index 115cb57..dbfc35e 100644 (file)
--- a/source/unres/src_MD/energy_p_new_barrier.F
+++ b/source/unres/src_MD/energy_p_new_barrier.F
@@ -4361,6 +4361,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
@@ -4396,6 +4397,16 @@ 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
+          if (energy_dec) write (iout,'(a6,2i5,f15.6,2f8.3)') 
+     &     "edisl",ii,jj,
+     &     fordepth(i)**4.0d0*rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)),
+     &     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
@@ -4413,7 +4424,8 @@ C
 C Evaluate gradient.
 C
             fac=waga*rdis/dd
-          endif  
+          endif
+          endif
           do j=1,3
             ggg(j)=fac*(c(j,jj)-c(j,ii))
           enddo
@@ -4429,6 +4441,18 @@ C
 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.0d0
+     &           *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i))
+            fac=fordepth(i)**4.0d0
+     &           *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd
+          if (energy_dec) write (iout,'(a6,2i5,f15.6,2f8.3)') 
+     7     "edisl",ii,jj,
+     &     fordepth(i)**4.0d0*rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)),
+     &     fordepth(i),dd
+c          if (energy_dec)
+c     &      write (iout,*) fac
+         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
@@ -4446,6 +4470,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
@@ -4471,7 +4496,10 @@ cgrad        enddo
           enddo
         endif
       enddo
-      ehpb=0.5D0*ehpb
+      if (constr_dist.ne.11) ehpb=0.5D0*ehpb
+c      do i=1,nres
+c        write (iout,*) "ghpbc",i,(ghpbc(j,i),j=1,3)
+c      enddo
       return
       end
 C--------------------------------------------------------------------------

diff --git a/source/unres/src_MD/gnmr1.f b/source/unres/src_MD/gnmr1.f
index 905e746..8bfc43a 100644 (file)
--- a/source/unres/src_MD/gnmr1.f
+++ b/source/unres/src_MD/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/unres/src_MD/gradient_p.F b/source/unres/src_MD/gradient_p.F
index a29f83a..5a47462 100644 (file)
--- a/source/unres/src_MD/gradient_p.F
+++ b/source/unres/src_MD/gradient_p.F
@@ -288,9 +288,9 @@ c     write (iout,*) "in cartgrad before sum_: duscdiff and duscdiffx"
       do i=1,nres
          write (iout,*) i,(duscdiff(j,i),j=1,3)
          write (iout,*) i,(duscdiffx(j,i),j=1,3)
-c        write (iout,*) "nphi+i",nphi+i," gloc",gloc(nphi+i,icg)
-c        write (iout,*) "i",i," gradc",(gradc(j,i,icg),j=1,3)
-c        write (iout,*) "i",i," gradx",(gradx(j,i,icg),j=1,3)
+         write (iout,*) "nphi+i",nphi+i," gloc",gloc(nphi+i,icg)
+         write (iout,*) "i",i," gradc",(gradc(j,i,icg),j=1,3)
+         write (iout,*) "i",i," gradx",(gradx(j,i,icg),j=1,3)
       enddo
 #endif
       call sum_gradient
@@ -369,11 +369,13 @@ cd     &    (gxcart(j,i),j=1,3),gloc(i,icg)
 #ifdef TIMING
       time_inttocart=time_inttocart+MPI_Wtime()-time01
 #endif
-cd      write (iout,*) "gcart and gxcart after int_to_cart"
-cd      do i=0,nres-1
-cd        write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
-cd     &      (gxcart(j,i),j=1,3)
-cd      enddo
+#ifdef DEBUG
+      write (iout,*) "gcart and gxcart after int_to_cart"
+      do i=0,nres-1
+        write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
+     &      (gxcart(j,i),j=1,3)
+      enddo
+#endif
 #ifdef TIMING
       time_cartgrad=time_cartgrad+MPI_Wtime()-time00
 #endif

diff --git a/source/unres/src_MD/readrtns.F b/source/unres/src_MD/readrtns.F
index 1d75b9a..8e73ba0 100644 (file)
--- a/source/unres/src_MD/readrtns.F
+++ b/source/unres/src_MD/readrtns.F
@@ -2681,8 +2681,14 @@ c            write (iout,*) "j",j," k",k
         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
@@ -2697,8 +2703,14 @@ c            write (iout,*) "j",j," k",k
       if (.not.out1file .or. me.eq.king) then
 #endif
       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
       call flush(iout)
 #ifdef MPI