Commit Adam 6/29/2014

[unres.git] / source / unres / src_MD-M-newcorr / energy_p_new-sep_barrier.F

diff --git a/source/unres/src_MD-M-newcorr/energy_p_new-sep_barrier.F b/source/unres/src_MD-M-newcorr/energy_p_new-sep_barrier.F
index 6592ace..5f9a128 100644 (file)
--- a/source/unres/src_MD-M-newcorr/energy_p_new-sep_barrier.F
+++ b/source/unres/src_MD-M-newcorr/energy_p_new-sep_barrier.F
@@ -1352,7 +1352,8 @@ C-------------------------------------------------------------------------------
       dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
      &          erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
       double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
-     &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
+     &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),gmuij1(4),gmuji1(4),
+     &    gmuij2(4),gmuji2(4)
       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
      &    dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
      &    num_conti,j1,j2
@@ -1582,6 +1583,14 @@ C
             do l=1,2
               kkk=kkk+1
               muij(kkk)=mu(k,i)*mu(l,j)
+#ifdef NEWCORR
+             gmuij1(kkk)=gtb1(k,i+1)*mu(l,j)
+c             write(iout,*) 'kkk=', gtb1(k,i)*mu(l,j),gtb1(k,i),k,i
+             gmuij2(kkk)=gUb2(k,i)*mu(l,j)
+             gmuji1(kkk)=mu(k,i)*gtb1(l,j+1)
+c             write(iout,*) 'kkk=', gtb1(k,i)*mu(l,j),gtb1(l,j),l,j
+             gmuji2(kkk)=mu(k,i)*gUb2(l,j)
+#endif
             enddo
           enddo  
 cd         write (iout,*) 'EELEC: i',i,' j',j
@@ -1748,7 +1757,40 @@ C Contribution to the local-electrostatic energy coming from the i-j pair
           eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
      &     +a33*muij(4)
 cd          write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
+C Calculate patrial derivative for theta angle
+#ifdef NEWCORR
+         geel_loc_ij=a22*gmuij1(1)
+     &     +a23*gmuij1(2)
+     &     +a32*gmuij1(3)
+     &     +a33*gmuij1(4)
+c         write(iout,*) "derivative over thatai"
+c         write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
+c     &   a33*gmuij1(4)
+         gloc(nphi+i,icg)=gloc(nphi+i,icg)+
+     &      geel_loc_ij*wel_loc
+c         write(iout,*) "derivative over thatai-1"
+c         write(iout,*) a22*gmuij2(1), a23*gmuij2(2) ,a32*gmuij2(3),
+c     &   a33*gmuij2(4)
+         geel_loc_ij=a22*gmuij2(1)+a23*gmuij2(2)+a32*gmuij2(3)
+     &     +a33*gmuij2(4)
+         gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
+     &      geel_loc_ij*wel_loc
+         geel_loc_ji=a22*gmuji1(1)+a23*gmuji1(2)+a32*gmuji1(3)
+     &     +a33*gmuji1(4)
+c         write(iout,*) "derivative over thataj"
+c         write(iout,*) a22*gmuji1(1), a23*gmuji1(2) ,a32*gmuji1(3),
+c     &   a33*gmuji1(4)
 
+         gloc(nphi+j,icg)=gloc(nphi+j,icg)+
+     &      geel_loc_ji*wel_loc
+         geel_loc_ji=a22*gmuji2(1)+a23*gmuji2(2)+a32*gmuji2(3)
+     &     +a33*gmuji2(4)
+c         write(iout,*) "derivative over thataj-1"
+c         write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3),
+c     &   a33*gmuji2(4)
+         gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+
+     &      geel_loc_ji*wel_loc
+#endif
           if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
      &            'eelloc',i,j,eel_loc_ij