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

@@ -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;
 }
 
 ////////////////////////////////////////////////////////////////////////////////

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_; }