+ 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=itype(i)
+ do j=i+kkk,ien
+ itj=itype(j)
+ 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(it1,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