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