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