Merge branch 'lipid' of mmka.chem.univ.gda.pl:unres into lipid

[unres.git] / source / wham / src-NEWSC-NEWCORR / contfunc.f

diff --git a/source/wham/src-NEWSC-NEWCORR/contfunc.f b/source/wham/src-NEWSC-NEWCORR/contfunc.f
new file mode 100644 (file)
index 0000000..7aed575
--- /dev/null
+++ b/source/wham/src-NEWSC-NEWCORR/contfunc.f
@@ -0,0 +1,96 @@
+      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