Commit 7d2b7219 by Francois Gygi

modif for ElectricEnthalpy


git-svn-id: http://qboxcode.org/svn/qb/branches/efield@1586 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 260a4beb
......@@ -27,6 +27,8 @@
#include "XCOperator.h"
#include "NonLocalPotential.h"
#include "ConfinementPotential.h"
#include "D3vector.h"
#include "ElectricEnthalpy.h"
#include "Timer.h"
#include "blas.h"
......@@ -158,6 +160,11 @@ EnergyFunctional::EnergyFunctional( Sample& s, const ChargeDensity& cd)
wf.sd(0,ikp)->basis());
}
// Electric enthalpy
el_enth_ = 0;
if ( norm(s_.ctrl.e_field) != 0.0 )
el_enth_ = new ElectricEnthalpy(s_);
sf.init(tau0,*vbasis_);
cell_moved();
......@@ -169,6 +176,7 @@ EnergyFunctional::EnergyFunctional( Sample& s, const ChargeDensity& cd)
////////////////////////////////////////////////////////////////////////////////
EnergyFunctional::~EnergyFunctional(void)
{
delete el_enth_;
delete xco;
for ( int ikp = 0; ikp < s_.wf.nkp(); ikp++ )
{
......@@ -301,6 +309,9 @@ void EnergyFunctional::update_vhxc(void)
v_r[1][i] += vloc;
}
}
if ( el_enth_ )
el_enth_->update();
}
////////////////////////////////////////////////////////////////////////////////
......@@ -612,7 +623,32 @@ double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
const double boltz = 1.0 / ( 11605.0 * 2.0 * 13.6058 );
ets_ = - wf_entropy * s_.ctrl.fermi_temp * boltz;
}
etotal_ = ekin_ + econf_ + eps_ + enl_ + ecoul_ + exc_ + ets_ + eexf_ + exhf_;
// Electric enthalpy
eefield_ = 0.0;
if ( el_enth_ )
{
tmap["el_enth_energy"].start();
eefield_ = el_enth_->energy(dwf,compute_hpsi);
tmap["el_enth_energy"].stop();
// add energy of ionic dipole in e_field
eefield_ += s_.atoms.dipole() * s_.ctrl.e_field;
if ( compute_forces )
{
for ( int is = 0; is < nsp_; is++ )
for ( int ia = 0; ia < na_[is]; ia++ )
{
D3vector f = zv_[is] * s_.ctrl.e_field;
fion[is][3*ia] += f.x;
fion[is][3*ia+1] += f.y;
fion[is][3*ia+2] += f.z;
}
}
}
etotal_ = ekin_ + econf_ + eps_ + enl_ + ecoul_ + exc_ +
ets_ + eexf_ + exhf_ + eefield_;
if ( compute_hpsi )
{
......@@ -1049,6 +1085,7 @@ void EnergyFunctional::print(ostream& os) const
<< " <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;
}
......
......@@ -26,8 +26,10 @@
#include <string>
#include "ChargeDensity.h"
#include "StructureFactor.h"
#include "ElectricEnthalpy.h"
#include "Timer.h"
class D3vector;
class Sample;
class Basis;
class AtomSet;
......@@ -51,6 +53,7 @@ class EnergyFunctional
std::vector<FourierTransform*> ft;
StructureFactor sf;
XCOperator* xco;
ElectricEnthalpy* el_enth_;
std::vector<std::vector<NonLocalPotential*> > nlp; // nlp[ispin][ikp]
std::vector<ConfinementPotential*> cfp; // cfp[ikp]
......@@ -65,7 +68,7 @@ class EnergyFunctional
std::vector<int> na_;
int namax_;
int nsp_;
double ekin_, econf_, eps_, enl_, ehart_,
double ekin_, econf_, eps_, enl_, ehart_, eefield_,
ecoul_, exc_, esr_, eself_, ets_, eexf_, exhf_, etotal_;
std::valarray<double> sigma_ekin,sigma_econf,sigma_eps,sigma_ehart,sigma_exc,
sigma_enl, sigma_esr, sigma;
......@@ -92,6 +95,9 @@ class EnergyFunctional
double eself(void) const { return eself_; }
double ets(void) const { return ets_; }
double eexf(void) const { return eexf_; }
double eefield(void) const { return eefield_; }
ElectricEnthalpy* el_enth() { return el_enth_; }
const ConfinementPotential *confpot(int ikp) const { return cfp[ikp]; }
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment