NucleationExternal.h

#ifndef SIMP_NUCLEATION_EXTERNAL_H
#define SIMP_NUCLEATION_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 NucleationExternal : public ParamComposite 
   {
   
   public:
   
      NucleationExternal();

      NucleationExternal(const NucleationExternal& other);

      NucleationExternal& operator = (const NucleationExternal& 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_;

      Vector shift_;

      double C_;

      int periodicity_;

      double interfaceWidth_;

      double nucleationClip_;

      double bias_;

      Boundary *boundaryPtr_;
   
      int    nAtomType_; 

      bool  isInitialized_;

   };
  
   // inline methods 

   /* 
   * Calculate external potential energy for a single atom.
   */
   inline double NucleationExternal::energy(const Vector& position, int type) const
   {
      const Vector cellLengths = boundaryPtr_->lengths();
      double clipParameter = 1.0/(2.0*M_PI*periodicity_*interfaceWidth_);
      double e;
      
      Vector r = position;
      r -= shift_;
      double cosine = 0.0;
      for (int i = 0; i < nWaveVectors_; ++i) {
         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 += cos(arg);
      }

      cosine *= clipParameter;
      e = prefactor_[type]*externalParameter_*tanh(C_+cosine) * 
      (tanh(nucleationClip_*(-bias_+cos(2.0*M_PI*position[0]/cellLengths[0]+acos(bias_))))+1)/2 *
      (tanh(nucleationClip_*(-bias_+cos(2.0*M_PI*position[1]/cellLengths[1]+acos(bias_))))+1)/2 * 
      (tanh(nucleationClip_*(-bias_+cos(2.0*M_PI*position[2]/cellLengths[2]+acos(bias_))))+1)/2 ;

      return e;
   }

   /* 
   * Calculate external force for a single atom.
   */
   inline 
   void NucleationExternal::getForce(const Vector& position, int type, 
                                     Vector& force) const
   {
      const Vector cellLengths = boundaryPtr_->lengths();
      double clipParameter = 1.0/(2.0*M_PI*periodicity_*interfaceWidth_);
      double ne;

      ne = 
      (tanh(nucleationClip_*(-bias_+cos(2.0*M_PI*position[0]/cellLengths[0]+acos(bias_))))+1)/2 * 
      (tanh(nucleationClip_*(-bias_+cos(2.0*M_PI*position[1]/cellLengths[1]+acos(bias_))))+1)/2 * 
      (tanh(nucleationClip_*(-bias_+cos(2.0*M_PI*position[2]/cellLengths[2]+acos(bias_))))+1)/2 ;

      Vector r = position;
      r -= shift_;
      double cosine = 0.0;
      Vector deriv;
      deriv.zero();
      for (int i = 0; i < nWaveVectors_; ++i) {
         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 += cos(arg);
         double sine = -1.0*sin(arg);
         q *= sine;
         deriv += q;
      }
      cosine *= clipParameter;

      deriv *= clipParameter;
      double tanH = tanh(C_+cosine);
      double sechSq = (1.0 - tanH*tanH);
      double f = prefactor_[type]*externalParameter_*sechSq;
      deriv *= -1.0*f;
      deriv *= ne;
      deriv[0] = deriv[0] + nucleationClip_*(M_PI/cellLengths[0]*sin(2.0*M_PI*position[0]/cellLengths[0]+acos(bias_)))/
                            cosh(2.0*M_PI*position[0]/cellLengths[0]+acos(bias_))/
                            cosh(2.0*M_PI*position[0]/cellLengths[0]+acos(bias_))*
                            (tanh(nucleationClip_*(-bias_+cos(2.0*M_PI*position[1]/cellLengths[1]+acos(bias_))))+1)/2*
                            (tanh(nucleationClip_*(-bias_+cos(2.0*M_PI*position[2]/cellLengths[2]+acos(bias_))))+1)/2;

      deriv[1] = deriv[1] + nucleationClip_*(M_PI/cellLengths[1]*sin(2.0*M_PI*position[1]/cellLengths[1]+acos(bias_)))/
                            cosh(2.0*M_PI*position[1]/cellLengths[1]+acos(bias_))/
                            cosh(2.0*M_PI*position[1]/cellLengths[1]+acos(bias_))*
                            (tanh(nucleationClip_*(-bias_+cos(2.0*M_PI*position[0]/cellLengths[0]+acos(bias_))))+1)/2*
                            (tanh(nucleationClip_*(-bias_+cos(2.0*M_PI*position[2]/cellLengths[2]+acos(bias_))))+1)/2;

      deriv[2] = deriv[2] + nucleationClip_*(M_PI/cellLengths[2]*sin(2.0*M_PI*position[2]/cellLengths[2]+acos(bias_)))/
                            cosh(2.0*M_PI*position[2]/cellLengths[2]+acos(bias_))/
                            cosh(2.0*M_PI*position[2]/cellLengths[2]+acos(bias_))*
                            (tanh(nucleationClip_*(-bias_+cos(2.0*M_PI*position[0]/cellLengths[0]+acos(bias_))))+1)/2*
                            (tanh(nucleationClip_*(-bias_+cos(2.0*M_PI*position[1]/cellLengths[1]+acos(bias_))))+1)/2;
      force = deriv;
   }
 
}
#endif