LocalLamellarOrderingExternal.h

#ifndef SIMP_LOCAL_LAMELLAR_ORDERING_EXTERNAL_H
#define SIMP_LOCAL_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 LocalLamellarOrderingExternal : public ParamComposite 
   {
   
   public:
   
      LocalLamellarOrderingExternal();

      LocalLamellarOrderingExternal(const LocalLamellarOrderingExternal& other);

      LocalLamellarOrderingExternal& operator = (const LocalLamellarOrderingExternal& 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 = 2;

      int perpDirection_;

      int parallelDirection_;

      double fraction_;

      double width_;

      DArray<double> prefactor_;   

      double externalParameter_;

      int periodicity_;

      Boundary *boundaryPtr_;
   
      int    nAtomType_; 

      bool  isInitialized_;

   };
  
   // inline methods 
 
   /* 
   inline double LocalLamellarOrderingExternal::energy(double dPerp, double dParallel, int type) const
   {
      double perpLength, parallelLength, q, clipParameter, arg, clipcos, parallelFactor;
      Vector lengths;
      lengths = boundaryPtr_->lengths();
      perpLength = lengths[perpDirection_];
      q = (2.0*M_PI*periodicity_)/perpLength;
      clipParameter   = 1.0/(q*width_*perpLength);
      arg = q*dPerp;
      clipcos = clipParameter*cos(arg);
      parallelLength = lengths[parallelDirection_];
      parallelFactor = tanh( 2.0*M_PI*((dParallel - (((1.0-fraction_)/2.0)*parallelLength))/(fraction_*parallelLength)) );
      parallelFactor *= tanh( 2.0*M_PI*(((((1.0+fraction_)/2.0)*parallelLength) - dParallel)/(fraction_*parallelLength)) );
      parallelFactor += 1.0;
      parallelFactor /= 2.0;
      if (parallelFactor < 0.0)
         UTIL_THROW("Negative base factor");
      return prefactor_[type]*externalParameter_*tanh(clipcos)*parallelFactor;
   }
   */

   /* 
   * Calculate external potential energy for a single atom.
   */
   inline 
   double LocalLamellarOrderingExternal::energy(const Vector& position, int type) const
   {
      double dPerp, dParallel, perpLength, parallelLength, q, clipParameter, arg, clipcos;
      double parallelFactor1, parallelFactor2, parallelFactor;
      dPerp = position[perpDirection_];
      dParallel = position[parallelDirection_];

      Vector lengths;
      lengths = boundaryPtr_->lengths();
      perpLength = lengths[perpDirection_];
      parallelLength = lengths[parallelDirection_];

      q = (2.0*M_PI*periodicity_)/perpLength;
      clipParameter   = 1.0/(q*width_*perpLength);
      arg = q*dPerp;
      clipcos = clipParameter*cos(arg);
      parallelFactor1 = tanh( 2.0*M_PI*((dParallel - (((1.0-fraction_)/2.0)*parallelLength))/(fraction_*parallelLength)) );
      parallelFactor2 = tanh( 2.0*M_PI*(((((1+fraction_)/2.0)*parallelLength) - dParallel)/(fraction_*parallelLength)) );
      parallelFactor = 1.0 + (parallelFactor1*parallelFactor2);
      parallelFactor /= 2.0;

      if (parallelFactor < 0.0)
         UTIL_THROW("Negative base factor");

      return prefactor_[type]*externalParameter_*tanh(clipcos)*parallelFactor;
   }

   /* 
   * Calculate external force for a single atom.
   */
   inline 
   void LocalLamellarOrderingExternal::getForce(const Vector& position, int type, 
                                     Vector& force) const
   {
      double dPerp, dParallel, perpLength, parallelLength, q, clipParameter, arg, clipcos, tanH, sechSq;
      double parallelFactor1, parallelFactor2, parallelFactor, forceParallelSech1, forceParallelSech2;
      dPerp = position[perpDirection_];
      dParallel = position[parallelDirection_];
      force.zero();

      Vector lengths;
      lengths = boundaryPtr_->lengths();
      perpLength = lengths[perpDirection_];
      parallelLength = lengths[parallelDirection_];

      q = (2.0*M_PI*periodicity_)/perpLength;
      clipParameter  = 1.0/(q*width_*perpLength);
      arg = q*dPerp;
      clipcos = clipParameter*cos(arg);
      tanH = tanh(clipcos);
      sechSq = (1.0 - tanH*tanH);
      parallelFactor1 = tanh( 2.0*M_PI*((dParallel - (((1.0-fraction_)/2.0)*parallelLength))/(fraction_*parallelLength)) );
      parallelFactor2 = tanh( 2.0*M_PI*(((((1+fraction_)/2.0)*parallelLength) - dParallel)/(fraction_*parallelLength)) );
      parallelFactor = 1.0 + (parallelFactor1*parallelFactor2);
      parallelFactor /= 2.0;

      if (parallelFactor < 0.0) {
         UTIL_THROW("Negative base factor");
      } else {
         forceParallelSech1 = (1.0 - (parallelFactor1*parallelFactor1));
         forceParallelSech1 *= 1.0*M_PI*(1/(fraction_*parallelLength))*parallelFactor2;
         forceParallelSech2 = (1.0 - (parallelFactor2*parallelFactor2));
         forceParallelSech2 *= -1.0*M_PI*(1/(fraction_*parallelLength))*parallelFactor1;
         force[parallelDirection_] = -1.0*prefactor_[type]*externalParameter_*tanh(clipcos);
         force[parallelDirection_] *= (forceParallelSech1 + forceParallelSech2);
         force[perpDirection_] = prefactor_[type]*externalParameter_*sechSq*clipParameter*sin(arg)*q*parallelFactor;
      }

   }
 
}
#endif