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