Tutorial__GS__T05_GS_DFT - fmauger1/QMol-grid GitHub Wiki
Tutorial 5: Density-functional theory ground state
This tutorial walks through the process of setting up and calculating the density-functional theory (DFT) ground state of a close-shell one-dimensional molecular model with 3 atomic centers and 6 active electrons.
Table of Contents
Introduction
The DFT ground-state corresponds to the lowest energy solution to the nonlinear eigen-value system
$$ {\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) $$
where ${\hat{\mathcal{H}} }_{{\mathrm{D}\mathrm{F}\mathrm{T}}}$ is the DFT hamiltonian operator. ${\hat{\mathcal{H}} }_{{\mathrm{D}\mathrm{F}\mathrm{T}}}$ is functionally dependent on the Kohn-Sham orbitals $\phi_k$ via the one-body density $\rho$ . 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{DFT} = -\frac{\Delta }{2} + \hat{\mathcal{V}} _\mathrm{ext} + \hat{\mathcal{V}} _\mathrm{H} + \hat{\mathcal{V}} _\mathrm{XC} .~~~~~~(2)$$
In this tutorial, we illustrate how to use the QMol-grid package to calculate the DFT ground-state orbitals of a close-shell one-dimensional molecular model with 3 atomic centers and 6 active electrons (thus 3 Kohn-Sham orbitals). The general workflow is similar to the ground-state Schrödinger equation calculation of tutorial 1 and we only highlight steps that are specific to the DFT calculation.
Setting up the DFT system
Molecular model
We build the molecular model out of three one-dimensional atomic models, each contributing 2 electrons to the molecule, using soft-Coulomb potentials $V_k (x)=-Z_k /\sqrt{(x-x_k )^2 +1}$ [Javanainen 1988]
V_1 = QMol_Va_softCoulomb('atom','X_1','charge',2,'position',-3);
V_2 = QMol_Va_softCoulomb('atom','X_2','charge',2,'position', 0);
V_3 = QMol_Va_softCoulomb('atom','X_3','charge',2,'position', 3);
with the charges $Z_k =2$ and centered around the positions $x_k =-3$ , $0$ , and $3$ .
Then, we obtain the molecular potential -- the so-called external potential ${\hat{\mathcal{V}} }_{{\mathrm{e}\mathrm{x}\mathrm{t}}}$ in DFT -- as
Vext = QMol_DFT_Vext('atom',{V_1,V_2,V_3});
DFT potentials
To complete our DFT model, we define the Hartree ${\hat{\mathcal{V}} }_{{\mathrm{H}}}$ and exchange\-correlation ${\hat{\mathcal{V}} }_{{\mathrm{X}\mathrm{C}}}$ potential objects. For the latter, we use a local-density approximation for the exchange and correlation functionals.
Vh = QMol_DFT_Vh_conv;
Vxc = {QMol_DFT_Vx_LDA_soft,QMol_DFT_Vc_LDA_soft};
Here we use a local-density approximation (LDA) Slater-exchange and correlation functionals, which we put together in a cell array for the exchange-correlation potential Vxc.
DFT model
We now have all the elements to define a DFT model object with
DFT = QMol_DFT_spinRes( ...
'xspan', -50:.1:50, ...
'occupation', [2 2 2], ...
'externalPotential', Vext, ...
'HartreePotential', Vh, ...
'exchangeCorrelationPotential', Vxc, ...
'selfInteractionCorrection', 'ADSIC' );
The domain xspan ranges from -50 to 50 atomic units with a discretization step of 0.1 a.u. The occupation parameter specifies the number of electrons in each of the Kohn-Sham orbitals, here with all orbitals fully occupied with 2 electrons, yielding a net total of 6 electrons. Finally with last line command, we impose an average-density self-interaction correction (ADSIC).
Like for Schrödinger-equation object, we can print out the configuration for the DFT object by calling its run-time documentation, which requires for object to be initialized first.
DFT.initialize;
DFT.showDocumentation;
yielding (removing the license and funding disclaimers)
(...)
=== External components ==================================================
* convertUnit
V-01.02.000 (09/30/2022) F. Mauger
* fourierTool
V-01.02.002 (07/18/2023) F. Mauger
=== Theory ===============================================================
* Density-functional theory (DFT) spin restricted (SR)
SR-DFT level of theory, using the Kohn-Sham formalism [Kohn 1695].
V-01.21.000 (06/17/2024) F. Mauger
* The variational DFT model matches the canonical Hamiltonian formalism
describing the time evolution of the system [Mauger 2024].
V-01.21.000 (06/17/2024) F. Mauger
=== Discretization =======================================================
* Domain discretization Cartesian grid
Grid = -50:0.1:50
Size = 1001 (7 x 11 x 13) points
V-01.21.000 (06/17/2024) F. Mauger
=== External potential ===================================================
* External-potential functional
V-01.21.000 (06/17/2024) F. Mauger
* Atom-like center(s)
> X_1, parameterized as (soft Coulomb)
Z = 2.00 | a = 1.00 | X0 = -3.00
> X_2, parameterized as (soft Coulomb)
Z = 2.00 | a = 1.00 | X0 = 0.00
> X_3, parameterized as (soft Coulomb)
Z = 2.00 | a = 1.00 | X0 = 3.00
* Soft-Coulomb potential [Javanainen 1988] (soft Coulomb)
Parameterized as V(x) = -Z ./ sqrt( (x-X0).^2 + a^2 ).
V-01.21.000 (06/17/2024) F. Mauger
=== DFT potential(s) =====================================================
* Hartree-potential functional explicit convolution
with average-density self-interaction correct. (ADSIC) [Legrand 2002].
Interaction pot. = @(x)1./sqrt(x.^2+2) (elec.-elec.)
V-01.21.000 (06/17/2024) F. Mauger
* Slater-exchange functional local-density approximation (LDA)
for a one-dimensional (1D) soft-Coulomb electron-electron interaction
potential, of the form
Vee(x) = -1 / sqrt( x^2 + 1.41^2 )
The LDA Slater exchange is approximated by the scaled energy-per-
particle, parameterized as [Mauger 2024]
eps_x(r_s) = -1/2 * (1+alpa*r_s) / (beta*r+2*alpa*m*r_s^2)
* log(1+beta*r+gamma*r.^m),
where r_s = 1/(2*rho) is the 1D Wigner-Seitz radius, and the
parameters alpha, beta, gamma and m are obtained by fit against the
exact LDA exchange potential, with average-density self-interaction
correction (ADSIC) [Legrand 2002].
V-01.21.001 (07/10/2024) F. Mauger
* Correlation functional local-density approximation (LDA)
for a one-dimensional (1D) soft-Coulomb electron-electron interaction
potential of the form Vee(x) = 1 / sqrt( x^2 + 1 ) [Helbig 2011], with
an average-density self-interaction correction (ADSIC) [Legrand 2002].
V-01.21.001 (07/10/2024) F. Mauger
=== System model =========================================================
* Electronic structure Kohn-Sham orbitals
Occupation = 2.00 | 2.00 | 2.00
Total charge = 6.00 (electrons)
=== References ===========================================================
[Helbig 2011] N. Helbig, J.I. Fuks, M. Casula, M.J. Verstraete,
M.A.L. Marques, I.V. Tokatly, and A. Rubio, "Density functional theory
beyond the linear regime: Validating an adiabatic local density
approximation," Physical Review A 83, 032503 (2011).
[Javanainen 1988] J. Javanainen, J.H. Eberly, and Q. Su, "Numerical
simulations of multiphoton ionization and above-threshold electron
spectra," Physical Review A 38, 3430 (1988).
[Kohn 1965] W. Kohn and L.J. Sham, "Self-Consistent Equations Including
Exchange and Correlation Effects," Phys. Rev. 140, A1133 (1965).
[Legrand 2002] C. Legrand, E. Suraud, and P.-G. Reinhard, "Comparison of
self-interaction-corrections for metal clusters," J. Phys. B 35, 1115
(2002).
[Mauger 2024] F. Mauger, C. Chandre, M.B. Gaarde, K. Lopata, and K.J.
Schafer, "Hamiltonian formulation and symplectic split-operator
schemes for time-dependent density-functional-theory equations of
electron dynamics in molecules," Communications in Nonlinear Science
and Numerical Simulation 129, 107685 (2024).
[Mauger 2024b] F. Mauger and C. Chandre, "QMol-grid: A MATLAB package
for quantum-mechanical simulations in atomic and molecular systems,"
arXiv:2406.17938 (2024).
(...)
The first four and last blocks ("External components", "Theory", "Discretization", "Potential", and "References") are the same as for Schrödinger-equation models -- see tutorial 1.
The fifth block ("DFT potential(s)") summarizes the Hartree and various exchange-correlation potentials included in the DFT model.
The sixth block ("System model") shows the electronic structure for the DFT model, including the number of electrons (charge) associated with each Kohn-Sham orbitals.
Ground-State Calculation
With the DFT object defined above, we move to calculating its associated ground-state (Kohn-Sham) molecular orbitals and energy using the two commands
SCF = QMol_DFT_SCF_Anderson;
SCF.solveSCF(DFT);
Unlike the Schrödinger-equation, the DF ground state of equation (1) is the solution of a nonlinear eigen-state problem, which is solved using self-consistent field (SCF) iterations. For our system, the SCF solver prints out (again, removing the disclaimers)
(...)
=== External components ==================================================
* convertUnit
V-01.02.000 (09/30/2022) F. Mauger
* fourierTool
V-01.02.002 (07/18/2023) F. Mauger
=== Build density-functional-theory (DFT) model ==========================
* Discretization OK
* Kohn-Sham orbitals OK
* External potential OK
* Hartree potential OK
* Exchange-correlation potential OK
=== Self-consistent-field (SCF) methods ==================================
* Eigen-state solver for DFT Hamiltonians MATLAB eigs function
Tolerance = 1e-12
Max. iter. = 300
Basis dim. = 100
V-01.21.000 (06/17/2024) F. Mauger
* Self-consistent-field (SCF) Anderson's mixing
DFT-SCF solver using an Anderson's mixing scheme [Anderson 1965], as
described in [Johnson 1988].
Tolerance = 1e-10
Max. iter. = 100
Mix. coeff. = 0.5
Mix. mode = density
Conv. test = density
V-01.21.000 (06/17/2024) F. Mauger
=== References ===========================================================
[Anderson 1965] D.G. Anderson, "Iterative procedures for nonlinear
integral equations," Journal of the ACM 12, 547 (1965).
[Johnson 1988] D. D. Johnson, "Modified Broyden's method for
accelerating convergence in self-consistent calculations," Physical
Review B 38, 12807 (1988).
[Mauger 2024b] F. Mauger and C. Chandre, "QMol-grid: A MATLAB package
for quantum-mechanical simulations in atomic and molecular systems,"
arXiv:2406.17938 (2024).
(...)
=== SCF iterations =======================================================
Iter. Orbital energies (-eV) Error
----- ----------------------------------------------------- --------
0 24.01 | 20.68 | 16.28 1.670e-01
1 24.73 | 20.87 | 16.51 9.203e-03
2 24.71 | 20.86 | 16.50 3.480e-03
3 24.71 | 20.86 | 16.50 9.523e-04
(...)
17 24.71 | 20.86 | 16.50 3.075e-10
18 24.71 | 20.86 | 16.50 1.178e-10
19 24.71 | 20.86 | 16.50 4.673e-11
------------------------------------------------------------------------
Kohn-Sham orbitals converged to tolerance
=== Orbital energies =====================================================
Orbital Occ. (elec.) Energy (-eV) Error(a.u.)
------- ------------ ------------ -----------
1 2.00 24.713 2.625e-12
2 2.00 20.859 2.124e-12
3 2.00 16.502 2.736e-12
------------------------------------------------------------------------
=== DFT-component energies ===============================================
Component Energy (a.u.) Energy (eV)
----------- ------------- -------------
Kinetic 1.058 28.801
External -16.466 -448.059
Hartree 5.583 151.923
Exch.-corr. -0.564 -15.358
----------- ------------- -------------
Total -10.389 -282.692
----------------------------------------------
##########################################################################
The main difference with ground-state Schrödinger-equation calculations is the "SCF iterations" block, which shows the progress in SCF iterations as they are being performed. These can be used to gauge whether the DFT ground-state calculation is converging while it is going on.
Plotting the result
At the end of the calculation, the ground-state Kohn-Sham orbitals are stored in the input DFT object, together with relevant information such as the domain discretization. For instance, solely relying on DFT, we can plot the ground-state Kohn-Sham molecular orbitals with
figure, hold on
plot(DFT.xspan,DFT.orbital.orbital(:,1),'-','LineWidth',2,'DisplayName','HOMO-2')
plot(DFT.xspan,DFT.orbital.orbital(:,2),'-','LineWidth',2,'DisplayName','HOMO-1')
plot(DFT.xspan,DFT.orbital.orbital(:,3),'-','LineWidth',2,'DisplayName','HOMO')
set(gca,'box','on','FontSize',12,'LineWidth',2)
xlabel('x (a.u.)')
ylabel('Kohn-Sham orbital (a.u.)')
xlim([-15 15])
legend show
producing (note that the ground-state calculation starts from a random seed and thus the resulting Kohn-Sham molecular orbitals are defined with an arbitrary sign that can change from calculation to calculation)
With this example, we see that the molecular orbitals are nested inside the object DFT.orbital, with the orbitals stored in the successive rows of DFT.orbital.orbital, from the deepest to the highest-energy one. We further explain how to interface input parameters and output results from calculations in the next tutorial.
References
[Javanainen 1988] J. Javanainen, J.H. Eberly, and Q. Su, "Numerical simulations of multiphoton ionization and above-threshold electron spectra," Physical Review A 38, 3430 (1988).
Notes
The results displayed in this tutorial were generated using version 01.21 of the QMol-grid package.