C
C $Id: pme_direct.F,v 1.4 1998/07/16 16:40:17 jjv5 Exp $
C
C------------------------------------------------------------------------
      subroutine pme_direct2(edir,eself,eclmb,dedir)

      implicit double precision(a-h,o-z)
#include "divcon.dim"
#include "divcon.h"
#include "constants.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      parameter (eta=0.22d0)

C
C     this routine returns the direct energy, the selfenergy and
C     the short-range (classical) Coulomb energy and gradients.
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     parameters:
      dimension dedir(maxpar)

CC----------------------------------------------------------CC
        
      call etimer(t1)
c        eta = betapme/6.0d0
C     units of eclmb, edir and eself are eV
      eclmb = 0.0
      edir = 0.0
      eself = 0.0
         do i=1,3*natoms
          dedir(i)=0.d0
         enddo




C     unit of gradient (dedir) is kcal/A

c      tokcal = 332.0704
c      tokcal = 332.0806806d0
c      bpisqrt = betapme/pisqrt
        eta=betapme                        
      do iat=1,natoms
         iat3 = 3*iat-2
         xi = xyz(1,iat)
         yi = xyz(2,iat)
         zi = xyz(3,iat)
         chi = chgpme(iat)
    

C     .  selfenergy (is a constant, it doesn't contribute to the gradient)
c         eself = eself - bpisqrt*chi*chi

         do jat=1,natoms
       
            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
      do nn1=-1,1
       do nn2=-1,1
        do nn3=-1,1
       if ((iat.eq.jat).and.((abs(nn1)+abs(nn2)+abs(nn3)).eq.0))then
            edir=edir - 2.0d0*chi*chi*eta/dsqrt(pi)
c           goto 128
c          edir = 0.0
       else

            x = xi-xj - nn1*dbox(1)
            y = yi-yj - nn2*dbox(1)
            z = zi-zj - nn3*dbox(1)


            rsqr = x*x + y*y + z*z

            rr = dsqrt(rsqr)
            rr3 = rr**3
            ch = chi*chgpme(jat)


C     .     short range Coulomb energy
c            eclmb = eclmb + ch/rr

C     .     calculate contribution to gradient of short range Coulomb
            dx = (ch*x)/rr3
            dy = (ch*y)/rr3
            dz = (ch*z)/rr3
            etar = eta*rr
            expar2 = exp(-etar*etar)
            qt     = 1.0/(1.0+qp*etar)
            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 = ch*(-twopisqrt*etar*expar2+erf-1)/(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
            dedir(iat3+1) = dedir(iat3+1) + y*dfe
            dedir(iat3+2) = dedir(iat3+2) + z*dfe

            dedir(jat3) = dedir(jat3) - x*dfe
            dedir(jat3+1) = dedir(jat3+1) - y*dfe
            dedir(jat3+2) = dedir(jat3+2) - z*dfe

           if ((abs(nn1)+abs(nn2)+abs(nn3)).eq.0) then
            dedir(iat3) = dedir(iat3) + dx
            dedir(iat3+1) = dedir(iat3+1) + dy
            dedir(iat3+2) = dedir(iat3+2) + dz

            dedir(jat3) = dedir(jat3) - dx
            dedir(jat3+1) = dedir(jat3+1) - dy
            dedir(jat3+2) = dedir(jat3+2) - dz    

            edir = edir - ch/rr
            eclmb= eclmb + ch/rr
           endif
    
       endif
           enddo
          enddo
         enddo
  
          enddo
         enddo
c            dedir(4) = - x*dfe - dx
c            dedir(5) = - y*dfe - dy
c            dedir(6) = - z*dfe - dz


 

C     energy conversion: until now energies are in e*e/A
C     multiply this by 1/(4 pi eps0) and then convert to V

         edir = 14.39964521*edir*0.5d0
         eclmb = 14.39964521*eclmb*0.5d0
          do i=1,3*natoms
            dedir(i) = 14.39964521*dedir(i)*eV2kcal*0.5d0
          enddo
c      edir = 14.39964521*edir
c 128     eself = 14.39964521*eself
c      eclmb = 14.39964521*eclmb
c       do i=1,3*natoms
c         dedir(i)=14.39964521* dedir(i)
c       enddo
      call etimer(t2)
      tpmedir = t2-t1

      end