C C $Id: pme_directmc.F,v 1.3 1998/07/16 16:40:17 jjv5 Exp $ C C------------------------------------------------------------------------ subroutine pme_directmc(edir,eself,eclmb,dedir) implicit double precision(a-h,o-z) #include "divcon.dim" #include "divcon.h" parameter (qp=0.3275911d0) parameter (a1=0.254829592d0) parameter (a2=-0.284496736d0) parameter (a3=1.421413741d00) parameter (a4=-1.453152027d0) parameter (a5=1.061405429d00) C C this routine returns the direct energy, the selfenergy and C the short-range (classical) Coulomb energy and gradient. C it will skip all the intramolecular terms to the gradient. C C C Error Function ( erf(x)= 1-erfc(x)) { erfc(x) is obtained from the C polynomial expression}. The term expar2 is the expotential appearing C in the Integral appearing in the exact expression for erf. C { erf(x) = (2/sqrt(pi))*Integ(exp(-t*t)dt)} C The polynomial expression for erfc(x) makes use of the same exponental. C However, the derivative is obtained from the actual expression, not the C polynomial. C C d(erf(ax)/x) = -erf(ax) + (2/sqrt(pi)) {ax*exp(-ax*ax)} C ------------ ---------------------------------------- C dx x*x C C C C parameters: dimension dedir(maxpar) CC----------------------------------------------------------CC call etimer(t1) C units of eclmb, edir and eself are eV eclmb = 0.0 edir = 0.0 eself = 0.0 C unit of gradient (dedir) is kcal/A tokcal = 332.0704 bpisqrt = betapme/pisqrt C intermolecular terms: include gradient do ires=1,nres ib = irpnt(ires) ie = irpnt(ires+1)-1 do iat=ib,ie iat3 = 3*iat-2 xi = xyz(1,iat) yi = xyz(2,iat) zi = xyz(3,iat) chi = atchg(iat) C . selfenergy (a constant, it doesn't contribute to the gradient) eself = eself - bpisqrt*chi*chi do 90 jres=1,ires-1 jb = irpnt(jres) je = irpnt(jres+1)-1 do jat=jb,je jat3 = 3*jat-2 if (pbc) then call pbcxyz(iat,jat,xj,yj,zj) else xj = xyz(1,jat) yj = xyz(2,jat) zj = xyz(3,jat) endif x = xi-xj y = yi-yj z = zi-zj rsqr = x*x + y*y + z*z rr = sqrt(rsqr) rr3 = rr**3 ch = chi*atchg(jat) chk = tokcal*ch C . short range Coulomb energy eclmb = eclmb + ch/rr C . calculate contribution to gradient of short range Coulomb dx = (chk*x)/rr3 dy = (chk*y)/rr3 dz = (chk*z)/rr3 betapmer = betapme*rr expar2 = exp(-betapmer*betapmer) qt = 1.0/(1.0+qp*betapmer) erfc = ((((a5*qt+a4)*qt+a3)*qt+a2)*qt+a1) & *qt*expar2 erf = 1.00D0 - erfc C . direct energy edir = edir + ch*erfc/rr C . calculate contribution to gradient of direct energy; C . dfe is gradient/|r| of direct energy dfe = chk*(-twopisqrt*betapmer*expar2+erf) & /(rr*rsqr) C . add direct and subtract coulomb contributions to gradient: C . dEclmb/dr_i = SUM(j) -q_i*q_j*r_ij/|r_ij|^3 C . but here dEclmb/dr_i is used as a correction, C . so change the sign. dedir(iat3) = dedir(iat3) + x*dfe + dx dedir(iat3+1) = dedir(iat3+1) + y*dfe + dy dedir(iat3+2) = dedir(iat3+2) + z*dfe + dz dedir(jat3) = dedir(jat3) - x*dfe - dx dedir(jat3+1) = dedir(jat3+1) - y*dfe - dy dedir(jat3+2) = dedir(jat3+2) - z*dfe - dz enddo 90 enddo enddo enddo c all intramolecular terms: only energy contribution, skip gradient do ires=1,nres ib = irpnt(ires) ie = irpnt(ires+1)-1 do iat=ib,ie iat3 = 3*iat-2 xi = xyz(1,iat) yi = xyz(2,iat) zi = xyz(3,iat) chi = atchg(iat) do jat=ib,iat-1 jat3 = 3*jat-2 if (pbc) then call pbcxyz(iat,jat,xj,yj,zj) else xj = xyz(1,jat) yj = xyz(2,jat) zj = xyz(3,jat) endif x = xi-xj y = yi-yj z = zi-zj rsqr = x*x + y*y + z*z rr = sqrt(rsqr) rr3 = rr**3 ch = chi*atchg(jat) C . short range Coulomb energy eclmb = eclmb + ch/rr betapmer = betapme*rr expar2 = exp(-betapmer*betapmer) qt = 1.0/(1.0+qp*betapmer) erfc = ((((a5*qt+a4)*qt+a3)*qt+a2)*qt+a1) & *qt*expar2 C . direct energy edir = edir + ch*erfc/rr enddo enddo enddo C energy conversion: until now energy is in e*e/A C multiply this by 1/(4 pi eps0) and convert to eV edir = 14.39964521*edir eself = 14.39964521*eself eclmb = 14.39964521*eclmb call etimer(t2) tpmedir = t2-t1 end