EwaldInteraction.h

#ifndef SIMP_EWALD_INTERACTION_H
#define SIMP_EWALD_INTERACTION_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 <util/param/ParamComposite.h>
#include <util/global.h>

#include <math.h>

namespace Simp
{

   using namespace Util;

   class EwaldInteraction : public ParamComposite
   {
   public:

      EwaldInteraction();

      EwaldInteraction(const EwaldInteraction& other);

      EwaldInteraction& operator = (const EwaldInteraction& other);

      ~EwaldInteraction() 
      {}



      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 rSpaceEnergy(double rSq, double qProduct) const;

      double rSpaceForceOverR(double rSq, double qProduct) const;
  
      double kSpacePotential(double kSq) const;



      double epsilon() const;
 
      double alpha() const;
 
      double rSpaceCutoff() const;

      double rSpaceCutoffSq() const;

      double get(std::string name) const;


   private:

      double epsilon_;          

      double alpha_; 

      double rSpaceCutoff_;  

      // Derived constants
      double rSpaceCutoffSq_; 
      double ce_;
      double cf_;
      double cg_;

      bool isInitialized_;

      void setDerivedConstants();

   };

   // Inline methods 

   /* 
   * Return medium dielectric permittivity.
   */
   inline double EwaldInteraction::epsilon() const
   {  return epsilon_; }

   /* 
   * Return Ewald smearing parameter alpha.
   */
   inline double EwaldInteraction::alpha() const
   {  return alpha_; }

   /* 
   * Return real space cutoff distance.
   */
   inline 
   double EwaldInteraction::rSpaceCutoff() const
   {  return rSpaceCutoff_; }

  /* 
   * Return real space cutoff distance squared.
   */
   inline 
   double EwaldInteraction::rSpaceCutoffSq() const
   {  return rSpaceCutoffSq_; }

   /* 
   * Compute and return r-space energy for a pair of charges.
   */
   inline 
   double EwaldInteraction::rSpaceEnergy(double rSq, double qProduct) 
   const 
   {
      double r = sqrt(rSq);
      return ce_*qProduct*erfc(alpha_*r)/r;
   }

   /*
   * Compute and return (r-space force) / distance for a pair of charges.
   */
   inline 
   double EwaldInteraction::rSpaceForceOverR(double rSq, double qProduct) 
   const 
   {
      double r = sqrt(rSq);
      double x = alpha_*r;
      return ce_*qProduct*(erfc(x) + cf_*r*exp(-x*x))/(r*rSq); 
   }

   /* 
   * Calculate k-space potential from squared wavenumber kSq.
   */
   inline
   double EwaldInteraction::kSpacePotential(double kSq) const
   {  return exp(cg_*kSq)/(kSq*epsilon_); }

}
#endif