subroutine contact(lprint,ncont,icont,ist,ien)
      implicit none
      include 'DIMENSIONS'
      include 'DIMENSIONS.ZSCOPT'
      include 'COMMON.CONTROL'
      include 'COMMON.IOUNITS'
      include 'COMMON.CHAIN'
      include 'COMMON.INTERACT'
      include 'COMMON.FFIELD'
      include 'COMMON.NAMES'
      include 'COMMON.CALC'
      include 'COMMON.CONTPAR'
      include 'COMMON.LOCAL'
      integer ist,ien,kkk,iti,itj,itypi,itypj,i1,i2,it1,it2
      real*8 csc,dist
      real*8 cscore(maxcont),omt1(maxcont),omt2(maxcont),omt12(maxcont),
     & ddsc(maxcont),ddla(maxcont),ddlb(maxcont)
      integer ncont,icont(2,maxcont)
      real*8 u,v,a(3),b(3),dla,dlb
      logical lprint
      ncont=0
      kkk=3
      if (lprint) then
      do i=1,nres
        write (iout,110) restyp(itype(i)),i,c(1,i),c(2,i),
     &    c(3,i),dc(1,nres+i),dc(2,nres+i),dc(3,nres+i),
     &    dc_norm(1,nres+i),dc_norm(2,nres+i),dc_norm(3,nres+i)
      enddo
      endif
  110 format (a,'(',i3,')',9f8.3)
      do i=ist,ien-kkk
        iti=iabs(itype(i))
        if (iti.le.0 .or. iti.gt.ntyp) cycle
        do j=i+kkk,ien
          itj=iabs(itype(j))
          if (itj.le.0 .or. itj.gt.ntyp) cycle
          itypi=iti
          itypj=itj
          xj = c(1,nres+j)-c(1,nres+i)    
          yj = c(2,nres+j)-c(2,nres+i)    
          zj = c(3,nres+j)-c(3,nres+i)    
          dxi = dc_norm(1,nres+i)
          dyi = dc_norm(2,nres+i)
          dzi = dc_norm(3,nres+i)
          dxj = dc_norm(1,nres+j)
          dyj = dc_norm(2,nres+j)
          dzj = dc_norm(3,nres+j)
          do k=1,3
            a(k)=dc(k,nres+i)
            b(k)=dc(k,nres+j)
          enddo
c          write (iout,*) (a(k),k=1,3),(b(k),k=1,3)
          if (icomparfunc.eq.1) then
            call contfunc(csc,iti,itj)
          else if (icomparfunc.eq.2) then
            call scdist(csc,iti,itj)
          else if (icomparfunc.eq.3 .or. icomparfunc.eq.5) then
            csc = dist(nres+i,nres+j)
          else if (icomparfunc.eq.4) then
            call odlodc(c(1,i),c(1,j),a,b,u,v,dla,dlb,csc)
          else
            write (*,*) "Error - Unknown sidechain contact function"
            write (iout,*) "Error - Unknown sidechain contact function"
          endif
          if (csc.lt.sc_cutoff(iti,itj)) then
c            write(iout,*) "i",i," j",j," dla",dla,dsc(iti),
c     &      " dlb",dlb,dsc(itj)," csc",csc,sc_cutoff(iti,itj),
c     &      dxi,dyi,dzi,dxi**2+dyi**2+dzi**2,
c     &      dxj,dyj,dzj,dxj**2+dyj**2+dzj**2,om1,om2,om12,
c     &      xj,yj,zj
c            write(iout,*)'egb',itypi,itypj,chi1,chi2,chip1,chip2,
c     &       sig0ij,rij,rrij,om1,om2,om12,chiom1,chiom2,chiom12,
c     &       chipom1,chipom2,chipom12,sig,eps2rt,rij_shift,e2,evdw,
c     &       csc
            ncont=ncont+1
            cscore(ncont)=csc
            icont(1,ncont)=i
            icont(2,ncont)=j
            omt1(ncont)=om1
            omt2(ncont)=om2
            omt12(ncont)=om12
            ddsc(ncont)=1.0d0/rij
            ddla(ncont)=dla
            ddlb(ncont)=dlb
          endif
        enddo
      enddo
      if (lprint) then
        write (iout,'(a)') 'Contact map:'
        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,5f8.3,3f10.5)') 
     &     i,restyp(it1),i1,restyp(it2),i2,cscore(i),
     &     sc_cutoff(iabs(it1),iabs(it2)),ddsc(i),ddla(i),ddlb(i),
     &     omt1(i),omt2(i),omt12(i)
        enddo
      endif
      return
      end
c----------------------------------------------------------------------------
      double precision function contact_fract(ncont,ncont_ref,
     &                                     icont,icont_ref)
      implicit none
      include 'DIMENSIONS'
      include 'COMMON.IOUNITS'
      integer i,j,nmatch
      integer ncont,ncont_ref,icont(2,maxcont),icont_ref(2,maxcont)
      nmatch=0
c     print *,'ncont=',ncont,' ncont_ref=',ncont_ref 
c     write (iout,'(20i4)') (icont_ref(1,i),i=1,ncont_ref)
c     write (iout,'(20i4)') (icont_ref(2,i),i=1,ncont_ref)
c     write (iout,'(20i4)') (icont(1,i),i=1,ncont)
c     write (iout,'(20i4)') (icont(2,i),i=1,ncont)
      do i=1,ncont
        do j=1,ncont_ref
          if (icont(1,i).eq.icont_ref(1,j) .and. 
     &        icont(2,i).eq.icont_ref(2,j)) nmatch=nmatch+1
        enddo
      enddo
c     print *,' nmatch=',nmatch
c     contact_fract=dfloat(nmatch)/dfloat(max0(ncont,ncont_ref))
      contact_fract=dfloat(nmatch)/dfloat(ncont_ref)
      return
      end
c------------------------------------------------------------------------------
      subroutine pept_cont(lprint,ncont,icont)
      implicit none
      include 'DIMENSIONS'
      include 'DIMENSIONS.ZSCOPT'
      include 'COMMON.IOUNITS'
      include 'COMMON.CHAIN'
      include 'COMMON.INTERACT'
      include 'COMMON.FFIELD'
      include 'COMMON.NAMES'
      integer ncont,icont(2,maxcont)
      integer i,j,k,kkk,i1,i2,it1,it2
      logical lprint
      real*8 dist
      real*8 rcomp /5.5d0/
      ncont=0
      kkk=0
      print *,'Entering pept_cont: nnt=',nnt,' nct=',nct
      do i=nnt,nct-3
        do k=1,3
          c(k,2*nres+1)=0.5d0*(c(k,i)+c(k,i+1))
        enddo
        do j=i+2,nct-1
          do k=1,3
            c(k,2*nres+2)=0.5d0*(c(k,j)+c(k,j+1))
          enddo
          if (dist(2*nres+1,2*nres+2).lt.rcomp) then
            ncont=ncont+1
            icont(1,ncont)=i
            icont(2,ncont)=j
          endif
        enddo
      enddo
      if (lprint) then
        write (iout,'(a)') 'PP contact map:'
        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)')
     &     i,restyp(it1),i1,restyp(it2),i2
        enddo
      endif
      return
      end