1 subroutine contfunc(cscore,itypi,itypj)
3 C This subroutine calculates the contact function based on
4 C the Gay-Berne potential of interaction.
8 include 'COMMON.CONTPAR'
12 double precision cscore
13 double precision sig0ij,rrij,sig,rij_shift,evdw,e2
15 sig0ij=sig_comp(itypi,itypj)
16 chi1=chi_comp(itypi,itypj)
17 chi2=chi_comp(itypj,itypi)
19 chip1=chip_comp(itypi,itypj)
20 chip2=chip_comp(itypj,itypi)
22 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
24 C Calculate angle-dependent terms of the contact function
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
32 c print *,'egb',itypi,itypj,chi1,chi2,chip1,chip2,
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)
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.
51 facp=1.0D0-om12*chipom12
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
57 sig=sig0ij*dsqrt(sigsq)
58 rij_shift=1.0D0/rij-sig+sig0ij
59 if (rij_shift.le.0.0D0) then
61 cscore = -dlog(evdw+1.0d-6)
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)
71 c------------------------------------------------------------------------------
72 subroutine scdist(cscore,itypi,itypj)
74 C This subroutine calculates the contact distance
78 include 'COMMON.CONTPAR'
81 double precision cscore
84 chi1=chi_comp(itypi,itypj)
85 chi2=chi_comp(itypj,itypi)
87 rrij=xj*xj+yj*yj+zj*zj
89 C Calculate angle-dependent terms of the contact function
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
100 cscore=dsqrt(rrij+chi1**2+chi2**2+2*rij*(chiom2-chiom1)-2*chiom12)