HOMOL wham/cluster new dihed constrain cos function

[unres.git] / source / wham / src / energy_p_new.F

diff --git a/source/wham/src/energy_p_new.F b/source/wham/src/energy_p_new.F
index 13267da..1ac587d 100644 (file)
--- a/source/wham/src/energy_p_new.F
+++ b/source/wham/src/energy_p_new.F
@@ -2948,12 +2948,14 @@ C Evaluate bridge-strain energy and its gradient in virtual-bond and SC vectors.
 C
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
+      include 'DIMENSIONS.FREE'
       include 'COMMON.SBRIDGE'
       include 'COMMON.CHAIN'
       include 'COMMON.DERIV'
       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
@@ -2986,6 +2988,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
@@ -3003,7 +3011,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
@@ -3019,6 +3028,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
@@ -3036,6 +3059,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
@@ -3056,7 +3080,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--------------------------------------------------------------------------
@@ -3200,8 +3224,12 @@ c     write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
          j = jres_homo(ii)
          dij=dist(i,j)
 c        write (iout,*) "dij(",i,j,") =",dij
+         nexl=0
          do k=1,constr_homology
-           if(.not.l_homo(k,ii)) cycle
+           if(.not.l_homo(k,ii)) then
+              nexl=nexl+1
+              cycle
+           endif
            distance(k)=odl(k,ii)-dij
 c          write (iout,*) "distance(",k,") =",distance(k)
 c
@@ -3239,7 +3267,15 @@ c        write (iout,* )"min_odl",min_odl
          write (iout,*) "distancek",(distancek(k),k=1,constr_homology)
          write (iout,* )"min_odl",min_odl
 #endif
+#ifdef OLDRESTR
          odleg2=0.0d0
+#else
+         if (waga_dist.ge.0.0d0) then
+           odleg2=nexl
+         else
+           odleg2=0.0d0
+         endif
+#endif
          do k=1,constr_homology
 c Nie wiem po co to liczycie jeszcze raz!
 c            odleg3=-waga_dist(iset)*((distance(i,j,k)**2)/ 
@@ -3384,9 +3420,13 @@ c      write (iout,*) idihconstr_start_homo,idihconstr_end_homo
       enddo
 #endif
       do i=idihconstr_start_homo,idihconstr_end_homo
+#ifdef OLDRESTR
         kat2=0.0d0
+#else
+        kat2=nexl
+#endif
 c        betai=beta(i,i+1,i+2,i+3)
-        betai = phi(i+3)
+        betai = phi(i)
 c       write (iout,*) "betai =",betai
         do k=1,constr_homology
           dih_diff(k)=pinorm(dih(k,i)-betai)
@@ -3395,8 +3435,11 @@ c          if (dih_diff(i,k).gt.3.14159) dih_diff(i,k)=
 c     &                                   -(6.28318-dih_diff(i,k))
 c          if (dih_diff(i,k).lt.-3.14159) dih_diff(i,k)=
 c     &                                   6.28318+dih_diff(i,k)
-
+#ifdef OLD_DIHED
           kat3=-0.5d0*dih_diff(k)**2*sigma_dih(k,i) ! waga_angle rmvd from Gaussian argument
+#else
+          kat3=(dcos(dih_diff(k))-1)*sigma_dih(k,i)
+#endif
 c         kat3=-0.5d0*waga_angle*dih_diff(k)**2*sigma_dih(k,i)
           gdih(k)=dexp(kat3)
           kat2=kat2+gdih(k)
@@ -3424,7 +3467,11 @@ c ----------------------------------------------------------------------
         sum_gdih=kat2
         sum_sgdih=0.0d0
         do k=1,constr_homology
+#ifdef OLD_DIHED
           sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i)  ! waga_angle rmvd
+#else
+          sgdih=-gdih(k)*dsin(dih_diff(k))*sigma_dih(k,i)
+#endif
 c         sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i)*waga_angle
           sum_sgdih=sum_sgdih+sgdih
         enddo
@@ -3492,7 +3539,11 @@ c
 c Deviation of theta angles wrt constr_homology ref structures
 c
         utheta_i=0.0d0 ! argument of Gaussian for single k
+#ifdef OLDRESTR
         gutheta_i=0.0d0 ! Sum of Gaussians over constr_homology ref structures
+#else
+        gutheta_i=nexl
+#endif
 c       do j=ifrag_back(1,i,iset)+2,ifrag_back(2,i,iset) ! original loop
 c       over residues in a fragment
 c       write (iout,*) "theta(",i,")=",theta(i)
@@ -3560,7 +3611,11 @@ c     write (iout,*) "waga_d",waga_d
 #endif
       do i=loc_start,loc_end
         usc_diff_i=0.0d0 ! argument of Gaussian for single k
+#ifdef OLDRESTR
         guscdiff(i)=0.0d0 ! Sum of Gaussians over constr_homology ref structures
+#else
+        guscdiff(i)=nexl
+#endif
 c       do j=ifrag_back(1,i,iset)+1,ifrag_back(2,i,iset)-1 ! Econstr_back legacy
 c       write(iout,*) "xxtab, yytab, zztab"
 c       write(iout,'(i5,3f8.2)') i,xxtab(i),yytab(i),zztab(i)