introduction of shielding to cluster DEBUG mode

[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 6c4d182..dcd4d23 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
+       
+      write(iout,*) "tu,", i,j,aa_aq(itypi,itypj)-aa,
+     & bb_aq(itypi,itypj)-bb,
+     & 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
+#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
+#undef DEBUG
 c            endif
             if (calc_grad) then
 C Calculate gradient components.
@@ -2962,6 +3008,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)
@@ -3121,6 +3169,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))