LamellarOrderingExternal.h

#ifndef SIMP_LAMELLAR_ORDERING_EXTERNAL_H
#define SIMP_LAMELLAR_ORDERING_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 LamellarOrderingExternal : public ParamComposite 
   {
   
   public:
   
      LamellarOrderingExternal();

      LamellarOrderingExternal(const LamellarOrderingExternal& other);

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

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

      int perpDirection_;

      double interfaceWidth_;

      DArray<double> prefactor_;   

      double externalParameter_;

      int periodicity_;

      Boundary *boundaryPtr_;
   
      int    nAtomType_; 

      bool  isInitialized_;

   };
  
   // inline methods 
 
    
   inline double LamellarOrderingExternal::energy(double d, int type) const
   {
      double perpLength, q, clipParameter, arg, clipcos;
      Vector lengths;
      lengths = boundaryPtr_->lengths();
      perpLength = lengths[perpDirection_];

      q = (2.0*M_PI*periodicity_)/perpLength;
      clipParameter   = 1.0/(q*interfaceWidth_*perpLength);
      arg = q*d;
      clipcos = clipParameter*cos(arg);
      
      return prefactor_[type]*externalParameter_*tanh(clipcos);
   }

   /* 
   * Calculate external potential energy for a single atom.
   */
   inline 
   double LamellarOrderingExternal::energy(const Vector& position, int type) const
   {
      double d, totalEnergy;
      totalEnergy = 0.0;
      d = position[perpDirection_];
      totalEnergy += energy(d, type);
       
      return totalEnergy;
   }

   /* 
   * Calculate force for a particle as a function of distance to boundary.
   */
   inline double LamellarOrderingExternal::forceScalar(double d, int type) const
   {
      double perpLength, q, clipParameter, arg, clipcos, tanH, sechSq;
      Vector lengths;
      lengths = boundaryPtr_->lengths();
      perpLength = lengths[perpDirection_];

      q = (2.0*M_PI*periodicity_)/perpLength;
      clipParameter = 1.0/(q*interfaceWidth_*perpLength);
      arg = q*d;
      clipcos = clipParameter*cos(arg);
      tanH = tanh(clipcos);
      sechSq = (1.0 - tanH*tanH);
      return prefactor_[type]*externalParameter_*sechSq*clipParameter*sin(arg)*q;
   }

   /* 
   * Calculate external force for a single atom.
   */
   inline 
   void LamellarOrderingExternal::getForce(const Vector& position, int type, 
                                     Vector& force) const
   {
      double d = position[perpDirection_];
      double scalarf;
      force.zero();
      scalarf = forceScalar(d, type);
      force[perpDirection_] = scalarf;
   }
 
}
#endif