Page 396 - Handbook Of Multiphase Flow Assurance

Computer code (Makogon, 1994, 1997) 395

xx=xn(natom2)
yy=yn(natom2)
zz=zn(natom2)
rsq=(x-xx)**2+(y-yy)**2+(z-zz)**2
C a quick test for close contacts (1.5 A)
if (rsq.lt.2.25)then
Ener=1.D8
return
endif
r=dsqrt(rsq)
Ehb = 0.D0
Elj = 0.D0
C Coulombic interaction
Ecou = Co*atomdata(natom1,4)*atomdata(natom2,4)/R
C Lennard-Jones interaction
IF (R .LE. Rcut) THEN
Itype=iatomtype(natom2,1)
R6red=(RR(Itype,natom1)*RR(Itype,natom1)/Rsq)**3
Elj = DD(Itype,natom1)*R6red*(R6red-2.D0)
C intramolecular hydrogen-bond interaction, if applicable
IF (Rsq.LT.Rsqoff) then
ndon=0
IF(rhb(natom1).eq.2.8D0.and.rhb(natom2).eq.0.01D0) then
ndon=natom1
nacc=natom2
endif
IF(rhb(natom2).eq.2.8D0.and.rhb(natom1).eq.0.01D0) then
ndon=natom2
nacc=natom1
endif
IF(ndon.ne.0) then
C ... check if within switching distance
Swtchr = 1.D0
IF (Rsq .GT. Rsqon) CALL SWITCH(Rsqon,Rsqoff,Rsq,Swtchr)
C ... get hydrogen atom vector
Xh = xn(ibonddata(nacc,1)) - xn(nacc)
Yh = yn(ibonddata(nacc,1)) - yn(nacc)
Zh = zn(ibonddata(nacc,1)) - zn(nacc)
Rhsq = Xh*Xh + Yh*Yh + Zh*Zh
C ... get donor-acceptor distance
Xah= xn(ndon) - xn(nacc)
Yah= yn(ndon) - yn(nacc)
Zah= zn(ndon) - zn(nacc)
Rahsq = Xah*Xah + Yah*Yah + Zah*Zah
C ... cos**2 of HB angle