bug fix after Ana and cluster lipid (still in progress)

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

diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F
index 2b99394..fb18913 100644 (file)
--- a/source/wham/src-M/energy_p_new.F
+++ b/source/wham/src-M/energy_p_new.F
@@ -224,6 +224,11 @@ c detecting NaNQ
 #ifdef MPL
 c     endif
 #endif
+#define DEBUG
+#ifdef DEBUG
+      call enerprint(energia,fact)
+#endif
+#undef DEBUG
       if (calc_grad) then
 C
 C Sum up the components of the Cartesian gradient.
@@ -265,12 +270,29 @@ C
      &                wturn6*fact(5)*gcorr6_turn(j,i)+
      &                wsccor*fact(2)*gsccorc(j,i)
      &               +wliptran*gliptranc(j,i)
+     &                 +welec*gshieldc(j,i)
+     &                 +welec*gshieldc_loc(j,i)
+     &                 +wcorr*gshieldc_ec(j,i)
+     &                 +wcorr*gshieldc_loc_ec(j,i)
+     &                 +wturn3*gshieldc_t3(j,i)
+     &                 +wturn3*gshieldc_loc_t3(j,i)
+     &                 +wturn4*gshieldc_t4(j,i)
+     &                 +wturn4*gshieldc_loc_t4(j,i)
+     &                 +wel_loc*gshieldc_ll(j,i)
+     &                 +wel_loc*gshieldc_loc_ll(j,i)
+
           gradx(j,i,icg)=fact(1)*wsc*gvdwx(j,i)
      &                 +fact(1)*wscp*gradx_scp(j,i)+
      &                  wbond*gradbx(j,i)+
      &                  wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+
      &                  wsccor*fact(2)*gsccorx(j,i)
      &                 +wliptran*gliptranx(j,i)
+     &                 +welec*gshieldx(j,i)
+     &                 +wcorr*gshieldx_ec(j,i)
+     &                 +wturn3*gshieldx_t3(j,i)
+     &                 +wturn4*gshieldx_t4(j,i)
+     &                 +wel_loc*gshieldx_ll(j,i)
+
 
         endif
         enddo
@@ -309,12 +331,29 @@ C
      &                wturn6*fact(5)*gcorr6_turn(j,i)+
      &                wsccor*fact(2)*gsccorc(j,i)
      &               +wliptran*gliptranc(j,i)
+     &                 +welec*gshieldc(j,i)
+     &                 +welec*gshieldc_loc(j,i)
+     &                 +wcorr*gshieldc_ec(j,i)
+     &                 +wcorr*gshieldc_loc_ec(j,i)
+     &                 +wturn3*gshieldc_t3(j,i)
+     &                 +wturn3*gshieldc_loc_t3(j,i)
+     &                 +wturn4*gshieldc_t4(j,i)
+     &                 +wturn4*gshieldc_loc_t4(j,i)
+     &                 +wel_loc*gshieldc_ll(j,i)
+     &                 +wel_loc*gshieldc_loc_ll(j,i)
+
           gradx(j,i,icg)=fact(1)*wsc*gvdwx(j,i)+
      &                  fact(1)*wscp*gradx_scp(j,i)+
      &                  wbond*gradbx(j,i)+
      &                  wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+
      &                  wsccor*fact(1)*gsccorx(j,i)
      &                 +wliptran*gliptranx(j,i)
+     &                 +welec*gshieldx(j,i)
+     &                 +wcorr*gshieldx_ec(j,i)
+     &                 +wturn3*gshieldx_t3(j,i)
+     &                 +wturn4*gshieldx_t4(j,i)
+     &                 +wel_loc*gshieldx_ll(j,i)
+
          endif
         enddo
 #endif
@@ -1034,7 +1073,12 @@ C lipbufthick is thickenes of lipid buffore
      &  +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
       bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
      &  +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
-C      write(iout,*) "tu,", i,j,aa_lip(itypi,itypj),bb_lip(itypi,itypj)
+C       if (aa.ne.aa_aq(itypi,itypj)) then
+       
+C      write(iout,*) "tu,", i,j,aa_aq(itypi,itypj)-aa,
+C     & bb_aq(itypi,itypj)-bb,
+C     & sslipi,sslipj
+C         endif
 
 C        write(iout,*),aa,aa_lip(itypi,itypj),aa_aq(itypi,itypj)
 C checking the distance
@@ -1118,6 +1162,7 @@ c     &         aux*e2/eps(itypi,itypj)
 c            if (lprn) then
             sigm=dabs(aa/bb)**(1.0D0/6.0D0)
             epsi=bb**2/aa
+C#define DEBUG
 #ifdef DEBUG
             write (iout,'(2(a3,i3,2x),17(0pf7.3))')
      &        restyp(itypi),i,restyp(itypj),j,
@@ -1127,6 +1172,7 @@ c            if (lprn) then
      &        evdwij
              write (iout,*) "partial sum", evdw, evdw_t
 #endif
+C#undef DEBUG
 c            endif
             if (calc_grad) then
 C Calculate gradient components.
@@ -2227,6 +2273,12 @@ c          write (iout,*) "i",i,iteli," j",j,itelj," eesij",eesij
 C 12/26/95 - for the evaluation of multi-body H-bonding interactions
           ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
           if (shield_mode.gt.0) then
+C#define DEBUG
+#ifdef DEBUG
+          write(iout,*) "ees_compon",i,j,el1,el2,
+     &    fac_shield(i),fac_shield(j)
+#endif
+C#undef DEBUG
 C          fac_shield(i)=0.4
 C          fac_shield(j)=0.6
           el1=el1*fac_shield(i)**2*fac_shield(j)**2
@@ -2962,6 +3014,8 @@ C Derivatives due to the contact function
      &          *fac_shield(i)*fac_shield(j)
 
                   gacontp_hb3(k,num_conti,i)=gggp(k)
+     &          *fac_shield(i)*fac_shield(j)
+
                   gacontm_hb1(k,num_conti,i)=ghalfm
      &              +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
      &              + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
@@ -3040,7 +3094,7 @@ C end of changes suggested by Ana
 C     &    .or. itype(i+5).eq.ntyp1
 C     &    .or. itype(i).eq.ntyp1
 C     &    .or. itype(i-1).eq.ntyp1
-     &    ) goto 178
+     &    ) goto 179
 
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 C
@@ -3121,6 +3175,7 @@ C Derivatives in gamma(i)
         call transpose2(auxmat2(1,1),pizda(1,1))
         call matmat2(a_temp(1,1),pizda(1,1),pizda(1,1))
         gel_loc_turn3(i)=gel_loc_turn3(i)+0.5d0*(pizda(1,1)+pizda(2,2))
+     &   *fac_shield(i)*fac_shield(j)
 C Derivatives in gamma(i+1)
         call matmat2(EUg(1,1,i+1),EUgder(1,1,i+2),auxmat2(1,1))
         call transpose2(auxmat2(1,1),pizda(1,1))
@@ -3169,6 +3224,7 @@ C Cartesian derivatives
 
         enddo
         endif
+  179 continue
       else if (j.eq.i+3 .and. itype(i+2).ne.ntyp1) then
       if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
 C changes suggested by Ana to avoid out of bounds