From 6e2b9157f27e5f80e2f1b58fb999e1c9ab64b354 Mon Sep 17 00:00:00 2001
From: Francois Gygi <fgygi@ucdavis.edu>
Date: Wed, 31 Dec 2014 01:50:21 +0000
Subject: [PATCH] Redefined enthalpy to include electric enthalpy and p*V

git-svn-id: http://qboxcode.org/svn/qb/branches/efield@1616 cba15fb0-1239-40c8-b417-11db7ca47a34
---
 src/EnergyFunctional.C | 51 ++++++++++++++++++++++++++++++++-------------------
 src/EnergyFunctional.h |  5 ++++-
 2 files changed, 36 insertions(+), 20 deletions(-)

diff --git a/src/EnergyFunctional.C b/src/EnergyFunctional.C
index 9337ec4..1718290 100644
--- a/src/EnergyFunctional.C
+++ b/src/EnergyFunctional.C
@@ -629,11 +629,9 @@ double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
   if ( el_enth_ )
   {
     tmap["el_enth_energy"].start();
-    eefield_ = el_enth_->energy(dwf,compute_hpsi);
+    eefield_ = el_enth_->enthalpy(dwf,compute_hpsi);
     tmap["el_enth_energy"].stop();
-
-    // add energy of ionic dipole in e_field
-    eefield_ += s_.atoms.dipole() * s_.ctrl.e_field;
+    enthalpy_ += eefield_;
 
     if ( compute_forces )
     {
@@ -648,7 +646,19 @@ double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
     }
   }
   etotal_ = ekin_ + econf_ + eps_ + enl_ + ecoul_ + exc_ +
-            ets_ + eexf_ + exhf_ + eefield_;
+            ets_ + eexf_ + exhf_;
+  enthalpy_ = etotal_;
+
+  if ( el_enth_ )
+    enthalpy_ += eefield_;
+
+  if ( compute_stress )
+  {
+    valarray<double> sigma_ext(s_.ctrl.ext_stress,6);
+    const double pext = (sigma_ext[0]+sigma_ext[1]+sigma_ext[2])/3.0;
+    epv_ = pext * omega;
+    enthalpy_ += epv_;
+  }
 
   if ( compute_hpsi )
   {
@@ -789,8 +799,10 @@ double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
   }
 
   if ( compute_stress )
+  {
     sigma = sigma_ekin + sigma_econf + sigma_eps + sigma_enl +
             sigma_ehart + sigma_exc + sigma_esr;
+  }
 
   if ( debug_stress && s_.ctxt_.onpe0() )
   {
@@ -1074,20 +1086,21 @@ void EnergyFunctional::print(ostream& os) const
   os.setf(ios::fixed,ios::floatfield);
   os.setf(ios::right,ios::adjustfield);
   os << setprecision(8);
-  os << "  <ekin>   " << setw(15) << ekin()   << " </ekin>\n"
-     << "  <econf>  " << setw(15) << econf()  << " </econf>\n"
-     << "  <eps>    " << setw(15) << eps()    << " </eps>\n"
-     << "  <enl>    " << setw(15) << enl()    << " </enl>\n"
-     << "  <ecoul>  " << setw(15) << ecoul()  << " </ecoul>\n"
-     << "  <exc>    " << setw(15) << exc()    << " </exc>\n"
-     << "  <exhf>   " << setw(15) << exhf()   << " </exhf>\n"
-     << "  <esr>    " << setw(15) << esr()    << " </esr>\n"
-     << "  <eself>  " << setw(15) << eself()  << " </eself>\n"
-     << "  <ets>    " << setw(15) << ets()    << " </ets>\n"
-     << "  <eexf>   " << setw(15) << eexf()   << " </eexf>\n"
-     << "  <eefield>" << setw(15) << eefield()<< " </eefield>\n"
-     << "  <etotal> " << setw(15) << etotal() << " </etotal>\n"
-     << flush;
+  os << "  <ekin>    " << setw(15) << ekin()   << " </ekin>\n"
+     << "  <econf>   " << setw(15) << econf()  << " </econf>\n"
+     << "  <eps>     " << setw(15) << eps()    << " </eps>\n"
+     << "  <enl>     " << setw(15) << enl()    << " </enl>\n"
+     << "  <ecoul>   " << setw(15) << ecoul()  << " </ecoul>\n"
+     << "  <exc>     " << setw(15) << exc()    << " </exc>\n"
+     << "  <exhf>    " << setw(15) << exhf()   << " </exhf>\n"
+     << "  <esr>     " << setw(15) << esr()    << " </esr>\n"
+     << "  <eself>   " << setw(15) << eself()  << " </eself>\n"
+     << "  <ets>     " << setw(15) << ets()    << " </ets>\n"
+     << "  <eexf>    " << setw(15) << eexf()   << " </eexf>\n"
+     << "  <etotal>  " << setw(15) << etotal() << " </etotal>\n"
+     << "  <epv>     " << setw(15) << epv() << " </epv>\n"
+     << "  <eefield> " << setw(15) << eefield() << " </eefield>\n"
+     << "  <enthalpy>" << setw(15) << enthalpy() << " </enthalpy>" << endl;
 }
 
 ////////////////////////////////////////////////////////////////////////////////
diff --git a/src/EnergyFunctional.h b/src/EnergyFunctional.h
index bea82a3..c70b0de 100644
--- a/src/EnergyFunctional.h
+++ b/src/EnergyFunctional.h
@@ -68,8 +68,9 @@ class EnergyFunctional
   std::vector<int> na_;
   int namax_;
   int nsp_;
-  double ekin_, econf_, eps_, enl_, ehart_, eefield_,
+  double ekin_, econf_, eps_, enl_, ehart_,
          ecoul_, exc_, esr_, eself_, ets_, eexf_, exhf_, etotal_;
+  double epv_, eefield_, enthalpy_;
   std::valarray<double> sigma_ekin,sigma_econf,sigma_eps,sigma_ehart,sigma_exc,
     sigma_enl, sigma_esr, sigma;
 
@@ -96,6 +97,8 @@ class EnergyFunctional
   double ets(void) const { return ets_; }
   double eexf(void) const { return eexf_; }
   double eefield(void) const { return eefield_; }
+  double epv(void) const { return epv_; }
+  double enthalpy(void) const { return enthalpy_; }
 
   ElectricEnthalpy* el_enth() { return el_enth_; }
 
--
libgit2 0.26.0