HarmonicAngle.h

#ifndef SIMP_HARMONIC_ANGLE_H
#define SIMP_HARMONIC_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/param/ParamComposite.h>  // base class
#include <simp/interaction/angle/BendForce.h>      // used in inline function
#include <util/global.h>
#include <cmath>

namespace Util {class Random;}

namespace Simp
{

   using namespace Util;

   class HarmonicAngle : public ParamComposite
   {
   
   public:

      HarmonicAngle();
   
      HarmonicAngle(const HarmonicAngle& other);
   
      HarmonicAngle& operator = (const HarmonicAngle& 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& b1, const Vector& b2,
                       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];   
      int  nAngleType_;                 

   };
   
   // Inline method definitions
   
   /*
   * Return angle energy.
   */
   inline double HarmonicAngle::energy(double cosTheta, int type) const
   {
      double dTheta = std::acos(cosTheta) - theta0_[type];
      return 0.5*kappa_[type]*dTheta*dTheta;
   }

   /* 
   * Return:
   *    F1 = d energy / d(b1)
   *    F2 = d energy / d(b2)
   * for use in MD and stress calculation.
   */
   inline
   void HarmonicAngle::force(const Vector& b1, const Vector& b2,
                             Vector& F1, Vector& F2, int type) const
   {
      BendForce bend;
      bend.computeDerivatives(b1, b2);
      double s = bend.sinTheta();
      if (s > 1.0E-10) {
         double factor = kappa_[type]*(theta0_[type] - bend.theta())/s;
         F1.multiply(bend.d1, factor);
         F2.multiply(bend.d2, factor);
      } else {
         F1.zero();
         F2.zero();
      }
   }

}
 
#endif