MdSpmePotential.h

#ifndef MCMD_MD_SPME_POTENTIAL_H
#define MCMD_MD_SPME_POTENTIAL_H

/*
* Simpatico - Simulation Package for Polymeric and Molecular Liquids
*
* Copyright 2010 - 2012, David Morse (morse012@umn.edu)
* Distributed under the terms of the GNU General Public License.
*/

#include <mcMd/potentials/coulomb/MdCoulombPotential.h>      // base class
#include <mcMd/potentials/coulomb/EwaldRSpaceAccumulator.h>  // member
#include <mcMd/chemistry/AtomType.h>     // member template parameter

#include <simp/interaction/coulomb/EwaldInteraction.h>       // member
#include <simp/boundary/Boundary.h>      // typedef

#include <util/space/IntVector.h>        // member template parameter
#include <util/space/Vector.h>           // member template parameter
#include <util/space/Tensor.h>           // member template parameter
#include <util/containers/Pair.h>        // member template parameter
#include <util/containers/GArray.h>      // member template
#include <util/containers/GridArray.h>   // member template
#include <util/misc/Setable.h>           // member template
#include <util/containers/Array.h>       // member class template

#include <complex>
#include <fftw3.h>

namespace McMd
{

   class Simulation;
   class System;

   typedef std::complex<double> DCMPLX;

   using namespace Util;
   using namespace Simp;

   class MdSpmePotential : public MdCoulombPotential
   {

   public:

      MdSpmePotential(System& system);

      virtual ~MdSpmePotential();



      virtual void readParameters(std::istream& in);

      virtual void loadParameters(Serializable::IArchive &ar);

      virtual void save(Serializable::OArchive &ar);



      void set(std::string name, double value);

      double get(std::string name) const;



      virtual void makeWaves();

      virtual int nWave() const;

      virtual void addForces();

      virtual void computeEnergy();

      virtual void computeStress();



      EwaldRSpaceAccumulator& rSpaceAccumulator()
      {  return rSpaceAccumulator_; }

      EwaldInteraction& ewaldInteraction()
      {  return ewaldInteraction_; }


   private:

      // Ewald Interaction - core Ewald computations
      EwaldInteraction ewaldInteraction_;

      // Pointer to parent Simulation
      Simulation* simulationPtr_;

      // Pointer to parent System
      System* systemPtr_;

      // Pointer to boundary of associated System.
      Boundary* boundaryPtr_;

      // Pointer to array of atom types
      const Array<AtomType>* atomTypesPtr_;

      IntVector gridDimensions_;

      GridArray<DCMPLX> rhoR_;

      GridArray<DCMPLX> rhoK_;

      GridArray<double> g_;
      
      GridArray<double> sqWaves_;
      
      GridArray<Vector> vecWaves_;
     
      GridArray<DCMPLX> xfield_;

      GridArray<DCMPLX> yfield_;

      GridArray<DCMPLX> zfield_;

      int order_;
      
      fftw_plan forward_plan;

      fftw_plan xfield_backward_plan, yfield_backward_plan, zfield_backward_plan;

      template<class T> 
      void setGridToZero(GridArray<T>& grid);

      void setGridDimensions();

      void computeWaves();

      double bfactor(double m , int dim);

      void assignCharges();

      double basisSpline(double x);

   };

}
#endif