Page 389 - Handbook Of Multiphase Flow Assurance

388 10. Research methods in flow assurance

natom2 = monodata(ii,jj)
x =xn(natom1)
y =yn(natom1)
z =zn(natom1)
xx=xn(natom2)
yy=yn(natom2)
zz=zn(natom2)
rsq=(x-xx)**2+(y-yy)**2+(z-zz)**2
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 hydrogen-bond interaction, if applicable
ndon=0
IF(rhb(natom1).eq.2.8D0.and.rhb(natom2).eq.0.01D0) then
ndon=monodata(i,j)
nacc=monodata(ii,jj)
endif
IF(rhb(natom2).eq.2.8D0.and.rhb(natom1).eq.0.01D0) then
ndon=monodata(ii,jj)
nacc=monodata(i,j)
endif
IF(ndon.gt.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
Chbsq = (Xh*Xah + Yh*Yah + Zh*Zah)**2/(Rhsq*Rahsq)
IF (Chbsq . LT. Casqof) goto 10