Commit bea47293 by Francois Gygi

Harris-Foulkes estimates for etotal_int


git-svn-id: http://qboxcode.org/svn/qb/trunk@1780 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 063f16b8
......@@ -161,6 +161,7 @@ void BOSampleStepper::initialize_density(void)
////////////////////////////////////////////////////////////////////////////////
void BOSampleStepper::step(int niter)
{
const double boltz = 1.0 / ( 11605.0 * 2.0 * 13.6058 );
const Context& ctxt = s_.ctxt_;
const bool onpe0 = ctxt.onpe0();
......@@ -812,9 +813,8 @@ void BOSampleStepper::step(int niter)
if ( nite_ > 0 ) wf_stepper->preprocess();
// non-self-consistent loop
// repeat until the change in etotal_int or in the
// eigenvalue sum is smaller than a fraction of the change in
// Hartree energy in the last scf iteration
// repeat until the change in eigenvalue_sum is smaller than a
// fraction of the change in Hartree energy in the last scf iteration
bool nonscf_converged = false;
if ( itscf > 0 )
ehart_m = ehart;
......@@ -825,23 +825,16 @@ void BOSampleStepper::step(int niter)
// if ( onpe0 && nite_ > 0 )
// cout << " delta_ehart = " << delta_ehart << endl;
int ite = 0;
double etotal_int, etotal_int_m = 0.0;
double energy, etotal_int, etotal_int_m = 0.0;
double etotal, etotal_m = 0.0;
double eigenvalue_sum, eigenvalue_sum_m = 0.0;
// if nite == 0: do 1 iteration, no screening in charge mixing
// if nite > 0: do nite iterations, use screening in charge mixing
//
double energy;
while ( !nonscf_converged && ite < max(nite_,1) )
{
energy = ef_.energy(true,dwf,false,fion,false,sigma_eks);
double enthalpy = ef_.enthalpy();
if ( ite > 0 )
etotal_int_m = etotal_int;
etotal_int = energy;
// compute the sum of eigenvalues (with fixed weight)
// to measure convergence of the subspace update
......@@ -860,33 +853,9 @@ void BOSampleStepper::step(int niter)
wf_stepper->update(dwf);
if ( onpe0 )
{
cout.setf(ios::fixed,ios::floatfield);
cout.setf(ios::right,ios::adjustfield);
cout << " <etotal_int> " << setprecision(8) << setw(15)
<< energy << " </etotal_int>\n";
if ( compute_stress || ef_.el_enth() )
{
cout << " <enthalpy_int>" << setw(15)
<< enthalpy << " </enthalpy_int>\n" << flush;
}
}
// compare delta_etotal_int only after first iteration
// compare delta_eig_sum only after first iteration
if ( ite > 0 )
{
#if 0
double delta_etotal_int = fabs(etotal_int - etotal_int_m);
nonscf_converged |= (delta_etotal_int < 0.01 * delta_ehart);
if ( onpe0 )
{
cout << " BOSampleStepper::step: delta_e_int: "
<< delta_etotal_int << endl;
cout << " BOSampleStepper::step: delta_ehart: "
<< delta_ehart << endl;
}
#else
double delta_eig_sum = fabs(eigenvalue_sum - eigenvalue_sum_m);
nonscf_converged |= (delta_eig_sum < 0.01 * delta_ehart);
#ifdef DEBUG
......@@ -898,20 +867,10 @@ void BOSampleStepper::step(int niter)
<< delta_ehart << endl;
}
#endif
#endif
}
ite++;
}
if ( itscf > 0 )
etotal_m = etotal;
etotal = energy;
// if ( onpe0 && nite_ > 0 && ite >= nite_ )
// cout << " BOSampleStepper::step: nscf loop not converged after "
// << nite_ << " iterations" << endl;
// subspace diagonalization
if ( compute_eigvec || s_.ctrl.wf_diag == "EIGVAL" )
{
......@@ -950,18 +909,71 @@ void BOSampleStepper::step(int niter)
}
// update occupation numbers if fractionally occupied states
// compute weighted sum of eigenvalues and entropy term
// default value if no fractional occupation
double w_eigenvalue_sum = 2.0 * eigenvalue_sum;
const double wf_entropy = wf.entropy();
if ( fractional_occ )
{
wf.update_occ(s_.ctrl.fermi_temp);
const double wf_entropy = wf.entropy();
if ( onpe0 )
{
cout << " Wavefunction entropy: " << wf_entropy << endl;
const double boltz = 1.0 / ( 11605.0 * 2.0 * 13.6058 );
cout << " Entropy contribution to free energy: "
<< - wf_entropy * s_.ctrl.fermi_temp * boltz << endl;
}
w_eigenvalue_sum = 0.0;
for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
{
for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
{
const int nst = wf.sd(ispin,ikp)->nst();
const double wkp = wf.weight(ikp);
for ( int n = 0; n < nst; n++ )
{
const double occ = wf.sd(ispin,ikp)->occ(n);
w_eigenvalue_sum += wkp * occ * wf.sd(ispin,ikp)->eig(n);
}
}
}
}
// Harris-Foulkes estimate of the total energy
etotal_int = w_eigenvalue_sum - ef_.ehart_e() + ef_.ehart_p() +
ef_.esr() - ef_.eself() + ef_.dxc();
double enthalpy = ef_.enthalpy() -
wf_entropy * s_.ctrl.fermi_temp * boltz;
#ifdef DEBUG
if ( onpe0 )
{
cout << setprecision(8);
cout << "w_eigenvalue_sum = " << setw(15) << w_eigenvalue_sum << endl;
cout << "ef.dxc() = " << setw(15) << ef_.dxc() << endl;
cout << "ef.ehart() = " << setw(15) << ef_.ehart() << endl;
cout << "ef.ehart_e() = " << setw(15) << ef_.ehart_e() << endl;
cout << "ef.ehart_ep() = " << setw(15) << ef_.ehart_ep() << endl;
cout << "ef.esr() = " << setw(15) << ef_.esr() << endl;
}
#endif
if ( onpe0 )
{
cout.setf(ios::fixed,ios::floatfield);
cout.setf(ios::right,ios::adjustfield);
cout << " <etotal_int> " << setprecision(8) << setw(15)
<< etotal_int << " </etotal_int>\n";
if ( compute_stress || ef_.el_enth() )
{
cout << " <enthalpy_int>" << setw(15)
<< enthalpy << " </enthalpy_int>\n" << flush;
}
}
if ( itscf > 0 )
etotal_m = etotal;
etotal = etotal_int;
if ( nite_ > 0 && onpe0 )
cout << " BOSampleStepper: end scf iteration" << endl;
......
......@@ -232,7 +232,7 @@ void EnergyFunctional::update_vhxc(bool compute_stress)
const double *const g2i = vbasis_->g2i_ptr();
const double fpi = 4.0 * M_PI;
const int ngloc = vbasis_->localsize();
double sum[2], tsum[2];
double sum[5], tsum[5];
// compute total electronic density: rhoelg = rho_up + rho_dn
if ( wf.nspin() == 1 )
......@@ -258,6 +258,7 @@ void EnergyFunctional::update_vhxc(bool compute_stress)
xco->update(v_r, compute_stress);
exc_ = xco->exc();
dxc_ = xco->dxc();
if ( compute_stress )
xco->compute_stress(sigma_exc);
......@@ -294,24 +295,36 @@ void EnergyFunctional::update_vhxc(bool compute_stress)
}
sum[0] *= omega; // sum[0] contains eps
// Hartree energy
ehart_ = 0.0;
// Hartree energy and electron-electron energy (without pseudocharges)
double ehsum = 0.0;
double ehesum = 0.0;
double ehepsum = 0.0;
double ehpsum = 0.0;
for ( int ig = 0; ig < ngloc; ig++ )
{
const complex<double> tmp = rhoelg[ig] + rhopst[ig];
ehsum += norm(tmp) * g2i[ig];
rhogt[ig] = tmp;
const complex<double> r = rhoelg[ig];
const complex<double> rp = rhopst[ig];
ehsum += norm(r+rp) * g2i[ig];
ehesum += norm(r) * g2i[ig];
ehepsum += 2.0*real(conj(r)*rp * g2i[ig]);
ehpsum += norm(rp) * g2i[ig];
rhogt[ig] = r+rp;
}
// factor 1/2 from definition of Ehart cancels with half sum over G
// Note: rhogt[ig] includes a factor 1/Omega
// Factor omega in next line yields prefactor 4 pi / omega in
// Factor omega in next line yields prefactor 4 pi / omega in Ehart
sum[1] = omega * fpi * ehsum;
// tsum[1] contains ehart
sum[2] = omega * fpi * ehesum;
sum[3] = omega * fpi * ehepsum;
sum[4] = omega * fpi * ehpsum;
MPI_Allreduce(sum,tsum,2,MPI_DOUBLE,MPI_SUM,vbasis_->comm());
eps_ = tsum[0];
ehart_ = tsum[1];
MPI_Allreduce(sum,tsum,5,MPI_DOUBLE,MPI_SUM,vbasis_->comm());
eps_ = tsum[0];
ehart_ = tsum[1];
ehart_e_ = tsum[2];
ehart_ep_ = tsum[3];
ehart_p_ = tsum[4];
// compute vlocal_g = vion_local_g + vhart_g
// where vhart_g = 4 * pi * (rhoelg + rhopst) * g2i
......
......@@ -67,8 +67,9 @@ class EnergyFunctional
std::vector<double> zv_, rcps_;
std::vector<int> na_;
int nsp_;
double ekin_, econf_, eps_, enl_, ehart_,
double ekin_, econf_, eps_, enl_, ehart_, ehart_e_, ehart_ep_, ehart_p_,
ecoul_, exc_, esr_, eself_, ets_, eexf_, etotal_;
double dxc_;
double epv_, eefield_, enthalpy_;
std::valarray<double> sigma_ekin,sigma_econf,sigma_eps,sigma_ehart,sigma_exc,
sigma_enl, sigma_esr, sigma;
......@@ -90,8 +91,12 @@ class EnergyFunctional
double eps(void) const { return eps_; }
double enl(void) const { return enl_; }
double ehart(void) const { return ehart_; }
double ehart_e(void) const { return ehart_e_; }
double ehart_ep(void) const { return ehart_ep_; }
double ehart_p(void) const { return ehart_p_; }
double ecoul(void) const { return ecoul_; }
double exc(void) const { return exc_; }
double dxc(void) const { return dxc_; }
double esr(void) const { return esr_; }
double eself(void) const { return eself_; }
double ets(void) const { return ets_; }
......
......@@ -28,6 +28,7 @@ XCOperator::XCOperator(Sample& s, const ChargeDensity& cd) :cd_(cd)
xcp_ = 0;
xop_ = 0;
exc_ = 0.0 ;
dxc_ = 0.0 ;
sigma_exc_.resize(6);
......@@ -105,6 +106,7 @@ void XCOperator::update(std::vector<std::vector<double> >& vr, bool compute_stre
// LDA/GGA exchange energy
exc_ = xcp_->exc();
dxc_ = xcp_->dxc();
if ( compute_stress )
xcp_->compute_stress(sigma_exc_);
......@@ -112,12 +114,15 @@ void XCOperator::update(std::vector<std::vector<double> >& vr, bool compute_stre
else
{
exc_ = 0.0;
dxc_ = 0.0;
sigma_exc_ = 0.0;
}
if ( hasHF() )
{
exc_ += xop_->update_operator(compute_stress);
double ex_hf = xop_->update_operator(compute_stress);
exc_ += ex_hf;
dxc_ -= ex_hf;
if ( compute_stress )
xop_->add_stress(sigma_exc_);
}
......
......@@ -35,6 +35,7 @@ class XCOperator
const ChargeDensity& cd_;
double HFmixCoeff_ ;
double exc_; // XC energy: includes local and HF terms
double dxc_;
std::valarray<double> sigma_exc_;
......@@ -64,6 +65,7 @@ class XCOperator
void compute_stress(std::valarray<double>& sigma);
void cell_moved(void);
double exc(void) { return exc_ ; };
double dxc(void) { return dxc_ ; };
};
class XCOperatorException
......
......@@ -103,15 +103,15 @@ void XCPotential::update(vector<vector<double> >& vr)
// Output: (through member function xcf())
//
// exc_, dxc, dxc0_, dxc1_, dxc2_
// exc_, dxc_
//
// LDA Functional:
// exc_, dxc
// exc_, dxc_
// spin unpolarized: xcf()->exc, xcf()->vxc1
// spin polarized: xcf()->exc, xcf()->vxc1_up, xcf()->vxc1_dn
//
// GGA Functional: (through member function xcf())
// exc_, dxc, dxc0_, dxc1_, dxc2_
// exc_, dxc_
// spin unpolarized: xcf()->exc, xcf()->vxc1, xcf()->vxc2
// spin polarized: xcf()->exc_up, xcf()->exc_dn,
// xcf()->vxc1_up, xcf()->vxc1_dn
......@@ -125,6 +125,7 @@ void XCPotential::update(vector<vector<double> >& vr)
xcf_->setxc();
exc_ = 0.0;
dxc_ = 0.0;
const double *const e = xcf_->exc;
const int size = xcf_->np();
......@@ -135,8 +136,12 @@ void XCPotential::update(vector<vector<double> >& vr)
const double *const v = xcf_->vxc1;
for ( int i = 0; i < size; i++ )
{
exc_ += rh[i] * e[i];
vr[0][i] += v[i];
const double e_i = e[i];
const double v_i = v[i];
const double rh_i = rh[i];
exc_ += rh_i * e_i;
dxc_ += rh_i * ( e_i - v_i );
vr[0][i] += v_i;
}
}
else
......@@ -148,15 +153,19 @@ void XCPotential::update(vector<vector<double> >& vr)
const double *const v_dn = xcf_->vxc1_dn;
for ( int i = 0; i < size; i++ )
{
exc_ += (rh_up[i] + rh_dn[i]) * e[i];
const double r_i = rh_up[i] + rh_dn[i];
exc_ += r_i * e[i];
dxc_ += r_i * e[i] - rh_up[i] * v_up[i] - rh_dn[i] * v_dn[i];
vr[0][i] += v_up[i];
vr[1][i] += v_dn[i];
}
}
double sum = exc_ * vbasis_.cell().volume() / vft_.np012();
double tsum;
MPI_Allreduce(&sum,&tsum,1,MPI_DOUBLE,MPI_SUM,vbasis_.comm());
exc_ = tsum;
double sum[2],tsum[2];
sum[0] = exc_ * vbasis_.cell().volume() / vft_.np012();
sum[1] = dxc_ * vbasis_.cell().volume() / vft_.np012();
MPI_Allreduce(&sum,&tsum,2,MPI_DOUBLE,MPI_SUM,vbasis_.comm());
exc_ = tsum[0];
dxc_ = tsum[1];
}
else
{
......@@ -294,6 +303,7 @@ void XCPotential::update(vector<vector<double> >& vr)
// div(vxc2*grad_rho) is stored in vxctmp[ispin][ir]
double esum=0.0;
double dsum=0.0;
if ( nspin_ == 1 )
{
const double *const e = xcf_->exc;
......@@ -302,8 +312,12 @@ void XCPotential::update(vector<vector<double> >& vr)
{
for ( int ir = 0; ir < np012loc_; ir++ )
{
esum += rh[ir] * e[ir];
vr[0][ir] += v1[ir] + vxctmp[0][ir];
const double e_i = e[ir];
const double rh_i = rh[ir];
const double v_i = v1[ir] + vxctmp[0][ir];
esum += rh_i * e_i;
dsum += rh_i * ( e_i - v_i );
vr[0][ir] += v_i;
}
}
}
......@@ -317,17 +331,22 @@ void XCPotential::update(vector<vector<double> >& vr)
const double *const rh_dn = xcf_->rho_dn;
for ( int ir = 0; ir < np012loc_; ir++ )
{
const double r_up = rh_up[ir];
const double r_dn = rh_dn[ir];
esum += r_up * eup[ir] + r_dn * edn[ir];
vr[0][ir] += v1_up[ir] + vxctmp[0][ir];
vr[1][ir] += v1_dn[ir] + vxctmp[1][ir];
const double r_up_i = rh_up[ir];
const double r_dn_i = rh_dn[ir];
esum += r_up_i * eup[ir] + r_dn_i * edn[ir];
const double v_up = v1_up[ir] + vxctmp[0][ir];
const double v_dn = v1_dn[ir] + vxctmp[1][ir];
dsum += r_up_i * ( eup[ir] - v_up ) + r_dn_i * ( edn[ir] - v_dn );
vr[0][ir] += v_up;
vr[1][ir] += v_dn;
}
}
double sum = esum * vbasis_.cell().volume() / vft_.np012();
double tsum;
MPI_Allreduce(&sum,&tsum,1,MPI_DOUBLE,MPI_SUM,vbasis_.comm());
exc_ = tsum;
double sum[2], tsum[2];
sum[0] = esum * vbasis_.cell().volume() / vft_.np012();
sum[1] = dsum * vbasis_.cell().volume() / vft_.np012();
MPI_Allreduce(&sum,&tsum,2,MPI_DOUBLE,MPI_SUM,vbasis_.comm());
exc_ = tsum[0];
dxc_ = tsum[1];
}
}
////////////////////////////////////////////////////////////////////////////////
......@@ -339,7 +358,7 @@ void XCPotential::compute_stress(valarray<double>& sigma_exc)
{
// LDA functional
dxc_ = 0.0;
double dsum = 0.0;
const double *const e = xcf_->exc;
const int size = xcf_->np();
......@@ -350,7 +369,7 @@ void XCPotential::compute_stress(valarray<double>& sigma_exc)
const double *const v = xcf_->vxc1;
for ( int i = 0; i < size; i++ )
{
dxc_ += rh[i] * (e[i] - v[i]);
dsum += rh[i] * (e[i] - v[i]);
}
}
else
......@@ -363,14 +382,14 @@ void XCPotential::compute_stress(valarray<double>& sigma_exc)
for ( int i = 0; i < size; i++ )
{
const double rh = rh_up[i] + rh_dn[i];
dxc_ += rh * e[i] - rh_up[i] * v_up[i] - rh_dn[i] * v_dn[i];
dsum += rh * e[i] - rh_up[i] * v_up[i] - rh_dn[i] * v_dn[i];
}
}
const double fac = 1.0 / vft_.np012();
double sum, tsum;
// Next line: factor omega in volume element cancels 1/omega in
// definition of sigma_exc
sum = - fac * dxc_;
sum = - fac * dsum;
MPI_Allreduce(&sum,&tsum,1,MPI_DOUBLE,MPI_SUM,vbasis_.comm());
// Note: contribution to sigma_exc is a multiple of the identity
......
......@@ -40,7 +40,7 @@ class XCPotential
std::vector<std::complex<double> > tmpr; // tmpr[ir]
std::vector<std::complex<double> > tmp1, tmp2; // tmp1[ig], tmp2[ig]
double exc_, dxc_, dxc0_, dxc1_, dxc2_;
double exc_, dxc_;
int nspin_;
int ngloc_;
int np012loc_;
......@@ -58,6 +58,7 @@ class XCPotential
void update(std::vector<std::vector<double> >& vr);
void compute_stress(std::valarray<double>& sigma_exc);
double exc(void) { return exc_; }
double dxc(void) { return dxc_; }
};
class XCPotentialException
......
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