Adasko's dir

[unres.git] / source / cluster / wham / src-NEWSC / contact.f

diff --git a/source/cluster/wham/src-NEWSC/contact.f b/source/cluster/wham/src-NEWSC/contact.f
new file mode 100644 (file)
index 0000000..b17f153
--- /dev/null
+++ b/source/cluster/wham/src-NEWSC/contact.f
@@ -0,0 +1,69 @@
+      subroutine contact(lprint,ncont,icont)
+      include 'DIMENSIONS'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CHAIN'
+      include 'COMMON.INTERACT'
+      include 'COMMON.FFIELD'
+      include 'COMMON.NAMES'
+      real*8 facont /1.569D0/  ! facont = (2/(1-sqrt(1-1/4)))**(1/6)
+      integer ncont,icont(2,maxcont)
+      logical lprint
+      ncont=0
+      kkk=3
+c     print *,'nnt=',nnt,' nct=',nct
+      do i=nnt+kkk,nct
+        iti=itype(i)
+        do j=nnt,i-kkk
+          itj=itype(j)
+          if (ipot.ne.4) then
+c           rcomp=sigmaii(iti,itj)+1.0D0
+            rcomp=facont*sigmaii(iti,itj)
+          else 
+c           rcomp=sigma(iti,itj)+1.0D0
+            rcomp=facont*sigma(iti,itj)
+          endif
+c         rcomp=6.5D0
+c         print *,'rcomp=',rcomp,' dist=',dist(nres+i,nres+j)
+         if (dist(nres+i,nres+j).lt.rcomp) then
+            ncont=ncont+1
+            icont(1,ncont)=i
+            icont(2,ncont)=j
+          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)') 
+     &     i,restyp(it1),i1,restyp(it2),i2 
+        enddo
+      endif
+      return
+      end
+c----------------------------------------------------------------------------
+      double precision function contact_fract(ncont,ncont_ref,
+     &                                     icont,icont_ref)
+      include 'DIMENSIONS'
+      include 'COMMON.IOUNITS'
+      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