Reciprocal Space Term

c **********************************************************************
c     recpTerm.F
c **********************************************************************
      subroutine CalcRecp
c **********************************************************************

      use global_parameters
      use scalars
      use serial_arrays
      use pmeVar

      implicit none

      integer, dimension(3) :: k 
      integer i1, i2, i3, ksq
      real Qsq

      pi = two*asin(one)
      rtPi = sqrt(pi)
      pisq = pi*pi
      alphasq = alpha*alpha

      recpTerm = zero
      Call MeshInterp   !use atom(),qatom() to calculate realQmesh() 
                        ![mesh=3d allocatable array]
      
      Call meshFFTforward   !transforms realQmesh() to recpQmesh(), 
                            ! includes FFTW setup?
      do i1 = 1, meshSize
         do i2 = 1, meshSize
            do i3 = 1, meshSize
               k(1) = i1
               k(2) = i2
               k(3) = i3
               ksq = dot_product(k,k)
               Qsq = dot_product(Qmesh(k), Qmesh(k))
               recpTerm = recpTerm+exp(-pisq*ksq/alphasq)/(two*pi*ksq)*
     &          Qsq
            end do
         end do
      end do
      print *, "recp term = ",recpTerm
      
      return

      end


c **********************************************************************
      subroutine MeshInterp
c **********************************************************************

      use global_parameters
      use scalars
      use serial_arrays
      use pmeVar

      implicit none
      integer i,j,k,m,n,even, k1,k2,k3,nMP,MP
      real sum1, spline
      integer, dimension(:,:), allocatable :: nearestMesh
      complex test
      Qmesh = zero
      allocate(nearestMesh(p,3))

      pi = two*asin(one)
      rtPi = sqrt(pi)
      do i=1,numa(0)
         do j = 1, 3
            nMP = anint(atom(j,i)*meshSize) 
            nearestMesh(1,j) = nMP
            if (nMP > meshSize) then
               print *, "error mesh point is ",nMP
               stop
            end if
            do k = 1, (p-1)/2
               if (nMP+k > meshSize) then
                  nearestMesh(1+k,j) = nMP+k - meshSize
                  nearestMesh(1+p/2+k,j) = nMP-k
               else
                  if (nMP-k < 0) then
                     nearestMesh(1+k,j) = nMP+k
                     nearestMesh(1+p/2+k,j) = nMP-k + meshSize
                  else
                     nearestMesh(1+k,j) = nMP+k
                     nearestMesh(1+p/2+k,j) = nMP-k
                  end if
               end if
            end do
            even = p/2 -(p-1)/2
            if (even .eq. 1) then
               if (nMP>meshSize) then
                  MP = nMP-meshSize
               else if (MP<0) then
                  MP = nMP+meshSize
               else
                  MP = nMP
               end if
               if (nMP>atom(j,i)) then
                  nearestMesh(p,j) = MP - k
                else
                  nearestMesh(p,j) = MP + k
               end if
            end if
         end do
         do n=1,p
            sum1 = zero
            do m=1,3
               sum1= sum1+qatom(i)*spline((atom(j,i)-
     .real(nearestMesh(n,m))/real(MeshSize))*meshSize,p)
            end do
            k1 = nearestMesh(n,1)
            k2 = nearestMesh(n,2)
            k3 = nearestMesh(n,3)
            Qmesh(k1,k2,k3) = Qmesh(k1,k2,k3)+cmplx(sum1,zero) 
         end do
      end do

      return

      end


c **********************************************************************
      subroutine meshFFTforward
c **********************************************************************

      use global_parameters
      use scalars
      use pmeVar
      use serial_arrays
      
      implicit none
      include 'fftw_f77.i'

      scaleG = one/meshSize**3
      
      call fftwnd_f77_one(meshForward,Qmesh,0)
      Qmesh = Qmesh*scaleG

      return

      end