CosineAngle.h

#ifndef SIMP_COSINE_ANGLE_H
#define SIMP_COSINE_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 CosineAngle : public ParamComposite 
   {
   
   public:

      CosineAngle();
   
      CosineAngle(const CosineAngle& other);
   
      ~CosineAngle();
   
      CosineAngle& operator = (const CosineAngle& other);
  
      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];  
      int    nAngleType_;            

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

   /* 
   * Return:
   *    F1 = d energy / d(R1)
   *    F2 = d energy / d(R2)
   * for use in MD and stress calculation.
   */
   inline
   void CosineAngle::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;

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