synch

jimlutsko · Oct 22, 2024 · 1f583f1 · 1f583f1
1 parent e3639d9
commit 1f583f1
Show file tree

Hide file tree

Showing 4 changed files with 34 additions and 7 deletions.
diff --git a/include/DFT.h b/include/DFT.h
@@ -57,7 +57,9 @@ class DFT : public Dynamical_Matrix
   const Lattice& get_lattice() const {return allSpecies_.front()->getLattice();}
   const Density& getDensity(int species) const {return allSpecies_[species]->getDensity();}
   const Density& get_density(int species) const {return allSpecies_[species]->getDensity();}
-  Interaction_Base* getInteraction(int which) { return Interactions_[which];}
+  Interaction_Base* getInteraction(int which) { return (Interactions_.size() == 0 ? NULL : Interactions_[which]);}
+  const int get_num_interactions() const { return Interactions_.size();}
+
 
   Species *getSpecies(int s) const { return allSpecies_[s];}
   int      getNumberOfSpecies() const {return allSpecies_.size();}
@@ -102,7 +104,7 @@ class DFT : public Dynamical_Matrix
   double calculateFreeEnergyAndDerivatives(bool onlyFex = false, bool H_dot_Force = false);
   double evaluate(bool onlyFex = false, bool H_dot_Force = false) { return calculateFreeEnergyAndDerivatives(onlyFex,H_dot_Force);}  
 
-  // Bulk Thermodynamics
+  // Bulk properties
 
   virtual double mu_times_beta(const vector<double> &densities, int species) const;
   virtual double omega_times_beta_over_volume(const vector<double> &x) const;

diff --git a/include/Species.h b/include/Species.h
@@ -631,10 +631,15 @@ class FMT_Species_EOS : public FMT_Species
 
   virtual bool is_eos() const { return true;}  
 
+  double get_eos_hsd() const { return hsd_*D_EOS_;}
+
   // pass through from the EOS object
   double get_bulk_dfex(double x, const void *param) const;
   double get_bulk_ddfex_deta(double x, const void *param) const;
   double get_bulk_d2dfex_deta2(double x, const void *param) const;
+
+  double get_bulk_ddfex_dx(double x, const void *param) const {return get_bulk_ddfex_deta(x,param)*(M_PI*pow(D_EOS_*hsd_,3))/6;}
+
 
   // Evaluate EOS free energy functional at point I
   double effDensity(long I);

diff --git a/src/DFT.cpp b/src/DFT.cpp
@@ -28,13 +28,26 @@ double DFT::real_space_dcf(double r, double x) const
     throw std:: runtime_error("DFT::real_space_dcf only implemented for single species");
 
   double dcf = 0;
-  double hsd = allSpecies_[0]->getHSD();
+  Species *s = allSpecies_[0];
+  double hsd = s->getHSD();
+  double dV = s->getLattice().dV();
 
+  // FMT contribution
   if(fmt_) dcf += fmt_->get_real_space_dcf(r,x,hsd);
 
+  // Interaction contribution
   for(auto &x: Interactions_)
-    dcf -= x->getW(r);
-
+    dcf -= x->getW(r)/dV; // spurious dV in definition of weights
+
+  // EOS correction contribution
+  FMT_Species_EOS *e = dynamic_cast<FMT_Species_EOS*>(s);
+  if(e)
+    {
+      double hsd1 = e->get_eos_hsd();
+      if(r < hsd1)
+	dcf -= (M_PI/12)*(r-hsd1)*(r-hsd1)*(r+2*hsd1)*e->get_bulk_ddfex_dx(x, fmt_);
+    }
+
   return dcf;
 }
 

diff --git a/src/FMT.cpp b/src/FMT.cpp
@@ -692,7 +692,9 @@ double FMT::d2Phi_dn_dn(int I[3], int si, int J[3], int sj, vector<Species*> &al
 
 double FMT::BulkMuex(const vector<double> &x, const vector<Species*> &allSpecies, int species) const
 {
-  if(allSpecies[0]->is_eos()) throw std::runtime_error("FMT::BulkMuex Not implemented for eos correction species");
+  //if(allSpecies[0]->is_eos()) throw std::runtime_error("FMT::BulkMuex Not implemented for eos correction species");
+
+  double mu = 0.0;
 
   FundamentalMeasures fm(0.0,1.0);
 
@@ -704,6 +706,11 @@ double FMT::BulkMuex(const vector<double> &x, const vector<Species*> &allSpecies
 	  FundamentalMeasures fm1(x[s],sfmt->getHSD());
 	  fm.add(fm1);
 	}
+
+      FMT_Species_EOS *eosfmt = dynamic_cast<FMT_Species_EOS*>(allSpecies[s]);
+      if(eosfmt)
+	  mu += eosfmt->get_bulk_ddfex_dx(x[s],this);
+
     }
   FundamentalMeasures dPhi;
   calculate_dPhi_wrt_fundamental_measures(fm, dPhi);
@@ -712,7 +719,7 @@ double FMT::BulkMuex(const vector<double> &x, const vector<Species*> &allSpecies
   FundamentalMeasures weights(1.0, allSpecies[species]->getHSD()); // a non-fmt species gives zero hsd so OK.
   vector<double> weights_v = weights.get_as_vector();
 
-  double mu = 0.0;
+
   for(int i=0;i<dPhi_v.size();i++) mu += dPhi_v[i]*weights_v[i];
 
   // I now assume there is only one species as I do not (now, yet) have the generalization