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