Commit 9c12d517 by Francois Gygi

Implemented Hartree-Fock stress


git-svn-id: http://qboxcode.org/svn/qb/trunk@1253 cba15fb0-1239-40c8-b417-11db7ca47a34
parent c7d77802
......@@ -366,7 +366,7 @@ void BOSampleStepper::step(int niter)
cd_.update_density();
tmap["charge"].stop();
ef_.update_vhxc();
ef_.update_vhxc(compute_stress);
const bool compute_forces = true;
double energy =
ef_.energy(false,dwf,compute_forces,fion,compute_stress,sigma_eks);
......@@ -838,7 +838,7 @@ void BOSampleStepper::step(int niter)
}
} // if nite_ > 1
ef_.update_vhxc();
ef_.update_vhxc(compute_stress);
// reset stepper only if multiple non-selfconsistent steps
if ( nite_ > 1 ) wf_stepper->preprocess();
......@@ -1048,7 +1048,7 @@ void BOSampleStepper::step(int niter)
cd_.update_density();
tmap["charge"].stop();
ef_.update_vhxc();
ef_.update_vhxc(compute_stress);
const bool compute_forces = true;
ef_.energy(false,dwf,compute_forces,fion,compute_stress,sigma_eks);
......@@ -1109,7 +1109,7 @@ void BOSampleStepper::step(int niter)
tmap["charge"].start();
cd_.update_density();
tmap["charge"].stop();
ef_.update_vhxc();
ef_.update_vhxc(compute_stress);
ef_.energy(true,dwf,false,fion,false,sigma_eks);
if ( onpe0 )
{
......@@ -1146,7 +1146,7 @@ void BOSampleStepper::step(int niter)
cd_.update_density();
tmap["charge"].stop();
ef_.update_vhxc();
ef_.update_vhxc(compute_stress);
const bool compute_forces = true;
double energy =
ef_.energy(false,dwf,compute_forces,fion,compute_stress,sigma_eks);
......@@ -1206,7 +1206,7 @@ void BOSampleStepper::step(int niter)
cd_.update_density();
tmap["charge"].stop();
ef_.update_vhxc();
ef_.update_vhxc(compute_stress);
const bool compute_forces = true;
double energy =
ef_.energy(false,dwf,compute_forces,fion,compute_stress,sigma_eks);
......
......@@ -132,7 +132,7 @@ void CPSampleStepper::step(int niter)
cd_.update_density();
tmap["charge"].stop();
ef_.update_vhxc();
ef_.update_vhxc(compute_stress);
double energy =
ef_.energy(compute_hpsi,dwf,compute_forces,fion,compute_stress,sigma_eks);
......@@ -273,7 +273,7 @@ void CPSampleStepper::step(int niter)
tmap["charge"].start();
cd_.update_density();
tmap["charge"].stop();
ef_.update_vhxc();
ef_.update_vhxc(compute_stress);
energy =
ef_.energy(compute_hpsi,dwf,compute_forces,fion,compute_stress,sigma_eks);
......
......@@ -196,7 +196,7 @@ EnergyFunctional::~EnergyFunctional(void)
}
////////////////////////////////////////////////////////////////////////////////
void EnergyFunctional::update_vhxc(void)
void EnergyFunctional::update_vhxc(bool compute_stress)
{
// called when the charge density has changed
// update Hartree and xc potentials using the charge density cd_
......@@ -227,20 +227,17 @@ void EnergyFunctional::update_vhxc(void)
}
}
// update XC energy and potential
// update XC operator
// compute xc energy, update self-energy operator and potential
tmap["exc"].start();
for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
fill(v_r[ispin].begin(),v_r[ispin].end(),0.0);
xco->update_v(v_r);
xco->update(v_r, compute_stress);
exc_ = xco->exc();
if ( compute_stress )
xco->compute_stress(sigma_exc);
tmap["exc"].stop();
// update self-energy operator
exhf_ = 0.0;
tmap["exhf"].start();
exhf_ = xco->update_sigma();
tmap["exhf"].stop();
// compute local potential energy:
// integral of el. charge times ionic local pot.
int len=2*ngloc,inc1=1;
......@@ -571,12 +568,6 @@ double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
sigma_ehart[5] = - 2.0 * fpi * tsum[5];
} // compute_stress
// Stress from exchange-correlation
if ( compute_stress )
{
xco->compute_stress(sigma_exc);
}
// Non local energy and forces
tmap["nonlocal"].start();
// modify next loop to span only local ikp
......@@ -612,7 +603,7 @@ 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_;
etotal_ = ekin_ + econf_ + eps_ + enl_ + ecoul_ + exc_ + ets_ + eexf_;
if ( compute_hpsi )
{
......@@ -660,7 +651,7 @@ double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
} //ispin
// apply self-energy operator
xco->apply_sigma(dwf);
xco->apply_self_energy(dwf);
tmap["hpsi"].stop();
} // if compute_hpsi
......@@ -1044,7 +1035,6 @@ void EnergyFunctional::print(ostream& os) const
<< " <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"
......
......@@ -66,7 +66,7 @@ class EnergyFunctional
int namax_;
int nsp_;
double ekin_, econf_, eps_, enl_, ehart_,
ecoul_, exc_, esr_, eself_, ets_, eexf_, exhf_, etotal_;
ecoul_, exc_, esr_, eself_, ets_, eexf_, etotal_;
std::valarray<double> sigma_ekin,sigma_econf,sigma_eps,sigma_ehart,sigma_exc,
sigma_enl, sigma_esr, sigma;
......@@ -87,7 +87,6 @@ class EnergyFunctional
double ehart(void) const { return ehart_; }
double ecoul(void) const { return ecoul_; }
double exc(void) const { return exc_; }
double exhf(void) const { return exhf_; }
double esr(void) const { return esr_; }
double eself(void) const { return eself_; }
double ets(void) const { return ets_; }
......@@ -95,7 +94,7 @@ class EnergyFunctional
const ConfinementPotential *confpot(int ikp) const { return cfp[ikp]; }
void update_vhxc(void);
void update_vhxc(bool compute_stress);
void atoms_moved(void);
void cell_moved(void);
......
......@@ -37,16 +37,21 @@ class ExchangeOperator
// mixing coefficient for exchange energy and dwf accumulation
double HFCoeff_;
// HF stress tensor
std::valarray<double> sigma_exhf_;
// reference wf and dwf for non scf iteration
Wavefunction wf0_;
Wavefunction dwf0_;
// copy of wf
Wavefunction wfc_;
double compute_exchange_for_general_case_(Sample* s, Wavefunction* dwf);
double compute_exchange_at_gamma_(const Wavefunction &wf, Wavefunction* dwf);
double compute_exchange_for_general_case_(Sample* s, Wavefunction* dwf,
bool compute_stress);
double compute_exchange_at_gamma_(const Wavefunction &wf, Wavefunction* dwf,
bool compute_stress);
void apply_VXC_(double mix, Wavefunction& wf_ref,
Wavefunction& dwf_ref, Wavefunction& dwf);
Wavefunction& dwf_ref, Wavefunction& dwf);
// Connectivity of the kpoint Grid
KPConnectivity KPGridPerm_;
......@@ -183,15 +188,14 @@ class ExchangeOperator
// parameters
void setmixCoeff(double value) { HFCoeff_ = value; };
void setIntegrationCoeff(double value) { rcut_ = value; };
double HFCoeff() { return HFCoeff_; };
double integrationCoeff() { return rcut_; };
// exchange energy and forces computation
double eex() { return eex_; };
double update_energy();
double update_sigma();
double apply_sigma(Wavefunction& dwf);
double update_energy(bool compute_stress);
double update_operator(bool compute_stress);
double apply_operator(Wavefunction& dwf);
void add_stress (std::valarray<double> & sigma_exc);
};
class ExchangeOperatorException
......
......@@ -27,7 +27,8 @@ XCOperator::XCOperator(Sample& s, const ChargeDensity& cd) :cd_(cd)
xcp_ = 0;
xop_ = 0;
exc_ = 0.0 ;
exhf_ = 0.0 ;
sigma_exc_.resize(6);
string functional_name = s.ctrl.xc;
......@@ -90,9 +91,10 @@ XCOperator::~XCOperator()
}
////////////////////////////////////////////////////////////////////////////////
void XCOperator::update_v(std::vector<std::vector<double> >& vr)
void XCOperator::update(std::vector<std::vector<double> >& vr, bool compute_stress)
{
// update: used whenever the charge density has changed
// update xc potential and self-energy
// used whenever the charge density and/or wave functions have changed
// compute vxc potential and energy
if ( hasPotential_ )
{
......@@ -101,45 +103,33 @@ void XCOperator::update_v(std::vector<std::vector<double> >& vr)
// LDA/GGA exchange energy
exc_ = xcp_->exc();
if ( compute_stress )
xcp_->compute_stress(sigma_exc_);
}
else
{
exc_ = 0.0;
sigma_exc_ = 0.0;
}
}
////////////////////////////////////////////////////////////////////////////////
double XCOperator::update_sigma(void)
{
if ( hasHF() )
{
return xop_->update_sigma();
exc_ += xop_->update_operator(compute_stress);
if ( compute_stress )
xop_->add_stress(sigma_exc_);
}
return 0.0;
}
////////////////////////////////////////////////////////////////////////////////
void XCOperator::apply_sigma(Wavefunction &dwf)
void XCOperator::apply_self_energy(Wavefunction &dwf)
{
if ( hasHF() )
{
xop_->apply_sigma(dwf);
}
xop_->apply_operator(dwf);
}
////////////////////////////////////////////////////////////////////////////////
void XCOperator::compute_stress(std::valarray<double>& sigma_exc)
void XCOperator::compute_stress(std::valarray<double>& sigma)
{
if ( hasPotential_ )
{
// LDA/GGA stress
xcp_->compute_stress( sigma_exc );
}
if ( hasHF_ )
{
throw XCOperatorException("stress not implemented with HF exchange");
}
sigma = sigma_exc_;
}
......@@ -33,8 +33,9 @@ class XCOperator
const ChargeDensity& cd_;
double HFmixCoeff_ ;
double exc_; // local XC energy
double exhf_; // Hartree-Fock exchange energy
double exc_; // XC energy: includes local and HF terms
std::valarray<double> sigma_exc_;
bool hasPotential_;
bool hasGGA_;
......@@ -57,12 +58,10 @@ class XCOperator
bool hasGGA(void) { return hasGGA_; };
bool hasHF(void) { return hasHF_; };
void update_v(std::vector<std::vector<double> >& vr);
double update_sigma(void);
void apply_sigma(Wavefunction &dwf);
void compute_stress(std::valarray<double>& sigma_exc);
double exc(void) { return exc_; };
double exhf(void) { return exhf_; };
void update(std::vector<std::vector<double> >& vr, bool compute_stress);
void apply_self_energy(Wavefunction &dwf);
void compute_stress(std::valarray<double>& sigma);
double exc(void) { return exc_ ; };
};
class XCOperatorException
......
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