synch

jimlutsko · Aug 29, 2024 · a739dcb · a739dcb
1 parent b89e339
commit a739dcb
Show file tree

Hide file tree

Showing 2 changed files with 58 additions and 5 deletions.
diff --git a/include/DFT_Factory.h b/include/DFT_Factory.h
@@ -188,6 +188,7 @@ class DFT_Factory
 
 	if(eos_correction_ == LJ_JZG_EOS)   eos_ = new LJ_JZG(kT_, rcut1_); // need to add no-shift option
 	if(eos_correction_ == EOS_NULL)     eos_ = new EOS_NULL_(kT_); // need to add no-shift option
+	if(eos_correction_ == HS_PY_EOS)    eos_ = new PY_EOS(kT_); // need to add no-shift option
 
 	double avdw = 0;
 	if(potential1_) avdw = potential1_->getVDW_Parameter(kT_);

diff --git a/include/EOS.h b/include/EOS.h
@@ -2,12 +2,14 @@
 #define __LUTSKO_EOS__
 
 #include <ostream>
+#include "Enskog.h"
 
 // Empirical equations of state.
 
 const string LJ_JZG_EOS("LJ_JZG_EOS");
-const string EOS_NULL("EOS_NULL");
 const string LJ_MECKE_EOS("LJ_MECKE_EOS");
+const string HS_PY_EOS("HS_PY_EOS");
+const string EOS_NULL("EOS_NULL");
 
 
 class EOS
@@ -60,7 +62,7 @@ class EOS_NULL_ : public EOS
 
   // Basic functionality
   virtual double get_energy_per_atom(  double density) {return 0.0;}  // F/(NkT) = FE per unit atom (phi)
-  virtual double get_energy_per_volume(double density) {return density*get_energy_per_atom(density);}
+  virtual double get_energy_per_volume(double density) {return 0.0;}
   virtual double get_pressure(double density) {return 0.0;}     // P/nkT
 
   // convenience helpers ...
@@ -69,9 +71,9 @@ class EOS_NULL_ : public EOS
   virtual double get_pressure(double density, double kT)          {kT_ = kT; return get_pressure(density);}
 
   // Excess F_ex/(kT*V) and derivatives
-  virtual double fex(double density)  {return density*phix(density);}
-  virtual double f1ex(double density) {return phix(density) + density*phi1x(density);}
-  virtual double f2ex(double density) {return 2*phi1x(density) + density*phi2x(density);}
+  virtual double fex(double density)  {return 0.0;} //density*phix(density);}
+  virtual double f1ex(double density) {return 0.0;} //phix(density) + density*phi1x(density);}
+  virtual double f2ex(double density) {return 0.0;} //2*phi1x(density) + density*phi2x(density);}
 
   // excess free energy per atom and density derivatives
   virtual double phix(double density)  const {return 0.0;} //{ throw std::runtime_error("phix not implemented in EOS object");}
@@ -88,6 +90,54 @@ class EOS_NULL_ : public EOS
   double kT_ = 1;
 };
 
+class PY_EOS : public EOS
+{
+ public:
+  PY_EOS(double kT, double rc = -1.0, int noShift = 0)
+   : EOS(kT) {}
+
+  virtual ~PY_EOS() {}
+
+  virtual char const * const getName() const { return "PY_EOS";}
+
+  virtual double get_energy_per_atom(double density)
+  {
+    Enskog e(density);
+    return e.freeEnergyPYC();
+  }
+  virtual double get_pressure(double density)
+  {
+    Enskog e(density);
+    return e.pressurePYC();
+  }
+  // excess free energy per atom and density derivatives
+  virtual double phix(double density)  const 
+  {
+    Enskog e(density);
+    return e.exFreeEnergyPYC();
+  }
+
+  virtual double phi1x(double density) const 
+  {
+    Enskog e(density);
+    return e.dexFreeEnergyPYCdRho();
+  }
+
+  virtual double phi2x(double density) const
+  {
+    Enskog e(density);
+    return e.d2exFreeEnergyPYCdRho2();    
+  }
+  virtual double phi3x(double density) const 
+  {
+    Enskog e(density);
+    return e.d3exFreeEnergyPYCdRho3();        
+  }
+
+ private:
+};
+
+
 class LJ_JZG : public EOS
 {
  public:
@@ -259,6 +309,8 @@ class LJ_JZG : public EOS
   double x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,x21,x22,x23,x24,x25,x26,x27,x28,x29,x30,x31,x32;
 };
 
+
+
 class LJ_Mecke : public EOS
 {
  public: