LJPair.h

#ifndef SIMP_LJ_PAIR_H
#define SIMP_LJ_PAIR_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>
#include <util/global.h>

#include <math.h>

namespace Simp
{

   using namespace Util;

   class LJPair : public ParamComposite 
   {
   
   public:
   
      LJPair();

      LJPair(const LJPair& other);

      LJPair& operator = (const LJPair& other);



      void setNAtomType(int nAtomType);

      void readParameters(std::istream &in);

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

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

      void setEpsilon(int i, int j, double epsilon);
 
      void setSigma(int i, int j, double sigma);

      void set(std::string name, int i, int j, double value);



      double energy(double rsq, int i, int j) const;
  
      double forceOverR(double rsq, int i, int j) const;
  
      double cutoffSq(int i, int j) const;
 
      double maxPairCutoff() const;

      double epsilon(int i, int j) const;
 
      double sigma(int i, int j) const;
 
      double get(std::string name, int i, int j) const;


   protected:
   
      static const int MaxAtomType = 4;
   
      // Lennard-Jones parameters for different types of atom pairs
      double epsilon_[MaxAtomType][MaxAtomType];  
      double sigma_[MaxAtomType][MaxAtomType];    
      double sigmaSq_[MaxAtomType][MaxAtomType];  
      double cutoff_[MaxAtomType][MaxAtomType];   
      double cutoffSq_[MaxAtomType][MaxAtomType]; 
      double ljShift_[MaxAtomType][MaxAtomType];  
      double eps48_[MaxAtomType][MaxAtomType];    
 
      double maxPairCutoff_;

      int    nAtomType_; 

      bool  isInitialized_;

   };
  
   // Inline methods 
 
   /* 
   * Calculate interaction energy for a pair, as function of squared distance.
   */
   inline double LJPair::energy(double rsq, int i, int j) const 
   {
      double r6i;
      if (rsq < cutoffSq_[i][j]) {
         if (rsq < 0.6*sigmaSq_[i][j]) {
             rsq = 0.6*sigmaSq_[i][j];
         }
         r6i = sigmaSq_[i][j]/rsq;
         r6i = r6i*r6i*r6i;
         return 4.0*epsilon_[i][j]*(r6i*r6i - r6i) + ljShift_[i][j];
      } else {
         return 0.0;
      }
   }
  
   /* 
   * Calculate force/distance for a pair as function of squared distance.
   */
   inline double LJPair::forceOverR(double rsq, int i, int j) const
   {
      double r2i, r6i;
      r2i = 1.0/rsq;
      r6i = sigmaSq_[i][j]*r2i;
      r6i = r6i*r6i*r6i;
      return eps48_[i][j]*(r6i - 0.5)*r6i*r2i;
   }

   /* 
   * Return cutoff parameter for a specific atom type pair.
   */
   inline double LJPair::cutoffSq(int i, int j) const
   {  return cutoffSq_[i][j]; }

}
#endif