C
C $Id: pme_recipmc.F,v 1.5 1998/07/16 16:40:20 jjv5 Exp $
C
C------------------------------------------------------------------------
      subroutine pme_recipmc(erec,derec)

      implicit double precision(a-h,o-z)
#include "divcon.dim"
#include "divcon.h"

C     this routine calculates the reciprocal enery and derivatives
C     by the PME method of J.Chem.Phys. 103 ('95) 8577.
C
C     written by Arjan van der Vaart, Dec. '97
C
C     parameters:
C     *) erec            the reciprocal energy
C     *) derec           the derivative of erec to the coordinates
C
C     compiler directives:
C     *) MEMORY_OVERLAP         if defined, the pme variables will
C     .                         share memory with the eigenvectors /
C     .                         eigenvalues
C
C

C     parameters:
      dimension derec(maxpar)

C     local:
      dimension ndimn(3), ch(maxatm)

#ifdef MEMORY_OVERLAP

C     now the following arrays share memory:
C     array1           common |  array2        common   | shared memory array
C     -----------------------------------------------------------------------
C     qpmec(maxkpme32) /pme/  |  ff(msorb2)    /work/   | ff(maxovlp1)
C     qpme(maxkpme3)   /pme/  |  ww(msorb2)    /work/   | ww(maxovlp2)
C

#define QPMEC_ ff
#define QPME_ ww

#else

C     no sharing of memory of the pme arrays with the /work/ arrays

#define QPMEC_ qpmec
#define QPME_ qpme

#endif

CC---------------------------------------------------------------------CC

      call etimer(t1)

C     initialization for the Fourier transforms
      ndim = 3
      ndimn(1) = k1pme
      ndimn(2) = k2pme
      ndimn(3) = k3pme

C     intra- and intermolecular contributions
      call pme_recip(erec, derec)

C     subtract intramolecular contributions
      do i=1,natoms
         ch(i) = 0.0
      enddo

      n = 1
      do igr=1,ngroup

         ir = ingroup(n)
         i1 = irpnt(ir)
         i2 = irpnt(ir+1)-1
         nat = i2-i1+1
         nat3 = 3*nat

C     .  average the charges
         do i=0,ingroupn(igr)-1
            irr = ingroup(n+i)
            ii1 = irpnt(irr)
            ii2 = irpnt(irr+1)-1
            k = i1
            do ii=ii1,ii2
               ch(k) = ch(k) + atchg(ii)
               k = k + 1
            enddo
         enddo
         xn = dble(ingroupn(igr))
         do k=i1,i2
            ch(k) = ch(k)/xn
         enddo

c     .  calculate the function Q
         call pme_calcqmc(ir, ch)

C     .  copy Q to QPMEC_, since QPMEC_ will be overwritten with F[Q]
         j = 1
         do 10 i=1,k123pme
            QPMEC_(j) = QPME_(i)
            j = j + 1
            QPMEC_(j) = 0.0
            j = j + 1
 10      continue

C     .  Fourier transform Q and multiply with theta
         call fourn(QPMEC_,ndimn,ndim,iminone)
         j = 1
         do 20 i=1,k123pme
            QPMEC_(j) = QPMEC_(j)*thetapme(i)
            j = j + 1
            QPMEC_(j) = QPMEC_(j)*thetapme(i)
            j = j + 1
 20      continue

C     .  after the following transformation, QPMEC_ will contain Conv(theta,Q)
         call fourn(QPMEC_,ndimn,ndim,ione)

C     .  calculate the gradient
         call pme_derecmc(ir, igr, n, ch, derec)

         n = n + ingroupn(igr)

      enddo

      call etimer(t2)
      tpmerec = t2-t1

      end

CC---------------------------------------------------------------------CC

#undef QPMEC_
#undef QPME_