introduction of shielding to cluster DEBUG mode

diff --git a/source/cluster/wham/src-M/COMMON.CLUSTER b/source/cluster/wham/src-M/COMMON.CLUSTER
index 4477d19..4fedc52 100644 (file)
--- a/source/cluster/wham/src-M/COMMON.CLUSTER
+++ b/source/cluster/wham/src-M/COMMON.CLUSTER
@@ -6,12 +6,12 @@
       integer ncut,ngr,licz,nconf,iass,icc,mult,list_conf,
      &  nss_all,ihpb_all,jhpb_all,iass_tot,iscore,nprop
       common /clu/ diss(maxdist),energy(0:maxconf),
-     &  enetb(max_ene,maxstr_proc),ecut,
+     &  enetb(max_ene,maxconf),ecut,
      &  entfac(maxconf),totfree(0:maxconf),totfree_gr(maxgr),
      &  rcutoff(max_cut+1),ncut,min_var,tree,plot_tree,lgrp
       common /clu1/ ngr,licz(maxgr),nconf(maxgr,maxingr),iass(maxgr),
      &  iass_tot(maxgr,max_cut),list_conf(maxconf)
-      common /alles/ allcart(3,maxres2,maxstr_proc),rmstb(maxconf),
+      common /alles/ allcart(3,maxres2,maxconf),rmstb(maxconf),
      & icc(maxconf),
-     & mult(maxres),nss_all(maxstr_proc),ihpb_all(maxss,maxstr_proc),
-     & jhpb_all(maxss,maxstr_proc),iscore(maxconf),nprop
+     & mult(maxres),nss_all(maxconf),ihpb_all(maxss,maxconf),
+     & jhpb_all(maxss,maxconf),iscore(maxconf),nprop

diff --git a/source/cluster/wham/src-M/COMMON.CONTROL b/source/cluster/wham/src-M/COMMON.CONTROL
index f339ae9..708fe2b 100644 (file)
--- a/source/cluster/wham/src-M/COMMON.CONTROL
+++ b/source/cluster/wham/src-M/COMMON.CONTROL
@@ -1,6 +1,6 @@
       double precision betaT
       integer iscode,indpdb,outpdb,outmol2,iopt,nstart,nend,symetr,
-     & constr_dist
+     & constr_dist,shield_mode,tor_mode
       logical refstr,pdbref,punch_dist,print_dist,caonly,lside,
      & lprint_cart,lprint_int,from_cart,efree,from_bx,from_cx,
      & with_dihed_constr,with_theta_constr
@@ -8,4 +8,5 @@
      & punch_dist,print_dist,caonly,lside,lprint_cart,lprint_int,
      & from_cart,from_bx,from_cx, with_dihed_constr,with_theta_constr,
      & efree,iopt,nstart,nend,symetr,
+     & tor_mode,shield_mode,
      & constr_dist

diff --git a/source/cluster/wham/src-M/COMMON.DERIV b/source/cluster/wham/src-M/COMMON.DERIV
index 79f8630..4e116e1 100644 (file)
--- a/source/cluster/wham/src-M/COMMON.DERIV
+++ b/source/cluster/wham/src-M/COMMON.DERIV
@@ -3,13 +3,31 @@
      & gradcorr5,gradcorr6,gel_loc,gcorr3_turn,gcorr4_turn,gcorr6_turn,
      & gel_loc_loc,gel_loc_turn3,gel_loc_turn4,gel_loc_turn6,gcorr_loc,
      & g_corr5_loc,g_corr6_loc,gradb,gradbx,gsccorc,gsccorx,gsccor_loc,
-     & gscloc,gsclocx
+     & gscloc,gsclocx,gshieldx,gradafm,
+     & gshieldc, gshieldc_loc, gshieldx_ec, gshieldc_ec,
+     & gshieldc_loc_ec, gshieldx_t3,gshieldc_t3,gshieldc_loc_t3,
+     & gshieldx_t4, gshieldc_t4,gshieldc_loc_t4,gshieldx_ll,
+     & gshieldc_ll, gshieldc_loc_ll
+
+
       integer nfl,icg
       logical calc_grad
       common /derivat/ dcdv(6,maxdim),dxdv(6,maxdim),dxds(6,maxres),
      & gradx(3,maxres,2),gradc(3,maxres,2),gvdwx(3,maxres),
      & gvdwc(3,maxres),gelc(3,maxres),gvdwpp(3,maxres),
      & gradx_scp(3,maxres),
+     & gliptranc(3,-1:maxres),
+     & gliptranx(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),
      & gvdwc_scp(3,maxres),ghpbx(3,maxres),ghpbc(3,maxres),
      & gloc(maxvar,2),gradcorr(3,maxres),gradxorr(3,maxres),
      & gradcorr5(3,maxres),gradcorr6(3,maxres),

diff --git a/source/cluster/wham/src-M/DIMENSIONS b/source/cluster/wham/src-M/DIMENSIONS
index 29e8f24..86f775a 100644 (file)
--- a/source/cluster/wham/src-M/DIMENSIONS
+++ b/source/cluster/wham/src-M/DIMENSIONS
@@ -9,7 +9,7 @@ C Max. number of processors.
       parameter (maxprocs=16)
 C Max. number of AA residues
       integer maxres,maxres2
-      parameter (maxres=650)
+      parameter (maxres=100)
 C Appr. max. number of interaction sites
       parameter (maxres2=2*maxres)
 C Max. number of variables
@@ -60,7 +60,7 @@ C Max number of symetric chains
       parameter (maxperm=120)
 C Max. number of energy components
       integer max_ene
-      parameter (max_ene=24)
+      parameter (max_ene=25)
 C Max. number of temperatures
       integer maxt
       parameter (maxT=5)

diff --git a/source/cluster/wham/src-M/energy_p_new.F b/source/cluster/wham/src-M/energy_p_new.F
index f640679..8203cc4 100644 (file)
--- a/source/cluster/wham/src-M/energy_p_new.F
+++ b/source/cluster/wham/src-M/energy_p_new.F
@@ -22,6 +22,8 @@ cMS$ATTRIBUTES C ::  proc_proc
       include 'COMMON.INTERACT'
       include 'COMMON.SBRIDGE'
       include 'COMMON.CHAIN'
+      include 'COMMON.SHIELD'
+      include 'COMMON.CONTROL'
       double precision fact(6)
 cd      write(iout, '(a,i2)')'Calling etotal ipot=',ipot
 cd    print *,'nnt=',nnt,' nct=',nct
@@ -47,7 +49,14 @@ C Gay-Berne-Vorobjev potential (shifted LJ, angular dependence).
 C
 C Calculate electrostatic (H-bonding) energy of the main chain.
 C
-  106 call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
+  106 continue
+      write(iout,*) "shield_mode",shield_mode,ethetacnstr 
+      if (shield_mode.eq.1) then
+       call set_shield_fac
+      else if  (shield_mode.eq.2) then
+       call set_shield_fac2
+      endif
+      call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
 C
 C Calculate excluded-volume interaction energy between peptide groups
 C and side chains.
@@ -102,8 +111,20 @@ c         print *,ecorr,ecorr5,ecorr6,eturn6
       if (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) then
          call multibody_hb(ecorr,ecorr5,ecorr6,n_corr,n_corr1)
       endif
-c      write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t
+      write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t
 #ifdef SPLITELE
+      if (shield_mode.gt.0) then
+      etot=fact(1)*wsc*(evdw+fact(6)*evdw_t)+fact(1)*wscp*evdw2
+     & +welec*fact(1)*ees
+     & +fact(1)*wvdwpp*evdw1
+     & +wang*ebe+wtor*fact(1)*etors+wscloc*escloc
+     & +wstrain*ehpb+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
+     & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4
+     & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
+     & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
+     & +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr
+C     & +wliptran*eliptran
+      else
       etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2+welec*fact(1)*ees
      & +wvdwpp*evdw1
      & +wang*ebe+wtor*fact(1)*etors+wscloc*escloc
@@ -112,7 +133,20 @@ c      write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t
      & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
      & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
      & +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr
+C     & +wliptran*eliptran
+      endif
 #else
+      if (shield_mode.gt.0) then
+      etot=fact(1)wsc*(evdw+fact(6)*evdw_t)+fact(1)*wscp*evdw2
+     & +welec*fact(1)*(ees+evdw1)
+     & +wang*ebe+wtor*fact(1)*etors+wscloc*escloc
+     & +wstrain*ehpb+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
+     & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4
+     & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
+     & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
+     & +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr
+C     & +wliptran*eliptran
+      else
       etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2
      & +welec*fact(1)*(ees+evdw1)
      & +wang*ebe+wtor*fact(1)*etors+wscloc*escloc
@@ -121,7 +155,10 @@ c      write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t
      & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
      & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
      & +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr
+C     & +wliptran*eliptran
+      endif
 #endif
+
       energia(0)=etot
       energia(1)=evdw
 #ifdef SCP14
@@ -181,6 +218,7 @@ C
 #ifdef SPLITELE
       do i=1,nct
         do j=1,3
+      if (shield_mode.eq.0) then
           gradc(j,i,icg)=wsc*gvdwc(j,i)+wscp*gvdwc_scp(j,i)+
      &                welec*fact(1)*gelc(j,i)+wvdwpp*gvdwpp(j,i)+
      &                wbond*gradb(j,i)+
@@ -193,14 +231,57 @@ C
      &                wcorr6*fact(5)*gradcorr6(j,i)+
      &                wturn6*fact(5)*gcorr6_turn(j,i)+
      &                wsccor*fact(2)*gsccorc(j,i)
+     &               +wliptran*gliptranc(j,i)
           gradx(j,i,icg)=wsc*gvdwx(j,i)+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)
+        else
+          gradc(j,i,icg)=fact(1)*wsc*gvdwc(j,i)
+     &                +fact(1)*wscp*gvdwc_scp(j,i)+
+     &               welec*fact(1)*gelc(j,i)+fact(1)*wvdwpp*gvdwpp(j,i)+
+     &                wbond*gradb(j,i)+
+     &                wstrain*ghpbc(j,i)+
+     &                wcorr*fact(3)*gradcorr(j,i)+
+     &                wel_loc*fact(2)*gel_loc(j,i)+
+     &                wturn3*fact(2)*gcorr3_turn(j,i)+
+     &                wturn4*fact(3)*gcorr4_turn(j,i)+
+     &                wcorr5*fact(4)*gradcorr5(j,i)+
+     &                wcorr6*fact(5)*gradcorr6(j,i)+
+     &                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
 #else
-      do i=1,nct
+       do i=1,nct
         do j=1,3
+                if (shield_mode.eq.0) then
           gradc(j,i,icg)=wsc*gvdwc(j,i)+wscp*gvdwc_scp(j,i)+
      &                welec*fact(1)*gelc(j,i)+wstrain*ghpbc(j,i)+
      &                wbond*gradb(j,i)+
@@ -212,11 +293,34 @@ C
      &                wcorr6*fact(5)*gradcorr6(j,i)+
      &                wturn6*fact(5)*gcorr6_turn(j,i)+
      &                wsccor*fact(2)*gsccorc(j,i)
+     &               +wliptran*gliptranc(j,i)
           gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+
      &                  wbond*gradbx(j,i)+
      &                  wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+
      &                  wsccor*fact(1)*gsccorx(j,i)
-        enddo
+     &                 +wliptran*gliptranx(j,i)
+              else
+          gradc(j,i,icg)=fact(1)*wsc*gvdwc(j,i)+
+     &                   fact(1)*wscp*gvdwc_scp(j,i)+
+     &                welec*fact(1)*gelc(j,i)+wstrain*ghpbc(j,i)+
+     &                wbond*gradb(j,i)+
+     &                wcorr*fact(3)*gradcorr(j,i)+
+     &                wel_loc*fact(2)*gel_loc(j,i)+
+     &                wturn3*fact(2)*gcorr3_turn(j,i)+
+     &                wturn4*fact(3)*gcorr4_turn(j,i)+
+     &                wcorr5*fact(4)*gradcorr5(j,i)+
+     &                wcorr6*fact(5)*gradcorr6(j,i)+
+     &                wturn6*fact(5)*gcorr6_turn(j,i)+
+     &                wsccor*fact(2)*gsccorc(j,i)
+     &               +wliptran*gliptranc(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)
+         endif
+        enddo     
 #endif
       enddo
 
@@ -774,7 +878,7 @@ C
       integer icant
       external icant
       integer xshift,yshift,zshift
-      logical energy_dec /.false./
+      logical energy_dec /.true./
 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
       evdw_t=0.0d0
@@ -942,13 +1046,13 @@ c     &         aux*e2/eps(itypi,itypj)
 c            if (lprn) then
             sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
             epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
-c            write (iout,'(2(a3,i3,2x),17(0pf7.3))')
-c     &        restyp(itypi),i,restyp(itypj),j,
-c     &        epsi,sigm,chi1,chi2,chip1,chip2,
-c     &        eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
-c     &        om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
-c     &        evdwij
-c             write (iout,*) "pratial sum", evdw,evdw_t
+            write (iout,'(2(a3,i3,2x),17(0pf7.3))')
+     &        restyp(itypi),i,restyp(itypj),j,
+     &        epsi,sigm,chi1,chi2,chip1,chip2,
+     &        eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
+     &        om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
+     &        evdwij
+             write (iout,*) "pratial sum", evdw,evdw_t
 c            endif
             if (calc_grad) then
 C Calculate gradient components.
@@ -1837,6 +1941,8 @@ C
       include 'COMMON.TORSION'
       include 'COMMON.VECTORS'
       include 'COMMON.FFIELD'
+      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),
@@ -1907,15 +2013,15 @@ cd      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)
@@ -1935,16 +2041,16 @@ cd      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
           if (itel(j).eq.0) goto 1216
           ind=ind+1
           iteli=itel(i)
@@ -1975,7 +2081,7 @@ C End diagnostics
           if (yj.lt.0) yj=yj+boxysize
           zj=mod(zj,boxzsize)
           if (zj.lt.0) zj=zj+boxzsize
-      dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+      dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
       xj_safe=xj
       yj_safe=yj
       zj_safe=zj
@@ -1986,7 +2092,7 @@ C End diagnostics
           xj=xj_safe+xshift*boxxsize
           yj=yj_safe+yshift*boxysize
           zj=zj_safe+zshift*boxzsize
-          dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+          dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
           if(dist_temp.lt.dist_init) then
             dist_init=dist_temp
             xj_temp=xj
@@ -2032,7 +2138,20 @@ c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
 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          fac_shield(i)=0.4
+C          fac_shield(j)=0.6
+          el1=el1*fac_shield(i)**2*fac_shield(j)**2
+          el2=el2*fac_shield(i)**2*fac_shield(j)**2
+          eesij=(el1+el2)
           ees=ees+eesij
+          else
+          fac_shield(i)=1.0
+          fac_shield(j)=1.0
+          eesij=(el1+el2)
+          ees=ees+eesij
+          endif
+C          ees=ees+eesij
           evdw1=evdw1+evdwij*sss
 cd          write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
 cd     &      iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
@@ -2055,6 +2174,63 @@ 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
+C           enddo
+C          enddo
+           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
+
           do k=1,3
             ghalf=0.5D0*ggg(k)
             gelc(k,i)=gelc(k,i)+ghalf
@@ -2138,15 +2314,19 @@ 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) 
+     &      *fac_shield(i)**2*fac_shield(j)**2
           enddo
           do k=1,3
             ghalf=0.5D0*ggg(k)
             gelc(k,i)=gelc(k,i)+ghalf
      &               +(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)+ghalf
      &               +(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
           enddo
           do k=i+1,j-1
             do l=1,3
@@ -2394,16 +2574,70 @@ C Contribution to the local-electrostatic energy coming from the i-j pair
      &     +a33*muij(4)
 cd          write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
 cd          write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3)
+          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
 C Partial derivatives in virtual-bond dihedral angles gamma
           if (calc_grad) then
+          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
           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)
+     &    *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)
+     &    *fac_shield(i)*fac_shield(j)
+
 cd          call checkint3(i,j,mu1,mu2,a22,a23,a32,a33,acipa,eel_loc_ij)
 cd          write(iout,*) 'agg  ',agg
 cd          write(iout,*) 'aggi ',aggi
@@ -2415,6 +2649,8 @@ 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)
+     &    *fac_shield(i)*fac_shield(j)
+
           enddo
           do k=i+2,j2
             do l=1,3
@@ -2425,12 +2661,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
           ENDIF
@@ -2536,9 +2780,21 @@ c               fac3=dsqrt(-ael6i)/r0ij**3
                 fac3=dsqrt(-ael6i)*r3ij
                 ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
                 ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
+                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
 c               ees0mij=0.0D0
                 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
@@ -2603,17 +2859,29 @@ C Derivatives due to the contact function
                   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
 C Diagnostics. Comment out or remove after debugging!
@@ -2659,6 +2927,8 @@ C Third- and fourth-order contributions from turns
       include 'COMMON.TORSION'
       include 'COMMON.VECTORS'
       include 'COMMON.FFIELD'
+      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),
@@ -2667,6 +2937,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
@@ -2681,22 +2963,80 @@ cd        call checkint_turn3(i,a_temp,eello_turn3_num)
         call matmat2(EUg(1,1,i+1),EUg(1,1,i+2),auxmat(1,1))
         call transpose2(auxmat(1,1),auxmat1(1,1))
         call matmat2(a_temp(1,1),auxmat1(1,1),pizda(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)
+
 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
         if (calc_grad) then
+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 Derivatives in gamma(i)
         call matmat2(EUgder(1,1,i+1),EUg(1,1,i+2),auxmat2(1,1))
         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))
         call matmat2(a_temp(1,1),pizda(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
           a_temp(1,1)=aggi(l,1)
@@ -2706,6 +3046,8 @@ C Cartesian derivatives
           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)
           a_temp(1,2)=aggi1(l,2)
           a_temp(2,1)=aggi1(l,3)
@@ -2713,6 +3055,8 @@ C Cartesian derivatives
           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)
           a_temp(1,2)=aggj(l,2)
           a_temp(2,1)=aggj(l,3)
@@ -2720,6 +3064,8 @@ C Cartesian derivatives
           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)
@@ -2727,9 +3073,24 @@ C Cartesian derivatives
           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
         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
+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
+C     &    .or. itype(i+5).eq.ntyp1
+     &    .or. itype(i).eq.ntyp1
+C     &    .or. itype(i-1).eq.ntyp1
+     &    ) goto 178
+
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 C
 C               Fourth-order contributions
@@ -2757,11 +3118,64 @@ cd        call checkint_turn4(i,a_temp,eello_turn4_num)
         call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
         call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
         s3=0.5d0*(pizda(1,1)+pizda(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.4
+C        fac_shield(j)=0.6
+        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)
+
 cd        write (2,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3),
 cd     &    ' eello_turn4_num',8*eello_turn4_num
 C Derivatives in gamma(i)
         if (calc_grad) then
+          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
+
         call transpose2(EUgder(1,1,i+1),e1tder(1,1))
         call matmat2(e1tder(1,1),a_temp(1,1),auxmat(1,1))
         call matvec2(auxmat(1,1),Ub2(1,i+3),auxvec(1))
@@ -2769,6 +3183,8 @@ 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)) 
@@ -2777,6 +3193,8 @@ 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))
@@ -2788,6 +3206,8 @@ C Derivatives in gamma(i+2)
         call matmat2(auxmat(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
@@ -2807,6 +3227,8 @@ 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
@@ -2825,6 +3247,8 @@ 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)
@@ -2839,6 +3263,8 @@ 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)
@@ -2853,6 +3279,8 @@ 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)
@@ -2867,8 +3295,11 @@ 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,j1)=gcorr4_turn(l,j1)-(s1+s2+s3)
+     &  *fac_shield(i)*fac_shield(j)
+
         enddo
         endif
+  178 continue
       endif          
       return
       end
@@ -5549,6 +5980,8 @@ c------------------------------------------------------------------------------
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+      include 'COMMON.SHIELD'
+
       double precision gx(3),gx1(3)
       logical lprn
       lprn=.false.
@@ -5606,7 +6039,85 @@ C Calculate multi-body contributions to the gradient.
      &     ekont*(coeffp*ees0pij*gacontp_hb3(ll,kk,k)+
      &     coeffm*ees0mij*gacontm_hb3(ll,kk,k))
         enddo
-      enddo 
+      enddo
+      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
       endif
       ehbcorr=ekont*ees
       return
@@ -7883,4 +8394,340 @@ C-----------------------------------------------------------------------
       return
       end
 C-----------------------------------------------------------------------
+C first for shielding is setting of function of side-chains
+       subroutine set_shield_fac2
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.CHAIN'
+      include 'COMMON.DERIV'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.SHIELD'
+      include 'COMMON.INTERACT'
+C this is the squar root 77 devided by 81 the epislion in lipid (in protein)
+      double precision div77_81/0.974996043d0/,
+     &div4_81/0.2222222222d0/,sh_frac_dist_grad(3)
+
+C the vector between center of side_chain and peptide group
+       double precision pep_side(3),long,side_calf(3),
+     &pept_group(3),costhet_grad(3),cosphi_grad_long(3),
+     &cosphi_grad_loc(3),pep_side_norm(3),side_calf_norm(3)
+C the line belowe needs to be changed for FGPROC>1
+      do i=1,nres-1
+      if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle
+      ishield_list(i)=0
+Cif there two consequtive dummy atoms there is no peptide group between them
+C the line below has to be changed for FGPROC>1
+      VolumeTotal=0.0
+      do k=1,nres
+       if ((itype(k).eq.ntyp1).or.(itype(k).eq.10)) cycle
+       dist_pep_side=0.0
+       dist_side_calf=0.0
+       do j=1,3
+C first lets set vector conecting the ithe side-chain with kth side-chain
+      pep_side(j)=c(j,k+nres)-(c(j,i)+c(j,i+1))/2.0d0
+C      pep_side(j)=2.0d0
+C and vector conecting the side-chain with its proper calfa
+      side_calf(j)=c(j,k+nres)-c(j,k)
+C      side_calf(j)=2.0d0
+      pept_group(j)=c(j,i)-c(j,i+1)
+C lets have their lenght
+      dist_pep_side=pep_side(j)**2+dist_pep_side
+      dist_side_calf=dist_side_calf+side_calf(j)**2
+      dist_pept_group=dist_pept_group+pept_group(j)**2
+      enddo
+       dist_pep_side=dsqrt(dist_pep_side)
+       dist_pept_group=dsqrt(dist_pept_group)
+       dist_side_calf=dsqrt(dist_side_calf)
+      do j=1,3
+        pep_side_norm(j)=pep_side(j)/dist_pep_side
+        side_calf_norm(j)=dist_side_calf
+      enddo
+C now sscale fraction
+       sh_frac_dist=-(dist_pep_side-rpp(1,1)-buff_shield)/buff_shield
+C       print *,buff_shield,"buff"
+C now sscale
+        if (sh_frac_dist.le.0.0) cycle
+C If we reach here it means that this side chain reaches the shielding sphere
+C Lets add him to the list for gradient       
+        ishield_list(i)=ishield_list(i)+1
+C ishield_list is a list of non 0 side-chain that contribute to factor gradient
+C this list is essential otherwise problem would be O3
+        shield_list(ishield_list(i),i)=k
+C Lets have the sscale value
+        if (sh_frac_dist.gt.1.0) then
+         scale_fac_dist=1.0d0
+         do j=1,3
+         sh_frac_dist_grad(j)=0.0d0
+         enddo
+        else
+         scale_fac_dist=-sh_frac_dist*sh_frac_dist
+     &                   *(2.0d0*sh_frac_dist-3.0d0)
+         fac_help_scale=6.0d0*(sh_frac_dist-sh_frac_dist**2)
+     &                  /dist_pep_side/buff_shield*0.5d0
+C remember for the final gradient multiply sh_frac_dist_grad(j) 
+C for side_chain by factor -2 ! 
+         do j=1,3
+         sh_frac_dist_grad(j)=fac_help_scale*pep_side(j)
+C         sh_frac_dist_grad(j)=0.0d0
+C         scale_fac_dist=1.0d0
+C         print *,"jestem",scale_fac_dist,fac_help_scale,
+C     &                    sh_frac_dist_grad(j)
+         enddo
+        endif
+C this is what is now we have the distance scaling now volume...
+      short=short_r_sidechain(itype(k))
+      long=long_r_sidechain(itype(k))
+      costhet=1.0d0/dsqrt(1.0d0+short**2/dist_pep_side**2)
+      sinthet=short/dist_pep_side*costhet
+C now costhet_grad
+C       costhet=0.6d0
+C       sinthet=0.8
+       costhet_fac=costhet**3*short**2*(-0.5d0)/dist_pep_side**4
+C       sinthet_fac=costhet**2*0.5d0*(short**3/dist_pep_side**4*costhet
+C     &             -short/dist_pep_side**2/costhet)
+C       costhet_fac=0.0d0
+       do j=1,3
+         costhet_grad(j)=costhet_fac*pep_side(j)
+       enddo
+C remember for the final gradient multiply costhet_grad(j) 
+C for side_chain by factor -2 !
+C fac alfa is angle between CB_k,CA_k, CA_i,CA_i+1
+C pep_side0pept_group is vector multiplication  
+      pep_side0pept_group=0.0d0
+      do j=1,3
+      pep_side0pept_group=pep_side0pept_group+pep_side(j)*side_calf(j)
+      enddo
+      cosalfa=(pep_side0pept_group/
+     & (dist_pep_side*dist_side_calf))
+      fac_alfa_sin=1.0d0-cosalfa**2
+      fac_alfa_sin=dsqrt(fac_alfa_sin)
+      rkprim=fac_alfa_sin*(long-short)+short
+C      rkprim=short
+
+C now costhet_grad
+       cosphi=1.0d0/dsqrt(1.0d0+rkprim**2/dist_pep_side**2)
+C       cosphi=0.6
+       cosphi_fac=cosphi**3*rkprim**2*(-0.5d0)/dist_pep_side**4
+       sinphi=rkprim/dist_pep_side/dsqrt(1.0d0+rkprim**2/
+     &      dist_pep_side**2)
+C       sinphi=0.8
+       do j=1,3
+         cosphi_grad_long(j)=cosphi_fac*pep_side(j)
+     &+cosphi**3*0.5d0/dist_pep_side**2*(-rkprim)
+     &*(long-short)/fac_alfa_sin*cosalfa/
+     &((dist_pep_side*dist_side_calf))*
+     &((side_calf(j))-cosalfa*
+     &((pep_side(j)/dist_pep_side)*dist_side_calf))
+C       cosphi_grad_long(j)=0.0d0
+        cosphi_grad_loc(j)=cosphi**3*0.5d0/dist_pep_side**2*(-rkprim)
+     &*(long-short)/fac_alfa_sin*cosalfa
+     &/((dist_pep_side*dist_side_calf))*
+     &(pep_side(j)-
+     &cosalfa*side_calf(j)/dist_side_calf*dist_pep_side)
+C       cosphi_grad_loc(j)=0.0d0
+       enddo
+C      print *,sinphi,sinthet
+      VofOverlap=VSolvSphere/2.0d0*(1.0d0-dsqrt(1.0d0-sinphi*sinthet))
+     &                    /VSolvSphere_div
+C     &                    *wshield
+C now the gradient...
+      do j=1,3
+      grad_shield(j,i)=grad_shield(j,i)
+C gradient po skalowaniu
+     &                +(sh_frac_dist_grad(j)*VofOverlap
+C  gradient po costhet
+     &       +scale_fac_dist*VSolvSphere/VSolvSphere_div/4.0d0*
+     &(1.0d0/(-dsqrt(1.0d0-sinphi*sinthet))*(
+     &       sinphi/sinthet*costhet*costhet_grad(j)
+     &      +sinthet/sinphi*cosphi*cosphi_grad_long(j)))
+     & )*wshield
+C grad_shield_side is Cbeta sidechain gradient
+      grad_shield_side(j,ishield_list(i),i)=
+     &        (sh_frac_dist_grad(j)*-2.0d0
+     &        *VofOverlap
+     &       -scale_fac_dist*VSolvSphere/VSolvSphere_div/2.0d0*
+     &(1.0d0/(-dsqrt(1.0d0-sinphi*sinthet))*(
+     &       sinphi/sinthet*costhet*costhet_grad(j)
+     &      +sinthet/sinphi*cosphi*cosphi_grad_long(j)))
+     &       )*wshield
+
+       grad_shield_loc(j,ishield_list(i),i)=
+     &       scale_fac_dist*VSolvSphere/VSolvSphere_div/2.0d0*
+     &(1.0d0/(dsqrt(1.0d0-sinphi*sinthet))*(
+     &       sinthet/sinphi*cosphi*cosphi_grad_loc(j)
+     &        ))
+     &        *wshield
+      enddo
+      VolumeTotal=VolumeTotal+VofOverlap*scale_fac_dist
+      enddo
+      fac_shield(i)=VolumeTotal*wshield+(1.0d0-wshield)
+C      write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i)
+      enddo
+      return
+      end
+C first for shielding is setting of function of side-chains
+       subroutine set_shield_fac
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.CHAIN'
+      include 'COMMON.DERIV'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.SHIELD'
+      include 'COMMON.INTERACT'
+C this is the squar root 77 devided by 81 the epislion in lipid (in protein)
+      double precision div77_81/0.974996043d0/,
+     &div4_81/0.2222222222d0/,sh_frac_dist_grad(3)
+
+C the vector between center of side_chain and peptide group
+       double precision pep_side(3),long,side_calf(3),
+     &pept_group(3),costhet_grad(3),cosphi_grad_long(3),
+     &cosphi_grad_loc(3),pep_side_norm(3),side_calf_norm(3)
+C the line belowe needs to be changed for FGPROC>1
+      do i=1,nres-1
+      if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle
+      ishield_list(i)=0
+Cif there two consequtive dummy atoms there is no peptide group between them
+C the line below has to be changed for FGPROC>1
+      VolumeTotal=0.0
+      do k=1,nres
+       if ((itype(k).eq.ntyp1).or.(itype(k).eq.10)) cycle
+       dist_pep_side=0.0
+       dist_side_calf=0.0
+       do j=1,3
+C first lets set vector conecting the ithe side-chain with kth side-chain
+      pep_side(j)=c(j,k+nres)-(c(j,i)+c(j,i+1))/2.0d0
+C      pep_side(j)=2.0d0
+C and vector conecting the side-chain with its proper calfa
+      side_calf(j)=c(j,k+nres)-c(j,k)
+C      side_calf(j)=2.0d0
+      pept_group(j)=c(j,i)-c(j,i+1)
+C lets have their lenght
+      dist_pep_side=pep_side(j)**2+dist_pep_side
+      dist_side_calf=dist_side_calf+side_calf(j)**2
+      dist_pept_group=dist_pept_group+pept_group(j)**2
+      enddo
+       dist_pep_side=dsqrt(dist_pep_side)
+       dist_pept_group=dsqrt(dist_pept_group)
+       dist_side_calf=dsqrt(dist_side_calf)
+      do j=1,3
+        pep_side_norm(j)=pep_side(j)/dist_pep_side
+        side_calf_norm(j)=dist_side_calf
+      enddo
+C now sscale fraction
+       sh_frac_dist=-(dist_pep_side-rpp(1,1)-buff_shield)/buff_shield
+C       print *,buff_shield,"buff"
+C now sscale
+        if (sh_frac_dist.le.0.0) cycle
+C If we reach here it means that this side chain reaches the shielding sphere
+C Lets add him to the list for gradient       
+        ishield_list(i)=ishield_list(i)+1
+C ishield_list is a list of non 0 side-chain that contribute to factor gradient
+C this list is essential otherwise problem would be O3
+        shield_list(ishield_list(i),i)=k
+C Lets have the sscale value
+        if (sh_frac_dist.gt.1.0) then
+         scale_fac_dist=1.0d0
+         do j=1,3
+         sh_frac_dist_grad(j)=0.0d0
+         enddo
+        else
+         scale_fac_dist=-sh_frac_dist*sh_frac_dist
+     &                   *(2.0*sh_frac_dist-3.0d0)
+         fac_help_scale=6.0*(sh_frac_dist-sh_frac_dist**2)
+     &                  /dist_pep_side/buff_shield*0.5
+C remember for the final gradient multiply sh_frac_dist_grad(j) 
+C for side_chain by factor -2 ! 
+         do j=1,3
+         sh_frac_dist_grad(j)=fac_help_scale*pep_side(j)
+C         print *,"jestem",scale_fac_dist,fac_help_scale,
+C     &                    sh_frac_dist_grad(j)
+         enddo
+        endif
+C        if ((i.eq.3).and.(k.eq.2)) then
+C        print *,i,sh_frac_dist,dist_pep,fac_help_scale,scale_fac_dist
+C     & ,"TU"
+C        endif
+
+C this is what is now we have the distance scaling now volume...
+      short=short_r_sidechain(itype(k))
+      long=long_r_sidechain(itype(k))
+      costhet=1.0d0/dsqrt(1.0+short**2/dist_pep_side**2)
+C now costhet_grad
+C       costhet=0.0d0
+       costhet_fac=costhet**3*short**2*(-0.5)/dist_pep_side**4
+C       costhet_fac=0.0d0
+       do j=1,3
+         costhet_grad(j)=costhet_fac*pep_side(j)
+       enddo
+C remember for the final gradient multiply costhet_grad(j) 
+C for side_chain by factor -2 !
+C fac alfa is angle between CB_k,CA_k, CA_i,CA_i+1
+C pep_side0pept_group is vector multiplication  
+      pep_side0pept_group=0.0
+      do j=1,3
+      pep_side0pept_group=pep_side0pept_group+pep_side(j)*side_calf(j)
+      enddo
+      cosalfa=(pep_side0pept_group/
+     & (dist_pep_side*dist_side_calf))
+      fac_alfa_sin=1.0-cosalfa**2
+      fac_alfa_sin=dsqrt(fac_alfa_sin)
+      rkprim=fac_alfa_sin*(long-short)+short
+C now costhet_grad
+       cosphi=1.0d0/dsqrt(1.0+rkprim**2/dist_pep_side**2)
+       cosphi_fac=cosphi**3*rkprim**2*(-0.5)/dist_pep_side**4
+
+       do j=1,3
+         cosphi_grad_long(j)=cosphi_fac*pep_side(j)
+     &+cosphi**3*0.5/dist_pep_side**2*(-rkprim)
+     &*(long-short)/fac_alfa_sin*cosalfa/
+     &((dist_pep_side*dist_side_calf))*
+     &((side_calf(j))-cosalfa*
+     &((pep_side(j)/dist_pep_side)*dist_side_calf))
+
+        cosphi_grad_loc(j)=cosphi**3*0.5/dist_pep_side**2*(-rkprim)
+     &*(long-short)/fac_alfa_sin*cosalfa
+     &/((dist_pep_side*dist_side_calf))*
+     &(pep_side(j)-
+     &cosalfa*side_calf(j)/dist_side_calf*dist_pep_side)
+       enddo
+
+      VofOverlap=VSolvSphere/2.0d0*(1.0-costhet)*(1.0-cosphi)
+     &                    /VSolvSphere_div
+     &                    *wshield
+C now the gradient...
+C grad_shield is gradient of Calfa for peptide groups
+C      write(iout,*) "shield_compon",i,k,VSolvSphere,scale_fac_dist,
+C     &               costhet,cosphi
+C       write(iout,*) "cosphi_compon",i,k,pep_side0pept_group,
+C     & dist_pep_side,dist_side_calf,c(1,k+nres),c(1,k),itype(k)
+      do j=1,3
+      grad_shield(j,i)=grad_shield(j,i)
+C gradient po skalowaniu
+     &                +(sh_frac_dist_grad(j)
+C  gradient po costhet
+     &-scale_fac_dist*costhet_grad(j)/(1.0-costhet)
+     &-scale_fac_dist*(cosphi_grad_long(j))
+     &/(1.0-cosphi) )*div77_81
+     &*VofOverlap
+C grad_shield_side is Cbeta sidechain gradient
+      grad_shield_side(j,ishield_list(i),i)=
+     &        (sh_frac_dist_grad(j)*-2.0d0
+     &       +scale_fac_dist*costhet_grad(j)*2.0d0/(1.0-costhet)
+     &       +scale_fac_dist*(cosphi_grad_long(j))
+     &        *2.0d0/(1.0-cosphi))
+     &        *div77_81*VofOverlap
+
+       grad_shield_loc(j,ishield_list(i),i)=
+     &   scale_fac_dist*cosphi_grad_loc(j)
+     &        *2.0d0/(1.0-cosphi)
+     &        *div77_81*VofOverlap
+      enddo
+      VolumeTotal=VolumeTotal+VofOverlap*scale_fac_dist
+      enddo
+      fac_shield(i)=VolumeTotal*div77_81+div4_81
+C      write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i)
+      enddo
+      return
+      end
+C--------------------------------------------------------------------------
 

diff --git a/source/cluster/wham/src-M/main_clust.F b/source/cluster/wham/src-M/main_clust.F
index 892a6c7..dc77f04 100644 (file)
--- a/source/cluster/wham/src-M/main_clust.F
+++ b/source/cluster/wham/src-M/main_clust.F
@@ -122,7 +122,7 @@ c      call flush(iout)
       ndis=ncon_work*(ncon_work-1)/2
       call work_partition(.true.,ndis)
 #endif
-
+      write(iout,*) "AFTET wort_part",NCON_work
       DO I=1,NCON_work
         ICC(I)=I
       ENDDO
@@ -132,6 +132,8 @@ C
 C CALCULATE DISTANCES
 C
       call daread_ccoords(1,ncon_work)
+      write (iout,*) "AM I HERE"
+      call flush(iout)
       ind1=0
       DO I=1,NCON_work-1
         if (mod(i,100).eq.0) print *,'Calculating RMS i=',i

diff --git a/source/cluster/wham/src-M/parmread.F b/source/cluster/wham/src-M/parmread.F
index 1368050..4292bd6 100644 (file)
--- a/source/cluster/wham/src-M/parmread.F
+++ b/source/cluster/wham/src-M/parmread.F
@@ -862,9 +862,13 @@ C----------------------- LJK potential --------------------------------
       endif
       goto 50
 C---------------------- GB or BP potential -----------------------------
-   30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
-     &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
-     &  (alp(i),i=1,ntyp)
+   30 do i=1,ntyp
+       read (isidep,*)(eps(i,j),j=i,ntyp)
+      enddo
+      read (isidep,*)(sigma0(i),i=1,ntyp)
+      read (isidep,*)(sigii(i),i=1,ntyp)
+      read (isidep,*)(chip0(i),i=1,ntyp)
+      read (isidep,*)(alp(i),i=1,ntyp)
 C For the GB potential convert sigma'**2 into chi'
       if (ipot.eq.4) then
        do i=1,ntyp

diff --git a/source/cluster/wham/src-M/probabl.F b/source/cluster/wham/src-M/probabl.F
index 0c45402..ec3fbff 100644 (file)
--- a/source/cluster/wham/src-M/probabl.F
+++ b/source/cluster/wham/src-M/probabl.F
@@ -42,7 +42,8 @@ c      enddo
       write (iout,*) me," indstart",indstart(me)," indend",indend(me)
       call daread_ccoords(indstart(me),indend(me))
 #endif
-c      write (iout,*) "ncon",ncon
+      write (iout,*) "ncon",ncon
+      call flush(iout)
       temper=1.0d0/(beta_h(ib)*1.987D-3)
 c      write (iout,*) "ib",ib," beta_h",beta_h(ib)," temper",temper
 c      quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
@@ -54,6 +55,7 @@ c        quotl=quotl*quot
 c        kfacl=kfacl*kfac
 c        fT(l)=kfacl/(kfacl-1.0d0+quotl)
 c      enddo
+#define DEBUG
             if (rescale_mode.eq.1) then
               quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
               quotl=1.0d0
@@ -88,8 +90,8 @@ c      enddo
                 fT(l)=1.12692801104297249644d0/
      &             dlog(dexp(quotl)+dexp(-quotl))
               enddo
-c              write (iout,*) 1.0d0/(beta_h(ib)*1.987D-3),ft
-c              call flush(iout)
+              write (iout,*) 1.0d0/(beta_h(ib)*1.987D-3),ft
+              call flush(iout)
 #if defined(FUNCTH)
               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
      &                  320.0d0
@@ -114,25 +116,33 @@ c              call flush(iout)
       do i=1,ncon
         ii=i
 #endif
-c        write (iout,*) "i",i," ii",ii
-c        call flush(iout)
+        write (iout,*) "i",i," ii",ii,"ib",ib,scount(me)
+        call flush(iout)
         if (ib.eq.1) then
           do j=1,nres
             do k=1,3
               c(k,j)=allcart(k,j,i)
               c(k,j+nres)=allcart(k,j+nres,i)
+              write(iout,*) "coord",i,j,k,allcart(k,j,i),c(k,j),
+     &        c(k,j+nres),allcart(k,j+nres,i)
             enddo
           enddo
+          write(iout,*) "out of j loop"
+          call flush(iout)
           do k=1,3
             c(k,nres+1)=c(k,1)
             c(k,nres+nres)=c(k,nres)
           enddo
+          write(iout,*) "after nres+nres",nss_all(i)
+          call flush(iout)
           nss=nss_all(i)
           do j=1,nss
             ihpb(j)=ihpb_all(j,i)
             jhpb(j)=jhpb_all(j,i)
           enddo 
           call int_from_cart1(.false.)
+          write(iout,*) "before etotal"
+          call flush(iout)
           call etotal(energia(0),fT)
           totfree(i)=energia(0)         
           totfree_buf(i)=totfree(i)
@@ -207,7 +217,7 @@ c#endif
         write (iout,*) "evdw2", wscp, evdw2
         write (iout,*) "welec", ft(1),welec,ees
         write (iout,*) "evdw1", wvdwpp,evdw1
-        write (iout,*) "ebe" ebe,wang
+        write (iout,*) "ebe", ebe,wang
 #endif        
         Fdimless(i)=beta_h(ib)*etot+entfac(ii)
         Fdimless_buf(i)=Fdimless(i)

diff --git a/source/cluster/wham/src-M/read_coords.F b/source/cluster/wham/src-M/read_coords.F
index 17f39a3..ac11ac0 100644 (file)
--- a/source/cluster/wham/src-M/read_coords.F
+++ b/source/cluster/wham/src-M/read_coords.F
@@ -178,13 +178,13 @@ c through a ring.
 #endif
         endif
 
-#define DEBUG
+C#define DEBUG
 #ifdef DEBUG
         write (iout,*) "Opening file ",intinname(:ilen(intinname))
         write (iout,*) "lenrec",lenrec_in
         call flush(iout)
 #endif
-#undef DEBUG
+C#undef DEBUG
 c        write (iout,*) "maxconf",maxconf
         i=0
         do while (.true.)
@@ -283,6 +283,7 @@ c            write (iout,*) "nss",nss
               enddo
             enddo
           endif
+#define DEBUG
 #ifdef DEBUG
           write (iout,'(5hREAD ,i5,3f15.4,i10)') 
      &     jj+1,energy(jj+1),entfac(jj+1),
@@ -292,6 +293,7 @@ c            write (iout,*) "nss",nss
           write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
           call flush(iout)
 #endif
+#undef DEBUG
           call add_new_cconf(jjj,jj,jj_old,icount,Next)
         enddo
   101   continue
@@ -640,10 +642,11 @@ c------------------------------------------------------------------------------
       integer i,j,ij,ii,iii
       integer len
       character*16 form,acc
-      character*32 nam
+      character*80 nam
 c
 c Read conformations off a DA scratchfile.
 c
+#define DEBUG
 #ifdef DEBUG
       write (iout,*) "DAREAD_COORDS"
       write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
@@ -680,7 +683,10 @@ c     &    nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss),
      &    jhpb_all(i,ij),i=1,nss)
         call flush(iout)
 #endif
+#undef DEBUG
       enddo
+      write (iout,*) "just before leave"
+      call flush(iout)
       return
       end
 c------------------------------------------------------------------------------

diff --git a/source/cluster/wham/src-M/readrtns.F b/source/cluster/wham/src-M/readrtns.F
index ea51dbf..6f8f5e4 100644 (file)
--- a/source/cluster/wham/src-M/readrtns.F
+++ b/source/cluster/wham/src-M/readrtns.F
@@ -16,12 +16,13 @@ C
       include 'COMMON.FREE'
       include 'COMMON.INTERACT'
       include "COMMON.SPLITELE"
+      include 'COMMON.SHIELD'
       character*320 controlcard,ucase
 #ifdef MPL
       include 'COMMON.INFO'
 #endif
       integer i,i1,i2,it1,it2
-
+      double precision pi
       read (INP,'(a80)') titel
       call card_concat(controlcard)
 
@@ -36,6 +37,26 @@ C Reading the dimensions of box in x,y,z coordinates
 c Cutoff range for interactions
       call reada(controlcard,"R_CUT",r_cut,15.0d0)
       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
+C Shielding mode
+      call readi(controlcard,'SHIELD',shield_mode,0)
+      write (iout,*) "SHIELD MODE",shield_mode
+      if (shield_mode.gt.0) then
+      pi=3.141592d0
+C VSolvSphere the volume of solving sphere
+C      print *,pi,"pi"
+C rpp(1,1) is the energy r0 for peptide group contact and will be used for it 
+C there will be no distinction between proline peptide group and normal peptide
+C group in case of shielding parameters
+      VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
+      VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
+      write (iout,*) VSolvSphere,VSolvSphere_div
+C long axis of side chain 
+      do i=1,ntyp
+      long_r_sidechain(i)=vbldsc0(1,i)
+      short_r_sidechain(i)=sigma0(i)
+      enddo
+      buff_shield=1.0d0
+      endif
       call readi(controlcard,'PDBOUT',outpdb,0)
       call readi(controlcard,'MOL2OUT',outmol2,0)
       refstr=(index(controlcard,'REFSTR').gt.0)
@@ -103,6 +124,7 @@ C
       include 'COMMON.CONTACTS'
       include 'COMMON.TIME1'
       include 'COMMON.TORCNSTR'
+      include 'COMMON.SHIELD'
 #ifdef MPL
       include 'COMMON.INFO'
 #endif
@@ -118,7 +140,9 @@ C Body
 C
 C Read weights of the subsequent energy terms.
       call card_concat(weightcard)
-      call reada(weightcard,'WSC',wsc,1.0d0)
+      write(iout,*) weightcard
+C      call reada(weightcard,'WSC',wsc,1.0d0)
+      write(iout,*) wsc
       call reada(weightcard,'WLONG',wsc,wsc)
       call reada(weightcard,'WSCP',wscp,1.0d0)
       call reada(weightcard,'WELEC',welec,1.0D0)
@@ -150,6 +174,8 @@ C Read weights of the subsequent energy terms.
       call reada(weightcard,"V2SS",v2ss,7.61d0)
       call reada(weightcard,"V3SS",v3ss,13.7d0)
       call reada(weightcard,"EBR",ebr,-5.50D0)
+      call reada(weightcard,'WSHIELD',wshield,1.0d0)
+      write(iout,*) 'WSHIELD',wshield
       call reada(weightcard,"ATRISS",atriss,0.301D0)
       call reada(weightcard,"BTRISS",btriss,0.021D0)
       call reada(weightcard,"CTRISS",ctriss,1.001D0)

diff --git a/source/cluster/wham/src-M/work_partition.F b/source/cluster/wham/src-M/work_partition.F
index e31db53..ba87bd2 100644 (file)
--- a/source/cluster/wham/src-M/work_partition.F
+++ b/source/cluster/wham/src-M/work_partition.F
@@ -74,12 +74,13 @@ c        print *,"N",n," NCON_WORK",ncon_work
 
       if (lprint) then
         write (iout,*) "Partition of work between processors"
-          do i=0,nprocs-1
-            write (iout,'(a,i5,a,i7,a,i7,a,i7)')
-     &        "Processor",i," indstart",indstart(i),
-     &        " indend",indend(i)," count",scount(i)
-          enddo
-      endif
+C          do i=0,nprocs-1
+C            write (iout,'(a,i5,a,i7,a,i7,a,i7)')
+C     &        "Processor",i," indstart",indstart(i),
+C     &        " indend",indend(i)," count",scount(i)
+C         enddo
+      endif 
+      write(iout,*) "just before leave"
       return
       end
 #endif

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 64288d6..3ad8f52 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -3661,6 +3661,7 @@ C 13-go grudnia roku pamietnego...
       double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
      &                   0.0d0,1.0d0,0.0d0,
      &                   0.0d0,0.0d0,1.0d0/
+       integer xshift,yshift,zshift
 c          time00=MPI_Wtime()
 cd      write (iout,*) "eelecij",i,j
 c          ind=ind+1

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))

diff --git a/source/wham/src-M/parmread.F b/source/wham/src-M/parmread.F
index 8db78b9..cd2bce9 100644 (file)
--- a/source/wham/src-M/parmread.F
+++ b/source/wham/src-M/parmread.F
@@ -1228,7 +1228,7 @@ c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
           endif
          if (lprint) then
             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') 
-     &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
+     &      restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
          endif
         enddo

diff --git a/source/wham/src-M/wham_calc1.F b/source/wham/src-M/wham_calc1.F
index eff98da..1d27235 100644 (file)
--- a/source/wham/src-M/wham_calc1.F
+++ b/source/wham/src-M/wham_calc1.F
@@ -1,4 +1,4 @@
-      subroutine WHAM_CALC(islice,*)
+             subroutine WHAM_CALC(islice,*)
 ! Weighed Histogram Analysis Method (WHAM) code
 ! Written by A. Liwo based on the work of Kumar et al., 
 ! J.Comput.Chem., 13, 1011 (1992)
@@ -36,6 +36,7 @@ c      parameter (MaxHdim=200000)
       include "COMMON.SBRIDGE"
       include "COMMON.PROT"
       include "COMMON.ENEPS"
+      include "COMMON.SHIELD"
       integer MaxPoint,MaxPointProc
       parameter (MaxPoint=MaxStr,
      & MaxPointProc=MaxStr_Proc)
@@ -113,6 +114,7 @@ c      parameter (MaxHdim=200000)
       do t=0,MaxN
         htot(t)=0
       enddo
+#define DEBUG
 #ifdef MPI
       do i=1,scount(me1)
 #else
@@ -820,7 +822,7 @@ C     &            +ftprim(6)*evdw_t
      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
      &            ftprim(1)*wsccor*esccor
             ebis=ftbis(1)*wsc*(evdw+ft(6)*evdw_t)
-     &            +ftbis(1)**wscp*evdw2+
+     &            +ftbis(1)*wscp*evdw2+
      &            ftbis(1)*welec*ees
      &            +ftbis(1)*wvdwpp*evdw
      &            +ftbis(1)*wtor*etors+
@@ -1284,5 +1286,5 @@ C     &            +ftprim(6)*evdw_t
 #endif
 
       return
-
+#undef DEBUG
       end