--- /dev/null
+ subroutine contfunc(cscore,itypi,itypj)
+C
+C This subroutine calculates the contact function based on
+C the Gay-Berne potential of interaction.
+C
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CONTPAR'
+ include 'COMMON.CALC'
+ integer expon /6/
+C
+ sig0ij=sig_comp(itypi,itypj)
+ chi1=chi_comp(itypi,itypj)
+ chi2=chi_comp(itypj,itypi)
+ chi12=chi1*chi2
+ chip1=chip_comp(itypi,itypj)
+ chip2=chip_comp(itypj,itypi)
+ chip12=chip1*chip2
+ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+ rij=dsqrt(rrij)
+C Calculate angle-dependent terms of the contact function
+ erij(1)=xj*rij
+ erij(2)=yj*rij
+ erij(3)=zj*rij
+ om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
+ om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
+ om12=dxi*dxj+dyi*dyj+dzi*dzj
+ chiom12=chi12*om12
+c print *,'egb',itypi,itypj,chi1,chi2,chip1,chip2,
+c & sig0ij,
+c & rij,rrij,om1,om2,om12
+C Calculate eps1(om12)
+ faceps1=1.0D0-om12*chiom12
+ faceps1_inv=1.0D0/faceps1
+ eps1=dsqrt(faceps1_inv)
+C Following variable is eps1*deps1/dom12
+ eps1_om12=faceps1_inv*chiom12
+C Calculate sigma(om1,om2,om12)
+ om1om2=om1*om2
+ chiom1=chi1*om1
+ chiom2=chi2*om2
+ facsig=om1*chiom1+om2*chiom2-2.0D0*om1om2*chiom12
+ sigsq=1.0D0-facsig*faceps1_inv
+C Calculate eps2 and its derivatives in om1, om2, and om12.
+ chipom1=chip1*om1
+ chipom2=chip2*om2
+ chipom12=chip12*om12
+ facp=1.0D0-om12*chipom12
+ facp_inv=1.0D0/facp
+ facp1=om1*chipom1+om2*chipom2-2.0D0*om1om2*chipom12
+C Following variable is the square root of eps2
+ eps2rt=1.0D0-facp1*facp_inv
+ sigsq=1.0D0/sigsq
+ sig=sig0ij*dsqrt(sigsq)
+ rij_shift=1.0D0/rij-sig+sig0ij
+ if (rij_shift.le.0.0D0) then
+ evdw=1.0D1
+ cscore = -dlog(evdw+1.0d-6)
+ return
+ endif
+ rij_shift=1.0D0/rij_shift
+ e2=(rij_shift*sig0ij)**expon
+ evdw=dabs(eps1*eps2rt**2*e2)
+ if (evdw.gt.1.0d1) evdw = 1.0d1
+ cscore = -dlog(evdw+1.0d-6)
+ return
+ end
+c------------------------------------------------------------------------------
+ subroutine scdist(cscore,itypi,itypj)
+C
+C This subroutine calculates the contact distance
+C
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CONTPAR'
+ include 'COMMON.CALC'
+C
+ chi1=chi_comp(itypi,itypj)
+ chi2=chi_comp(itypj,itypi)
+ chi12=chi1*chi2
+ rrij=xj*xj+yj*yj+zj*zj
+ rij=dsqrt(rrij)
+C Calculate angle-dependent terms of the contact function
+ erij(1)=xj/rij
+ erij(2)=yj/rij
+ erij(3)=zj/rij
+ om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
+ om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
+ om12=dxi*dxj+dyi*dyj+dzi*dzj
+ chiom12=chi12*om12
+ om1om2=om1*om2
+ chiom1=chi1*om1
+ chiom2=chi2*om2
+ cscore=dsqrt(rrij+chi1**2+chi2**2+2*rij*(chiom2-chiom1)-2*chiom12)
+ return
+ end