HarmonicBond.h

#ifndef SIMP_HARMONIC_BOND_H
#define SIMP_HARMONIC_BOND_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/param/ParamComposite.h>  // base class

#include <util/random/Random.h>           // Util namespace 

#include <cmath>

namespace Simp
{

   using namespace Util;

   //class Random;
   
   class HarmonicBond : public ParamComposite 
   {
   
   public:

      HarmonicBond();
   
      HarmonicBond(const HarmonicBond& other);
   
      HarmonicBond& operator = (const HarmonicBond& other);
  
      // Default C++ destructor.
 
      void setNBondType(int nBondType);

      void readParameters(std::istream &in);

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

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

      double energy(double rSq, int type) const;
   
      double forceOverR(double rSq, int type) const;

      double randomBondLength(Random *random, double beta, int type) const;
   
      void set(std::string name, int type, double value);

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

      double kappa(int type) const;

      double length(int type) const;

   private:
 
      static const int MaxNBondType = 8;
   
      double kappa_[MaxNBondType];  
      double length_[MaxNBondType]; 
      int nBondType_;               

   };

   // Inline method definitions

   /*
   * Return bond spring constant.
   */
   inline double HarmonicBond::kappa(int type) const
   {
      assert(type >= 0);
      assert(type < nBondType_);
      return kappa_[type];
   }

   /*
   * Return equilibrium bond length.
   */
   inline double HarmonicBond::length(int type) const
   {
      assert(type < nBondType_);
      return length_[type];
   }

   /* 
   * Return interaction energy for one bond
   */
   inline double HarmonicBond::energy(double rSq, int type) const
   {
      double dl = sqrt(rSq) - length_[type];
      return 0.5*kappa_[type]*dl*dl;
   }

   /* 
   * Return force / distance for one bond, for use in MD
   */
   inline 
   double HarmonicBond::forceOverR(double rSq, int type) const
   {
      double r  = sqrt(rSq);
      double dl = length_[type] - r;
      return kappa_[type]*dl/r;
   }
  
} 
#endif