SlitExternal.h

#ifndef SIMP_SLIT_EXTERNAL_H
#define SIMP_SLIT_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/param/ParamComposite.h>
#include <util/global.h>
#include <cmath>

namespace Simp
{

   using namespace Util;

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

      SlitExternal(const SlitExternal& other);

      SlitExternal& operator = (const SlitExternal& 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 SlitExternal::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 SlitExternal::energy(const Vector& position, int type) const
   {
      double d = position[2];
      double halfL = boundaryPtr_->lengths()[2] * 0.5;
      if (d > halfL) d = halfL + halfL - d;
      return energy(d, type);
   }
   
   /* 
   * Calculate force for a particle as a function of distance to boundary.
   */
   inline double SlitExternal::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 
   SlitExternal::getForce(const Vector& position, int type, Vector& force) const
   {
      double d = position[2];
      double halfL = boundaryPtr_->lengths()[2] * 0.5;
      double scalarf;
      if (d > halfL) {
         d = halfL + halfL - d;
         scalarf = forceScalar(d, type);
         force = Vector(0.0, 0.0, -scalarf);
      } else {
         scalarf = forceScalar(d, type);
         force = Vector(0.0, 0.0, scalarf);
      }
   }
 
}
#endif