SphericalTabulatedExternal.h

#ifndef SIMP_SPHERICAL_TABULATED_EXTERNAL_H
#define SIMP_SPHERICAL_TABULATED_EXTERNAL_H

/*
* Simpatico - Simulation Package for Polymeric and Molecular Liquids
*
* Copyright 2010 - 2017, The Regents of the University of Minnesota
* Distributed under the terms of the GNU General Public License.
*/

#include <simp/boundary/Boundary.h>
#include <util/space/Dimension.h>
#include <util/param/ParamComposite.h>
#include <util/global.h>

#include <cmath>
#include <string>

namespace Simp
{

   using namespace Util;

   class SphericalTabulatedExternal : public ParamComposite 
   {
   
   public:
   
      SphericalTabulatedExternal();

      SphericalTabulatedExternal(const SphericalTabulatedExternal& other);

      SphericalTabulatedExternal& 
      operator = (const SphericalTabulatedExternal& other);



      void setNAtomType(int nAtomType);

      void setBoundary(Boundary &boundary);

      void readParameters(std::istream &in);

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

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

      /*
      * Set a potential energy parameter, identified by a string.
      */
      void set(std::string name, double value);

      /*
      * Get a parameter value, identified by a string.
      */
      double get(std::string name) const;



      std::string className() const;

      double energy(const Vector& position, int i) const;

      void getForce(const Vector& position, int type, Vector& force) const;
 

   private:
   
      static const int MaxAtomType = 6;

      std::string filename_;

      double rMax_;

      double rMaxSq_;

      double dr_;
  
      Boundary *boundaryPtr_;
   
      int nr_; 

      int nAtomType_; 

      DArray< DArray<double> > potentials_;

      bool isInitialized_;

   };
  
   // inline methods 
 
   /* 
   * Calculate external potential energy for a single atom.
   */
   inline 
   double SphericalTabulatedExternal::energy(const Vector& position, int type) const
   {
      Vector origin = boundaryPtr_->lengths();
      origin *= 0.5;
      double rSq = boundaryPtr_->distanceSq(position, origin);
      if (rSq < rMaxSq_) {
         double r = sqrt(rSq);
         int k = (int) (r/dr_);
         double fp = r - k*dr_;
         UTIL_CHECK(k >= 0);
         UTIL_CHECK(k < nr_ -1);
         UTIL_CHECK(fp >= 0.0);
         UTIL_CHECK(fp <= 1.0);
         double fm = 1.0 - fp;
         return fm*potentials_[type][k] + fp*potentials_[type][k+1]; 
      } else {
         return 0.0;
      }
   }
   
   /* 
   * Calculate external force for a single atom.
   */
   inline void 
   SphericalTabulatedExternal::getForce(const Vector& position, int type, 
                                        Vector& force) 
   const
   {
      Vector origin = boundaryPtr_->lengths();
      origin *= 0.5;
      double rSq = boundaryPtr_->distanceSq(position, origin, force);
      if (rSq < rMaxSq_) {
         double r = sqrt(rSq);
         int k = (int) (r/dr_);
         UTIL_CHECK(k >= 0);
         UTIL_CHECK(k < nr_ -1);
         double fs = potentials_[type][k] - potentials_[type][k+1];
         fs /= dr_;
         force *= fs/r;
      } else {
         for (int i = 0; i < Dimension; ++i) {
            force[i] = 0.0;
         }
      }
   }
 
}
#endif