CosineSqAngle.h

#ifndef SIMP_COSINE_SQ_ANGLE_H
#define SIMP_COSINE_SQ_ANGLE_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 <util/global.h>
#include <util/space/Vector.h>
#include <util/param/ParamComposite.h>  // base class

#include <cmath>

namespace Util {class Random;}

namespace Simp
{

   using namespace Util;

   class CosineSqAngle : public ParamComposite 
   {
   
   public:

      CosineSqAngle();
   
      CosineSqAngle(const CosineSqAngle& other);
   
      CosineSqAngle& operator = (const CosineSqAngle& other);
  
      // Default C++ destructor.
 
      void setNAngleType(int nAngleType);
       
      void readParameters(std::istream &in);

      virtual void loadParameters(Serializable::IArchive &ar);

      virtual void save(Serializable::OArchive &ar);

      void set(std::string name, int type, double value);

      double energy(double cosTheta, int type) const;
 
      void force(const Vector& R1, const Vector& R2,
                       Vector& F1, Vector& F2, int type) const;

      double randomAngle(Random *random, double beta, int type) const;

      double randomCosineAngle(Random *random, double beta, int type) const;


      double get(std::string name, int type) const;

      std::string className() const;
 
   private:
   
      static const int MaxNAngleType = 4;
   
      double kappa_[MaxNAngleType];      
      double theta0_[MaxNAngleType];     
      double cosTheta0_[MaxNAngleType];  
      int    nAngleType_;                

   };
   
   // Inline method definitions
   
   /* 
   * Return angle energy.
   */
   inline double CosineSqAngle::energy(double cosTheta, int type) const
   {
      return ( kappa_[type]
              *(cosTheta-cosTheta0_[type])
              *(cosTheta-cosTheta0_[type])*0.5 );
   }

   /* 
   * Return:
   *    F1 = d energy / d(R1)
   *    F2 = d energy / d(R2)
   * for use in MD and stress calculation.
   */
   inline
   void CosineSqAngle::force(const Vector& R1, const Vector& R2,
                   Vector& F1, Vector& F2, int type) const
   {
      Vector u1 = R1;
      Vector u2 = R2;
      double r1 = R1.abs();
      double r2 = R2.abs();
      double cosTheta;
     
      u1 /= r1;
      u2 /= r2;
      cosTheta = u1.dot(u2);

      F1.multiply(u1, cosTheta);
      F1 -= u2;
      F1 *= kappa_[type]/r1*(cosTheta0_[type]-cosTheta);

      F2.multiply(u2, cosTheta);
      F2 -= u1;
      F2 *= kappa_[type]/r2*(cosTheta0_[type]-cosTheta);
   }
  
}
 
#endif