HF - fmauger1/QMol-grid GitHub Wiki
Ground-state Hartree Fock (HF)
Hartree Fock (HF) is obtained by running density-functional theory (DFT) computations with (i) an exact-exchange functional, (ii) no correlation functional, and (iii) no self-interaction correction (SIC) in the Hartree functional. For details about these, we refer to the DFT documentation. Here, we only review the key steps to define HF models and converge ground-state calculations.
Hartree-Fock model
Since HF models are implemented via QMol-grid's DFT feature, spin polarized HF is described with QMol_DFT_spinPol while spin-restricted HF is with QMol_DFT_spinRes. For example, a minimal example of a molecular system with three atomic centers and 7 electrons (4 in the up- and 3 in the down-spin channels) is
% Domain
x = -20:.1:15;
% Molecule
A1 = QMol_Va_softCoulomb('name','atom 1','charge',3,'position',-3);
A2 = QMol_Va_softCoulomb('name','atom 2','charge',3,'position', 0);
A3 = QMol_Va_softCoulomb('name','atom 3','charge',1,'position', 1.5);
Vmol= QMol_DFT_Vext('atom',{A1,A2,A3});
% HF model
HF = QMol_DFT_spinPol( ...
'xspan', x, ...
'occupation', {[1 1 1 1],[1 1 1]},...
'externalPotential', Vmol, ...
'HartreePotential', QMol_DFT_Vh_conv, ...
'exchangeCorrelationPotential', QMol_DFT_Vx_XX_conv);
Next, we initialize the HF model and display its run-time documentation
HF.initialize;
HF.showDocumentation;
yielding
=== Theory ===============================================================
* Density-functional theory (DFT) spin polarized (SP)
SP-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 = -20:0.1:15
Size = 351 (3 x 3 x 3 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)
> atom 1, parameterized as (soft Coulomb)
Z = 3.00 | a = 1.00 | X0 = -3.00
> atom 2, parameterized as (soft Coulomb)
Z = 3.00 | a = 1.00 | X0 = 0.00
> atom 3, parameterized as (soft Coulomb)
Z = 1.00 | a = 1.00 | X0 = 1.50
* 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
Interaction pot. = @(x)1./sqrt(x.^2+2) (elec.-elec.)
V-01.21.000 (06/17/2024) F. Mauger
* Exact-exchange functional explicit convolution
Interaction pot. = @(x)1./sqrt(x.^2+2) (elec.-elec.)
V-01.21.000 (06/17/2024) F. Mauger
=== System model =========================================================
* Electronic structure Kohn-Sham orbitals
Up-spin occ. = 1.00 | 1.00 | 1.00 | 1.00
Down-spin = 1.00 | 1.00 | 1.00
Total charge = 4.00 (up) + 3.00 (down) = 7.00 (electrons)
=== 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).
[Kohn 1965] W. Kohn and L.J. Sham, "Self-Consistent Equations Including
Exchange and Correlation Effects," Phys. Rev. 140, A1133 (1965).
[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).
Note that the run-time documentation specifies a spin-polarized DFT model with only exact exchange, which is HF.
Ground-state computation
The DFT ground-state Anderson scheme implemented in the QMol-grid package is ill-equipped for converging systems with implicit exchange functionals like Hartree-Fock. Thus, ground-state HF computations can thus be tedious to converge and brute-force calculations from scratch will often fail. Some strategies to compute ground-state HF include:
- Reduce the
mixingcoefficient in theQMol_DFT_SCF_Andersonsolver, - Start the ground-state HF computation from another (converged) ground state using an approximate exchange functional like LDA, and
- Include self-interaction correction in the preliminary approximate exchange ground-state computation, to ensure that the Kohn-Sham potential has the appropriate tail.
Here is an example using all three strategies
% Change the exchange functional
HF.set('exchangeCorrelationPotential',QMol_DFT_Vx_LDA_exp, ...
'selfInteractionCorrection','ADSIC');
% Run ground-state computation
SCF = QMol_DFT_SCF_Anderson;
SCF.solveSCF(HF);
that produces
=== 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).
=== Funding ==============================================================
The original development of the QMol-grid toolbox, and its (TD)DFT
features, was supported by the U.S. Department of Energy, Office of
Science, Office of Basic Energy Sciences, under Award No. DE-SC0012462.
Addition of the (TD)SE features was supported by the National Science
Foundation under Grant No. 2207656.
=== SCF iterations =======================================================
Iter. Orbital energies (-eV) Error
----- ----------------------------------------------------- --------
0 35.87 | 29.66 | 14.60 | 7.30 1.645e-01
35.63 | 29.48 | 14.42 8.036e-02
1 37.52 | 30.75 | 15.47 | 7.18 3.360e-02
37.23 | 30.53 | 15.25 1.372e-02
2 37.39 | 30.71 | 15.38 | 7.21 7.966e-03
37.10 | 30.48 | 15.16 8.327e-03
3 37.41 | 30.70 | 15.40 | 7.21 5.343e-03
37.12 | 30.47 | 15.17 2.423e-03
(...)
22 37.43 | 30.72 | 15.41 | 7.23 1.793e-10
37.14 | 30.49 | 15.18 7.634e-11
23 37.43 | 30.72 | 15.41 | 7.23 5.057e-11
37.14 | 30.49 | 15.18 3.614e-11
24 37.43 | 30.72 | 15.41 | 7.23 3.789e-11
37.14 | 30.49 | 15.18 1.667e-11
------------------------------------------------------------------------
Kohn-Sham orbitals converged to tolerance
=== Orbital energies =====================================================
Orbital Occ. (elec.) Energy (-eV) Error(a.u.)
------- ------------ ------------ -----------
1 1.00 37.429 7.420e-12
2 1.00 30.723 9.359e-12
3 1.00 15.405 1.516e-11
4 1.00 7.227 2.128e-11
------- ------------ ------------ -----------
1 1.00 37.137 7.654e-12
2 1.00 30.489 9.546e-12
3 1.00 15.182 1.518e-11
------------------------------------------------------------------------
=== DFT-component energies ===============================================
Component Energy (a.u.) Energy (eV) [ spin up/down (eV) ]
----------- ------------- ------------- ---------------------
Kinetic 2.187 59.498 [ 34.90/ 24.60 ]
External -24.844 -676.037 [ -370.91/ -305.12 ]
Hartree 8.306 226.029
Exch.-corr. -0.553 -15.054
----------- ------------- -------------
Total -14.904 -405.564
----------------------------------------------
This approximate-exchange ground state can then be used as an initial guess in a new computation of the ground state reverting back to the exact exchange (also, we need to remove the self-interaction correction)
% Reinstate exact exchange (HF)
HF.set('exchangeCorrelationPotential',QMol_DFT_Vx_XX_conv, ...
'selfInteractionCorrection','none');
% Run ground-state calculation from current orbitals
HF.initialize;
SCF.set('mixing',.2);
SCF.solveSCF(HF,[],'DFT');
yielding
=== 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.2
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).
=== Funding ==============================================================
The original development of the QMol-grid toolbox, and its (TD)DFT
features, was supported by the U.S. Department of Energy, Office of
Science, Office of Basic Energy Sciences, under Award No. DE-SC0012462.
Addition of the (TD)SE features was supported by the National Science
Foundation under Grant No. 2207656.
=== SCF iterations =======================================================
Iter. Orbital energies (-eV) Error
----- ----------------------------------------------------- --------
0 42.46 | 35.74 | 18.47 | 8.09 1.954e-02
42.23 | 35.10 | 16.78 1.279e-02
1 42.45 | 35.75 | 18.49 | 8.27 3.148e-02
42.23 | 35.09 | 16.87 1.904e-02
2 42.45 | 35.77 | 18.51 | 8.37 4.091e-02
42.24 | 35.09 | 16.91 2.355e-02
3 42.42 | 35.75 | 18.48 | 8.40 3.792e-02
42.21 | 35.06 | 16.89 2.078e-02
(...)
85 42.34 | 35.64 | 18.37 | 8.31 4.602e-11
42.13 | 34.95 | 16.81 3.630e-12
86 42.34 | 35.64 | 18.37 | 8.31 2.242e-11
42.13 | 34.95 | 16.81 4.320e-12
87 42.34 | 35.64 | 18.37 | 8.31 1.807e-11
42.13 | 34.95 | 16.81 1.176e-11
------------------------------------------------------------------------
Kohn-Sham orbitals converged to tolerance
=== Orbital energies =====================================================
Orbital Occ. (elec.) Energy (-eV) Error(a.u.)
------- ------------ ------------ -----------
1 1.00 42.343 4.318e-12
2 1.00 35.638 4.522e-12
3 1.00 18.372 8.431e-12
4 1.00 8.313 7.899e-12
------- ------------ ------------ -----------
1 1.00 42.133 4.086e-12
2 1.00 34.949 4.348e-12
3 1.00 16.807 6.178e-12
------------------------------------------------------------------------
=== DFT-component energies ===============================================
Component Energy (a.u.) Energy (eV) [ spin up/down (eV) ]
----------- ------------- ------------- ---------------------
Kinetic 2.248 61.161 [ 35.64/ 25.53 ]
External -24.957 -679.120 [ -372.72/ -306.40 ]
Hartree 9.765 265.714
Exch.-corr. -2.058 -56.011
----------- ------------- -------------
Total -15.003 -408.256
----------------------------------------------
Finally, we plot the HF orbitals
figure; hold on
plot(HF.xspan,HF.KSO.KSOup(:,1),'-' ,'LineWidth',2,'DisplayName','MO #1 (up)')
plot(HF.xspan,HF.KSO.KSOup(:,2),'-' ,'LineWidth',2,'DisplayName','MO #2 (up)')
plot(HF.xspan,HF.KSO.KSOup(:,3),'-' ,'LineWidth',2,'DisplayName','MO #3 (up)')
plot(HF.xspan,HF.KSO.KSOup(:,4),'-' ,'LineWidth',2,'DisplayName','MO #4 (up)')
plot(HF.xspan,HF.KSO.KSOdw(:,1),'--','LineWidth',2,'DisplayName','MO #1 (down)')
plot(HF.xspan,HF.KSO.KSOdw(:,2),'--','LineWidth',2,'DisplayName','MO #2 (down)')
plot(HF.xspan,HF.KSO.KSOdw(:,3),'--','LineWidth',2,'DisplayName','MO #3 (down)')
xlabel('position (a.u.)'); xlim(HF.xspan([1 end]));
ylabel('orbital')
legend show
Basis-set discretization
We briefly illustrate the use of a basis-set dicretization for HF simulations. First we define the domain and set of vectors to use for building the basis set
% Domain
x = -20:.1:15;
% Atomic centers
A1 = QMol_Va_softCoulomb('name','atom 1','charge',3,'position',-3);
A2 = QMol_Va_softCoulomb('name','atom 2','charge',3,'position', 0);
A3 = QMol_Va_softCoulomb('name','atom 3','charge',1,'position', 1.5);
% Atomic orbital vectors
AO = @(s,n,x0,x) (x(:)-x0).^n .*exp(-(x(:)-x0).^2 * .5/s^2);
V = [AO(1.3,0,A1.position,x),AO(1.7,1,A1.position,x),AO(0.8,2,A1.position,x), ...
AO(1.3,0,A2.position,x),AO(1.7,1,A2.position,x),AO(0.8,2,A2.position,x), ...
AO(1.3,0,A3.position,x),AO(1.7,1,A3.position,x),AO(0.8,2,A3.position,x)];
At this point, the family of vectors in V are not orthonormal to each other. We correct for this using the built-in facilities of QMol_disc_basis
% Discretization domain
disc= QMol_disc_basis('x',x,'basis',V);
disc.initialize;
% Orthonormalize the basis
disc.set('basis',disc.DFT_normalizeOrbital(disc.basis)*sqrt(disc.dx));
disc.orthonormalizeBasis('overlap');
The first DFT_normalizeOrbital individually normalizes the vectors in V and is optional -- note that we multiply the result by the square root of the grid stencil. This is done to ensure that all vectors are treated on an equal footing when we orthonormalize the basis.
We can now define the basis-set HF model and look at its run-time documentation
% HF model
HF = QMol_DFT_spinPol( ...
'discretization', disc, ...
'occupation', {[1 1 1 1],[1 1 1]}, ...
'externalPotential', QMol_DFT_Vext('atom',{A1,A2,A3}), ...
'HartreePotential', QMol_DFT_Vh_conv, ...
'exchangeCorrelationPotential', QMol_DFT_Vx_XX_conv);
% Run-time documentation
HF.initialize;
HF.showDocumentation;
yielding (we only show the sections that are different from the grid discretization above)
(...)
=== Discretization =======================================================
* Domain discretization basis set on a Cartesian grid
Grid = -20:0.1:15
Size = 351 (3 x 3 x 3 x 13) points
Basis = 9 vectors
V-01.21.000 (06/17/2024) F. Mauger
(...)
We can also look at the memory profiling for the basis-set model
HF.getMemoryProfile(true);
* Domain discretization
> domain 2.7 KB
> gradient 5.5 KB
> Laplacian 2.7 KB
> Laplacian (velocity gauge) 5.5 KB
> basis 24.7 KB
> gradient matrix 648 B
> Laplacian matrix 648 B
> Laplacian matrix (velocity gauge) 1.3 KB
> reconstructed orbitals 38.4 KB
* Kohn-Sham orbitals (basis set)
> spin up 576 B
> spin down 432 B
* One-body density
> density 5.5 KB
> gradient 5.5 KB
* Kohn-Sham potential 1.3 KB
* External functional
> potential 2.7 KB
> potential gradient 2.7 KB
* Hartree functional 5.5 KB
* Exact-exchange functional (conv.)
> interaction potential 5.5 KB
> spin up kernel 21.9 KB
> spin down kernel 16.5 KB
Finally, we compute the basis-set HF ground state (here we can compute the HF ground state from scratch)
SCF = QMol_DFT_SCF_Anderson;
SCF.solveSCF(HF);
yielding
=== 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 eig function
using a direct diagonalization of the DFT Hamiltonian matrix.
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).
=== Funding ==============================================================
The original development of the QMol-grid toolbox, and its (TD)DFT
features, was supported by the U.S. Department of Energy, Office of
Science, Office of Basic Energy Sciences, under Award No. DE-SC0012462.
Addition of the (TD)SE features was supported by the National Science
Foundation under Grant No. 2207656.
=== SCF iterations =======================================================
Iter. Orbital energies (-eV) Error
----- ----------------------------------------------------- --------
0 41.52 | 34.87 | 17.18 | 8.34 2.101e-01
41.25 | 34.12 | 15.14 1.236e-01
1 43.35 | 35.91 | 18.24 | 8.05 4.128e-02
43.17 | 35.21 | 16.69 5.651e-02
2 43.07 | 35.85 | 18.37 | 7.85 4.669e-02
42.87 | 35.19 | 16.61 4.050e-02
3 42.97 | 35.83 | 18.43 | 7.76 7.246e-02
42.74 | 35.18 | 16.55 1.049e-02
(...)
68 42.31 | 35.43 | 18.07 | 7.71 7.190e-11
42.08 | 34.79 | 16.11 4.544e-11
69 42.31 | 35.43 | 18.07 | 7.71 1.536e-10
42.08 | 34.79 | 16.11 3.309e-10
70 42.31 | 35.43 | 18.07 | 7.71 2.395e-11
42.08 | 34.79 | 16.11 3.479e-11
------------------------------------------------------------------------
Kohn-Sham orbitals converged to tolerance
=== Orbital energies =====================================================
Orbital Occ. (elec.) Energy (-eV) Error(a.u.)
------- ------------ ------------ -----------
1 1.00 42.308 1.418e-11
2 1.00 35.433 1.434e-11
3 1.00 18.069 2.562e-11
4 1.00 7.708 2.046e-11
------- ------------ ------------ -----------
1 1.00 42.080 1.432e-11
2 1.00 34.791 1.519e-11
3 1.00 16.107 2.626e-11
------------------------------------------------------------------------
=== DFT-component energies ===============================================
Component Energy (a.u.) Energy (eV) [ spin up/down (eV) ]
----------- ------------- ------------- ---------------------
Kinetic 2.246 61.123 [ 36.04/ 25.09 ]
External -24.938 -678.609 [ -373.13/ -305.48 ]
Hartree 9.801 266.685
Exch.-corr. -2.065 -56.190
----------- ------------- -------------
Total -14.957 -406.991
----------------------------------------------
Finally, we plot the result
% Reconstruct orbitals on grid
KSO = disc.DFT_reconstructOrbital(HF.KSO);
% Plot results
figure; hold on
plot(HF.xspan,KSO.KSOup(:,1),'-' ,'LineWidth',2,'DisplayName','MO #1 (up)')
plot(HF.xspan,KSO.KSOup(:,2),'-' ,'LineWidth',2,'DisplayName','MO #2 (up)')
plot(HF.xspan,KSO.KSOup(:,3),'-' ,'LineWidth',2,'DisplayName','MO #3 (up)')
plot(HF.xspan,KSO.KSOup(:,4),'-' ,'LineWidth',2,'DisplayName','MO #4 (up)')
plot(HF.xspan,KSO.KSOdw(:,1),'--','LineWidth',2,'DisplayName','MO #1 (down)')
plot(HF.xspan,KSO.KSOdw(:,2),'--','LineWidth',2,'DisplayName','MO #2 (down)')
plot(HF.xspan,KSO.KSOdw(:,3),'--','LineWidth',2,'DisplayName','MO #3 (down)')
xlabel('position (a.u.)'); xlim(HF.xspan([1 end]));
ylabel('orbital')
legend show
Notes
The results displayed in this documentation page were generated using version 01.21 of the QMol-grid package.