working cluser DEBUG OFF
[unres.git] / source / wham / src-M / energy_p_new.F
index eb6146a..730500e 100644 (file)
@@ -224,6 +224,11 @@ c detecting NaNQ
 #ifdef MPL
 c     endif
 #endif
+C#define DEBUG
+#ifdef DEBUG
+      call enerprint(energia,fact)
+#endif
+C#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
@@ -895,7 +934,7 @@ C
       include 'COMMON.SBRIDGE'
       logical lprn
       common /srutu/icall
-      integer icant
+      integer icant,xshift,yshift,zshift
       external icant
       do i=1,210
         do j=1,2
@@ -1034,6 +1073,13 @@ 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       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
       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
@@ -1116,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,
@@ -1125,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.
@@ -2098,15 +2146,15 @@ C      write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
         gcorr_loc(i)=0.0d0
       enddo
       do i=iatel_s,iatel_e
-          if (i.eq.1) then 
+C          if (i.eq.1) then 
            if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
-     &  .or. itype(i+2).eq.ntyp1) cycle
-          else
-        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
-     &  .or. itype(i+2).eq.ntyp1
-     &  .or. itype(i-1).eq.ntyp1
+C     &  .or. itype(i+2).eq.ntyp1) cycle
+C          else
+C        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
+C     &  .or. itype(i+2).eq.ntyp1
+C     &  .or. itype(i-1).eq.ntyp1
      &) cycle
-         endif
+C         endif
         if (itel(i).eq.0) goto 1215
         dxi=dc(1,i)
         dyi=dc(2,i)
@@ -2126,16 +2174,16 @@ C      write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
         num_conti=0
 C        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
         do j=ielstart(i),ielend(i)
-          if (j.eq.1) then
-           if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
-     & .or.itype(j+2).eq.ntyp1
-     &) cycle  
-          else     
+          if (j.le.1) cycle
+C           if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
+C     & .or.itype(j+2).eq.ntyp1
+C     &) cycle  
+C          else     
           if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
-     & .or.itype(j+2).eq.ntyp1
-     & .or.itype(j-1).eq.ntyp1
+C     & .or.itype(j+2).eq.ntyp1
+C     & .or.itype(j-1).eq.ntyp1
      &) cycle
-         endif
+C         endif
 C
 C) cycle
           if (itel(j).eq.0) goto 1216
@@ -2225,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
@@ -2960,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)
@@ -3028,6 +3084,18 @@ C Third- and fourth-order contributions from turns
      &    aggj(3,4),aggj1(3,4),a_temp(2,2)
       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,j1,j2
       if (j.eq.i+2) then
+      if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
+C changes suggested by Ana to avoid out of bounds
+C     & .or.((i+5).gt.nres)
+C     & .or.((i-1).le.0)
+C end of changes suggested by Ana
+     &    .or. itype(i+2).eq.ntyp1
+     &    .or. itype(i+3).eq.ntyp1
+C     &    .or. itype(i+5).eq.ntyp1
+C     &    .or. itype(i).eq.ntyp1
+C     &    .or. itype(i-1).eq.ntyp1
+     &    ) goto 179
+
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 C
 C               Third-order contributions
@@ -3107,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))
@@ -3155,17 +3224,19 @@ 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
-     & .or.((i+5).gt.nres)
-     & .or.((i-1).le.0)
+C     & .or.((i+5).gt.nres)
+C     & .or.((i-1).le.0)
 C end of changes suggested by Ana
      &    .or. itype(i+3).eq.ntyp1
      &    .or. itype(i+4).eq.ntyp1
-     &    .or. itype(i+5).eq.ntyp1
+C     &    .or. itype(i+5).eq.ntyp1
      &    .or. itype(i).eq.ntyp1
-     &    .or. itype(i-1).eq.ntyp1) goto 178
+C     &    .or. itype(i-1).eq.ntyp1
+     &    ) goto 178
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 C
 C               Fourth-order contributions
@@ -3851,6 +3922,7 @@ C     &       gnmr1(vbld(i),-1.0d0,distchainmax)
 C        else
          if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then
         diff = vbld(i)-vbldpDUM
+C         write(iout,*) i,diff
          else
           diff = vbld(i)-vbldp0
 c          write (iout,*) i,vbld(i),vbldp0,diff,AKP*diff*diff