C
C $Id: pme_recip.F,v 1.8 1998/07/16 16:40:20 jjv5 Exp $
C
C------------------------------------------------------------------------
      subroutine pme_recip(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)

#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     calculate the function Q
      call pme_calcq

C     copy Q to QPME_, 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_derec(derec)

C     calculate the reciprocal energy in kcal/mol
      erec = 0.0
      kk = 1
      do 30 k=1,k123pme
         erec = erec + QPMEC_(kk)*QPME_(k)
         kk = kk + 2
 30   continue
     
      erec = 7.199822605*erec
c      do i=1,3*natoms
c         derec(i)=derec(i)*7.199822605
c      enddo


      call etimer(t2)
      tpmerec = t2-t1

      end

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

#undef QPMEC_
#undef QPME_