source/wham/src/contact.f

   1       subroutine contact(lprint,ncont,icont,ist,ien)
   2       implicit none
   3       include 'DIMENSIONS'
   4       include 'DIMENSIONS.ZSCOPT'
   5       include 'COMMON.CONTROL'
   6       include 'COMMON.IOUNITS'
   7       include 'COMMON.CHAIN'
   8       include 'COMMON.INTERACT'
   9       include 'COMMON.FFIELD'
  10       include 'COMMON.NAMES'
  11       include 'COMMON.CALC'
  12       include 'COMMON.CONTPAR'
  13       include 'COMMON.LOCAL'
  14       integer ist,ien,kkk,iti,itj,itypi,itypj,i1,i2,it1,it2
  15       real*8 csc,dist
  16       real*8 cscore(maxcont),omt1(maxcont),omt2(maxcont),omt12(maxcont),
  17      & ddsc(maxcont),ddla(maxcont),ddlb(maxcont)
  18       integer ncont,icont(2,maxcont)
  19       real*8 u,v,a(3),b(3),dla,dlb
  20       logical lprint
  21       ncont=0
  22       kkk=3
  23       if (lprint) then
  24       do i=1,nres
  25         write (iout,110) restyp(itype(i)),i,c(1,i),c(2,i),
  26      &    c(3,i),dc(1,nres+i),dc(2,nres+i),dc(3,nres+i),
  27      &    dc_norm(1,nres+i),dc_norm(2,nres+i),dc_norm(3,nres+i)
  28       enddo
  29       endif
  30   110 format (a,'(',i3,')',9f8.3)
  31       do i=ist,ien-kkk
  32         iti=itype(i)
  33         do j=i+kkk,ien
  34           itj=itype(j)
  35           itypi=iti
  36           itypj=itj
  37           xj = c(1,nres+j)-c(1,nres+i)
  38           yj = c(2,nres+j)-c(2,nres+i)
  39           zj = c(3,nres+j)-c(3,nres+i)
  40           dxi = dc_norm(1,nres+i)
  41           dyi = dc_norm(2,nres+i)
  42           dzi = dc_norm(3,nres+i)
  43           dxj = dc_norm(1,nres+j)
  44           dyj = dc_norm(2,nres+j)
  45           dzj = dc_norm(3,nres+j)
  46           do k=1,3
  47             a(k)=dc(k,nres+i)
  48             b(k)=dc(k,nres+j)
  49           enddo
  50 c          write (iout,*) (a(k),k=1,3),(b(k),k=1,3)
  51           if (icomparfunc.eq.1) then
  52             call contfunc(csc,iti,itj)
  53           else if (icomparfunc.eq.2) then
  54             call scdist(csc,iti,itj)
  55           else if (icomparfunc.eq.3 .or. icomparfunc.eq.5) then
  56             csc = dist(nres+i,nres+j)
  57           else if (icomparfunc.eq.4) then
  58             call odlodc(c(1,i),c(1,j),a,b,u,v,dla,dlb,csc)
  59           else
  60             write (*,*) "Error - Unknown sidechain contact function"
  61             write (iout,*) "Error - Unknown sidechain contact function"
  62           endif
  63           if (csc.lt.sc_cutoff(iti,itj)) then
  64 c            write(iout,*) "i",i," j",j," dla",dla,dsc(iti),
  65 c     &      " dlb",dlb,dsc(itj)," csc",csc,sc_cutoff(iti,itj),
  66 c     &      dxi,dyi,dzi,dxi**2+dyi**2+dzi**2,
  67 c     &      dxj,dyj,dzj,dxj**2+dyj**2+dzj**2,om1,om2,om12,
  68 c     &      xj,yj,zj
  69 c            write(iout,*)'egb',itypi,itypj,chi1,chi2,chip1,chip2,
  70 c     &       sig0ij,rij,rrij,om1,om2,om12,chiom1,chiom2,chiom12,
  71 c     &       chipom1,chipom2,chipom12,sig,eps2rt,rij_shift,e2,evdw,
  72 c     &       csc
  73             ncont=ncont+1
  74             cscore(ncont)=csc
  75             icont(1,ncont)=i
  76             icont(2,ncont)=j
  77             omt1(ncont)=om1
  78             omt2(ncont)=om2
  79             omt12(ncont)=om12
  80             ddsc(ncont)=1.0d0/rij
  81             ddla(ncont)=dla
  82             ddlb(ncont)=dlb
  83           endif
  84         enddo
  85       enddo
  86       if (lprint) then
  87         write (iout,'(a)') 'Contact map:'
  88         do i=1,ncont
  89           i1=icont(1,i)
  90           i2=icont(2,i)
  91           it1=itype(i1)
  92           it2=itype(i2)
  93           write (iout,'(i3,2x,a,i4,2x,a,i4,5f8.3,3f10.5)')
  94      &     i,restyp(it1),i1,restyp(it2),i2,cscore(i),
  95      &     sc_cutoff(it1,it2),ddsc(i),ddla(i),ddlb(i),
  96      &     omt1(i),omt2(i),omt12(i)
  97         enddo
  98       endif
  99       return
 100       end
 101 c----------------------------------------------------------------------------
 102       double precision function contact_fract(ncont,ncont_ref,
 103      &                                     icont,icont_ref)
 104       implicit none
 105       include 'DIMENSIONS'
 106       include 'COMMON.IOUNITS'
 107       integer i,j,nmatch
 108       integer ncont,ncont_ref,icont(2,maxcont),icont_ref(2,maxcont)
 109       nmatch=0
 110 c     print *,'ncont=',ncont,' ncont_ref=',ncont_ref
 111 c     write (iout,'(20i4)') (icont_ref(1,i),i=1,ncont_ref)
 112 c     write (iout,'(20i4)') (icont_ref(2,i),i=1,ncont_ref)
 113 c     write (iout,'(20i4)') (icont(1,i),i=1,ncont)
 114 c     write (iout,'(20i4)') (icont(2,i),i=1,ncont)
 115       do i=1,ncont
 116         do j=1,ncont_ref
 117           if (icont(1,i).eq.icont_ref(1,j) .and.
 118      &        icont(2,i).eq.icont_ref(2,j)) nmatch=nmatch+1
 119         enddo
 120       enddo
 121 c     print *,' nmatch=',nmatch
 122 c     contact_fract=dfloat(nmatch)/dfloat(max0(ncont,ncont_ref))
 123       contact_fract=dfloat(nmatch)/dfloat(ncont_ref)
 124       return
 125       end
 126 c------------------------------------------------------------------------------
 127       subroutine pept_cont(lprint,ncont,icont)
 128       implicit none
 129       include 'DIMENSIONS'
 130       include 'DIMENSIONS.ZSCOPT'
 131       include 'COMMON.IOUNITS'
 132       include 'COMMON.CHAIN'
 133       include 'COMMON.INTERACT'
 134       include 'COMMON.FFIELD'
 135       include 'COMMON.NAMES'
 136       integer ncont,icont(2,maxcont)
 137       integer i,j,k,kkk,i1,i2,it1,it2
 138       logical lprint
 139       real*8 dist
 140       real*8 rcomp /5.5d0/
 141       ncont=0
 142       kkk=0
 143       print *,'Entering pept_cont: nnt=',nnt,' nct=',nct
 144       do i=nnt,nct-3
 145         do k=1,3
 146           c(k,2*nres+1)=0.5d0*(c(k,i)+c(k,i+1))
 147         enddo
 148         do j=i+2,nct-1
 149           do k=1,3
 150             c(k,2*nres+2)=0.5d0*(c(k,j)+c(k,j+1))
 151           enddo
 152           if (dist(2*nres+1,2*nres+2).lt.rcomp) then
 153             ncont=ncont+1
 154             icont(1,ncont)=i
 155             icont(2,ncont)=j
 156           endif
 157         enddo
 158       enddo
 159       if (lprint) then
 160         write (iout,'(a)') 'PP contact map:'
 161         do i=1,ncont
 162           i1=icont(1,i)
 163           i2=icont(2,i)
 164           it1=itype(i1)
 165           it2=itype(i2)
 166           write (iout,'(i3,2x,a,i4,2x,a,i4)')
 167      &     i,restyp(it1),i1,restyp(it2),i2
 168         enddo
 169       endif
 170       return
 171       end