Commit 5117e207 by Martin Schlipf

add special G=0 treatement for HSE functional

Although G=0 is not divergent for HSE, the large changes near the
origin cause problems in the numerical calculation of the nonlocal
exchange energy. To correct this a "divergence" correction analogously
to PBE0 is introduced.

subtract exp(-rcut^2*G^2)*V(G) to make Fourier transform stable near G=0
set all energy contributions of G=0 in main loop to 0
add divergence correction at the end modifying the PBE0 one appropriately

new function in HSEFunctional is used to calculate the scaling of the
HSE "divergence" relative to the PBE0 one
InteractionPotential.h contains this additional function

the integrals ~k^2 are not evaluated because they are of higher order
(in PBE0 the potential is ~k^{-2} near G=0, so that the k^2 terms are
important, in HSE the potential is ~const. near G=0, so that the k^2
terms do not contribute a lot

git-svn-id: http://qboxcode.org/svn/qb/branches/hse-dev@1404 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 7876ddd8
......@@ -70,12 +70,14 @@ class HSEFunctional : public XCFunctional
static double interaction_potential(const double& g2);
// derivative of interaction potential
static double derivative_interaction_potential(const double& g2);
// scaling of the divergence correction relative to the Coulomb potential
static double divergence_scaling(const double& rcut);
// construct interaction potential class
static const InteractionPotential make_interaction_potential()
{
return InteractionPotential(&interaction_potential,
&derivative_interaction_potential);
&derivative_interaction_potential,&divergence_scaling);
}
};
......
......@@ -34,27 +34,29 @@ class InteractionPotential
public:
// default constructor = Coulomb potential
InteractionPotential():
InteractionPotential() :
coulomb_(true)
{
{
}
// constructor - define function and derivative
InteractionPotential(double(*V)(const double&), double(*dV)(const double&)) :
V_(V), dV_(dV), coulomb_(false)
InteractionPotential(double(*V)(const double&), double(*dV)(const double&),
double(*div_scale)(const double&)) :
V_(V), dV_(dV), div_scale_(div_scale), coulomb_(false)
{
}
// is the interaction potential a coulomb potential?
inline bool coulomb() const {
inline bool coulomb() const
{
return coulomb_;
}
// evaluate the interaction potential for given norm of G vector
inline double operator()(const double G2) const
{
// the current implementation expects that coulomb potential is treated externaly
assert( not coulomb_ );
assert(not coulomb_);
return V_(G2);
}
......@@ -62,15 +64,23 @@ class InteractionPotential
inline double derivative(const double G2) const
{
// the current implementation expects that coulomb potential is treated externaly
assert( not coulomb_ );
assert(not coulomb_);
return dV_(G2);
}
inline double divergence_scaling(const double rcut) const
{
// the current implementation expects that coulomb potential is treated externaly
assert(not coulomb_);
return div_scale_(rcut);
}
private:
const bool coulomb_;
double (*V_)(const double&);
double (*dV_)(const double&);
double (*div_scale_)(const double&);
};
......
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