source/maxlik/src_MD_T_maxlik-NEWCORR-PMF-PDB/contfunc.f

   1       subroutine contfunc(cscore,itypi,itypj)
   2 C
   3 C This subroutine calculates the contact function based on
   4 C the Gay-Berne potential of interaction.
   5 C
   6       implicit none
   7       include 'DIMENSIONS'
   8       include 'COMMON.CONTPAR'
   9       include 'COMMON.CALC'
  10       integer expon /6/
  11       integer itypi,itypj
  12       double precision cscore
  13       double precision sig0ij,rrij,sig,rij_shift,evdw,e2
  14 C
  15       sig0ij=sig_comp(itypi,itypj)
  16       chi1=chi_comp(itypi,itypj)
  17       chi2=chi_comp(itypj,itypi)
  18       chi12=chi1*chi2
  19       chip1=chip_comp(itypi,itypj)
  20       chip2=chip_comp(itypj,itypi)
  21       chip12=chip1*chip2
  22       rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
  23       rij=dsqrt(rrij)
  24 C Calculate angle-dependent terms of the contact function
  25       erij(1)=xj*rij
  26       erij(2)=yj*rij
  27       erij(3)=zj*rij
  28       om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
  29       om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
  30       om12=dxi*dxj+dyi*dyj+dzi*dzj
  31       chiom12=chi12*om12
  32 c      print *,'egb',itypi,itypj,chi1,chi2,chip1,chip2,
  33 c     &  sig0ij,
  34 c     &  rij,rrij,om1,om2,om12
  35 C Calculate eps1(om12)
  36       faceps1=1.0D0-om12*chiom12
  37       faceps1_inv=1.0D0/faceps1
  38       eps1=dsqrt(faceps1_inv)
  39 C Following variable is eps1*deps1/dom12
  40       eps1_om12=faceps1_inv*chiom12
  41 C Calculate sigma(om1,om2,om12)
  42       om1om2=om1*om2
  43       chiom1=chi1*om1
  44       chiom2=chi2*om2
  45       facsig=om1*chiom1+om2*chiom2-2.0D0*om1om2*chiom12
  46       sigsq=1.0D0-facsig*faceps1_inv
  47 C Calculate eps2 and its derivatives in om1, om2, and om12.
  48       chipom1=chip1*om1
  49       chipom2=chip2*om2
  50       chipom12=chip12*om12
  51       facp=1.0D0-om12*chipom12
  52       facp_inv=1.0D0/facp
  53       facp1=om1*chipom1+om2*chipom2-2.0D0*om1om2*chipom12
  54 C Following variable is the square root of eps2
  55       eps2rt=1.0D0-facp1*facp_inv
  56       sigsq=1.0D0/sigsq
  57       sig=sig0ij*dsqrt(sigsq)
  58       rij_shift=1.0D0/rij-sig+sig0ij
  59       if (rij_shift.le.0.0D0) then
  60         evdw=1.0D1
  61         cscore = -dlog(evdw+1.0d-6)
  62         return
  63       endif
  64       rij_shift=1.0D0/rij_shift
  65       e2=(rij_shift*sig0ij)**expon
  66       evdw=dabs(eps1*eps2rt**2*e2)
  67       if (evdw.gt.1.0d1) evdw = 1.0d1
  68       cscore = -dlog(evdw+1.0d-6)
  69       return
  70       end
  71 c------------------------------------------------------------------------------
  72       subroutine scdist(cscore,itypi,itypj)
  73 C
  74 C This subroutine calculates the contact distance
  75 C
  76       implicit none
  77       include 'DIMENSIONS'
  78       include 'COMMON.CONTPAR'
  79       include 'COMMON.CALC'
  80       integer itypi,itypj
  81       double precision cscore
  82       double precision rrij
  83 C
  84       chi1=chi_comp(itypi,itypj)
  85       chi2=chi_comp(itypj,itypi)
  86       chi12=chi1*chi2
  87       rrij=xj*xj+yj*yj+zj*zj
  88       rij=dsqrt(rrij)
  89 C Calculate angle-dependent terms of the contact function
  90       erij(1)=xj/rij
  91       erij(2)=yj/rij
  92       erij(3)=zj/rij
  93       om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
  94       om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
  95       om12=dxi*dxj+dyi*dyj+dzi*dzj
  96       chiom12=chi12*om12
  97       om1om2=om1*om2
  98       chiom1=chi1*om1
  99       chiom2=chi2*om2
 100       cscore=dsqrt(rrij+chi1**2+chi2**2+2*rij*(chiom2-chiom1)-2*chiom12)
 101       return
 102       end