working gradient for shielding - even the fine grain mode

diff --git a/source/unres/src_MD-M/COMMON.DERIV b/source/unres/src_MD-M/COMMON.DERIV
index 17e3184..9b3b081 100644 (file)
--- a/source/unres/src_MD-M/COMMON.DERIV
+++ b/source/unres/src_MD-M/COMMON.DERIV
@@ -10,9 +10,16 @@
      & gvdwpp(3,-1:maxres),gvdwc_scpp(3,-1:maxres),
      & gliptranc(3,-1:maxres),
      & gliptranx(3,-1:maxres),
-     & gshieldx(3,-1:maxres),
-     & gshieldc(3,-1:maxres),
+     & gshieldx(3,-1:maxres), gshieldc(3,-1:maxres),
      & gshieldc_loc(3,-1:maxres),
+     & gshieldx_ec(3,-1:maxres), gshieldc_ec(3,-1:maxres),
+     & gshieldc_loc_ec(3,-1:maxres),
+     & gshieldx_t3(3,-1:maxres), gshieldc_t3(3,-1:maxres),
+     & gshieldc_loc_t3(3,-1:maxres),
+     & gshieldx_t4(3,-1:maxres), gshieldc_t4(3,-1:maxres),
+     & gshieldc_loc_t4(3,-1:maxres),
+     & gshieldx_ll(3,-1:maxres), gshieldc_ll(3,-1:maxres),
+     & gshieldc_loc_ll(3,-1:maxres),
      & gradafm(3,-1:maxres),
      & gradx_scp(3,-1:maxres),gvdwc_scp(3,-1:maxres),
      & ghpbx(3,-1:maxres),

diff --git a/source/unres/src_MD-M/COMMON.SHIELD b/source/unres/src_MD-M/COMMON.SHIELD
index 65d7891..77126a2 100644 (file)
--- a/source/unres/src_MD-M/COMMON.SHIELD
+++ b/source/unres/src_MD-M/COMMON.SHIELD
@@ -1,12 +1,14 @@
        double precision VSolvSphere,VSolvSphere_div,long_r_sidechain,
      & short_r_sidechain,fac_shield,grad_shield_side,grad_shield,
      & buff_shield            
-       integer  ishield_list,shield_list
+       integer  ishield_list,shield_list,ees0plist
        common /shield/ VSolvSphere,VSolvSphere_div,buff_shield,
      & long_r_sidechain(ntyp),
      & short_r_sidechain(ntyp),fac_shield(maxres),
      & grad_shield_side(3,15,-1:maxres),grad_shield(3,-1:maxres),
      &  grad_shield_loc(3,15,-1:maxres),
-     & ishield_list(maxres),shield_list(15,maxres)
+     & ishield_list(maxres),shield_list(15,maxres),
+     & ees0plist(maxcont,maxres)
+
 
        

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 5cb9a2a..a889c75 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -547,6 +547,11 @@ c      enddo
      &                +wliptran*gliptranc(j,i)
      &                +gradafm(j,i)
      &                 +welec*gshieldc(j,i)
+     &                 +wcorr*gshieldc_ec(j,i)
+     &                 +wturn3*gshieldc_t3(j,i)
+     &                 +wturn4*gshieldc_t4(j,i)
+     &                 +wel_loc*gshieldc_ll(j,i)
+
 
         enddo
       enddo 
@@ -566,6 +571,10 @@ c      enddo
      &                +wliptran*gliptranc(j,i)
      &                +gradafm(j,i)
      &                 +welec*gshieldc(j,i)
+     &                 +wcorr*gshieldc_ec(j,i)
+     &                 +wturn4*gshieldc_t4(j,i)
+     &                 +wel_loc*gshieldc_ll(j,i)
+
 
         enddo
       enddo 
@@ -713,13 +722,25 @@ C          print *,gradafm(1,13),"AFM"
      &                +gradafm(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)
+
+
+
+
 
 
 #else
           gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+
      &                wel_loc*gel_loc(j,i)+
      &                0.5d0*(wscp*gvdwc_scpp(j,i)+
-     &                welec*gelc_long(j,i)
+     &                welec*gelc_long(j,i)+
      &                wel_loc*gel_loc_long(j,i)+
      &                wcorr*gcorr_long(j,i)+
      &                wcorr5*gradcorr5_long(j,i)+
@@ -738,6 +759,17 @@ C          print *,gradafm(1,13),"AFM"
      &                +gradafm(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)
+
+
+
 
 
 #endif
@@ -748,6 +780,13 @@ C          print *,gradafm(1,13),"AFM"
      &                 +wscloc*gsclocx(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)
+
+
+
         enddo
       enddo 
 #ifdef DEBUG
@@ -4164,15 +4203,71 @@ cgrad            endif
 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)
+          if (shield_mode.eq.0) then 
+           fac_shield(i)=1.0
+           fac_shield(j)=1.0
+C          else
+C           fac_shield(i)=0.4
+C           fac_shield(j)=0.6
+          endif
+          eel_loc_ij=eel_loc_ij
+     &    *fac_shield(i)*fac_shield(j)
+C Now derivative over eel_loc
+          if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
+     &  (shield_mode.gt.0)) then
+C          print *,i,j     
+
+          do ilist=1,ishield_list(i)
+           iresshield=shield_list(ilist,i)
+           do k=1,3
+           rlocshield=grad_shield_side(k,ilist,i)*eel_loc_ij
+     &                                          /fac_shield(i)
+C     &      *2.0
+           gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
+     &              rlocshield
+     & +grad_shield_loc(k,ilist,i)*eel_loc_ij/fac_shield(i)
+            gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
+     &      +rlocshield
+           enddo
+          enddo
+          do ilist=1,ishield_list(j)
+           iresshield=shield_list(ilist,j)
+           do k=1,3
+           rlocshield=grad_shield_side(k,ilist,j)*eel_loc_ij
+     &                                       /fac_shield(j)
+C     &     *2.0
+           gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
+     &              rlocshield
+     & +grad_shield_loc(k,ilist,j)*eel_loc_ij/fac_shield(j)
+           gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
+     &             +rlocshield
+
+           enddo
+          enddo
+
+          do k=1,3
+            gshieldc_ll(k,i)=gshieldc_ll(k,i)+
+     &              grad_shield(k,i)*eel_loc_ij/fac_shield(i)
+            gshieldc_ll(k,j)=gshieldc_ll(k,j)+
+     &              grad_shield(k,j)*eel_loc_ij/fac_shield(j)
+            gshieldc_ll(k,i-1)=gshieldc_ll(k,i-1)+
+     &              grad_shield(k,i)*eel_loc_ij/fac_shield(i)
+            gshieldc_ll(k,j-1)=gshieldc_ll(k,j-1)+
+     &              grad_shield(k,j)*eel_loc_ij/fac_shield(j)
+           enddo
+           endif
+
+
 c          write (iout,*) 'i',i,' j',j,itype(i),itype(j),
 c     &                     ' eel_loc_ij',eel_loc_ij
 C          write(iout,*) 'muije=',i,j,muij(1),muij(2),muij(3),muij(4)
 C Calculate patrial derivative for theta angle
 #ifdef NEWCORR
-         geel_loc_ij=a22*gmuij1(1)
+         geel_loc_ij=(a22*gmuij1(1)
      &     +a23*gmuij1(2)
      &     +a32*gmuij1(3)
-     &     +a33*gmuij1(4)         
+     &     +a33*gmuij1(4))
+     &    *fac_shield(i)*fac_shield(j)
 c         write(iout,*) "derivative over thatai"
 c         write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
 c     &   a33*gmuij1(4) 
@@ -4188,6 +4283,8 @@ c     &   a33*gmuij2(4)
      &     +a33*gmuij2(4)
          gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
      &      geel_loc_ij*wel_loc
+     &    *fac_shield(i)*fac_shield(j)
+
 c  Derivative over j residue
          geel_loc_ji=a22*gmuji1(1)
      &     +a23*gmuji1(2)
@@ -4199,6 +4296,8 @@ c     &   a33*gmuji1(4)
 
         gloc(nphi+j,icg)=gloc(nphi+j,icg)+
      &      geel_loc_ji*wel_loc
+     &    *fac_shield(i)*fac_shield(j)
+
          geel_loc_ji=
      &     +a22*gmuji2(1)
      &     +a23*gmuji2(2)
@@ -4209,6 +4308,7 @@ 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
+     &    *fac_shield(i)*fac_shield(j)
 #endif
 cd          write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
 
@@ -4222,15 +4322,19 @@ c     &     i,j,a22,muij(1),a23,muij(2),a32,muij(3),a33,muij(4)
 C Partial derivatives in virtual-bond dihedral angles gamma
           if (i.gt.1)
      &    gel_loc_loc(i-1)=gel_loc_loc(i-1)+ 
-     &            a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
-     &           +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j)
+     &            (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
+     &           +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j))
+     &    *fac_shield(i)*fac_shield(j)
+
           gel_loc_loc(j-1)=gel_loc_loc(j-1)+ 
-     &            a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
-     &           +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)
+     &           (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
+     &           +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j))
+     &    *fac_shield(i)*fac_shield(j)
 C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
           do l=1,3
-            ggg(l)=agg(l,1)*muij(1)+
-     &          agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)
+            ggg(l)=(agg(l,1)*muij(1)+
+     &          agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4))
+     &    *fac_shield(i)*fac_shield(j)
             gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
             gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
 cgrad            ghalf=0.5d0*ggg(l)
@@ -4246,12 +4350,20 @@ C Remaining derivatives of eello
           do l=1,3
             gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+
      &        aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4))
+     &    *fac_shield(i)*fac_shield(j)
+
             gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+
      &     aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4))
+     &    *fac_shield(i)*fac_shield(j)
+
             gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+
      &       aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4))
+     &    *fac_shield(i)*fac_shield(j)
+
             gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+
      &     aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4))
+     &    *fac_shield(i)*fac_shield(j)
+
           enddo
           ENDIF
 C Change 12/26/95 to calculate four-body contributions to H-bonding energy
@@ -4330,8 +4442,18 @@ c                ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
                   ees0mij=0
                 endif
 c               ees0mij=0.0D0
+                if (shield_mode.eq.0) then
+                fac_shield(i)=1.0d0
+                fac_shield(j)=1.0d0
+                else
+                ees0plist(num_conti,i)=j
+C                fac_shield(i)=0.4d0
+C                fac_shield(j)=0.6d0
+                endif
                 ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
+     &          *fac_shield(i)*fac_shield(j) 
                 ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
+     &          *fac_shield(i)*fac_shield(j)
 C Diagnostics. Comment out or remove after debugging!
 c               ees0p(num_conti,i)=0.5D0*fac3*ees0pij
 c               ees0m(num_conti,i)=0.5D0*fac3*ees0mij
@@ -4399,17 +4521,29 @@ cgrad                  ghalfm=0.5D0*gggm(k)
                   gacontp_hb1(k,num_conti,i)=!ghalfp
      &              +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
      &              + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
+     &          *fac_shield(i)*fac_shield(j)
+
                   gacontp_hb2(k,num_conti,i)=!ghalfp
      &              +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
      &              + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
+     &          *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)
+     &          *fac_shield(i)*fac_shield(j)
+
                   gacontm_hb2(k,num_conti,i)=!ghalfm
      &              +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
      &              + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
+     &          *fac_shield(i)*fac_shield(j)
+
                   gacontm_hb3(k,num_conti,i)=gggm(k)
+     &          *fac_shield(i)*fac_shield(j)
+
                 enddo
 C Diagnostics. Comment out or remove after debugging!
 cdiag           do k=1,3
@@ -4461,6 +4595,7 @@ C Third- and fourth-order contributions from turns
       include 'COMMON.VECTORS'
       include 'COMMON.FFIELD'
       include 'COMMON.CONTROL'
+      include 'COMMON.SHIELD'
       dimension ggg(3)
       double precision auxmat(2,2),auxmat1(2,2),auxmat2(2,2),pizda(2,2),
      &  e1t(2,2),e2t(2,2),e3t(2,2),e1tder(2,2),e2tder(2,2),e3tder(2,2),
@@ -4501,15 +4636,71 @@ c auxalary matrix for i+2 and constant i+1
         call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1))
         call matmat2(a_temp(1,1),auxgmatt1(1,1),gpizda1(1,1))
         call matmat2(a_temp(1,1),auxgmatt2(1,1),gpizda2(1,1))
+        if (shield_mode.eq.0) then
+        fac_shield(i)=1.0
+        fac_shield(j)=1.0
+C        else
+C        fac_shield(i)=0.4
+C        fac_shield(j)=0.6
+        endif
         eello_turn3=eello_turn3+0.5d0*(pizda(1,1)+pizda(2,2))
+     &  *fac_shield(i)*fac_shield(j)
+        eello_t3=0.5d0*(pizda(1,1)+pizda(2,2))
+     &  *fac_shield(i)*fac_shield(j)
 C Derivatives in theta
         gloc(nphi+i,icg)=gloc(nphi+i,icg)
      &  +0.5d0*(gpizda1(1,1)+gpizda1(2,2))*wturn3
+     &   *fac_shield(i)*fac_shield(j)
         gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)
      &  +0.5d0*(gpizda2(1,1)+gpizda2(2,2))*wturn3
+     &   *fac_shield(i)*fac_shield(j)
 
-        if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
-     &          'eturn3',i,j,0.5d0*(pizda(1,1)+pizda(2,2))
+
+C        if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
+C Derivatives in shield mode
+          if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
+     &  (shield_mode.gt.0)) then
+C          print *,i,j     
+
+          do ilist=1,ishield_list(i)
+           iresshield=shield_list(ilist,i)
+           do k=1,3
+           rlocshield=grad_shield_side(k,ilist,i)*eello_t3/fac_shield(i)
+C     &      *2.0
+           gshieldx_t3(k,iresshield)=gshieldx_t3(k,iresshield)+
+     &              rlocshield
+     & +grad_shield_loc(k,ilist,i)*eello_t3/fac_shield(i)
+            gshieldc_t3(k,iresshield-1)=gshieldc_t3(k,iresshield-1)
+     &      +rlocshield
+           enddo
+          enddo
+          do ilist=1,ishield_list(j)
+           iresshield=shield_list(ilist,j)
+           do k=1,3
+           rlocshield=grad_shield_side(k,ilist,j)*eello_t3/fac_shield(j)
+C     &     *2.0
+           gshieldx_t3(k,iresshield)=gshieldx_t3(k,iresshield)+
+     &              rlocshield
+     & +grad_shield_loc(k,ilist,j)*eello_t3/fac_shield(j)
+           gshieldc_t3(k,iresshield-1)=gshieldc_t3(k,iresshield-1)
+     &             +rlocshield
+
+           enddo
+          enddo
+
+          do k=1,3
+            gshieldc_t3(k,i)=gshieldc_t3(k,i)+
+     &              grad_shield(k,i)*eello_t3/fac_shield(i)
+            gshieldc_t3(k,j)=gshieldc_t3(k,j)+
+     &              grad_shield(k,j)*eello_t3/fac_shield(j)
+            gshieldc_t3(k,i-1)=gshieldc_t3(k,i-1)+
+     &              grad_shield(k,i)*eello_t3/fac_shield(i)
+            gshieldc_t3(k,j-1)=gshieldc_t3(k,j-1)+
+     &              grad_shield(k,j)*eello_t3/fac_shield(j)
+           enddo
+           endif
+
+C        if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
 cd        write (2,*) 'i,',i,' j',j,'eello_turn3',
 cd     &    0.5d0*(pizda(1,1)+pizda(2,2)),
 cd     &    ' eello_turn3_num',4*eello_turn3_num
@@ -4518,12 +4709,14 @@ C Derivatives in gamma(i)
         call transpose2(auxmat2(1,1),auxmat3(1,1))
         call matmat2(a_temp(1,1),auxmat3(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),auxmat3(1,1))
         call matmat2(a_temp(1,1),auxmat3(1,1),pizda(1,1))
         gel_loc_turn3(i+1)=gel_loc_turn3(i+1)
      &    +0.5d0*(pizda(1,1)+pizda(2,2))
+     &   *fac_shield(i)*fac_shield(j)
 C Cartesian derivatives
         do l=1,3
 c            ghalf1=0.5d0*agg(l,1)
@@ -4537,6 +4730,8 @@ c            ghalf4=0.5d0*agg(l,4)
           call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1))
           gcorr3_turn(l,i)=gcorr3_turn(l,i)
      &      +0.5d0*(pizda(1,1)+pizda(2,2))
+     &   *fac_shield(i)*fac_shield(j)
+
           a_temp(1,1)=aggi1(l,1)!+agg(l,1)
           a_temp(1,2)=aggi1(l,2)!+agg(l,2)
           a_temp(2,1)=aggi1(l,3)!+agg(l,3)
@@ -4544,6 +4739,7 @@ c            ghalf4=0.5d0*agg(l,4)
           call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1))
           gcorr3_turn(l,i+1)=gcorr3_turn(l,i+1)
      &      +0.5d0*(pizda(1,1)+pizda(2,2))
+     &   *fac_shield(i)*fac_shield(j)
           a_temp(1,1)=aggj(l,1)!+ghalf1
           a_temp(1,2)=aggj(l,2)!+ghalf2
           a_temp(2,1)=aggj(l,3)!+ghalf3
@@ -4551,6 +4747,7 @@ c            ghalf4=0.5d0*agg(l,4)
           call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1))
           gcorr3_turn(l,j)=gcorr3_turn(l,j)
      &      +0.5d0*(pizda(1,1)+pizda(2,2))
+     &   *fac_shield(i)*fac_shield(j)
           a_temp(1,1)=aggj1(l,1)
           a_temp(1,2)=aggj1(l,2)
           a_temp(2,1)=aggj1(l,3)
@@ -4558,6 +4755,7 @@ c            ghalf4=0.5d0*agg(l,4)
           call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1))
           gcorr3_turn(l,j1)=gcorr3_turn(l,j1)
      &      +0.5d0*(pizda(1,1)+pizda(2,2))
+     &   *fac_shield(i)*fac_shield(j)
         enddo
       return
       end
@@ -4578,6 +4776,7 @@ C Third- and fourth-order contributions from turns
       include 'COMMON.VECTORS'
       include 'COMMON.FFIELD'
       include 'COMMON.CONTROL'
+      include 'COMMON.SHIELD'
       dimension ggg(3)
       double precision auxmat(2,2),auxmat1(2,2),auxmat2(2,2),pizda(2,2),
      &  e1t(2,2),e2t(2,2),e3t(2,2),e1tder(2,2),e2tder(2,2),e3tder(2,2),
@@ -4678,20 +4877,82 @@ c i+3
         gsEE1=0.5d0*(gtEpizda1(1,1)+gtEpizda1(2,2))
         gsEE2=0.5d0*(gtEpizda2(1,1)+gtEpizda2(2,2))
         gsEE3=0.5d0*(gtEpizda3(1,1)+gtEpizda3(2,2))
-
+        if (shield_mode.eq.0) then
+        fac_shield(i)=1.0
+        fac_shield(j)=1.0
+C        else
+C        fac_shield(i)=0.6
+C        fac_shield(j)=0.4
+        endif
         eello_turn4=eello_turn4-(s1+s2+s3)
+     &  *fac_shield(i)*fac_shield(j)
+        eello_t4=-(s1+s2+s3)
+     &  *fac_shield(i)*fac_shield(j)
 c             write(iout,*)'chujOWO', auxvec(1),b1(1,iti2)
         if (energy_dec) write (iout,'(a6,2i5,0pf7.3,3f7.3)')
      &      'eturn4',i,j,-(s1+s2+s3),s1,s2,s3
+C Now derivative over shield:
+          if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
+     &  (shield_mode.gt.0)) then
+C          print *,i,j     
+
+          do ilist=1,ishield_list(i)
+           iresshield=shield_list(ilist,i)
+           do k=1,3
+           rlocshield=grad_shield_side(k,ilist,i)*eello_t4/fac_shield(i)
+C     &      *2.0
+           gshieldx_t4(k,iresshield)=gshieldx_t4(k,iresshield)+
+     &              rlocshield
+     & +grad_shield_loc(k,ilist,i)*eello_t4/fac_shield(i)
+            gshieldc_t4(k,iresshield-1)=gshieldc_t4(k,iresshield-1)
+     &      +rlocshield
+           enddo
+          enddo
+          do ilist=1,ishield_list(j)
+           iresshield=shield_list(ilist,j)
+           do k=1,3
+           rlocshield=grad_shield_side(k,ilist,j)*eello_t4/fac_shield(j)
+C     &     *2.0
+           gshieldx_t4(k,iresshield)=gshieldx_t4(k,iresshield)+
+     &              rlocshield
+     & +grad_shield_loc(k,ilist,j)*eello_t4/fac_shield(j)
+           gshieldc_t4(k,iresshield-1)=gshieldc_t4(k,iresshield-1)
+     &             +rlocshield
+
+           enddo
+          enddo
+
+          do k=1,3
+            gshieldc_t4(k,i)=gshieldc_t4(k,i)+
+     &              grad_shield(k,i)*eello_t4/fac_shield(i)
+            gshieldc_t4(k,j)=gshieldc_t4(k,j)+
+     &              grad_shield(k,j)*eello_t4/fac_shield(j)
+            gshieldc_t4(k,i-1)=gshieldc_t4(k,i-1)+
+     &              grad_shield(k,i)*eello_t4/fac_shield(i)
+            gshieldc_t4(k,j-1)=gshieldc_t4(k,j-1)+
+     &              grad_shield(k,j)*eello_t4/fac_shield(j)
+           enddo
+           endif
+
+
+
+
+
+
 cd        write (2,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3),
 cd     &    ' eello_turn4_num',8*eello_turn4_num
 #ifdef NEWCORR
         gloc(nphi+i,icg)=gloc(nphi+i,icg)
      &                  -(gs13+gsE13+gsEE1)*wturn4
+     &  *fac_shield(i)*fac_shield(j)
         gloc(nphi+i+1,icg)= gloc(nphi+i+1,icg)
      &                    -(gs23+gs21+gsEE2)*wturn4
+     &  *fac_shield(i)*fac_shield(j)
+
         gloc(nphi+i+2,icg)= gloc(nphi+i+2,icg)
      &                    -(gs32+gsE31+gsEE3)*wturn4
+     &  *fac_shield(i)*fac_shield(j)
+
 c         gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)-
 c     &   gs2
 #endif
@@ -4707,6 +4968,7 @@ C Derivatives in gamma(i)
         call matmat2(ae3e2(1,1),e1tder(1,1),pizda(1,1))
         s3=0.5d0*(pizda(1,1)+pizda(2,2))
         gel_loc_turn4(i)=gel_loc_turn4(i)-(s1+s3)
+     &  *fac_shield(i)*fac_shield(j)
 C Derivatives in gamma(i+1)
         call transpose2(EUgder(1,1,i+2),e2tder(1,1))
         call matvec2(ae3(1,1),Ub2der(1,i+2),auxvec(1)) 
@@ -4715,6 +4977,7 @@ C Derivatives in gamma(i+1)
         call matmat2(auxmat(1,1),e1t(1,1),pizda(1,1))
         s3=0.5d0*(pizda(1,1)+pizda(2,2))
         gel_loc_turn4(i+1)=gel_loc_turn4(i+1)-(s2+s3)
+     &  *fac_shield(i)*fac_shield(j)
 C Derivatives in gamma(i+2)
         call transpose2(EUgder(1,1,i+3),e3tder(1,1))
         call matvec2(e1a(1,1),Ub2der(1,i+3),auxvec(1))
@@ -4726,6 +4989,7 @@ C Derivatives in gamma(i+2)
         call matmat2(auxmat3(1,1),e1t(1,1),pizda(1,1))
         s3=0.5d0*(pizda(1,1)+pizda(2,2))
         gel_loc_turn4(i+2)=gel_loc_turn4(i+2)-(s1+s2+s3)
+     &  *fac_shield(i)*fac_shield(j)
 C Cartesian derivatives
 C Derivatives of this turn contributions in DC(i+2)
         if (j.lt.nres-1) then
@@ -4745,6 +5009,7 @@ C Derivatives of this turn contributions in DC(i+2)
             s3=0.5d0*(pizda(1,1)+pizda(2,2))
             ggg(l)=-(s1+s2+s3)
             gcorr4_turn(l,i+2)=gcorr4_turn(l,i+2)-(s1+s2+s3)
+     &  *fac_shield(i)*fac_shield(j)
           enddo
         endif
 C Remaining derivatives of this turn contribution
@@ -4763,6 +5028,7 @@ C Remaining derivatives of this turn contribution
           call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
           s3=0.5d0*(pizda(1,1)+pizda(2,2))
           gcorr4_turn(l,i)=gcorr4_turn(l,i)-(s1+s2+s3)
+     &  *fac_shield(i)*fac_shield(j)
           a_temp(1,1)=aggi1(l,1)
           a_temp(1,2)=aggi1(l,2)
           a_temp(2,1)=aggi1(l,3)
@@ -4777,6 +5043,7 @@ C Remaining derivatives of this turn contribution
           call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
           s3=0.5d0*(pizda(1,1)+pizda(2,2))
           gcorr4_turn(l,i+1)=gcorr4_turn(l,i+1)-(s1+s2+s3)
+     &  *fac_shield(i)*fac_shield(j)
           a_temp(1,1)=aggj(l,1)
           a_temp(1,2)=aggj(l,2)
           a_temp(2,1)=aggj(l,3)
@@ -4791,6 +5058,7 @@ C Remaining derivatives of this turn contribution
           call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
           s3=0.5d0*(pizda(1,1)+pizda(2,2))
           gcorr4_turn(l,j)=gcorr4_turn(l,j)-(s1+s2+s3)
+     &  *fac_shield(i)*fac_shield(j)
           a_temp(1,1)=aggj1(l,1)
           a_temp(1,2)=aggj1(l,2)
           a_temp(2,1)=aggj1(l,3)
@@ -4806,6 +5074,7 @@ C Remaining derivatives of this turn contribution
           s3=0.5d0*(pizda(1,1)+pizda(2,2))
 c          write (iout,*) "s1",s1," s2",s2," s3",s3," s1+s2+s3",s1+s2+s3
           gcorr4_turn(l,j1)=gcorr4_turn(l,j1)-(s1+s2+s3)
+     &  *fac_shield(i)*fac_shield(j)
         enddo
       return
       end
@@ -8147,9 +8416,11 @@ c------------------------------------------------------------------------------
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
       include 'COMMON.SHIELD'
+      include 'COMMON.CONTROL'
       double precision gx(3),gx1(3)
       logical lprn
       lprn=.false.
+C      print *,"wchodze",fac_shield(i),shield_mode
       eij=facont_hb(jj,i)
       ekl=facont_hb(kk,k)
       ees0pij=ees0p(jj,i)
@@ -8157,8 +8428,9 @@ c------------------------------------------------------------------------------
       ees0mij=ees0m(jj,i)
       ees0mkl=ees0m(kk,k)
       ekont=eij*ekl
-      ees=-(coeffp*ees0pij*ees0pkl+coeffm*ees0mij*ees0mkl)*
-     & fac_shield(i)**2*fac_shield(j)**2
+      ees=-(coeffp*ees0pij*ees0pkl+coeffm*ees0mij*ees0mkl)
+C*
+C     & fac_shield(i)**2*fac_shield(j)**2
 cd    ees=-(coeffp*ees0pkl+coeffm*ees0mkl)
 C Following 4 lines for diagnostics.
 cd    ees0pkl=0.0D0
@@ -8222,7 +8494,89 @@ cgrad     &     coeffm*ees0mij*gacontm_hb3(ll,kk,k))
 cgrad        enddo
 cgrad      enddo 
 c      write (iout,*) "ehbcorr",ekont*ees
+C      print *,ekont,ees,i,k
       ehbcorr=ekont*ees
+C now gradient over shielding
+C      return
+      if (shield_mode.gt.0) then
+       j=ees0plist(jj,i)
+       l=ees0plist(kk,k)
+C        print *,i,j,fac_shield(i),fac_shield(j),
+C     &fac_shield(k),fac_shield(l)
+        if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
+     &      (fac_shield(k).gt.0).and.(fac_shield(l).gt.0)) then
+          do ilist=1,ishield_list(i)
+           iresshield=shield_list(ilist,i)
+           do m=1,3
+           rlocshield=grad_shield_side(m,ilist,i)*ehbcorr/fac_shield(i)
+C     &      *2.0
+           gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+
+     &              rlocshield
+     & +grad_shield_loc(m,ilist,i)*ehbcorr/fac_shield(i)
+            gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1)
+     &+rlocshield
+           enddo
+          enddo
+          do ilist=1,ishield_list(j)
+           iresshield=shield_list(ilist,j)
+           do m=1,3
+           rlocshield=grad_shield_side(m,ilist,j)*ehbcorr/fac_shield(j)
+C     &     *2.0
+           gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+
+     &              rlocshield
+     & +grad_shield_loc(m,ilist,j)*ehbcorr/fac_shield(j)
+           gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1)
+     &     +rlocshield
+           enddo
+          enddo
+
+          do ilist=1,ishield_list(k)
+           iresshield=shield_list(ilist,k)
+           do m=1,3
+           rlocshield=grad_shield_side(m,ilist,k)*ehbcorr/fac_shield(k)
+C     &     *2.0
+           gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+
+     &              rlocshield
+     & +grad_shield_loc(m,ilist,k)*ehbcorr/fac_shield(k)
+           gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1)
+     &     +rlocshield
+           enddo
+          enddo
+          do ilist=1,ishield_list(l)
+           iresshield=shield_list(ilist,l)
+           do m=1,3
+           rlocshield=grad_shield_side(m,ilist,l)*ehbcorr/fac_shield(l)
+C     &     *2.0
+           gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+
+     &              rlocshield
+     & +grad_shield_loc(m,ilist,l)*ehbcorr/fac_shield(l)
+           gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1)
+     &     +rlocshield
+           enddo
+          enddo
+C          print *,gshieldx(m,iresshield)
+          do m=1,3
+            gshieldc_ec(m,i)=gshieldc_ec(m,i)+
+     &              grad_shield(m,i)*ehbcorr/fac_shield(i)
+            gshieldc_ec(m,j)=gshieldc_ec(m,j)+
+     &              grad_shield(m,j)*ehbcorr/fac_shield(j)
+            gshieldc_ec(m,i-1)=gshieldc_ec(m,i-1)+
+     &              grad_shield(m,i)*ehbcorr/fac_shield(i)
+            gshieldc_ec(m,j-1)=gshieldc_ec(m,j-1)+
+     &              grad_shield(m,j)*ehbcorr/fac_shield(j)
+
+            gshieldc_ec(m,k)=gshieldc_ec(m,k)+
+     &              grad_shield(m,k)*ehbcorr/fac_shield(k)
+            gshieldc_ec(m,l)=gshieldc_ec(m,l)+
+     &              grad_shield(m,l)*ehbcorr/fac_shield(l)
+            gshieldc_ec(m,k-1)=gshieldc_ec(m,k-1)+
+     &              grad_shield(m,k)*ehbcorr/fac_shield(k)
+            gshieldc_ec(m,l-1)=gshieldc_ec(m,l-1)+
+     &              grad_shield(m,l)*ehbcorr/fac_shield(l)
+
+           enddo       
+      endif
+      endif
       return
       end
 #ifdef MOMENT

diff --git a/source/unres/src_MD-M/gradient_p.F b/source/unres/src_MD-M/gradient_p.F
index edbf1bc..4861790 100644 (file)
--- a/source/unres/src_MD-M/gradient_p.F
+++ b/source/unres/src_MD-M/gradient_p.F
@@ -358,9 +358,24 @@ C
           gvdwc_scp(j,i)=0.0D0
           gvdwc_scpp(j,i)=0.0d0
          gelc (j,i)=0.0D0
+C below is zero grad for shielding in order: ees (p-p)
+C ecorr4, eturn3, eturn4, eel_loc, c denotes calfa,x is side-chain
           gshieldx(j,i)=0.0d0
           gshieldc(j,i)=0.0d0
           gshieldc_loc(j,i)=0.0d0
+          gshieldx_ec(j,i)=0.0d0
+          gshieldc_ec(j,i)=0.0d0
+          gshieldc_loc_ec(j,i)=0.0d0
+          gshieldx_t3(j,i)=0.0d0
+          gshieldc_t3(j,i)=0.0d0
+          gshieldc_loc_t3(j,i)=0.0d0
+          gshieldx_t4(j,i)=0.0d0
+          gshieldc_t4(j,i)=0.0d0
+          gshieldc_loc_t4(j,i)=0.0d0
+          gshieldx_ll(j,i)=0.0d0
+          gshieldc_ll(j,i)=0.0d0
+          gshieldc_loc_ll(j,i)=0.0d0
+C end of zero grad for shielding
          gelc_long(j,i)=0.0D0
           gradb(j,i)=0.0d0
           gradbx(j,i)=0.0d0