small fix in .gitingnore
[unres.git] / 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