From 727d240a3a9446f4169e127e8d1e6062ae4c9b86 Mon Sep 17 00:00:00 2001
From: Jim Lutsko <jim@lutsko.com>
Date: Mon, 21 Oct 2024 16:27:09 +0200
Subject: [PATCH] synch

---
 include/DFT.h           |  5 ++++-
 src/DFT.cpp             | 21 ++++++++++++++++++++-
 src/FMT.cpp             | 11 ++++++-----
 src/FMT_Species_EOS.cpp |  2 --
 4 files changed, 30 insertions(+), 9 deletions(-)

diff --git a/include/DFT.h b/include/DFT.h
index a696af1..0897523 100644
--- a/include/DFT.h
+++ b/include/DFT.h
@@ -134,6 +134,7 @@ class DFT : public Dynamical_Matrix
   double get_f_ext() const { return F_ext_;} //external field contribution to free energy including chemical potential
   double get_f_hs()  const { return F_hs_;}  //Hard-sphere contribution to free energy
   double get_f_mf()  const { return F_mf_;} //mean field contribution to free energy
+  double get_f_eos()  const { return F_eos_;} 
 
   // Implement Dynamical_Matrix interface.
   // Second derivatives contracted into arbitrary vector
@@ -167,7 +168,8 @@ class DFT : public Dynamical_Matrix
     ar & F_id_;
     ar & F_ext_;
     ar & F_hs_;
-    ar & F_mf_;	
+    ar & F_mf_;
+    ar & F_eos_;
   }
 
  protected:
@@ -180,6 +182,7 @@ class DFT : public Dynamical_Matrix
   double F_ext_ = 0.0; ///< External field contribution to free energy (including chemical potential)
   double F_hs_  = 0.0; ///< Hard-sphere contribution to free energy
   double F_mf_  = 0.0; ///< Mean-field contribution to free energy
+  double F_eos_ = 0.0;
   
   mutable bool full_hessian_ = true; // used in matrix_holder to distinguish full hessian from excess hessian.
   mutable double rms_force_ = 0.0;   // updated everytime the forces are calculated.
diff --git a/src/DFT.cpp b/src/DFT.cpp
index ccfd011..dcd5331 100644
--- a/src/DFT.cpp
+++ b/src/DFT.cpp
@@ -228,6 +228,26 @@ double DFT::calculateFreeEnergyAndDerivatives_internal_(bool onlyFex)
       species->doFFT();
   }
 
+  // EOS term
+  F_eos_ = 0;
+  for(auto &species : allSpecies_)
+    {
+      FMT_Species_EOS *eos_species = dynamic_cast<FMT_Species_EOS*>(species);
+      if(eos_species)
+	{
+	  const Density& density = species->getDensity();	  
+	  double dV = density.dV();
+	  long Ntot = density.Ntot();
+
+	  double F = 0;	  
+	  for(long I=0;I<Ntot;I++)
+	    F += eos_species->dfex(I, fmt_);
+	  F_eos_ += F*dV;      	  
+	}
+    }
+  F += F_eos_;
+
+  
   //< Mean field contribution to F and dF
   F_mf_ = 0;
   for(auto &interaction: DFT::Interactions_)    
@@ -242,7 +262,6 @@ double DFT::calculateFreeEnergyAndDerivatives_internal_(bool onlyFex)
   
   for(auto &species: allSpecies_)
     F_ext_ += species->evaluate_contribution_chemical_potential();
-
   
   for(auto &field : external_fields_)
     F_ext_ += allSpecies_[field->get_species()]->evaluate_external_field(*field);
diff --git a/src/FMT.cpp b/src/FMT.cpp
index 46f898b..d72cf19 100644
--- a/src/FMT.cpp
+++ b/src/FMT.cpp
@@ -354,14 +354,15 @@ double FMT::calculateFreeEnergy(vector<Species*> &allSpecies)
     throw Eta_Too_Large_Exception();
 
   // Add in contributions from any EOS corrections
-  double FEOS = 0;
+  /*  double FEOS = 0;
   for(auto s: allSpecies)
     {
       FMT_Species_EOS *eos_species = dynamic_cast<FMT_Species_EOS*>(s);
       if(eos_species)
 	FEOS += EOS_Correction(*eos_species);
     }
-    
+
+  */
   // For the AO species, there is additional work to do for both the free energy and the forces. 
   // Do FFT of density and compute the fundamental measures by convolution
   double FAO = 0;
@@ -372,9 +373,9 @@ double FMT::calculateFreeEnergy(vector<Species*> &allSpecies)
 	  FAO += fao_species->free_energy_post_process(needsTensor());
     }
 
-  return F.sum()+FEOS+FAO;
+  return F.sum()+FAO; //FEOS+FAO;
 }
-
+/*
 double FMT::EOS_Correction(FMT_Species_EOS &eos_species)
 {
   double F = 0;  
@@ -384,7 +385,7 @@ double FMT::EOS_Correction(FMT_Species_EOS &eos_species)
       F += eos_species.dfex(I, this);  
   return F;
 }
-
+*/
 // Calculate dF[i] = dPhi/drho(i)
 //                 = sum_j dV * dPhi(j)/drho(i)
 //                 = sum_j dV * sum_alpha dPhi(j)/dn_{alpha}(j)  dn_{alpha}(j)/drho(i)
diff --git a/src/FMT_Species_EOS.cpp b/src/FMT_Species_EOS.cpp
index 6b88276..3cc86a1 100644
--- a/src/FMT_Species_EOS.cpp
+++ b/src/FMT_Species_EOS.cpp
@@ -169,8 +169,6 @@ void FMT_Species_EOS::add_second_derivative(const DFT_FFT &v, DFT_Vec &d2F, cons
   // N.B. the fmt weights have all necessary normalization factors built in, including dV
   bool bConjugate = false;
   convolute_eos_eta_weight_with(v, Psi, bConjugate);
-  for(int i=0;i<1;i++)
-    cout << "i = " << i << " r = " << i*0.05 << " Psi[" << i << "] = " << Psi.Real().get(i) << endl; 
 
   // Get Lambda: Lambda(K) = (d2dfex(K)/deta(K) deta(K))psi(K)  
   DFT_FFT Lambda(Nx,Ny,Nz);