Ground-state density-functional theory (DFT)

The QMol-grid package uses the Kohn-Sham formulation of density functional theory (DFT) [Kohn 1965], where the electronic structure of a polyelectronic atom or molecule is described, in atomic units, by the nonlinear system eigen-state problem

$$ {\hat{\mathcal{H}} }_{{\mathrm{D}\mathrm{F}\mathrm{T}}} \left\lbrack \lbrace \phi_k \rbrace_k \right\rbrack ~\phi_k (x)=E_k ~\phi_k (x)~~{\mathrm{f}\mathrm{o}\mathrm{r}}~~k=1,2,\ldots~~~~~~(1) $$

often referred to as the self-consistent field (SCF) equation. Formally, the DFT Hamiltonian operator ${\hat{\mathcal{H}} }_{{\mathrm{D}\mathrm{F}\mathrm{T}}}$ is only functionally dependent on the one-body density $\rho$ but some DFT functionals still involve the Kohn-Sham orbitals in their definition. Traditionally, the DFT Hamiltonian is decomposed between the kinetic, external, Hartree, and exchange-correlation operators, each associated with an energy functional

$$ {\hat{\mathcal{H}} } _{{\mathrm{D}\mathrm{F}\mathrm{T}}} =-\frac{\Delta }{2}+{\hat{\mathcal{V}} } _{{\mathrm{e}\mathrm{x}\mathrm{t}}} +{\hat{\mathcal{V}} } _{{\mathrm{H}}} +{\hat{\mathcal{V}} } _{{\mathrm{X}\mathrm{C}}} .~~~~~~(2) $$

The QMol-grid package provides support both for spin-polarized, where the up- and down-spin electronic components may be different, and spin-polarized, where they are forced to be identical DFT models.

• QMol_DFT_spinPol specifically describes spin polarized DFT systems.
• QMol_DFT_spinRes specifically describes spin restricted DFT systems.

The two classes share very similar interfaces to provide a consistent end-user experience -- see their respective documentation pages for a detailed desciption of the classes and features they support.

DFT components

For clarity and streamlined future development of DFT capabilities in the QMol-grid package, we sort DFT component (classes) into different topical groups. Each component is equally compatible with both spin-polarized and -restricted DFT models

Model discretization

The following classes are used to discretize DFT systems

• QMol_DFT_density is the Cartesian-grid discretization of a one-body density. It also provides support for the computation of density gradient.
• QMol_DFT_orbital is the Cartesian-grid discretization of the Kohn-Sham orbitals.
• QMol_DFT_Vks is the Cartesian-grid discretization of the explicit part of the Kohn-Sham potential. Implicit components are also supported and provide an interface for the computation of the action of the potential operator on arbitrary wave functions.
• QMol_DFT_Vks_grad is the Cartesian-grid discretization of the explicit part of the Kohn-Sham potential gradient. Similar to the Kohn-Sham potential, it also handles implicit components and the computation of the action of the potential gradient operator on arbitrary wave functions.

For basis-set discretization models, with QMol_disc_basis, some of the classes above are redefined (overloaded)

• QMol_DFT_orbital_basis for Kohn-Sham orbitals.
• QMol_DFT_Vks_basis for the Kohn-Sham potential.

Model components

The following functionals are currently available for DFT models and simulations in the QMol-grid package

• QMol_DFT_Vext describes the external functional and potential ${\hat{\mathcal{V}} }_{{\mathrm{e}\mathrm{x}\mathrm{t}}}$ of Eq. (2).
• QMol_DFT_Vh_conv describes the Hartree functional and potential ${\hat{\mathcal{V}} }_{{\mathrm{H}}}$ of Eq. (2), using an explicit-convolution scheme.
• QMol_DFT_Vh_fft describes the Hartree functional and potential ${\hat{\mathcal{V}} }_{{\mathrm{H}}}$ of Eq. (2), using a fast-Fourier-transform convolution scheme.
• QMol_DFT_Vx_LDA_exp describes local-density-approximation Slater exchange ${\hat{\mathcal{V}} }_{{\mathrm{X}\mathrm{C}}}$ of Eq. (2), for exponential electron-electron interaction potentials.
• QMol_DFT_Vx_LDA_soft describes local-density-approximation Slater exchange ${\hat{\mathcal{V}} }_{{\mathrm{X}\mathrm{C}}}$ of Eq. (2), for soft-Coulomb electron-electron interaction potentials.
• QMol_DFT_Vx_XX_conv describes the exact-exchange functional and potential of Eq. (2), using an explicit-convolution scheme.
• QMol_DFT_Vx_XX_fft describes the exact-exchange functional and potential of Eq. (2), using a fast-Fourier-transform convolution scheme.
• QMol_DFT_Vc_LDA_soft describes local-density-approximation correlation ${\hat{\mathcal{V}} }_{{\mathrm{X}\mathrm{C}}}$ of Eq. (2), for soft-Coulomb electron-electron interaction potentials.

All DFT functionals can be used both for regular grid-based and basis-set deiscretization models, using the same end-user interfaces.

Memory and execution-time profiler

The choice of the domain discretization and DFT components can greatly affect the resources, both in terms of memory and execution time, required to carry out simulations.

• QMol_DFT_profiler provides estimates of the memory and run-time requirements for a given DFT model or components.

Keep in mind that differentiation and (some) convolutions are performed via fast-Fourier transforms. Thus, as a rule of thumb, domain discretization with small prime factor number of grid points tend to produce faster results -- for best performance, prime factors not greater than 7.

Ground-state computation

The DFT ground state, or more generally some DFT eigen state, is computed by numerically solving the SCF problem of Eq. (1). In the QMol-grid package, this can be achieved with

• QMol_DFT_SCF_Anderson performs SCF iterations for DFT ground/eigen-state calculations using an Anderson's mixing scheme.
• QMol_DFT_eigs is a linear eigen solver, to be used in combination with SCF iterations.
• QMol_DFT_eig_basis is a linear eigen solver for basis-set DFT models, to be used in combination with SCF iterations.

Self-interaction correction (SIC)

Various approaches have been proposed to provide self-interaction correction to the DFT Eq. (1). Currently, the QMol-grid package supports:

No self-interaction correction

In this case, self-interaction is ignored and no SIC is provided. This is the default options for all DFT models in the QMol-grid package.

Average-density self-interaction correction (ADSIC)

ADSIC implements an average-density level of SIC [Legrand 2002] that is easy to implement in the general DFT framework of Eq. (1). Unless otherwise specified, all Hartree and exchange-correlation DFT functionals in the QMol-grid package support ADSIC.

Notes for developers

DFT model implementation

The QMol_DFT abstract class defines the common interface and run-time documentation for both spin-polarized (QMol_DFT_spinPol) and spin-restricted (QMol_DFT_spinRes) DFT models.

Many DFT components and methods require convergence tolerance or thresholds in their implementation. Throughout, we parameterize these through class properties (generally named tol). The default configuration of the QMol-grid package roughly aims for $10^{-10}$ as the target tolerance/threshold. Changing the values for the tol parameters may have unintended side effects -- see the respective documentation pages for more details on these tolerance/threshold parameters.

Exchange-correlation functionals

To be compatible with other DFT components of the QMol-grid package, exchange-correlation functionals must fulfill the following requirements:

Class definition

Define exchange-correlation classes by overloading the QMol_suite class. This will facilitate name-value pair assignment for property values in the constructor, set, and clear methods -- see the QMol_suite documentation page for more details.

Class properties

Exchange-correlation classes should define the following member properties (with the specified attributes)

• DFT: with recommended attributes Transient,Hidden,GetAccess=?QMol_suite,SetAccess={?QMol_DFT_###,?QMol_DFT}, where QMol_DFT_### is the name of the exchange-correlation class being implemented, and at minimum such that QMol_DFT can edit the value of DFT. This feature is used by QMol_DFT_profiler to estimate the memory footprint of DFT models without allocating its components.
• type: as a constant character array with public access. This is used by TDDFT propagators to check that the DFT functionals used in the model are supported. Possible values are 'GGA_C' (for GGA exchange), 'GGA_X' (for GGA exchange), 'GGA_XC' (for GGA exchange correlation), 'LDA_C' (for LDA exchange), 'LDA_X' (for LDA exchange), 'LDA_XC' (for LDA exchange correlation), 'XX' (for exact exchange).

Run-time documentation

Implement run-time documentation with the follwing interface

methods (Access=public)
function ref = showDocumentation(obj)
%showDocumentation displays the documentation reflecting member property 
%   values
    
    % Add run-time documentation

end
end

Each exchange-correlation functional class must implement a run-time documentation as it is automatically called by QMol_DFT in its own run-time documentation.

• If no reference is cited, return ref = {} otherwise ref = {'ref 1','ref 2',___}, where each of the ref corresponds to the bibliographic code for the reference (generally of the form 'name_of_first_author publication_year') and is case insensitive -- see the documentation page for a list of citable references.

Additionally, exchange-correlation functional class should implement an estimator of its memory footprint with the following interface

methods (Access=public)
function mem = getMemoryProfile(obj,opt)
%getMemoryProfile computes and returns an estimate of the total memory
%   footprint (in bytes) of the DFT object with all its components
%   initialized and used.
    
    % Add memory footprint estimate

end
end

• getMemoryProfile must be able to return the memory footprint estimate even if the object has not been initialized. To do so, it can rely on its internal copy of the parent DFT object, that is properly assigned by QMol_DFT_profiler before calling getMemoryProfile.
• Only properties with large memory footprint -- e.g., the discretization of internal properties on the domain discretization grid -- need to be listed. Other small properties (scalar, object handle, flags, etc.) can be ignored; if the class only uses such small properties, it may return mem = 0.
• Use the QMol_DFT_profiler.getMemoryFootprint(N,'real') and QMol_DFT_profiler.getMemoryFootprint(N,'imag') methods to convert the size (number of elements N) of the memory components to the number of bytes they take in memory.
• opt = false requests the memory footprint with no details printed out.
• opt = true should additionally display the details of the memory estimate. Use QMol_DFT_profiler.showMemoryFootprint(msg,mem,lvl) to do this, with msg (character array) the description of the component, mem the estimated number of bytes associated with it, and lvl the indentation level (optional, lvl = 1 and lvl = 2 are the only two levels supported).

Class initialization

Initialize the class with the following interface

methods (Access=public)
function initialize(obj,DFT,SIC)
%initialize initializes the object
    
    % Initialize the functional object

end
end

• the input DFT is a QMol_DFT object containing the DFT model associated with the exchange-correlation functional object to be initialized. For run-time considerations it is often practical to keep a local copy of DFT in the object in the DFT member method.
• The input SIC specifies the flavor of SIC to be used with the exchange-correlation functional. 'none' notify to disable SIC while 'ADSIC' specifies the ADSIC method [Legrand 2002].
• Both DFT and SIC information are systematically passed along by the parent DFT model object when it initializes its exchange-correlation functional components.

Exchange-correlation functional kernel

Implicit functionals may require specific initialization of the exchange-correlation functional object before they can compute the energy, potential, or potential gradient. This is done with

methods (Access=public)
function setPotentialKernel(obj)
%setPotentialKernel sets the kernel to be used for the computation of 
%   exact-exchange potentials

    % Proper initialization
end
end

This typically involves storing a local copy of the Kohn-Sham orbitals of the parent DFT model. If no kernel initialization is required, include the following block in the class definition

methods (Access=public)
function setPotentialKernel(~)
%setPotentialKernel no kernel to set
    
end
end

Exchange-correlation energy

Compute and return the exchange-correlation energy with

methods (Access=public)
function E = getEnergy(obj,rho)
%getEnergy returns the exchange-correlation energy
    
    % Exchange-correlation energy computation
end
end

The getEnergy method should support the following computation modes -- see one of the available exchange-correlation functional documentation page for more details.

Get the exchange-correlation energy for the parent DFT object

E = obj.getEnergy;

• This computes the exchange energy associated with the Kohn-Sham orbitals in the parent DFT model. To do so, it should first compute the one-body density associated with the Kohn-Sham orbitals and their occupation parameters.
• This should be equivalent to, but more efficient than, obj.getEnergy(DFT.getDensity).

Get the exchange-correlation energy for a specific one-body density

E = obj.getEnergy(rho);

• Same as above but with using the input one-body-density object rho in the numerical computation of the energy.
• rho should be a QMol_DFT_density object defined over the same domain discretization as the parent DFT model.

Exchange-correlation potential

Compute and return the exchange-correlation potential with

methods (Access=public)
function E = getPotential(obj,rho)
%getPotential returns the exchange-correlation potential
    
    % Exchange-correlation potential computation
end
end

The getPotential method should support the following computation modes -- see one of the available exchange-correlation functional documentation page for more details.

Get the exchange-correlation potential for the parent DFT object with either

V = obj.getPotential;
V = obj.getPotential([]);

• This computes the exchange-correlation potential associated with the Kohn-Sham orbitals in the parent DFT model. To do so, it should first compute the one-body density $\rho$ associated with the Kohn-Sham orbitals and their occupation parameters.
• The output Kohn-Sham potential object contains the numerical evaluation of the exchange-correlation potential.
• For implicit exchange-correlation functionals, it provides the proper function handle to the Kohn-Sham potential class.
• This creates a new Kohn-Sham potential object V in which the exchange-correlation potential is stored.
• Note that getPotential should not initialize the output potential object V.

Get the exchange-correlation potential for a specific one-body density

V = obj.getPotential(rho);

• Same as above but with using the input one-body-density object rho in the numerical computation of the potential.
• rho should be a QMol_DFT_density object defined over the same domain discretization as the parent DFT model.

Overwrite the exchange-correlation potential in an existing Kohn-Sham potential object with any of

obj.getPotential([],V);         % use parent DFT density
obj.getPotential([],V,false);
obj.getPotential(rho,V);        % supply the density
obj.getPotential(rho,V,false);

• This is similar to above without creating a new Kohn-Sham potential object.
• Any content in the input object V is erased before assigning the external potential to it.

Add the exchange-correlation potential to an existing Kohn-Sham potential object

obj.getPotential([],V,true);    % use parent DFT density
obj.getPotential(rho,V,true);   % supply the density

• This is formally equivalent to the in-place addition $\mathcal{V}\gets \mathcal{V}+{\mathcal{V}}_{{\mathrm{X}\mathrm{C}}}$ .

Exchange-correlation potential gradient

Compute and return the exchange-correlation potential gradient with

methods (Access=public)
function E = getPotentialDerivative(obj,rho)
%getPotentialDerivative returns the exchange-correlation potential
%   gradient
    
    % Exchange-correlation potential computation
end
end

The getPotentialDerivative method should support the following computation modes -- see one of the available exchange-correlation functional documentation page for more details.

Get the exchange-correlation potential gradient for the parent DFT object with either

DV = obj.getPotentialDerivative(dim);
DV = obj.getPotentialDerivative(dim,[]);

• This computes the exchange-correlation potential gradient, along the dimension dim (1 for x, 2 for y, and 3 for z), associated with the Kohn-Sham orbitals in the parent DFT model. To do so, it should first compute the one-body density $\rho$ associated with the Kohn-Sham orbitals and their occupation parameters.
• The output Kohn-Sham potential gradient object contains the numerical evaluation of the exchange-correlation potential gradient.
• For implicit exchange-correlation functionals, it provides the proper function handle to the Kohn-Sham potential gradient class.
• This creates a new Kohn-Sham potential gradient object DV in which the exchange potential gradient is stored.
• Note that getPotentialDerivative does not initialize the output potential gradient object DV.

Get the exchange-correlation potential for a specific one-body density

V = obj.getPotentialDerivative(1,rho);

• Same as above but with using the input one-body-density object rho in the numerical computation of the potential gradient.
• rho should be a QMol_DFT_density object defined over the same domain discretization as the parent DFT model.

Overwrite the exchange-correlation potential gradient in an existing Kohn-Sham potential gradient object with any of

obj.getPotentialDerivative(1,[],DV);        % use parent DFT density
obj.getPotentialDerivative(1,[],DV,false);
obj.getPotentialDerivative(1,rho,DV);       % supply the density
obj.getPotentialDerivative(1,rho,DV,false);

• This is similar to above without creating a new Kohn-Sham potential gradient object.
• Any content in the input object DV is ereased before assigning the exchange-correlation potential gradient to it.

Add the exchange-correlation potential gradient to an existing Kohn-Sham potential gradient object with either

obj.getPotentialDerivative(1,[],DV,true);   % use parent DFT density
obj.getPotentialDerivative(1,rho,DV,true);  % supply the density

• This is formally equivalent to the in-place addition $\nabla \mathcal{V}\gets \nabla \mathcal{V}+\nabla {\mathcal{V}}_{{\mathrm{X}\mathrm{C}}}$ .

References

[Kohn 1965] W. Kohn and L.J. Sham, "Self-Consistent Equations Including Exchange and Correlation Effects," Physical Review 140, A1133 (1965).

[Legrand 2002] C. Legrand, E. Suraud, and P.-G. Reinhard, "Comparison of self-interaction-corrections for metal clusters," Journal of Physics B: Atomic, Molecular and Optical Physics 35, 1115 (2002).