FeneBond.h

#ifndef SIMP_FENE_BOND_H
#define SIMP_FENE_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>
#include <sstream>

namespace Simp
{

   using namespace Util;

   class FeneBond : public ParamComposite 
   {
   
   public:

      FeneBond();
   
      FeneBond(const FeneBond& other);
   
      FeneBond& operator = (const FeneBond& other);
  
      // Default destructor.
 
      void setNBondType(int nBondType);
       
      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 rSq, int type) const;
   
      double forceOverR(double rSq, int type) const;

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

   private:
   
      static const int MaxNBondType = 2;
  
      double kappa_[MaxNBondType];

      double r0_[MaxNBondType];  

      double r0Sq_[MaxNBondType];  

      double r0SqInv_[MaxNBondType];  

      double ce_[MaxNBondType];  

      double rlSq_[MaxNBondType];

      double rl_[MaxNBondType];
      
      double  energyCutoff_[MaxNBondType];  
      
      // double a_[MaxNBondType];
      
      double  forceCutoff_;
      
      int    nBondType_;

   };
   
   // Inline method definitions
   
   /* 
   * Return interaction energy for one bond
   */
   inline double FeneBond::energy(double rSq, int type) const
   {
     if (rSq < rlSq_[type]){     
       double g = 1.0 - rSq*r0SqInv_[type];
       return ce_[type]*log(g);
     } else {
       //double b = rl_[type]*kappa_[type]/a_[type];
       return energyCutoff_[type] + forceCutoff_*(sqrt(rSq) - rl_[type]);  
     } 
   }

   /* 
   * Return force / distance for one bond, for use in MD
   */
   inline 
   double FeneBond::forceOverR(double rSq, int type) const
   {
     if (rSq < rlSq_[type]){ 
       double g = 1.0 - rSq*r0SqInv_[type];
       return -kappa_[type]/g;
     } else {
       if (rSq > r0Sq_[type]) {
          std::cout << "Warning: Long FENE bond: Rsq = " << rSq << std::endl;
       }
       //return -rl_[type]*kappa_[type]/(a_[type]*sqrt(rSq));  
       return -forceCutoff_/sqrt(rSq);  
     }              
   }

} 
#endif