projects / unres.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
respa cleaning (my bug) and shielding MPI and previous bug cleaning
[unres.git] / source / unres / src_MD-M / energy_p_new-sep_barrier.F
diff --git a/source/unres/src_MD-M/energy_p_new-sep_barrier.F b/source/unres/src_MD-M/energy_p_new-sep_barrier.F
index 7d2d27f..3815c13 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new-sep_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new-sep_barrier.F
@@ -1,4 +1,33 @@
 C-----------------------------------------------------------------------
+      double precision function sscalelip(r)
+      double precision r,gamm
+      include "COMMON.SPLITELE"
+C      if(r.lt.r_cut-rlamb) then
+C        sscale=1.0d0
+C      else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
+C        gamm=(r-(r_cut-rlamb))/rlamb
+        sscalelip=1.0d0+r*r*(2*r-3.0d0)
+C      else
+C        sscale=0d0
+C      endif
+      return
+      end
+C-----------------------------------------------------------------------
+      double precision function sscagradlip(r)
+      double precision r,gamm
+      include "COMMON.SPLITELE"
+C     if(r.lt.r_cut-rlamb) then
+C        sscagrad=0.0d0
+C      else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
+C        gamm=(r-(r_cut-rlamb))/rlamb
+        sscagradlip=r*(6*r-6.0d0)
+C      else
+C        sscagrad=0.0d0
+C      endif
+      return
+      end
+
+C-----------------------------------------------------------------------
       double precision function sscale(r)
       double precision r,gamm
       include "COMMON.SPLITELE"
@@ -75,8 +104,8 @@ cd   &                  'iend=',iend(i,iint)
               rrij=1.0D0/rij
               eps0ij=eps(itypi,itypj)
               fac=rrij**expon2
-              e1=fac*fac*aa(itypi,itypj)
-              e2=fac*bb(itypi,itypj)
+              e1=fac*fac*aa
+              e2=fac*bb
               evdwij=e1+e2
               evdw=evdw+(1.0d0-sss)*evdwij
 C 
@@ -164,8 +193,8 @@ C Change 12/1/95 to calculate four-body interactions
               rrij=1.0D0/rij
               eps0ij=eps(itypi,itypj)
               fac=rrij**expon2
-              e1=fac*fac*aa(itypi,itypj)
-              e2=fac*bb(itypi,itypj)
+              e1=fac*fac*aa
+              e2=fac*bb
               evdwij=e1+e2
               evdw=evdw+sss*evdwij
 C 
@@ -248,8 +277,8 @@ C
             if (sss.lt.1.0d0) then
               r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
               fac=r_shift_inv**expon
-              e1=fac*fac*aa(itypi,itypj)
-              e2=fac*bb(itypi,itypj)
+              e1=fac*fac*aa
+              e2=fac*bb
               evdwij=e_augm+e1+e2
 cd            sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
 cd            epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
@@ -331,8 +360,8 @@ C
             if (sss.gt.0.0d0) then
               r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
               fac=r_shift_inv**expon
-              e1=fac*fac*aa(itypi,itypj)
-              e2=fac*bb(itypi,itypj)
+              e1=fac*fac*aa
+              e2=fac*bb
               evdwij=e_augm+e1+e2
 cd            sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
 cd            epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
@@ -445,16 +474,16 @@ C Calculate the angle-dependent terms of energy & contributions to derivatives.
 C Calculate whole angle-dependent part of epsilon and contributions
 C to its derivatives
               fac=(rrij*sigsq)**expon2
-              e1=fac*fac*aa(itypi,itypj)
-              e2=fac*bb(itypi,itypj)
+              e1=fac*fac*aa
+              e2=fac*bb
               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
               eps2der=evdwij*eps3rt
               eps3der=evdwij*eps2rt
               evdwij=evdwij*eps2rt*eps3rt
               evdw=evdw+evdwij*(1.0d0-sss)
               if (lprn) then
-              sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
-              epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+              sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+              epsi=bb**2/aa
 cd              write (iout,'(2(a3,i3,2x),15(0pf7.3))')
 cd     &          restyp(itypi),i,restyp(itypj),j,
 cd     &          epsi,sigm,chi1,chi2,chip1,chip2,
@@ -558,16 +587,16 @@ C Calculate the angle-dependent terms of energy & contributions to derivatives.
 C Calculate whole angle-dependent part of epsilon and contributions
 C to its derivatives
               fac=(rrij*sigsq)**expon2
-              e1=fac*fac*aa(itypi,itypj)
-              e2=fac*bb(itypi,itypj)
+              e1=fac*fac*aa
+              e2=fac*bb
               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
               eps2der=evdwij*eps3rt
               eps3der=evdwij*eps2rt
               evdwij=evdwij*eps2rt*eps3rt
               evdw=evdw+evdwij*sss
               if (lprn) then
-              sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
-              epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+              sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+              epsi=bb**2/aa
 cd              write (iout,'(2(a3,i3,2x),15(0pf7.3))')
 cd     &          restyp(itypi),i,restyp(itypj),j,
 cd     &          epsi,sigm,chi1,chi2,chip1,chip2,
@@ -672,6 +701,33 @@ c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
           if (yj.lt.0) yj=yj+boxysize
           zj=mod(zj,boxzsize)
           if (zj.lt.0) zj=zj+boxzsize
+       if ((zj.gt.bordlipbot)
+     &.and.(zj.lt.bordliptop)) then
+C the energy transfer exist
+        if (zj.lt.buflipbot) then
+C what fraction I am in
+         fracinbuf=1.0d0-
+     &        ((positi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+         sslipj=sscalelip(fracinbuf)
+         ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
+        elseif (zi.gt.bufliptop) then
+         fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
+         sslipj=sscalelip(fracinbuf)
+         ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
+        else
+         sslipj=1.0d0
+         ssgradlipj=0.0
+        endif
+       else
+         sslipj=0.0d0
+         ssgradlipj=0.0
+       endif
+      aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+     &  +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
+
       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
       xj_safe=xj
       yj_safe=yj
@@ -732,8 +788,8 @@ cd     &          rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
 c---------------------------------------------------------------
               rij_shift=1.0D0/rij_shift 
               fac=rij_shift**expon
-              e1=fac*fac*aa(itypi,itypj)
-              e2=fac*bb(itypi,itypj)
+              e1=fac*fac*aa
+              e2=fac*bb
               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
               eps2der=evdwij*eps3rt
               eps3der=evdwij*eps2rt
@@ -742,8 +798,8 @@ c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
               evdwij=evdwij*eps2rt*eps3rt
               evdw=evdw+evdwij*(1.0d0-sss)
               if (lprn) then
-              sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
-              epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+              sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+              epsi=bb**2/aa
               write (iout,'(2(a3,i3,2x),17(0pf7.3))')
      &          restyp(itypi),i,restyp(itypj),j,
      &          epsi,sigm,chi1,chi2,chip1,chip2,
@@ -766,6 +822,8 @@ C Calculate the radial part of the gradient
               gg(1)=xj*fac
               gg(2)=yj*fac
               gg(3)=zj*fac
+              gg_lipi(3)=ssgradlipi*evdwij
+              gg_lipj(3)=ssgradlipj*evdwij
 C Calculate angular part of the gradient.
               call sc_grad_scale(1.0d0-sss)
             endif
@@ -854,6 +912,32 @@ c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
           if (yj.lt.0) yj=yj+boxysize
           zj=mod(zj,boxzsize)
           if (zj.lt.0) zj=zj+boxzsize
+       if ((zj.gt.bordlipbot)
+     &.and.(zj.lt.bordliptop)) then
+C the energy transfer exist
+        if (zj.lt.buflipbot) then
+C what fraction I am in
+         fracinbuf=1.0d0-
+     &        ((positi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+         sslipj=sscalelip(fracinbuf)
+         ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
+        elseif (zi.gt.bufliptop) then
+         fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
+         sslipj=sscalelip(fracinbuf)
+         ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
+        else
+         sslipj=1.0d0
+         ssgradlipj=0.0
+        endif
+       else
+         sslipj=0.0d0
+         ssgradlipj=0.0
+       endif
+      aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+     &  +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
       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
       xj_safe=xj
       yj_safe=yj
@@ -914,8 +998,8 @@ cd     &          rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
 c---------------------------------------------------------------
               rij_shift=1.0D0/rij_shift 
               fac=rij_shift**expon
-              e1=fac*fac*aa(itypi,itypj)
-              e2=fac*bb(itypi,itypj)
+              e1=fac*fac*aa
+              e2=fac*bb
               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
               eps2der=evdwij*eps3rt
               eps3der=evdwij*eps2rt
@@ -924,8 +1008,8 @@ c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
               evdwij=evdwij*eps2rt*eps3rt
               evdw=evdw+evdwij*sss
               if (lprn) then
-              sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
-              epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+              sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+              epsi=bb**2/aa
               write (iout,'(2(a3,i3,2x),17(0pf7.3))')
      &          restyp(itypi),i,restyp(itypj),j,
      &          epsi,sigm,chi1,chi2,chip1,chip2,
@@ -948,6 +1032,8 @@ C Calculate the radial part of the gradient
               gg(1)=xj*fac
               gg(2)=yj*fac
               gg(3)=zj*fac
+              gg_lipi(3)=ssgradlipi*evdwij
+              gg_lipj(3)=ssgradlipj*evdwij
 C Calculate angular part of the gradient.
               call sc_grad_scale(sss)
             endif
@@ -1044,8 +1130,8 @@ C I hate to put IF's in the loops, but here don't have another choice!!!!
 c---------------------------------------------------------------
               rij_shift=1.0D0/rij_shift 
               fac=rij_shift**expon
-              e1=fac*fac*aa(itypi,itypj)
-              e2=fac*bb(itypi,itypj)
+              e1=fac*fac*aa
+              e2=fac*bb
               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
               eps2der=evdwij*eps3rt
               eps3der=evdwij*eps2rt
@@ -1054,8 +1140,8 @@ c---------------------------------------------------------------
               evdwij=evdwij*eps2rt*eps3rt
               evdw=evdw+(evdwij+e_augm)*(1.0d0-sss)
               if (lprn) then
-              sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
-              epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+              sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+              epsi=bb**2/aa
               write (iout,'(2(a3,i3,2x),17(0pf7.3))')
      &          restyp(itypi),i,restyp(itypj),j,
      &          epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
@@ -1166,8 +1252,8 @@ C I hate to put IF's in the loops, but here don't have another choice!!!!
 c---------------------------------------------------------------
               rij_shift=1.0D0/rij_shift 
               fac=rij_shift**expon
-              e1=fac*fac*aa(itypi,itypj)
-              e2=fac*bb(itypi,itypj)
+              e1=fac*fac*aa
+              e2=fac*bb
               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
               eps2der=evdwij*eps3rt
               eps3der=evdwij*eps2rt
@@ -1176,8 +1262,8 @@ c---------------------------------------------------------------
               evdwij=evdwij*eps2rt*eps3rt
               evdw=evdw+(evdwij+e_augm)*sss
               if (lprn) then
-              sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
-              epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+              sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+              epsi=bb**2/aa
               write (iout,'(2(a3,i3,2x),17(0pf7.3))')
      &          restyp(itypi),i,restyp(itypj),j,
      &          epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
@@ -1234,10 +1320,10 @@ c      write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
       enddo 
 c      write (iout,*) "gg",(gg(k),k=1,3)
       do k=1,3
-        gvdwx(k,i)=gvdwx(k,i)-gg(k)
+        gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
      &            +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
      &          +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*scalfac
-        gvdwx(k,j)=gvdwx(k,j)+gg(k)
+        gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
      &            +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
      &          +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac
 c        write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
@@ -1249,8 +1335,8 @@ C
 C Calculate the components of the gradient in DC and X
 C
       do l=1,3
-        gvdwc(l,i)=gvdwc(l,i)-gg(l)
-        gvdwc(l,j)=gvdwc(l,j)+gg(l)
+        gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(l)
+        gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l)
       enddo
       return
       end
@@ -1282,6 +1368,7 @@ C
       include 'COMMON.VECTORS'
       include 'COMMON.FFIELD'
       include 'COMMON.TIME1'
+      include 'COMMON.SHIELD'
       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),
@@ -1368,8 +1455,8 @@ C
       do i=iturn3_start,iturn3_end
         if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
      &  .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1
-     &  .or. itype(i-1).eq.ntyp1
-     &  .or. itype(i+4).eq.ntyp1
+C     &  .or. itype(i-1).eq.ntyp1
+C     &  .or. itype(i+4).eq.ntyp1
      &   ) cycle
         dxi=dc(1,i)
         dyi=dc(2,i)
@@ -1395,8 +1482,8 @@ C
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
      &    .or. itype(i+3).eq.ntyp1
      &    .or. itype(i+4).eq.ntyp1
-     &    .or. itype(i+5).eq.ntyp1
-     &    .or. itype(i-1).eq.ntyp1
+C     &    .or. itype(i+5).eq.ntyp1
+C     &    .or. itype(i-1).eq.ntyp1
      &    ) cycle
         dxi=dc(1,i)
         dyi=dc(2,i)
@@ -1424,8 +1511,8 @@ c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
 c
       do i=iatel_s,iatel_e
         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
+C     &  .or. itype(i-1).eq.ntyp1
      &) cycle
         dxi=dc(1,i)
         dyi=dc(2,i)
@@ -1446,8 +1533,8 @@ c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
         num_conti=num_cont_hb(i)
         do j=ielstart(i),ielend(i)
           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
           call eelecij_scale(i,j,ees,evdw1,eel_loc)
         enddo ! j
@@ -1483,6 +1570,7 @@ C-------------------------------------------------------------------------------
       include 'COMMON.VECTORS'
       include 'COMMON.FFIELD'
       include 'COMMON.TIME1'
+      include 'COMMON.SHIELD'
       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),
@@ -1579,8 +1667,14 @@ c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
           fac3=ael6i*r6ij
           fac4=ael3i*r3ij
           evdwij=ev1+ev2
+          if (shield_mode.eq.0) then
+          fac_shield(i)=1.0
+          fac_shield(j)=1.0
+          endif
           el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
           el2=fac4*fac       
+          el1=el1*fac_shield(i)**2*fac_shield(j)**2
+          el2=el2*fac_shield(i)**2*fac_shield(j)**2
           eesij=el1+el2
 C 12/26/95 - for the evaluation of multi-body H-bonding interactions
           ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
@@ -1612,6 +1706,60 @@ C
           ggg(1)=facel*xj
           ggg(2)=facel*yj
           ggg(3)=facel*zj
+          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)*eesij/fac_shield(i)
+     &      *2.0
+           gshieldx(k,iresshield)=gshieldx(k,iresshield)+
+     &              rlocshield
+     & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)*2.0
+            gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
+C           gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield)
+C     & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
+C             if (iresshield.gt.i) then
+C               do ishi=i+1,iresshield-1
+C                gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield
+C     & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
+C
+C              enddo
+C             else
+C               do ishi=iresshield,i
+C                gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield
+C     & -grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
+C
+C               enddo
+C              endif
+           enddo
+          enddo
+          do ilist=1,ishield_list(j)
+           iresshield=shield_list(ilist,j)
+           do k=1,3
+           rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j)
+     &     *2.0
+           gshieldx(k,iresshield)=gshieldx(k,iresshield)+
+     &              rlocshield
+     & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0
+           gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
+           enddo
+          enddo
+
+          do k=1,3
+            gshieldc(k,i)=gshieldc(k,i)+
+     &              grad_shield(k,i)*eesij/fac_shield(i)*2.0
+            gshieldc(k,j)=gshieldc(k,j)+
+     &              grad_shield(k,j)*eesij/fac_shield(j)*2.0
+            gshieldc(k,i-1)=gshieldc(k,i-1)+
+     &              grad_shield(k,i)*eesij/fac_shield(i)*2.0
+            gshieldc(k,j-1)=gshieldc(k,j-1)+
+     &              grad_shield(k,j)*eesij/fac_shield(j)*2.0
+
+           enddo
+           endif
+
 c          do k=1,3
 c            ghalf=0.5D0*ggg(k)
 c            gelc(k,i)=gelc(k,i)+ghalf
@@ -1710,7 +1858,9 @@ C          ggg(3)=facvdw*zj
 cd        print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
 cd   &          (dcosg(k),k=1,3)
           do k=1,3
-            ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k) 
+            ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*
+     &      fac_shield(i)**2*fac_shield(j)**2
+
           enddo
 c          do k=1,3
 c            ghalf=0.5D0*ggg(k)
@@ -1728,11 +1878,14 @@ cgrad            enddo
 cgrad          enddo
           do k=1,3
             gelc(k,i)=gelc(k,i)
-     &               +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
-     &               + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
+     &               +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
+     &               + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))
+     &           *fac_shield(i)**2*fac_shield(j)**2
+
             gelc(k,j)=gelc(k,j)
-     &               +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
-     &               + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
+     &               +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
+     &               + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))
+     &           *fac_shield(i)**2*fac_shield(j)**2
             gelc_long(k,j)=gelc_long(k,j)+ggg(k)
             gelc_long(k,i)=gelc_long(k,i)-ggg(k)
           enddo
@@ -1928,19 +2081,80 @@ cd          write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
           if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
      &            'eelloc',i,j,eel_loc_ij
 
+
+          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)
           eel_loc=eel_loc+eel_loc_ij
+
+          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 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)
@@ -1954,14 +2168,22 @@ cgrad            enddo
 cgrad          enddo
 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)
-            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)
-            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)
-            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)
+            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
@@ -2040,8 +2262,19 @@ 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
@@ -2109,17 +2342,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
               ENDIF ! wcorr
               endif  ! num_conti.le.maxconts
UNRESPACK 3.x