MultiHarmonicDihedral.h

#ifndef SIMP_MULTI_HARMONIC_DIHEDRAL_H
#define SIMP_MULTI_HARMONIC_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/param/ParamComposite.h>     // base class
#include <util/containers/DArray.h>        // member
#include <util/space/Vector.h>             // inline function
#include <simp/interaction/dihedral/Torsion.h>        // inline function
#include <simp/interaction/dihedral/TorsionForce.h>   // inline function
#include <util/global.h>

namespace Simp
{

   using namespace Util;

   class MultiHarmonicDihedral : public ParamComposite 
   {
   
   public:

      MultiHarmonicDihedral();
   
      MultiHarmonicDihedral(const MultiHarmonicDihedral& other);
   
      // Default C++ destructor.
 
      MultiHarmonicDihedral& operator = (const MultiHarmonicDihedral& other);
  
      void setNDihedralType(int nDihedralType);
       
      void readParameters(std::istream &in);

      void writeParam(std::ostream& out);

      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:
   
      struct Parameter 
      {

      public:

         // Coefficients in expansion V(phi) = \sum_m k_m cos(m phi)
         double k0; 
         double k1; 
         double k2;
         double k3; 
         double k4;

         // Coefficients in expansion V(phi) = \sum_m a_m cos^m(phi)
         double a0;
         double a1; 
         double a2; 
         double a3; 
         double a4; 

         double g2; 
         double g3; 
         double g4; 

         void init()
         {
            a0 = k0 - k2 + k4;
            a1 = k1 - 3.0*k3;
            a2 = 2.0*k2 - 8.0*k4;
            a3 = 4.0*k3;
            a4 = 8.0*k4;

            g2 = 2.0*a2;
            g3 = 3.0*a3;
            g4 = 4.0*a4;
         }

      };

      DArray<MultiHarmonicDihedral::Parameter> parameters_;

      int nDihedralType_;        

   };
 
   // Inline method definitions
   
   /* 
   * Return dihedral energy.
   */
   inline
   double MultiHarmonicDihedral::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) { 
         double c = torsion.cosPhi;
         const Parameter* p = &parameters_[type];
         return (p->a0 + c*(p->a1 + c*(p->a2 + c*(p->a3 + c*p->a4))));
      } else {
         return 0.0;
      }
   }

   /* 
   * Compute and return vectors:
   *
   *    f1 = d energy / d(b1)
   *    f2 = d energy / d(b2)
   *    f3 = d energy / d(b3)
   *
   * for use in MD and stress calculation.
   */
   inline
   void MultiHarmonicDihedral::force(const Vector& b1, const Vector& b2,
        const Vector& b3, Vector& f1, Vector& f2, Vector& f3, int type) const
   {
      TorsionForce torsion;
      bool status;
      status = torsion.computeDerivatives(b1, b2, b3);
      if (!status) {
         double c = torsion.cosPhi;
         const Parameter* p = &parameters_[type];
         double dEdC = p->a1 + c*(p->g2 + c*(p->g3 + c*p->g4));
         f1.multiply(torsion.d1, dEdC);
         f2.multiply(torsion.d2, dEdC);
         f3.multiply(torsion.d3, dEdC);
      } else {
         f1.zero();
         f2.zero();
         f3.zero();
      }
   }

} 
#endif