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