From c45326004a77063de99e919c46db1a8bba1d3eba Mon Sep 17 00:00:00 2001
From: Adam Sieradzan <adasko@piasek4.chem.univ.gda.pl>
Date: Fri, 10 Jul 2015 13:46:08 +0200
Subject: [PATCH] introduction of quartic an Lonrentz constrains to WHAM

---
 source/wham/src-M/energy_p_new.F     |   58 ++++++++++++++++++++++++++++------
 source/wham/src-M/read_dist_constr.F |    6 ++++
 2 files changed, 54 insertions(+), 10 deletions(-)

diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F
index 83adc13..6e19387 100644
--- a/source/wham/src-M/energy_p_new.F
+++ b/source/wham/src-M/energy_p_new.F
@@ -2976,21 +2976,59 @@ C     & iabs(itype(jjj)).eq.1) then
           call ssbond_ene(iii,jjj,eij)
           ehpb=ehpb+2*eij
            endif
+        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
+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)
-        rdis=dd-dhpb(i)
+            rdis=dd-dhpb(i)
 C Get the force constant corresponding to this distance.
-        waga=forcon(i)
+            waga=forcon(i)
 C Calculate the contribution to energy.
-        ehpb=ehpb+waga*rdis*rdis
+            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
+            fac=waga*rdis/dd
+          endif
+          endif
+
         do j=1,3
           ggg(j)=fac*(c(j,jj)-c(j,ii))
         enddo
@@ -3010,7 +3048,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/wham/src-M/read_dist_constr.F b/source/wham/src-M/read_dist_constr.F
index 3d803bb..2e9ba3b 100644
--- a/source/wham/src-M/read_dist_constr.F
+++ b/source/wham/src-M/read_dist_constr.F
@@ -96,6 +96,11 @@ c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
         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),forcon(nhpb+1)
         if (forcon(nhpb+1).gt.0.0d0) then
           nhpb=nhpb+1
@@ -103,6 +108,7 @@ c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
           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