DpdPair.h

#ifndef SIMP_DPD_PAIR_H
#define SIMP_DPD_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 DpdPair : public ParamComposite 
   {
   
   public:
   
      DpdPair();

      DpdPair(const DpdPair& other);

      DpdPair& operator = (const DpdPair& 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 get(std::string name, int i, int j) const;



      double epsilon(int i, int j) const;
 
      double sigma(int i, int j) const;
 

   private:
   
      static const int MaxAtomType = 4;
   
      // Parameters for different types of particle pairs
      double epsilon_[MaxAtomType][MaxAtomType];   
      double sigma_[MaxAtomType][MaxAtomType];     
      double sigmaSq_[MaxAtomType][MaxAtomType];   
      double ce_[MaxAtomType][MaxAtomType];        
      double cf_[MaxAtomType][MaxAtomType];        
 
      double maxPairCutoff_;

      int    nAtomType_; 

      bool  isInitialized_;

   };
  
   // inline methods 
 
   /* 
   * Calculate interaction energy for a pair, as function of squared distance.
   */
   inline double DpdPair::energy(double rsq, int i, int j) const 
   {
      double dr;
      if (rsq < sigmaSq_[i][j]) {
         dr = sqrt(rsq) - sigma_[i][j];
         return ce_[i][j]*dr*dr;
      } else {
         return 0.0;
      }
   }
   
   /* 
   * Calculate force/distance for a pair as function of squared distance.
   */
   inline double DpdPair::forceOverR(double rsq, int i, int j) const
   {
      if (rsq < sigmaSq_[i][j]) {
         return cf_[i][j]*(sigma_[i][j]/sqrt(rsq) - 1.0);
      } else {
         return 0.0;
      }
   }

   /* 
   * Calculate force/distance for a pair as function of squared distance.
   */
   inline double DpdPair::cutoffSq(int i, int j) const
   {  return sigmaSq_[i][j]; }

}
#endif