Commit 4e6d0f58 by Francois Gygi

Fix non-linear core correction

parent c59340b5
......@@ -106,8 +106,6 @@ EnergyFunctional::EnergyFunctional(Sample& s, ChargeDensity& cd)
rhops[is].resize(ngloc);
}
xco = new XCOperator(s_, cd_);
nlp.resize(wf.nspin());
for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
{
......@@ -158,10 +156,8 @@ EnergyFunctional::EnergyFunctional(Sample& s, ChargeDensity& cd)
rhocore_sp_g.resize(nsp_);
for ( int is = 0; is < nsp_; is++ )
rhocore_sp_g[is].resize(ngloc);
rhocore_g.resize(ngloc);
rhocore_r.resize(vft->np012loc());
// set rhocore_r ptr in ChargeDensity
cd_.rhocore_r = &rhocore_r[0];
cd_.rhocore_g.resize(ngloc);
cd_.rhocore_r.resize(vft->np012loc());
}
// FT for interpolation of wavefunctions on the fine grid
......@@ -189,11 +185,9 @@ EnergyFunctional::EnergyFunctional(Sample& s, ChargeDensity& cd)
el_enth_ = new ElectricEnthalpy(s_);
sf.init(tau0,*vbasis_);
xco = new XCOperator(s_, cd_);
cell_moved();
atoms_moved();
}
////////////////////////////////////////////////////////////////////////////////
......@@ -342,7 +336,7 @@ void EnergyFunctional::update_vhxc(bool compute_stress, bool update_vh,
vlocal_g[ig] = vion_local_g[ig] + fpi * rhogt[ig] * g2i[ig];
}
// FT to tmpr_r
// FT to tmp_r
vft->backward(&vlocal_g[0],&tmp_r[0]);
// add external potential vext to tmp_r
......@@ -992,6 +986,7 @@ void EnergyFunctional::atoms_moved(void)
{
const AtomSet& atoms = s_.atoms;
const Wavefunction& wf = s_.wf;
const UnitCell& cell = s_.wf.cell();
int ngloc = vbasis_->localsize();
// fill tau0 with values in atom_list
......@@ -1007,23 +1002,22 @@ void EnergyFunctional::atoms_moved(void)
if ( core_charge_ )
{
// recalculate Fourier coefficients of the core charge
assert(rhocore_g.size()==ngloc);
memset( (void*)&rhocore_g[0], 0, 2*ngloc*sizeof(double) );
const double spin_fac = wf.cell().volume() / wf.nspin();
assert(cd_.rhocore_g.size()==ngloc);
memset( (void*)&cd_.rhocore_g[0], 0, 2*ngloc*sizeof(double) );
// divide core charge by two if two spins
const double spin_fac = cell.volume() / wf.nspin();
for ( int is = 0; is < atoms.nsp(); is++ )
{
complex<double> *s = &sf.sfac[is][0];
double *r = &rhocore_sp_g[is][0];
for ( int ig = 0; ig < ngloc; ig++ )
{
const complex<double> sg = s[ig];
// divide core charge by two if two spins
rhocore_g[ig] += spin_fac * sg * rhocore_sp_g[is][ig];
}
cd_.rhocore_g[ig] += spin_fac * s[ig] * r[ig];
}
// recalculate real-space core charge
vft->backward(&rhocore_g[0],&tmp_r[0]);
vft->backward(&cd_.rhocore_g[0],&tmp_r[0]);
const double omega_inv = 1.0 / cell.volume();
for ( int i = 0; i < tmp_r.size(); i++ )
rhocore_r[i] = tmp_r[i].real();
cd_.rhocore_r[i] = tmp_r[i].real() * omega_inv;
}
for ( int is = 0; is < atoms.nsp(); is++ )
......@@ -1039,7 +1033,6 @@ void EnergyFunctional::atoms_moved(void)
}
// compute esr: pseudocharge repulsion energy
const UnitCell& cell = s_.wf.cell();
const double omega_inv = 1.0 / cell.volume();
// distance between opposite planes of the cell
const double d0 = 2.0 * M_PI / length(cell.b(0));
......
......@@ -35,43 +35,46 @@ using namespace std;
XCPotential::XCPotential(const ChargeDensity& cd, const string functional_name,
const Control& ctrl): cd_(cd), vft_(*cd_.vft()), vbasis_(*cd_.vbasis())
{
// copy arrays to resize rhototal_r_ and rhototal_g_
rhototal_r_ = cd_.rhor;
rhototal_g_ = cd_.rhog;
if ( functional_name == "LDA" )
{
xcf_ = new LDAFunctional(cd_.rhor);
xcf_ = new LDAFunctional(rhototal_r_);
}
else if ( functional_name == "VWN" )
{
xcf_ = new VWNFunctional(cd_.rhor);
xcf_ = new VWNFunctional(rhototal_r_);
}
else if ( functional_name == "PBE" )
{
xcf_ = new PBEFunctional(cd_.rhor);
xcf_ = new PBEFunctional(rhototal_r_);
}
else if ( functional_name == "BLYP" )
{
xcf_ = new BLYPFunctional(cd_.rhor);
xcf_ = new BLYPFunctional(rhototal_r_);
}
else if ( functional_name == "PBE0" )
{
const double x_coeff = 1.0 - ctrl.alpha_PBE0;
const double c_coeff = 1.0;
xcf_ = new PBEFunctional(cd_.rhor,x_coeff,c_coeff);
xcf_ = new PBEFunctional(rhototal_r_,x_coeff,c_coeff);
}
else if ( functional_name == "HSE" )
{
xcf_ = new HSEFunctional(cd_.rhor);
xcf_ = new HSEFunctional(rhototal_r_);
}
else if ( functional_name == "RSH" )
{
xcf_ = new RSHFunctional(cd_.rhor,ctrl.alpha_RSH,ctrl.beta_RSH,ctrl.mu_RSH);
xcf_ = new RSHFunctional(rhototal_r_,ctrl.alpha_RSH,ctrl.beta_RSH,ctrl.mu_RSH);
}
else if ( functional_name == "B3LYP" )
{
xcf_ = new B3LYPFunctional(cd_.rhor);
xcf_ = new B3LYPFunctional(rhototal_r_);
}
else if ( functional_name == "BHandHLYP" )
{
xcf_ = new BHandHLYPFunctional(cd_.rhor);
xcf_ = new BHandHLYPFunctional(rhototal_r_);
}
else
{
......@@ -116,6 +119,10 @@ void XCPotential::update(vector<vector<double> >& vr)
// The array cd_.rhog is only used if xcf->isGGA() == true
// to compute the density gradients
// if a non-linear core correction is included,
// rhototal_r_ = cd_.rhor+cd_.rhocore_r. Otherwise
// rhototal_r_ = cd_.rhor
// Output: (through member function xcf())
//
// exc_, dxc_
......@@ -133,6 +140,24 @@ void XCPotential::update(vector<vector<double> >& vr)
// xcf()->vxc2_upup, xcf()->vxc2_dndn,
// xcf()->vxc2_updn, xcf()->vxc2_dnup
rhototal_r_ = cd_.rhor;
rhototal_g_ = cd_.rhog;
// test if a non-linear core correction is used
// if so, add core density to rhototal_r_ and rhototal_g_
if ( !cd_.rhocore_r.empty() )
{
// add core charge
// note: if nspin==2, the cd_.rhocore_{rg} vectors each
// contain half of the total core charge
for ( int ispin = 0; ispin < rhototal_r_.size(); ispin++ )
{
for ( int i = 0; i < rhototal_r_[ispin].size(); i++ )
rhototal_r_[ispin][i] += cd_.rhocore_r[i];
for ( int i = 0; i < rhototal_g_[ispin].size(); i++ )
rhototal_g_[ispin][i] += cd_.rhocore_g[i];
}
}
if ( !xcf_->isGGA() )
{
// LDA functional
......@@ -197,7 +222,7 @@ void XCPotential::update(vector<vector<double> >& vr)
for ( int ig = 0; ig < ngloc_; ig++ )
{
/* i*G_j*c(G) */
tmp1[ig] = complex<double>(0.0,omega_inv*gxj[ig]) * cd_.rhog[0][ig];
tmp1[ig] = complex<double>(0.0,omega_inv*gxj[ig])*rhototal_g_[0][ig];
}
vft_.backward(&tmp1[0],&tmpr[0]);
int inc2=2, inc1=1;
......@@ -210,8 +235,8 @@ void XCPotential::update(vector<vector<double> >& vr)
for ( int j = 0; j < 3; j++ )
{
const double *const gxj = vbasis_.gx_ptr(j);
const complex<double>* rhg0 = &cd_.rhog[0][0];
const complex<double>* rhg1 = &cd_.rhog[1][0];
const complex<double>* rhg0 = &rhototal_g_[0][0];
const complex<double>* rhg1 = &rhototal_g_[1][0];
for ( int ig = 0; ig < ngloc_; ig++ )
{
/* i*G_j*c(G) */
......
......@@ -37,6 +37,11 @@ class XCPotential
const ChargeDensity& cd_;
XCFunctional* xcf_;
// rhototal_r_: total density (valence+core)
std::vector<std::vector<double> > rhototal_r_;
// rhototal_g_: total density (valence+core)
std::vector<std::vector<std::complex<double> > > rhototal_g_;
std::vector<std::vector<double> > vxctmp; // vxctmp[ispin][ir]
std::vector<std::complex<double> > tmpr; // tmpr[ir]
std::vector<std::complex<double> > tmp1, tmp2; // tmp1[ig], tmp2[ig]
......
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