subroutine elecont(lprint,ncont,icont,ist,ien,ipermmin)
      implicit none
      include 'DIMENSIONS'
      include 'DIMENSIONS.ZSCOPT'
      include 'DIMENSIONS.COMPAR'
      include 'COMMON.IOUNITS'
      include 'COMMON.CHAIN'
      include 'COMMON.INTERACT'
      include 'COMMON.FFIELD'
      include 'COMMON.NAMES'
      include 'COMMON.LOCAL'
      logical lprint
      integer iperm,ipermmin,ii,jj
      integer i,j,k,ist,ien,iteli,itelj,ind,i1,i2,it1,it2,ic1,ic2
      double precision rri,xi,yi,zi,dxi,dyi,dzi,xmedi,ymedi,zmedi,
     &  xj,yj,zj,dxj,dyj,dzj,aaa,bbb,ael6i,ael3i,rrmij,rmij,r3ij,r6ij,
     &  vrmij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4,evdwij,el1,el2,
     &  eesij,ees,evdw,ene, rij,zj_temp,xj_temp,yj_temp,
     & sscale,sscagrad,dist_temp,xj_safe,yj_safe,zj_safe,dist_init
      double precision elpp6c(2,2),elpp3c(2,2),ael6c(2,2),ael3c(2,2),
     &  appc(2,2),bppc(2,2)
      double precision elcutoff,elecutoff_14
      integer ncont,icont(2,maxcont),xshift,yshift,zshift,isubchap
      double precision econt(maxcont)
*
* Load the constants of peptide bond - peptide bond interactions.
* Type 1 - ordinary peptide bond, type 2 - alkylated peptide bond (e.g.
* proline) - determined by averaging ECEPP energy.      
*
* as of 7/06/91.
*
c      data epp    / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/
c      data rpp    / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
      data elpp6c  /-0.2379d0,-0.2056d0,-0.2056d0,-0.0610d0/
      data elpp3c  / 0.0503d0, 0.0000d0, 0.0000d0, 0.0692d0/
      data elcutoff /-0.3d0/,elecutoff_14 /-0.5d0/
      ees=0.0d0
      evdw=0.0d0
      if (lprint) write (iout,'(a)') 
     &  "Constants of electrostatic interaction energy expression."
      do i=1,2
        do j=1,2
        rri=rpp(i,j)**6
        appc(i,j)=epp(i,j)*rri*rri 
        bppc(i,j)=-2.0*epp(i,j)*rri
        ael6c(i,j)=elpp6c(i,j)*4.2**6
        ael3c(i,j)=elpp3c(i,j)*4.2**3
        if (lprint)
     &  write (iout,'(2i2,4e15.4)') i,j,appc(i,j),bppc(i,j),ael6c(i,j),
     &                               ael3c(i,j)
        enddo
      enddo
      ncont=0
      do 1 i=ist,ien-2
        ii=iperm(i,ipermmin)
        xi=c(1,ii)
        yi=c(2,ii)
        zi=c(3,ii)
        dxi=c(1,ii+1)-c(1,ii)
        dyi=c(2,ii+1)-c(2,ii)
        dzi=c(3,ii+1)-c(3,ii)
        xmedi=xi+0.5*dxi
        ymedi=yi+0.5*dyi
        zmedi=zi+0.5*dzi
          xmedi=mod(xmedi,boxxsize)
          if (xmedi.lt.0) xmedi=xmedi+boxxsize
          ymedi=mod(ymedi,boxysize)
          if (ymedi.lt.0) ymedi=ymedi+boxysize
          zmedi=mod(zmedi,boxzsize)
          if (zmedi.lt.0) zmedi=zmedi+boxzsize
        do 4 j=i+2,ien-1
          jj=iperm(j,ipermmin)
          ind=ind+1
          iteli=itel(i)
          itelj=itel(j)
          if (j.eq.i+2 .and. itelj.eq.2) iteli=2
          if (iteli.eq.2 .and. itelj.eq.2 
     &      .or.iteli.eq.0 .or.itelj.eq.0) goto 4
          aaa=appc(iteli,itelj)
          bbb=bppc(iteli,itelj)
          ael6i=ael6c(iteli,itelj)
          ael3i=ael3c(iteli,itelj) 
          dxj=c(1,jj+1)-c(1,jj)
          dyj=c(2,jj+1)-c(2,jj)
          dzj=c(3,jj+1)-c(3,jj)
          xj=c(1,jj)+0.5*dxj
          yj=c(2,jj)+0.5*dyj
          zj=c(3,jj)+0.5*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
            sss=sscale(sqrt(rij))
            sssgrad=sscagrad(sqrt(rij))
          rrmij=1.0/(xj*xj+yj*yj+zj*zj)
          rmij=sqrt(rrmij)
          r3ij=rrmij*rmij
          r6ij=r3ij*r3ij  
          vrmij=vblinv*rmij
          cosa=(dxi*dxj+dyi*dyj+dzi*dzj)*vblinv2      
          cosb=(xj*dxi+yj*dyi+zj*dzi)*vrmij
          cosg=(xj*dxj+yj*dyj+zj*dzj)*vrmij
          fac=cosa-3.0*cosb*cosg
          ev1=aaa*r6ij*r6ij
          ev2=bbb*r6ij
          fac3=ael6i*r6ij
          fac4=ael3i*r3ij
          evdwij=ev1+ev2
          el1=fac3*(4.0+fac*fac-3.0*(cosb*cosb+cosg*cosg))
          el2=fac4*fac       
          eesij=el1+el2
          if (j.gt.i+2 .and. eesij.le.elcutoff .or.
     &        j.eq.i+2 .and. eesij.le.elecutoff_14) then
             ncont=ncont+1
             icont(1,ncont)=i
             icont(2,ncont)=j
	     econt(ncont)=eesij
          endif
          ees=ees+eesij
          evdw=evdw+evdwij*sss
    4   continue
    1 continue
      if (lprint) then
        write (iout,*) 'Total average electrostatic energy: ',ees
        write (iout,*) 'VDW energy between peptide-group centers: ',evdw
        write (iout,*)
        write (iout,*) 'Electrostatic contacts before pruning: '
        do i=1,ncont
          i1=icont(1,i)
          i2=icont(2,i)
          it1=itype(i1)
          it2=itype(i2)
          write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)')
     &     i,restyp(it1),i1,restyp(it2),i2,econt(i)
        enddo
      endif
c For given residues keep only the contacts with the greatest energy.
      i=0
      do while (i.lt.ncont)
        i=i+1
        ene=econt(i)
        ic1=icont(1,i)
        ic2=icont(2,i)
        j=i
        do while (j.lt.ncont)
          j=j+1
          if (ic1.eq.icont(1,j).and.iabs(icont(2,j)-ic2).le.2 .or.
     &        ic2.eq.icont(2,j).and.iabs(icont(1,j)-ic1).le.2) then
c            write (iout,*) "i",i," j",j," ic1",ic1," ic2",ic2,
c     &       " jc1",icont(1,j)," jc2",icont(2,j)," ncont",ncont
            if (econt(j).lt.ene .and. icont(2,j).ne.icont(1,j)+2) then
              if (ic1.eq.icont(1,j)) then
                do k=1,ncont
                  if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.icont(2,j)
     &               .and. iabs(icont(1,k)-ic1).le.2 .and. 
     &               econt(k).lt.econt(j) ) goto 21 
                enddo
              else if (ic2.eq.icont(2,j) ) then
                do k=1,ncont
                  if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.icont(1,j)
     &               .and. iabs(icont(2,k)-ic2).le.2 .and. 
     &               econt(k).lt.econt(j) ) goto 21 
                enddo
              endif
c Remove ith contact
              do k=i+1,ncont
                icont(1,k-1)=icont(1,k)
                icont(2,k-1)=icont(2,k)
                econt(k-1)=econt(k) 
              enddo
              i=i-1
              ncont=ncont-1
c              write (iout,*) "ncont",ncont
c              do k=1,ncont
c                write (iout,*) icont(1,k),icont(2,k)
c              enddo
              goto 20
            else if (econt(j).gt.ene .and. ic2.ne.ic1+2) 
     &      then
              if (ic1.eq.icont(1,j)) then
                do k=1,ncont
                  if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.ic2
     &               .and. iabs(icont(1,k)-icont(1,j)).le.2 .and. 
     &               econt(k).lt.econt(i) ) goto 21 
                enddo
              else if (ic2.eq.icont(2,j) ) then
                do k=1,ncont
                  if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.ic1
     &               .and. iabs(icont(2,k)-icont(2,j)).le.2 .and. 
     &               econt(k).lt.econt(i) ) goto 21 
                enddo
              endif
c Remove jth contact
              do k=j+1,ncont
                icont(1,k-1)=icont(1,k)
                icont(2,k-1)=icont(2,k)
                econt(k-1)=econt(k) 
              enddo
              ncont=ncont-1
c              write (iout,*) "ncont",ncont
c              do k=1,ncont
c                write (iout,*) icont(1,k),icont(2,k)
c              enddo
              j=j-1
            endif   
          endif
   21     continue
        enddo
   20   continue
      enddo
      if (lprint) then
        write (iout,*)
        write (iout,*) 'Electrostatic contacts after pruning: '
        do i=1,ncont
          i1=icont(1,i)
          i2=icont(2,i)
          it1=itype(i1)
          it2=itype(i2)
          write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)')
     &     i,restyp(it1),i1,restyp(it2),i2,econt(i)
        enddo
      endif
      return
      end