subroutine eelecij_scale(i,j,ees,evdw1,eel_loc)
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
#ifdef MPI
      include "mpif.h"
#endif
      include 'COMMON.CONTROL'
      include 'COMMON.IOUNITS'
      include 'COMMON.GEO'
      include 'COMMON.VAR'
      include 'COMMON.LOCAL'
      include 'COMMON.CHAIN'
      include 'COMMON.DERIV'
      include 'COMMON.INTERACT'
      include 'COMMON.CONTACTS'
      include 'COMMON.TORSION'
      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),
     &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
      common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
     &    dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
     &    num_conti,j1,j2
c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
#ifdef MOMENT
      double precision scal_el /1.0d0/
#else
      double precision scal_el /0.5d0/
#endif
C 12/13/98 
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/
c          time00=MPI_Wtime()
cd      write (iout,*) "eelecij",i,j
C      print *,"WCHODZE2"
          ind=ind+1
          iteli=itel(i)
          itelj=itel(j)
          if (j.eq.i+2 .and. itelj.eq.2) iteli=2
          aaa=app(iteli,itelj)
          bbb=bpp(iteli,itelj)
          ael6i=ael6(iteli,itelj)
          ael3i=ael3(iteli,itelj) 
          dxj=dc(1,j)
          dyj=dc(2,j)
          dzj=dc(3,j)
          dx_normj=dc_norm(1,j)
          dy_normj=dc_norm(2,j)
          dz_normj=dc_norm(3,j)
          xj=c(1,j)+0.5D0*dxj
          yj=c(2,j)+0.5D0*dyj
          zj=c(3,j)+0.5D0*dzj
         xj=mod(xj,boxxsize)
          if (xj.lt.0) xj=xj+boxxsize
          yj=mod(yj,boxysize)
          if (yj.lt.0) yj=yj+boxysize
          zj=mod(zj,boxzsize)
          if (zj.lt.0) zj=zj+boxzsize
      dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
      xj_safe=xj
      yj_safe=yj
      zj_safe=zj
      isubchap=0
      do xshift=-1,1
      do yshift=-1,1
      do zshift=-1,1
          xj=xj_safe+xshift*boxxsize
          yj=yj_safe+yshift*boxysize
          zj=zj_safe+zshift*boxzsize
          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
            yj_temp=yj
            zj_temp=zj
            isubchap=1
          endif
       enddo
       enddo
       enddo
       if (isubchap.eq.1) then
          xj=xj_temp-xmedi
          yj=yj_temp-ymedi
          zj=zj_temp-zmedi
       else
          xj=xj_safe-xmedi
          yj=yj_safe-ymedi
          zj=zj_safe-zmedi
       endif

          rij=xj*xj+yj*yj+zj*zj
          rrmij=1.0D0/rij
          rij=dsqrt(rij)
          rmij=1.0D0/rij
c For extracting the short-range part of Evdwpp
          sss=sscale(rij/rpp(iteli,itelj))
          sssgrad=sscagrad(rij/rpp(iteli,itelj))
          r3ij=rrmij*rmij
          r6ij=r3ij*r3ij  
          cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
          cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
          cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
          fac=cosa-3.0D0*cosb*cosg
          ev1=aaa*r6ij*r6ij
c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
          if (j.eq.i+2) ev1=scal_el*ev1
          ev2=bbb*r6ij
          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)
          ees=ees+eesij
          evdw1=evdw1+evdwij*(1.0d0-sss)
cd          write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
cd     &      iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
cd     &      1.0D0/dsqrt(rrmij),evdwij,eesij,
cd     &      xmedi,ymedi,zmedi,xj,yj,zj

          if (energy_dec) then 
              write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
              write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
          endif

C
C Calculate contributions to the Cartesian gradient.
C
#ifdef SPLITELE
          facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)
          facel=-3*rrmij*(el1+eesij)
          fac1=fac
          erij(1)=xj*rmij
          erij(2)=yj*rmij
          erij(3)=zj*rmij
*
* Radial derivatives. First process both termini of the fragment (i,j)
*
          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
c            gelc(k,j)=gelc(k,j)+ghalf
c          enddo
c 9/28/08 AL Gradient compotents will be summed only at the end
          do k=1,3
            gelc_long(k,j)=gelc_long(k,j)+ggg(k)
            gelc_long(k,i)=gelc_long(k,i)-ggg(k)
          enddo
*
* Loop over residues i+1 thru j-1.
*
cgrad          do k=i+1,j-1
cgrad            do l=1,3
cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
cgrad            enddo
cgrad          enddo
          ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
          ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
          ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
c          do k=1,3
c            ghalf=0.5D0*ggg(k)
c            gvdwpp(k,i)=gvdwpp(k,i)+ghalf
c            gvdwpp(k,j)=gvdwpp(k,j)+ghalf
c          enddo
c 9/28/08 AL Gradient compotents will be summed only at the end
          do k=1,3
            gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
            gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
          enddo
*
* Loop over residues i+1 thru j-1.
*
cgrad          do k=i+1,j-1
cgrad            do l=1,3
cgrad              gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
cgrad            enddo
cgrad          enddo
#else
          facvdw=ev1+evdwij*(1.0d0-sss) 
          facel=el1+eesij  
          fac1=fac
          fac=-3*rrmij*(facvdw+facvdw+facel)
          erij(1)=xj*rmij
          erij(2)=yj*rmij
          erij(3)=zj*rmij
*
* Radial derivatives. First process both termini of the fragment (i,j)
* 
          ggg(1)=fac*xj
          ggg(2)=fac*yj
          ggg(3)=fac*zj
c          do k=1,3
c            ghalf=0.5D0*ggg(k)
c            gelc(k,i)=gelc(k,i)+ghalf
c            gelc(k,j)=gelc(k,j)+ghalf
c          enddo
c 9/28/08 AL Gradient compotents will be summed only at the end
          do k=1,3
            gelc_long(k,j)=gelc(k,j)+ggg(k)
            gelc_long(k,i)=gelc(k,i)-ggg(k)
          enddo
*
* Loop over residues i+1 thru j-1.
*
cgrad          do k=i+1,j-1
cgrad            do l=1,3
cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
cgrad            enddo
cgrad          enddo
c 9/28/08 AL Gradient compotents will be summed only at the end
C          ggg(1)=facvdw*xj
C          ggg(2)=facvdw*yj
C          ggg(3)=facvdw*zj
          ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
          ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
          ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
          do k=1,3
            gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
            gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
          enddo
#endif
*
* Angular part
*          
          ecosa=2.0D0*fac3*fac1+fac4
          fac4=-3.0D0*fac4
          fac3=-6.0D0*fac3
          ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
          ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
          do k=1,3
            dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
            dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
          enddo
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
c          do k=1,3
c            ghalf=0.5D0*ggg(k)
c            gelc(k,i)=gelc(k,i)+ghalf
c     &               +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
c     &               + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
c            gelc(k,j)=gelc(k,j)+ghalf
c     &               +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
c     &               + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
c          enddo
cgrad          do k=i+1,j-1
cgrad            do l=1,3
cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
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))
     &           *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))
     &           *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
          IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
     &        .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 
     &        .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
C
C 9/25/99 Mixed third-order local-electrostatic terms. The local-interaction 
C   energy of a peptide unit is assumed in the form of a second-order 
C   Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al.
C   Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
C   are computed for EVERY pair of non-contiguous peptide groups.
C
          if (j.lt.nres-1) then
            j1=j+1
            j2=j-1
          else
            j1=j-1
            j2=j-2
          endif
          kkk=0
          do k=1,2
            do l=1,2
              kkk=kkk+1
              muij(kkk)=mu(k,i)*mu(l,j)
            enddo
          enddo  
cd         write (iout,*) 'EELEC: i',i,' j',j
cd          write (iout,*) 'j',j,' j1',j1,' j2',j2
cd          write(iout,*) 'muij',muij
          ury=scalar(uy(1,i),erij)
          urz=scalar(uz(1,i),erij)
          vry=scalar(uy(1,j),erij)
          vrz=scalar(uz(1,j),erij)
          a22=scalar(uy(1,i),uy(1,j))-3*ury*vry
          a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz
          a32=scalar(uz(1,i),uy(1,j))-3*urz*vry
          a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz
          fac=dsqrt(-ael6i)*r3ij
          a22=a22*fac
          a23=a23*fac
          a32=a32*fac
          a33=a33*fac
cd          write (iout,'(4i5,4f10.5)')
cd     &     i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
cd          write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
cd          write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
cd     &      uy(:,j),uz(:,j)
cd          write (iout,'(4f10.5)') 
cd     &      scalar(uy(1,i),uy(1,j)),scalar(uy(1,i),uz(1,j)),
cd     &      scalar(uz(1,i),uy(1,j)),scalar(uz(1,i),uz(1,j))
cd          write (iout,'(4f10.5)') ury,urz,vry,vrz
cd           write (iout,'(9f10.5/)') 
cd     &      fac22,a22,fac23,a23,fac32,a32,fac33,a33,eel_loc_ij
C Derivatives of the elements of A in virtual-bond vectors
          call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
          do k=1,3
            uryg(k,1)=scalar(erder(1,k),uy(1,i))
            uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1))
            uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1))
            urzg(k,1)=scalar(erder(1,k),uz(1,i))
            urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1))
            urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1))
            vryg(k,1)=scalar(erder(1,k),uy(1,j))
            vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1))
            vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1))
            vrzg(k,1)=scalar(erder(1,k),uz(1,j))
            vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1))
            vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1))
          enddo
C Compute radial contributions to the gradient
          facr=-3.0d0*rrmij
          a22der=a22*facr
          a23der=a23*facr
          a32der=a32*facr
          a33der=a33*facr
          agg(1,1)=a22der*xj
          agg(2,1)=a22der*yj
          agg(3,1)=a22der*zj
          agg(1,2)=a23der*xj
          agg(2,2)=a23der*yj
          agg(3,2)=a23der*zj
          agg(1,3)=a32der*xj
          agg(2,3)=a32der*yj
          agg(3,3)=a32der*zj
          agg(1,4)=a33der*xj
          agg(2,4)=a33der*yj
          agg(3,4)=a33der*zj
C Add the contributions coming from er
          fac3=-3.0d0*fac
          do k=1,3
            agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury)
            agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury)
            agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz)
            agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz)
          enddo
          do k=1,3
C Derivatives in DC(i) 
cgrad            ghalf1=0.5d0*agg(k,1)
cgrad            ghalf2=0.5d0*agg(k,2)
cgrad            ghalf3=0.5d0*agg(k,3)
cgrad            ghalf4=0.5d0*agg(k,4)
            aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j))
     &      -3.0d0*uryg(k,2)*vry)!+ghalf1
            aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j))
     &      -3.0d0*uryg(k,2)*vrz)!+ghalf2
            aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j))
     &      -3.0d0*urzg(k,2)*vry)!+ghalf3
            aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j))
     &      -3.0d0*urzg(k,2)*vrz)!+ghalf4
C Derivatives in DC(i+1)
            aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j))
     &      -3.0d0*uryg(k,3)*vry)!+agg(k,1)
            aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j))
     &      -3.0d0*uryg(k,3)*vrz)!+agg(k,2)
            aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j))
     &      -3.0d0*urzg(k,3)*vry)!+agg(k,3)
            aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j))
     &      -3.0d0*urzg(k,3)*vrz)!+agg(k,4)
C Derivatives in DC(j)
            aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i))
     &      -3.0d0*vryg(k,2)*ury)!+ghalf1
            aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i))
     &      -3.0d0*vrzg(k,2)*ury)!+ghalf2
            aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i))
     &      -3.0d0*vryg(k,2)*urz)!+ghalf3
            aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i)) 
     &      -3.0d0*vrzg(k,2)*urz)!+ghalf4
C Derivatives in DC(j+1) or DC(nres-1)
            aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i))
     &      -3.0d0*vryg(k,3)*ury)
            aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i))
     &      -3.0d0*vrzg(k,3)*ury)
            aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i))
     &      -3.0d0*vryg(k,3)*urz)
            aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i)) 
     &      -3.0d0*vrzg(k,3)*urz)
cgrad            if (j.eq.nres-1 .and. i.lt.j-2) then
cgrad              do l=1,4
cgrad                aggj1(k,l)=aggj1(k,l)+agg(k,l)
cgrad              enddo
cgrad            endif
          enddo
          acipa(1,1)=a22
          acipa(1,2)=a23
          acipa(2,1)=a32
          acipa(2,2)=a33
          a22=-a22
          a23=-a23
          do l=1,2
            do k=1,3
              agg(k,l)=-agg(k,l)
              aggi(k,l)=-aggi(k,l)
              aggi1(k,l)=-aggi1(k,l)
              aggj(k,l)=-aggj(k,l)
              aggj1(k,l)=-aggj1(k,l)
            enddo
          enddo
          if (j.lt.nres-1) then
            a22=-a22
            a32=-a32
            do l=1,3,2
              do k=1,3
                agg(k,l)=-agg(k,l)
                aggi(k,l)=-aggi(k,l)
                aggi1(k,l)=-aggi1(k,l)
                aggj(k,l)=-aggj(k,l)
                aggj1(k,l)=-aggj1(k,l)
              enddo
            enddo
          else
            a22=-a22
            a23=-a23
            a32=-a32
            a33=-a33
            do l=1,4
              do k=1,3
                agg(k,l)=-agg(k,l)
                aggi(k,l)=-aggi(k,l)
                aggi1(k,l)=-aggi1(k,l)
                aggj(k,l)=-aggj(k,l)
                aggj1(k,l)=-aggj1(k,l)
              enddo
            enddo 
          endif    
          ENDIF ! WCORR
          IF (wel_loc.gt.0.0d0) THEN
C Contribution to the local-electrostatic energy coming from the i-j pair
          eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
     &     +a33*muij(4)
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))
     &    *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)

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)

            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)
cgrad            gel_loc(l,i)=gel_loc(l,i)+ghalf
cgrad            gel_loc(l,j)=gel_loc(l,j)+ghalf
          enddo
cgrad          do k=i+1,j2
cgrad            do l=1,3
cgrad              gel_loc(l,k)=gel_loc(l,k)+ggg(l)
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))
     &    *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
c          if (j.gt.i+1 .and. num_conti.le.maxconts) then
          if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0
     &       .and. num_conti.le.maxconts) then
c            write (iout,*) i,j," entered corr"
C
C Calculate the contact function. The ith column of the array JCONT will 
C contain the numbers of atoms that make contacts with the atom I (of numbers
C greater than I). The arrays FACONT and GACONT will contain the values of
C the contact function and its derivative.
c           r0ij=1.02D0*rpp(iteli,itelj)
c           r0ij=1.11D0*rpp(iteli,itelj)
            r0ij=2.20D0*rpp(iteli,itelj)
c           r0ij=1.55D0*rpp(iteli,itelj)
            call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
            if (fcont.gt.0.0D0) then
              num_conti=num_conti+1
              if (num_conti.gt.maxconts) then
                write (iout,*) 'WARNING - max. # of contacts exceeded;',
     &                         ' will skip next contacts for this conf.'
              else
                jcont_hb(num_conti,i)=j
cd                write (iout,*) "i",i," j",j," num_conti",num_conti,
cd     &           " jcont_hb",jcont_hb(num_conti,i)
                IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. 
     &          wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
C 9/30/99 (AL) - store components necessary to evaluate higher-order loc-el
C  terms.
                d_cont(num_conti,i)=rij
cd                write (2,'(3e15.5)') rij,r0ij+0.2d0*r0ij,rij
C     --- Electrostatic-interaction matrix --- 
                a_chuj(1,1,num_conti,i)=a22
                a_chuj(1,2,num_conti,i)=a23
                a_chuj(2,1,num_conti,i)=a32
                a_chuj(2,2,num_conti,i)=a33
C     --- Gradient of rij
                do kkk=1,3
                  grij_hb_cont(kkk,num_conti,i)=erij(kkk)
                enddo
                kkll=0
                do k=1,2
                  do l=1,2
                    kkll=kkll+1
                    do m=1,3
                      a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll)
                      a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll)
                      a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll)
                      a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll)
                      a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll)
                    enddo
                  enddo
                enddo
                ENDIF
                IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN
C Calculate contact energies
                cosa4=4.0D0*cosa
                wij=cosa-3.0D0*cosb*cosg
                cosbg1=cosb+cosg
                cosbg2=cosb-cosg
c               fac3=dsqrt(-ael6i)/r0ij**3     
                fac3=dsqrt(-ael6i)*r3ij
c                 ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
                ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
                if (ees0tmp.gt.0) then
                  ees0pij=dsqrt(ees0tmp)
                else
                  ees0pij=0
                endif
c                ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
                ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
                if (ees0tmp.gt.0) then
                  ees0mij=dsqrt(ees0tmp)
                else
                  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
c               ees0m(num_conti,i)=0.0D0
C End diagnostics.
c               write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij,
c    & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont
C Angular derivatives of the contact function
                ees0pij1=fac3/ees0pij 
                ees0mij1=fac3/ees0mij
                fac3p=-3.0D0*fac3*rrmij
                ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
                ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
c               ees0mij1=0.0D0
                ecosa1=       ees0pij1*( 1.0D0+0.5D0*wij)
                ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
                ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
                ecosa2=       ees0mij1*(-1.0D0+0.5D0*wij)
                ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2) 
                ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
                ecosap=ecosa1+ecosa2
                ecosbp=ecosb1+ecosb2
                ecosgp=ecosg1+ecosg2
                ecosam=ecosa1-ecosa2
                ecosbm=ecosb1-ecosb2
                ecosgm=ecosg1-ecosg2
C Diagnostics
c               ecosap=ecosa1
c               ecosbp=ecosb1
c               ecosgp=ecosg1
c               ecosam=0.0D0
c               ecosbm=0.0D0
c               ecosgm=0.0D0
C End diagnostics
                facont_hb(num_conti,i)=fcont
                fprimcont=fprimcont/rij
cd              facont_hb(num_conti,i)=1.0D0
C Following line is for diagnostics.
cd              fprimcont=0.0D0
                do k=1,3
                  dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
                  dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
                enddo
                do k=1,3
                  gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
                  gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
                enddo
                gggp(1)=gggp(1)+ees0pijp*xj
                gggp(2)=gggp(2)+ees0pijp*yj
                gggp(3)=gggp(3)+ees0pijp*zj
                gggm(1)=gggm(1)+ees0mijp*xj
                gggm(2)=gggm(2)+ees0mijp*yj
                gggm(3)=gggm(3)+ees0mijp*zj
C Derivatives due to the contact function
                gacont_hbr(1,num_conti,i)=fprimcont*xj
                gacont_hbr(2,num_conti,i)=fprimcont*yj
                gacont_hbr(3,num_conti,i)=fprimcont*zj
                do k=1,3
c
c 10/24/08 cgrad and ! comments indicate the parts of the code removed 
c          following the change of gradient-summation algorithm.
c
cgrad                  ghalfp=0.5D0*gggp(k)
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
            endif  ! fcont.gt.0
          endif    ! j.gt.i+1
          if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
            do k=1,4
              do l=1,3
                ghalf=0.5d0*agg(l,k)
                aggi(l,k)=aggi(l,k)+ghalf
                aggi1(l,k)=aggi1(l,k)+agg(l,k)
                aggj(l,k)=aggj(l,k)+ghalf
              enddo
            enddo
            if (j.eq.nres-1 .and. i.lt.j-2) then
              do k=1,4
                do l=1,3
                  aggj1(l,k)=aggj1(l,k)+agg(l,k)
                enddo
              enddo
            endif
          endif
c          t_eelecij=t_eelecij+MPI_Wtime()-time00
      return
      end