Simplified interface to XC functionals. More...

#include <xcfunctional.h>

Public Types
enum	xc_arg { enum_rhoa =0 , enum_rhob =1 , enum_rho_pt =2 , enum_taua =3 , enum_taub =4 , enum_saa =10 , enum_sab =11 , enum_sbb =12 , enum_sigtot =13 , enum_sigma_pta_div_rho =14 , enum_sigma_ptb_div_rho =15 , enum_zetaa_x =16 , enum_zetaa_y =17 , enum_zetaa_z =18 , enum_zetab_x =19 , enum_zetab_y =20 , enum_zetab_z =21 , enum_ddens_ptx =25 , enum_ddens_pty =26 , enum_ddens_ptz =27 , enum_nemo_R2 =28 , enum_gradfa =29 , enum_gradfb =30 , enum_na =31 , enum_nb =32 , enum_Ga_x =33 , enum_Ga_y =34 , enum_Ga_z =35 , enum_Gb_x =36 , enum_Gb_y =37 , enum_Gb_z =38 , enum_u1_x =39 , enum_u1_y =40 , enum_u1_z =41 , enum_u1sq =42 }

Public Member Functions
	XCfunctional ()
	Default constructor is required.

	~XCfunctional ()
	Destructor.

madness::Tensor< double >	exc (const std::vector< madness::Tensor< double > > &t) const
	Computes the energy functional at given points.

std::vector< madness::Tensor< double > >	fxc_apply (const std::vector< madness::Tensor< double > > &t, const int ispin) const
	compute the second derivative of the XC energy wrt the density and apply

double	get_rhotol () const
	return the munging threshold for the density

double	get_tautol () const
	return the floor for the kinetic energy density

double	get_tauwmargin () const
	return the margin by which the von Weizsaecker clamp overshoots

bool	has_fxc () const
	Returns true if the second derivative of the functional is available (not yet supported)

bool	has_kxc () const
	Returns true if the third derivative of the functional is available (not yet supported)

double	hf_exchange_coefficient () const
	Returns the value of the hf exact exchange coefficient.

void	initialize (const std::string &input_line, bool polarized, World &world, const bool verbose=false)
	Initialize the object from the user input data.

bool	is_dft () const
	Returns true if there is a DFT functional (false probably means Hatree-Fock exchange only)

bool	is_gga () const
	Returns true if the potential is gga (needs first derivatives)

bool	is_lda () const
	Returns true if the potential is lda.

bool	is_meta () const
	Returns true if the potential is meta gga (needs the kinetic energy density)

bool	is_spin_polarized () const
	Returns true if the functional is spin_polarized.

bool	needs_sigma () const
	Returns true if the functional needs the reduced density gradients sigma.

bool	needs_tau () const
	Returns true if the functional needs the kinetic energy density tau.

void	plot () const
	Crude function to plot the energy and potential functionals.

double	tau_w_bound (const double sigma, const double rho) const
	the von Weizsaecker lower bound on tau, as the clamp actually applies it

bool	uses_libxc_backend () const
	Returns true when libxc is the active XC backend for this functional.

std::vector< madness::Tensor< double > >	vxc (const std::vector< madness::Tensor< double > > &t, const int ispin) const
	Computes components of the potential (derivative of the energy functional) at np points.

Static Public Member Functions
static double	munge_old (double rho)

Static Public Attributes
static const int	number_xc_args =43
	max number of intermediates

Protected Member Functions
void	make_libxc_args (const std::vector< madness::Tensor< double > > &t, madness::Tensor< double > &rho, madness::Tensor< double > &sigma, madness::Tensor< double > &tau, madness::Tensor< double > &rho_pt, madness::Tensor< double > &sigma_pt, std::vector< madness::Tensor< double > > &drho, std::vector< madness::Tensor< double > > &drho_pt, const bool need_response) const
	convert the raw density (gradient) data to be used by the xc operators

void	reset_screening_defaults ()
	put the configurable screening thresholds back to their defaults

Static Protected Member Functions
static void	polyn (const double x, double &p, double &dpdx)
	Smoothly switches between constant (x<xmin) and linear function (x>xmax)

Protected Attributes
double	hf_coeff =0.0
	Factor multiplying HF exchange (+1.0 gives HF)

int	nderiv =0
	Jacob's ladder rung; 0: lda, 1: gga, 2: mgga.

double	rhomin =default_rhomin
	what munge() puts in place of a density

double	rhotol =default_rhotol
	See initialize and munge*.

bool	spin_polarized =false
	True if the functional is spin polarized.

double	tautol =default_tautol
	floor for the kinetic energy density, see initialize

double	tauwmargin =default_tauwmargin
	von Weizsaecker clamp overshoot, see tau_w_bound

Static Protected Attributes
static constexpr double	default_rhomin =1.e-12
	our lda will divide by rho

static constexpr double	default_rhotol =1.e-7

static constexpr double	default_tautol =1.e-12

static constexpr double	default_tauwmargin =1.e-6

Private Member Functions
double	binary_munge (double rho, double refrho, const double thresh) const
	zero a quantity where the reference density is small

double	munge (double rho) const
	simple munging for the density only (LDA)

Detailed Description

Simplified interface to XC functionals.

Member Enumeration Documentation

◆ xc_arg

enum madness::XCfunctional::xc_arg

The ordering of the intermediates is fixed, but the code can handle non-initialized functions, so if e.g. no GGA is requested, all the corresponding vector components may be left empty.

Note the additional quantities $\zeta$ and $\chi$ , which are defined as

$\rho = \exp(\zeta)$

and thus the derivative of rho is given by

$\nabla_x\rho = \exp(\zeta)\nabla_x\zeta = \rho \nabla_x\zeta$

The reduced gradients \sigma may then be expressed as

$\sigma = |\nabla\rho|^2 = |\rho|^2 |\nabla\zeta|^2 = |\rho|^2 \chi$

Enumerator
enum_rhoa	alpha density $\rho_\alpha$
enum_rhob	beta density $\rho_\beta$
enum_rho_pt	perturbed density (CPHF, TDKS) $\rho_{pt}$
enum_taua	alpha kinetic energy density $\tau_\alpha = \frac{1}{2}\sum_i\|\nabla\psi_{i\alpha}\|^2$
enum_taub	beta kinetic energy density $\tau_\beta$
enum_saa	$\sigma_{aa} = \nabla \rho_{\alpha}.\nabla \rho_{\alpha}$
enum_sab	$\sigma_{ab} = \nabla \rho_{\alpha}.\nabla \rho_{\beta}$
enum_sbb	$\sigma_{bb} = \nabla \rho_{\beta}.\nabla \rho_{\beta}$
enum_sigtot	$\sigma = \nabla \rho.\nabla \rho$
enum_sigma_pta_div_rho	$\zeta_{\alpha}.\nabla\rho_{pt}$
enum_sigma_ptb_div_rho	$\zeta_{\beta}.\nabla\rho_{pt}$
enum_zetaa_x	$\zeta_{a,x}=\partial/{\partial x} \ln(\rho_a)$
enum_zetaa_y	$\zeta_{a,y}=\partial/{\partial y} \ln(\rho_a)$
enum_zetaa_z	$\zeta_{a,z}=\partial/{\partial z} \ln(\rho_a)$
enum_zetab_x	$\zeta_{b,x} = \partial/{\partial x} \ln(\rho_b)$
enum_zetab_y	$\zeta_{b,y} = \partial/{\partial y} \ln(\rho_b)$
enum_zetab_z	$\zeta_{b,z} = \partial/{\partial z} \ln(\rho_b)$
enum_ddens_ptx	$\nabla\rho_{pt}$
enum_ddens_pty	$\nabla\rho_{pt}$
enum_ddens_ptz	$\nabla\rho_{pt}$
enum_nemo_R2	, the ncf squared [MRA, smooth]
enum_gradfa	$\sum_i w_i\|\nabla F_{i\alpha}\|^2$ [MRA, smooth]
enum_gradfb	beta counterpart [MRA, smooth]
enum_na	$n_\alpha=\sum_i w_iF_{i\alpha}^2$ [MRA, smooth]
enum_nb	beta counterpart [MRA, smooth]
enum_Ga_x	$G_{\alpha,x}=\sum_i w_iF_i\partial_xF_i$ [MRA, smooth]
enum_Ga_y	$G_{\alpha,y}$ [MRA, smooth]
enum_Ga_z	$G_{\alpha,z}$ [MRA, smooth]
enum_Gb_x	beta counterpart [MRA, smooth]
enum_Gb_y	beta counterpart [MRA, smooth]
enum_Gb_z	beta counterpart [MRA, smooth]
enum_u1_x	$U_{1,x}$ [functor, never MRA]
enum_u1_y	$U_{1,y}$ [functor, never MRA]
enum_u1_z	$U_{1,z}$ [functor, never MRA]
enum_u1sq	$\|\mathbf U_1\|^2$ [functor, never MRA]

Constructor & Destructor Documentation

◆ XCfunctional()

madness::XCfunctional::XCfunctional ( )

inline

Default constructor is required.

◆ ~XCfunctional()

madness::XCfunctional::~XCfunctional ( )

Destructor.

Member Function Documentation

◆ binary_munge()

double madness::XCfunctional::binary_munge	(	double	rho,
		double	refrho,
		const double	thresh
	)		const

inlineprivate

zero a quantity where the reference density is small

Used for perturbed densities, which may be negative and much more diffuse than the ground state, and for screening outputs. Only where the reference density is large enough is DFT numerically well defined.

Substitutes zero, not rhomin. The argument is not necessarily a density – de/dtau and the semilocal response terms go through here too – so a density floor is the wrong thing to leave behind, and "screened" means "contributes nothing". rhomin stays what munge() puts in place of a density.

Parameters

[in]	rho	number to be munged
[in]	refrho	reference value for munging
[in]	thresh	threshold for munging

References madness::nonlinear_vector_solver(), and thresh.

◆ exc()

madness::Tensor< double > madness::XCfunctional::exc ( const std::vector< madness::Tensor< double > > & t ) const

Computes the energy functional at given points.

This uses the convention that the total energy is $E[\rho] = \int \epsilon[\rho(x)] dx$ Any HF exchange contribution must be separately computed. Items in the vector argument t are interpreted similarly to the xc_arg enum.

Parameters

[in] t The input densities and derivatives as required by the functional

Returns: The exchange-correlation energy functional

References madness::c_rks_vwn5__(), madness::c_uks_vwn5__(), madness::f, std::isnan(), MADNESS_PRAGMA_CLANG, munge(), madness::nonlinear_vector_solver(), madness::pop(), madness::print(), madness::Tensor< T >::ptr(), madness::BaseTensor::size(), spin_polarized, v, madness::x_rks_s__(), madness::x_uks_s__(), and madness::xf().

Referenced by madness::xc_functional::operator()(), and plot().

◆ fxc_apply()

std::vector< madness::Tensor< double > > madness::XCfunctional::fxc_apply	(	const std::vector< madness::Tensor< double > > &	t,
		const int	ispin
	)		const

compute the second derivative of the XC energy wrt the density and apply

Return the following quantities (RHF only) (see Yanai2005, Eq. (13))

$\begin{eqnarray*} \mbox{result[0]} &:& \qquad \frac{\partial^2 \epsilon}{\partial \rho^2} \rho_\mathrm{pt} + 2.0 * \frac{\partial^2 \epsilon}{\partial \rho\partial\sigma}\sigma_\mathrm{pt}\\ \mbox{result[1-3]} &:& \qquad 2.0 * \frac{\partial\epsilon}{\partial\sigma}\nabla\rho_\mathrm{pt} + 2.0 * \frac{\partial^2\epsilon}{\partial\rho\partial\sigma} \rho_\mathrm{pt}\nabla\rho + 4.0 * \frac{\partial^2\epsilon}{\partial^2\sigma} \sigma_\mathrm{pt}\nabla\rho \end{eqnarray*}$

Parameters

[in]	t	The input densities and derivatives as required by the functional, as in the xc_arg enum
[in]	ispin	not referenced since only RHF is implemented, always 0

Returns: a vector of Functions containing the contributions to the kernel apply

References MADNESS_EXCEPTION.

Referenced by madness::xc_kernel_apply::operator()().

◆ get_rhotol()

double madness::XCfunctional::get_rhotol ( ) const

inline

return the munging threshold for the density

References rhotol.

◆ get_tautol()

double madness::XCfunctional::get_tautol ( ) const

inline

return the floor for the kinetic energy density

meta-gga functionals build the iso-orbital indicators alpha and z with tau in the denominator, so tau needs a floor well above libxc's own 1e-20 default for a real-space code, where tau is a numerical derivative

References tautol.

◆ get_tauwmargin()

double madness::XCfunctional::get_tauwmargin ( ) const

inline

return the margin by which the von Weizsaecker clamp overshoots

References tauwmargin.

◆ has_fxc()

bool madness::XCfunctional::has_fxc ( ) const

Returns true if the second derivative of the functional is available (not yet supported)

◆ has_kxc()

bool madness::XCfunctional::has_kxc ( ) const

Returns true if the third derivative of the functional is available (not yet supported)

◆ hf_exchange_coefficient()

double madness::XCfunctional::hf_exchange_coefficient ( ) const

inline

Returns the value of the hf exact exchange coefficient.

References hf_coeff.

Referenced by madness::SCF::apply_potential().

◆ initialize()

void madness::XCfunctional::initialize	(	const std::string &	input_line,
		bool	polarized,
		World &	world,
		const bool	verbose = `false`
	)

Initialize the object from the user input data.

Parameters

[in]	input_line	User input line (without beginning XC keyword)
[in]	polarized	Boolean flag indicating if the calculation is spin-polarized

References hf_coeff, madness::nonlinear_vector_solver(), reset_screening_defaults(), rhomin, rhotol, and spin_polarized.

Referenced by MiniDFT::MiniDFT(), madness::SCF::SCF(), main(), test_hybrid_coefficients(), test_lda(), test_meta_gga_dedtau_polarized(), test_meta_gga_one_orbital_energy(), and test_xcfunctional().

◆ is_dft()

bool madness::XCfunctional::is_dft ( ) const

Returns true if there is a DFT functional (false probably means Hatree-Fock exchange only)

References is_gga(), is_lda(), and is_meta().

Referenced by madness::SCF::apply_potential().

◆ is_gga()

bool madness::XCfunctional::is_gga ( ) const

Returns true if the potential is gga (needs first derivatives)

Referenced by madness::XCOperator< T, NDIM >::apply_xc_kernel(), madness::xc_kernel_apply::get_result_size(), is_dft(), and madness::XCOperator< T, NDIM >::prep_xc_args_response().

◆ is_lda()

bool madness::XCfunctional::is_lda ( ) const

Returns true if the potential is lda.

References hf_coeff.

Referenced by madness::xc_potential::get_result_size(), and is_dft().

◆ is_meta()

bool madness::XCfunctional::is_meta ( ) const

Returns true if the potential is meta gga (needs the kinetic energy density)

Referenced by is_dft().

◆ is_spin_polarized()

bool madness::XCfunctional::is_spin_polarized ( ) const

inline

Returns true if the functional is spin_polarized.

References spin_polarized.

Referenced by madness::SCF::apply_potential(), madness::XCOperator< T, NDIM >::apply_xc_kernel(), madness::xc_potential::get_result_size(), madness::XCOperator< T, NDIM >::make_xc_potential_impl(), plot(), madness::XCOperator< T, NDIM >::prep_xc_args(), and madness::XCOperator< T, NDIM >::prep_xc_args_response().

◆ make_libxc_args()

void madness::XCfunctional::make_libxc_args	(	const std::vector< madness::Tensor< double > > &	xc_args,
		madness::Tensor< double > &	rho,
		madness::Tensor< double > &	sigma,
		madness::Tensor< double > &	tau,
		madness::Tensor< double > &	rho_pt,
		madness::Tensor< double > &	sigma_pt,
		std::vector< madness::Tensor< double > > &	drho,
		std::vector< madness::Tensor< double > > &	drho_pt,
		const bool	need_response
	)		const

protected

convert the raw density (gradient) data to be used by the xc operators

assemble the arguments libxc is evaluated on, munged and floored

Involves mainly munging of the densities and multiplying with 2 if the calculation is spin-restricted. Response densities and density gradients are munged based on the value of the ground state density, since they may become negative and may also be much more diffuse. dimensions of the output tensors are for spin-restricted and unrestricted (with np the number of grid points in the box): rho(np) or rho(2*np) sigma(np) sigma(3*np) rho_pt(np) sigma_pt(2*np)

Parameters

[in]	t	input density (gradients)
[out]	rho	ground state (spin) density, properly munged
[out]	sigma	ground state (spin) density gradients, properly munged
[out]	tau	ground state (spin) kinetic energy density, properly munged
[out]	rho_pt	response density, properly munged (no spin)
[out]	sigma_pt	response (spin) density gradients, properly munged
[out]	drho	density derivative, constructed from rho and zeta
[out]	drho_pt	response density derivative directly from xc_args
[in]	need_response	flag if rho_pt and sigma_pt need to be calculated

MUNGING, IN ONE PLACE

Every point handed to libxc must describe a state some wavefunction could be in. libxc relies on that: sigma's diagonals are non-negative, |sigma_ab| obeys Cauchy-Schwarz, the total sigma_aa + 2 sigma_ab + sigma_bb = |grad rho|^2 is non-negative (the correlation kernels return NaN for vsigma as soon as it is not), and tau >= tau_W = sigma/(8 rho) keeps the Fermi hole curvature positive and the iso-orbital indicators inside their domain (the clamp overshoots that bound by tauwmargin – see tau_w_bound). Screening one argument without its partners manufactures points no density can produce, and the functional's derivatives there answer no question.

One threshold decides all of it: rhotol (1e-7). It says whether there is any density at a point, and everything in that channel goes away together – and the outputs that need screening are screened on the same floor, so that a quantity is cut exactly where the density it came from was cut.

One mask per spin channel, decided once from the raw density and then applied to every argument of that channel:

m_s = (rho_s^raw > rhotol) rho_s = m_s ? rho_s^raw : rhomin grad = m_s * grad(rho_s)^raw sig_ss= m_s * max(1e-14, |grad(rho_s)|^2) sig_ab= clamped to +-sqrt(sig_aa sig_bb) – zero if either mask is zero tau_s = m_s ? max(tautol, tau_s^raw) : 0 tau_s = max(tau_s, tau_w_bound(sig_ss, rho_s)) – inert when m_s is zero

The floors (1e-14 on sigma, tautol on tau) are conditioning devices for where density exists, so the mask goes outside them: applied unconditionally they would resurrect exactly what the mask removed. sigma is contracted from the very gradients written into drho, by either the zeta or the split route, so it is the Gram matrix of those gradients and the three sigma bounds above hold for free.

Response densities are munged on the ground-state density instead, via binary_munge: they may be negative and much more diffuse, and DFT is only well defined where the reference density is there.

In the spin-restricted branch there is one channel and it is the total density, twice the alpha one; the rules above read the same with s dropped.

Outputs are screened elsewhere, and only where no factor already does the job. vxc() cuts de/dtau on its own spin's density, because it multiplies grad(psi) and so does not vanish with the density; fxc_apply() cuts its semilocal response terms the same way. Both use rhotol. Everything else needs no cut: the semilocal terms carry a factor of grad(rho_s), masked above, and de/drho is bare but stays finite once the arguments assembled here are consistent – which is what the tau rule is for. A channel left on the tautol floor at rho = 0 sends de/drho_s off like rho_s^(-8/3) through tau_unif ~ rho_s^(5/3).

exc() multiplies by the density, but by the TOTAL density in the spin-polarized branch rather than the channel's own, so it is damped only where every channel is empty. Where one spin is munged away and the other is not, it is again the consistency of the arguments, not a factor, that keeps it finite.

References MADNESS_EXCEPTION.

◆ munge()

double madness::XCfunctional::munge ( double rho ) const

inlineprivate

simple munging for the density only (LDA)

References rhomin, and rhotol.

Referenced by exc(), and vxc().

◆ munge_old()

static double madness::XCfunctional::munge_old ( double rho )

inlinestatic

References madness::nonlinear_vector_solver(), p(), and polyn().

◆ needs_sigma()

bool madness::XCfunctional::needs_sigma ( ) const

Returns true if the functional needs the reduced density gradients sigma.

True for gga AND meta-gga – a meta-gga needs the density gradients as well as the kinetic energy density. Use this, not is_gga(), to decide whether the gradient intermediates have to be computed.

Referenced by madness::XCOperator< T, NDIM >::is_weak_form(), madness::XCOperator< T, NDIM >::make_xc_potential_impl(), plot(), and madness::XCOperator< T, NDIM >::prep_xc_args().

◆ needs_tau()

bool madness::XCfunctional::needs_tau ( ) const

Returns true if the functional needs the kinetic energy density tau.

Referenced by madness::xc_potential::get_result_size(), madness::XCOperator< T, NDIM >::has_tau_term(), and plot().

◆ plot()

void madness::XCfunctional::plot ( ) const

inline

Crude function to plot the energy and potential functionals.

References e(), enum_rhoa, enum_rhob, enum_saa, enum_taua, exc(), madness::f, is_spin_polarized(), lo, needs_sigma(), needs_tau(), madness::nonlinear_vector_solver(), and vxc().

◆ polyn()

static void madness::XCfunctional::polyn	(	const double	x,
		double &	p,
		double &	dpdx
	)

inlinestaticprotected

Smoothly switches between constant (x<xmin) and linear function (x>xmax)

$f(x,x_{\mathrm{min}},x_{\mathrm{max}}) = \left\{ \begin{array}{ll} x_{\mathrm{min}} & x < x_{\mathrm{min}} \\ p(x,x_{\mathrm{min}},x_{\mathrm{max}}) & x_{\mathrm{min}} \leq x_{\mathrm{max}} \\ x & x_{\mathrm{max}} < x \end{array} \right.$

where $p(x)$ is the unique quintic polynomial that satisfies $p(x_{min})=x_{min}$ , $p(x_{max})=x_{max}$ , $dp(x_{max})/dx=1$ , and $dp(x_{min})/dx=d^2p(x_{min})/dx^2=d^2p(x_{max})/dx^2=0$ .

References a, b, c, d, e(), madness::nonlinear_vector_solver(), and p().

Referenced by munge_old().

◆ reset_screening_defaults()

void madness::XCfunctional::reset_screening_defaults ( )

inlineprotected

put the configurable screening thresholds back to their defaults

initialize() may be called more than once on the same object, and the xc input line can override any of these (RHOMIN/RHOTOL/TAUTOL), so they have to be restored before the next line is parsed – otherwise one functional's thresholds leak into the next one.

References default_rhomin, default_rhotol, default_tautol, default_tauwmargin, rhomin, rhotol, tautol, and tauwmargin.

Referenced by initialize().

◆ tau_w_bound()

double madness::XCfunctional::tau_w_bound	(	const double	sigma,
		const double	rho
	)		const

inline

the von Weizsaecker lower bound on tau, as the clamp actually applies it

tau >= tau_W = |grad rho_s|^2/(8 rho_s) is exact for any wavefunction, and meta-ggas need it: they are built on z = tau_W/tau, whose domain is [0,1], and outside it the Fermi hole curvature turns negative and libxc's correlation kernels return NaN.

Build it from SIGMA, not from chi. libxc forms z from the sigma it is handed, so the bound has to be built from that same sigma. The form rho*chi/8 equals sigma/(8 rho) only when sigma is literally rho^2 chi, which is false: sigma carries a positivity floor that chi does not.

And overshoot it. Clamping to exactly tau_W puts z at 1 to within one ulp – the endpoint of the domain, and the worst-conditioned point in it, where a one-ulp change in tau moves de/dtau by a factor of 30. Landing strictly inside costs nothing: the clamp only ever fires where tau was below a bound it should have satisfied anyway. Same device as r2SCAN's eta*tau_W.

Returns 0 where the density has been munged away, so the clamp is inert there: sigma's floor divided by a vanishing rho would grow like 1/rho.

References sigma, and tauwmargin.

◆ uses_libxc_backend()

bool madness::XCfunctional::uses_libxc_backend ( ) const

inline

Returns true when libxc is the active XC backend for this functional.

References madness::nonlinear_vector_solver().

◆ vxc()

std::vector< madness::Tensor< double > > madness::XCfunctional::vxc	(	const std::vector< madness::Tensor< double > > &	t,
		const int	ispin
	)		const

Computes components of the potential (derivative of the energy functional) at np points.

Any HF exchange contribution must be separately computed. Items in the vector argument t are interpreted similarly to the xc_arg enum.

We define $\sigma_{\mu \nu} = \nabla \rho_{\mu} . \nabla \rho_{\nu}$ with $\mu, \nu = \alpha$ or $\beta$ .

For unpolarized GGA, matrix elements of the potential are

$< \phi | \hat V | \psi > = \int \left( \frac{\partial \epsilon}{\partial \rho} \phi \psi + \left( 2 \frac{\partial \epsilon}{\partial \sigma} \right) \nabla \rho \cdot \nabla \left( \phi \psi \right) \right) dx$

For polarized GGA, matrix elements of the potential are

$< \phi_{\alpha} | \hat V | \psi_{\alpha} > = \int \left( \frac{\partial \epsilon}{\partial \rho_{\alpha}} \phi \psi + \left( 2 \frac{\partial \epsilon}{\partial \sigma_{\alpha \alpha}} \nabla \rho_{\alpha} + \frac{\partial \epsilon}{\partial \sigma_{\alpha \beta}} \nabla \rho_{\beta} \right) . \nabla \left( \phi \psi \right) \right) dx$

Integrating the above by parts and assuming free-space or periodic boundary conditions we obtain that the local multiplicative form of the GGA potential is

$V_{\alpha} = \frac{\partial \epsilon}{\partial \rho_{\alpha}} - \left(\nabla . \left(2 \frac{\partial \epsilon}{\partial \sigma_{\alpha \alpha}} \nabla \rho_{\alpha} + \frac{\partial \epsilon}{\partial \sigma_{\alpha \beta}} \nabla \rho_{\beta} \right) \right)$

Return the following quantities for RHF: (see Yanai2005, Eq. (12))

$\begin{eqnarray*} \mbox{result[0]} &:& \qquad \frac{\partial \epsilon}{\partial \rho} \\ \mbox{result[1-3]} &:& \qquad 2 \rho \frac{\partial \epsilon}{\partial \sigma} \nabla\rho \end{eqnarray*}$

and for UHF same-spin and other-spin quantities

$\begin{eqnarray*} \mbox{result[0]} &:& \qquad \frac{\partial \epsilon}{\partial \rho_{\alpha}} \\ \mbox{result[1-3]} &:& \qquad \rho_\alpha \frac{\partial \epsilon}{\partial \sigma_{\alpha \alpha}} \nabla\rho_\alpha\\ \mbox{result[4-6]} &:& \qquad \rho_\alpha \frac{\partial \epsilon}{\partial \sigma_{\alpha \beta}} \nabla\rho_\beta \end{eqnarray*}$

Parameters

[in]	t	The input densities and derivatives as required by the functional
[in]	ispin	Specifies which component of the potential is to be computed as described above

Returns: the requested quantity, based on ispin (0: same spin, 1: other spin)

References madness::c_rks_vwn5__(), madness::c_uks_vwn5__(), madness::f, std::isnan(), MADNESS_PRAGMA_CLANG, munge(), madness::nonlinear_vector_solver(), madness::pop(), madness::print(), madness::Tensor< T >::ptr(), q(), spin_polarized, v, madness::x_rks_s__(), madness::x_uks_s__(), and madness::xf().

Referenced by madness::xc_potential::operator()(), and plot().

Member Data Documentation

◆ default_rhomin

constexpr double madness::XCfunctional::default_rhomin =1.e-12

staticconstexprprotected

our lda will divide by rho

Referenced by reset_screening_defaults().

◆ default_rhotol

constexpr double madness::XCfunctional::default_rhotol =1.e-7

staticconstexprprotected

Referenced by reset_screening_defaults().

◆ default_tautol

constexpr double madness::XCfunctional::default_tautol =1.e-12

staticconstexprprotected

Referenced by reset_screening_defaults().

◆ default_tauwmargin

constexpr double madness::XCfunctional::default_tauwmargin =1.e-6

staticconstexprprotected

Referenced by reset_screening_defaults().

◆ hf_coeff

double madness::XCfunctional::hf_coeff =0.0

protected

Factor multiplying HF exchange (+1.0 gives HF)

Referenced by hf_exchange_coefficient(), initialize(), and is_lda().

◆ nderiv

int madness::XCfunctional::nderiv =0

protected

Jacob's ladder rung; 0: lda, 1: gga, 2: mgga.

◆ number_xc_args

const int madness::XCfunctional::number_xc_args =43

static

max number of intermediates

Referenced by madness::nemo_u1_functors::append(), madness::XCOperator< T, NDIM >::prep_xc_args(), and test_meta_gga_dedtau_polarized().

◆ rhomin

double madness::XCfunctional::rhomin =default_rhomin

protected

what munge() puts in place of a density

Referenced by initialize(), munge(), and reset_screening_defaults().

◆ rhotol

double madness::XCfunctional::rhotol =default_rhotol

protected

See initialize and munge*.

Referenced by get_rhotol(), initialize(), munge(), and reset_screening_defaults().

◆ spin_polarized

bool madness::XCfunctional::spin_polarized =false

protected

True if the functional is spin polarized.

Referenced by exc(), initialize(), is_spin_polarized(), and vxc().

◆ tautol

double madness::XCfunctional::tautol =default_tautol

protected

floor for the kinetic energy density, see initialize

Referenced by get_tautol(), and reset_screening_defaults().

◆ tauwmargin

double madness::XCfunctional::tauwmargin =default_tauwmargin

protected

von Weizsaecker clamp overshoot, see tau_w_bound

Referenced by get_tauwmargin(), reset_screening_defaults(), and tau_w_bound().

The documentation for this class was generated from the following files:

Public Types

Public Member Functions

Static Public Member Functions

Static Public Attributes

Protected Member Functions

Static Protected Member Functions

Protected Attributes

Static Protected Attributes

Private Member Functions

Detailed Description

Member Enumeration Documentation

◆ xc_arg

Constructor & Destructor Documentation

◆ XCfunctional()

◆ ~XCfunctional()

Member Function Documentation

◆ binary_munge()

◆ exc()

◆ fxc_apply()

◆ get_rhotol()

◆ get_tautol()

◆ get_tauwmargin()

◆ has_fxc()

◆ has_kxc()

◆ hf_exchange_coefficient()

◆ initialize()

◆ is_dft()

◆ is_gga()

◆ is_lda()

◆ is_meta()

◆ is_spin_polarized()

◆ make_libxc_args()

◆ munge()

◆ munge_old()

◆ needs_sigma()

◆ needs_tau()

◆ plot()

◆ polyn()

◆ reset_screening_defaults()

◆ tau_w_bound()

◆ uses_libxc_backend()

◆ vxc()

Member Data Documentation

◆ default_rhomin

◆ default_rhotol

◆ default_tautol

◆ default_tauwmargin

◆ hf_coeff

◆ nderiv

◆ number_xc_args

◆ rhomin

◆ rhotol

◆ spin_polarized

◆ tautol

◆ tauwmargin