GeneralPeriodicExternal.h

#ifndef SIMP_GENERAL_PERIODIC_EXTERNAL_H
#define SIMP_GENERAL_PERIODIC_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/space/Vector.h>
#include <util/param/ParamComposite.h>
#include <util/global.h>
#include <cmath>

namespace Simp
{

   using namespace Util;

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

      GeneralPeriodicExternal(const GeneralPeriodicExternal& other);

      GeneralPeriodicExternal& operator = (const GeneralPeriodicExternal& other);

      void setNAtomType(int nAtomType);

      void setExternalParameter(double externalParameter);

      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;

      double externalParameter() const;

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

      void getForce(const Vector& position, int type, Vector& force) const;
 
      std::string className() const;
 
   private:
   
      static const int MaxAtomType = 3;

      DArray<double> prefactor_;

      double externalParameter_;

      int  nWaveVectors_;

      DArray<Vector>  waveVectors_;

      DArray<double> phases_;

      DArray<double> shifts_;

      int periodicity_;

      double interfaceWidth_;

      Boundary *boundaryPtr_;
   
      int    nAtomType_; 

      bool  isInitialized_;

   };
  
   // inline methods 

   /* 
   * Calculate external potential energy for a single atom.
   */
   inline double GeneralPeriodicExternal::energy(const Vector& position, int type) const
   {
      const Vector cellLengths = boundaryPtr_->lengths();
      double clipParameter = 1.0/(2.0*M_PI*periodicity_*interfaceWidth_);
      
      Vector r;
      Vector a;
      r.zero();
      a.zero();
      double cosine = 0.0;
      for (int i = 0; i < nWaveVectors_; ++i) {
         r = position;
         for (int j = 0; j < Dimension; j++) {
            a[j] = shifts_[i*Dimension + j];
         }
         r -= a;
         Vector q;
         q[0] = 2.0*M_PI*periodicity_*waveVectors_[i][0]/cellLengths[0];
         q[1] = 2.0*M_PI*periodicity_*waveVectors_[i][1]/cellLengths[1]; 
         q[2] = 2.0*M_PI*periodicity_*waveVectors_[i][2]/cellLengths[2];
         double arg = q.dot(r)+phases_[i];
         cosine += prefactor_[i*nAtomType_ + type]*cos(arg);
      }
      cosine *= clipParameter;
      return externalParameter_*tanh(cosine);
   }

   /* 
   * Calculate external force for a single atom.
   */
   inline 
   void GeneralPeriodicExternal::getForce(const Vector& position, int type, 
                                     Vector& force) const
   {
      const Vector cellLengths = boundaryPtr_->lengths();
      double clipParameter = 1.0/(2.0*M_PI*periodicity_*interfaceWidth_);
 
      Vector r;
      Vector a;
      r.zero();
      a.zero();
      double cosine = 0.0;
      Vector deriv;
      deriv.zero();
      for (int i = 0; i < nWaveVectors_; ++i) {
         r = position;
         for (int j = 0; j < Dimension; j++) {
            a[j] = shifts_[i*Dimension + j];
         }
         r -= a;
         Vector q;
         q[0] = 2.0*M_PI*periodicity_*waveVectors_[i][0]/cellLengths[0];
         q[1] = 2.0*M_PI*periodicity_*waveVectors_[i][1]/cellLengths[1]; 
         q[2] = 2.0*M_PI*periodicity_*waveVectors_[i][2]/cellLengths[2];
         double arg = q.dot(r)+phases_[i];
         cosine += prefactor_[i*nAtomType_ + type]*cos(arg);
         double sine = -1.0*sin(arg);
         q *= prefactor_[i*nAtomType_ + type]*sine;
         deriv += q;
      }
      cosine *= clipParameter;
      deriv *= clipParameter;
      double tanH = tanh(cosine);
      double sechSq = (1.0 - tanH*tanH);
      double f = externalParameter_*sechSq;
      deriv *= -1.0*f;
      force = deriv;
   }
 
}
#endif