BoxExternal.h

#ifndef SIMP_BOX_EXTERNAL_H
#define SIMP_BOX_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>

namespace Simp
{

   using namespace Util;

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

      BoxExternal(const BoxExternal& other);

      BoxExternal& operator = (const BoxExternal& other);



      void setNAtomType(int nAtomType);

      void setBoundary(Boundary &boundary);

      void setExternalParameter(double externalParameter);

      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 externalParameter() const;
 
      double energy(double d, int i) const;
 
      double energy(const Vector& position, int i) const;

      double forceScalar(double d, int i) const;
 
      void getForce(const Vector& position, int type, Vector& force) const;
 

   private:
   
      static const int MaxAtomType = 2;

      double epsilon_;

      double sigma_;

      double sigmaCb_;
 
      double cutoff_;

      double coeff_;
  
      Boundary *boundaryPtr_;
   
      int    nAtomType_; 

      bool  isInitialized_;

   };
  
   // inline methods 
 
   /* 
   * Calculate repulsive potential energy for particle of type i.
   *
   * 7.905694 = 25/sqrt(10), which is the negative of the minimum of the
   * unshifted 9-3 potential.
   */
   inline double BoxExternal::energy(double d, int type) const
   {
      if (d < cutoff_) {
         double r3i;
         r3i = d * d * d;
         if (r3i < 0.129615*sigmaCb_) r3i = 0.129615 * sigmaCb_;
         r3i = sigmaCb_ / r3i;
         return coeff_ * (r3i*(r3i*r3i - 7.5) + 7.905694);
      } else {
         return 0.0;
      }
   }
   
   /* 
   * Calculate external potential energy for a single atom.
   */
   inline double BoxExternal::energy(const Vector& position, int type) const
   {
      double d, totalEnergy;
      Vector halfL = boundaryPtr_->lengths();
      halfL *= 0.5;
      totalEnergy = 0.0;
      for (int i = 0; i < Dimension; ++i) {
         d = position[i];
         if (d > halfL[i]) d = halfL[i] + halfL[i] - d;
            totalEnergy += energy(d, type);
      }
      return totalEnergy;
   }
   
   /* 
   * Calculate force for a particle as a function of distance to boundary.
   */
   inline double BoxExternal::forceScalar(double d, int i) const
   {
      if (d < cutoff_) {
         double r3i;
         r3i = d * d * d;
         if (r3i < 0.129615 * sigmaCb_) r3i = 0.129615 * sigmaCb_;
         r3i = sigmaCb_ / r3i;
         return coeff_ * (9.0*r3i*r3i - 22.5) * r3i / d;
      } else {
         return 0.0;
      }
   }

   /* 
   * Calculate external force for a single atom.
   */
   inline void 
   BoxExternal::getForce(const Vector& position, int type, Vector& force) const
   {
      double d = position[2];
      double scalarf;
      Vector halfL = boundaryPtr_->lengths();
      halfL *= 0.5;
      for (int i = 0; i < Dimension; ++i) {
         d = position[i];
         if (d > halfL[i]) {
            d = halfL[i] + halfL[i] - d;
            scalarf = forceScalar(d, type);
            force[i] = -scalarf;
         } else {
            scalarf = forceScalar(d, type);
            force[i] = scalarf;
         }
      }
   }
 
}
#endif