unres src_MD-M now HOMOL energy

[unres.git] / source / unres / src_MD-M / energy_p_new_barrier.F

diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F
index 526b4ab..ea2fe9e 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -6895,7 +6895,7 @@ c
 c
 
 
-      do i=1,19
+      do i=1,max_template
         distancek(i)=9999999.9
       enddo
 
@@ -6914,9 +6914,13 @@ 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
 c           write(iout,*) ii,k,i,j,l_homo(k,ii),dij,odl(k,ii)
-           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
@@ -6955,7 +6959,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)/ 
@@ -7110,8 +7122,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) ! waga_angle rmvd from Gaussian argument
+#endif
 c         kat3=-0.5d0*waga_angle*dih_diff(k)**2*sigma_dih(k,i)
           gdih(k)=dexp(kat3)
           kat2=kat2+gdih(k)
@@ -7138,7 +7153,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)  ! waga_angle rmvd
+#endif
 c         sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i)*waga_angle
           sum_sgdih=sum_sgdih+sgdih
         enddo
@@ -7220,7 +7239,7 @@ c
           utheta_i=-0.5d0*theta_diff(k)**2*sigma_theta(k,i) ! waga_theta rmvd from Gaussian argument
 c         utheta_i=-0.5d0*waga_theta*theta_diff(k)**2*sigma_theta(k,i) ! waga_theta?
           gtheta(k)=dexp(utheta_i) ! + min_utheta_i?
-          gutheta_i=gutheta_i+dexp(utheta_i)   ! Sum of Gaussians (pk)
+          gutheta_i=gutheta_i+gtheta(k)  ! Sum of Gaussians (pk)
 c         Gradient for single Gaussian restraint in subr Econstr_back
 c         dutheta(j-2)=dutheta(j-2)+wfrag_back(1,i,iset)*dtheta_i/(ii-1)
 c
@@ -7292,7 +7311,9 @@ c
 c         usc_diff(i)=-0.5d0*waga_d*(dxx**2+dyy**2+dzz**2)*sigma_d(k,i) ! waga_d?
 c         uscdiffk(k)=usc_diff(i)
           guscdiff2(k)=dexp(usc_diff_i) ! without min_scdiff
-          guscdiff(i)=guscdiff(i)+dexp(usc_diff_i)   !Sum of Gaussians (pk)
+c          write(iout,*) "i",i," k",k," sigma_d",sigma_d(k,i),
+c     &       " guscdiff2",guscdiff2(k)
+          guscdiff(i)=guscdiff(i)+guscdiff2(k)  !Sum of Gaussians (pk)
 c          write (iout,'(i5,6f10.5)') j,xxtab(j),yytab(j),zztab(j),
 c     &      xxref(j),yyref(j),zzref(j)
         enddo