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