Main Program

c **********************************************************************
c     PME.F
c **********************************************************************
c 
c  ProTeam's Particle Mesh Ewald (PME) code.  
c
c **********************************************************************

      program PME

      use global_parameters
      use scalars
      use serial_arrays
      use pmeVar
      implicit none

c     define constants
      pi = two*asin(one)
      rtPi = sqrt(pi)

c     hardware whoami, whomai is used in parallel code, but at the present
c     we our running our code on a serial machine
      whoami = 0

      alpha = two

      maxAtoms = 100
      numt     = 20
      maxt     = 30
     
      allocate(    numa(0:maxt)    )
      allocate(   itype(maxAtoms)      )
      allocate(  iatnum(maxAtoms)      )
      allocate(    atom(3,maxAtoms)    )
      allocate(   qatom(maxAtoms)      )
 
      call ReadInput 
      
c     Now that we know meshsize, allocate Qmesh

      allocate( Qmesh(0:meshSize,0:meshSize,0:meshSize) )

      call SetupFFT ! move outside timesteps in real MD simulation

      call CalcPME
     
      print *, "coulombic potential  = ", eCouloumb

      stop

      end program PME


c **********************************************************************
      subroutine CalcPME
c **********************************************************************

      use global_parameters
      use scalars
      use serial_arrays
      use pmeVar
      implicit none
      
      call CalcReal
      call CalcRecp
      call CalcSi
      
      
      eCouloumb = realTerm + recpTerm + selfTerm
      return
      end

c------------------------------------------------------------------------------
      subroutine SetupFFT
c     setup for 3D complex-to-complex FFT.
c------------------------------------------------------------------------------
c
      use scalars
      use serial_arrays
      use pmeVar
      implicit none
      include 'fftw_f77.i'

      call fftw3d_f77_create_plan(meshForward,meshSize,meshSize,
     .     meshSize,FFTW_FORWARD,FFTW_MEASURE + FFTW_IN_PLACE)

      call fftw3d_f77_create_plan(meshBackward,meshSize,meshSize, 
     .     meshSize,FFTW_BACKWARD,FFTW_MEASURE + FFTW_IN_PLACE)

      return
      end