CosineDihedral.h

#ifndef SIMP_COSINE_DIHEDRAL_H
#define SIMP_COSINE_DIHEDRAL_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 "Torsion.h"                    // used in in-line function
#include "TorsionForce.h"               // used in in-line function

//#include <cmath>

namespace Simp
{

   using namespace Util;

   class CosineDihedral : public ParamComposite 
   {
   
   public:

      CosineDihedral();
   
      CosineDihedral(const CosineDihedral& other);
   
      CosineDihedral& operator = (const CosineDihedral& other);
  
      // Default C++ destructor.
 
      void setNDihedralType(int nDihedralType);
       
      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(const Vector& b1, const Vector& b2, const Vector& b3,
          int type) const;
 
      void force(const Vector& b1, const Vector& b2, const Vector& b3,
                 Vector& f1, Vector& f2, Vector& f3, int type) const;

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

      std::string className() const;
 
   private:
   
      static const int MaxNDihedralType = 8;
   
      double kappa_[MaxNDihedralType];  
      int    nDihedralType_;            

   };
   
   // Inline method definitions
   
   /* 
   * Return dihedral energy.
   */
   inline
   double CosineDihedral::energy(const Vector& b1, const Vector& b2,
          const Vector& b3, int type) const
   {
      Torsion torsion;
      bool status;
      status = torsion.computeAngle(b1, b2, b3); // computes cosPhi
      if (!status) {
         return (kappa_[type] * (1.0 + torsion.cosPhi));
      } else {
         return 0.0;
      }
   }

   /* 
   * Return:
   *    f1 = d energy / d(b1)
   *    f2 = d energy / d(b2)
   *    f3 = d energy / d(b3)
   * for use in MD and stress calculation.
   */
   inline
   void CosineDihedral::force(const Vector& b1, const Vector& b2,
        const Vector& b3, Vector& f1, Vector& f2, Vector& f3, int type) const
   {
      TorsionForce torsion;
      bool status; // Error code, 0 is normal, 1 is error
      status = torsion.computeDerivatives(b1, b2, b3);

      if (!status) {
         double dEdCosPhi = kappa_[type];
         f1.multiply(torsion.d1, dEdCosPhi);
         f2.multiply(torsion.d2, dEdCosPhi);
         f3.multiply(torsion.d3, dEdCosPhi);
      } else {
         f1.zero();
         f2.zero();
         f3.zero();
      }
   }

} 
#endif