subroutine contact(lprint,ncont,icont,ist,ien)
      implicit none
      include 'DIMENSIONS'
      include 'DIMENSIONS.ZSCOPT'
      include 'DIMENSIONS.FREE'
      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