subroutine eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
C
C This subroutine calculates the average interaction energy and its gradient
C in the virtual-bond vectors between non-adjacent peptide groups, based on 
C the potential described in Liwo et al., Protein Sci., 1993, 2, 1715. 
C The potential depends both on the distance of peptide-group centers and on 
C the orientation of the CA-CA virtual bonds.
C 
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      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'
      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,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/
cd      write(iout,*) 'In EELEC'
cd      do i=1,nloctyp
cd        write(iout,*) 'Type',i
cd        write(iout,*) 'B1',B1(:,i)
cd        write(iout,*) 'B2',B2(:,i)
cd        write(iout,*) 'CC',CC(:,:,i)
cd        write(iout,*) 'DD',DD(:,:,i)
cd        write(iout,*) 'EE',EE(:,:,i)
cd      enddo
cd      call check_vecgrad
cd      stop
      if (icheckgrad.eq.1) then
        do i=1,nres-1
          fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i)))
          do k=1,3
            dc_norm(k,i)=dc(k,i)*fac
          enddo
c          write (iout,*) 'i',i,' fac',fac
        enddo
      endif
      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        call vec_and_deriv
        call set_matrices
      endif
cd      do i=1,nres-1
cd        write (iout,*) 'i=',i
cd        do k=1,3
cd        write (iout,'(i5,2f10.5)') k,uy(k,i),uz(k,i)
cd        enddo
cd        do k=1,3
cd          write (iout,'(f10.5,2x,3f10.5,2x,3f10.5)') 
cd     &     uz(k,i),(uzgrad(k,l,1,i),l=1,3),(uzgrad(k,l,2,i),l=1,3)
cd        enddo
cd      enddo
      num_conti_hb=0
      ees=0.0D0
      evdw1=0.0D0
      eel_loc=0.0d0 
      eello_turn3=0.0d0
      eello_turn4=0.0d0
      ind=0
      do i=1,nres
        num_cont_hb(i)=0
      enddo
cd      print '(a)','Enter EELEC'
cd      write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
      do i=1,nres
        gel_loc_loc(i)=0.0d0
        gcorr_loc(i)=0.0d0
      enddo
      do i=iatel_s,iatel_e
        dxi=dc(1,i)
        dyi=dc(2,i)
        dzi=dc(3,i)
        dx_normi=dc_norm(1,i)
        dy_normi=dc_norm(2,i)
        dz_normi=dc_norm(3,i)
        xmedi=c(1,i)+0.5d0*dxi
        ymedi=c(2,i)+0.5d0*dyi
        zmedi=c(3,i)+0.5d0*dzi
        num_conti=0
c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
        do j=ielstart(i),ielend(i)
          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) 
C Diagnostics only!!!
c         aaa=0.0D0
c         bbb=0.0D0
c         ael6i=0.0D0
c         ael3i=0.0D0
C End diagnostics
          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-xmedi
          yj=c(2,j)+0.5D0*dyj-ymedi
          zj=c(3,j)+0.5D0*dzj-zmedi
          rij=xj*xj+yj*yj+zj*zj
          rrmij=1.0D0/rij
          rij=dsqrt(rij)
          rmij=1.0D0/rij
          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
          el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
          el2=fac4*fac       
          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
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)') 'evdw1',i,j,evdwij
              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)
          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
          do k=1,3
            ghalf=0.5D0*ggg(k)
            gelc(k,i)=gelc(k,i)+ghalf
            gelc(k,j)=gelc(k,j)+ghalf
          enddo
*
* Loop over residues i+1 thru j-1.
*
caug8          do k=i+1,j-1
caug8            do l=1,3
caug8              gelc(l,k)=gelc(l,k)+ggg(l)
caug8            enddo
caug8          enddo
          ggg(1)=facvdw*xj
          ggg(2)=facvdw*yj
          ggg(3)=facvdw*zj
          do k=1,3
            ghalf=0.5D0*ggg(k)
            gvdwpp(k,i)=gvdwpp(k,i)+ghalf
            gvdwpp(k,j)=gvdwpp(k,j)+ghalf
          enddo
*
* Loop over residues i+1 thru j-1.
*
caug8          do k=i+1,j-1
caug8            do l=1,3
caug8              gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
caug8            enddo
caug8          enddo
#else
          facvdw=ev1+evdwij 
          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
          do k=1,3
            ghalf=0.5D0*ggg(k)
            gelc(k,i)=gelc(k,i)+ghalf
            gelc(k,j)=gelc(k,j)+ghalf
          enddo
*
* Loop over residues i+1 thru j-1.
*
caug8          do k=i+1,j-1
caug8            do l=1,3
caug8              gelc(l,k)=gelc(l,k)+ggg(l)
caug8            enddo
caug8          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) 
          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)
            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)
          enddo
caug8          do k=i+1,j-1
caug8            do l=1,3
caug8              gelc(l,k)=gelc(l,k)+ggg(l)
caug8            enddo
caug8          enddo

        enddo ! j
        num_cont_hb(i)=num_conti
      enddo   ! i
c      write (iout,*) "Number of loop steps in EELEC:",ind
cd      do i=1,nres
cd        write (iout,'(i3,3f10.5,5x,3f10.5)') 
cd     &     i,(gel_loc(k,i),k=1,3),gel_loc_loc(i)
cd      enddo
c 12/7/99 Adam eello_turn3 will be considered as a separate energy term
ccc      eel_loc=eel_loc+eello_turn3
      return
      end