Merge branch 'devel' into AFM
[unres.git] / source / wham / src-NEWSC-NEWCORR / contact.f
diff --git a/source/wham/src-NEWSC-NEWCORR/contact.f b/source/wham/src-NEWSC-NEWCORR/contact.f
new file mode 100644 (file)
index 0000000..5b05d57
--- /dev/null
@@ -0,0 +1,171 @@
+      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