7.1. Self-consistent field (SCF) methods

Modules: scf, pbc.scf, soscf

7.1.1. Introduction

In self-consistent field (SCF) methods, the electron interactions are treated in a mean-field way. Here, the SCF methods include both Hartree-Fock (HF) theory and Kohn-Sham (KS) density functional theory (DFT). This Chapter summarizes the general SCF capabilities of PySCF. For more details specific to DFT, see Section 7.2.

A minimal example of using the SCF module is as follows:

from pyscf import gto, scf
mol = gto.M(
    atom = 'H 0 0 0; F 0 0 1.1',  # in Angstrom
    basis = 'ccpvdz',
    symmetry = True,
mf = scf.HF(mol)

This will run a HF calculation for the hydrogen fluoride molecule using the default SCF settings.

7.1.2. Theory

In HF and KS-DFT methods, the ground-state wavefunction is approximated by a single Slater determinant,

\[\begin{split}|\Phi_0\rangle = \frac{1}{\sqrt{N!}} \begin{vmatrix} \psi_1(\mathbf{r}_1) &\psi_2(\mathbf{r}_1) &\dots &\psi_N(\mathbf{r}_1)\\ \psi_1(\mathbf{r}_2) &\psi_2(\mathbf{r}_2) &\dots &\psi_N(\mathbf{r}_2)\\ \vdots &\vdots &\ddots &\vdots\\ \psi_1(\mathbf{r}_N) &\psi_2(\mathbf{r}_N) &\dots &\psi_N(\mathbf{r}_N) \end{vmatrix} \;,\end{split}\]

where the molecular orbitals (MO) \(|\psi_i\rangle\) are obtained by solving the Hartree-Fock equtions:

\[\hat{F}|\psi_i\rangle = \varepsilon_i |\psi_i\rangle \;.\]

The Fock operator \(\hat{F}\), when represented within the atomic orbital (AO) basis, has the following form (in a spin unrestricted formalism):

\[F_{\mu\nu}^{\alpha} = h_{\mu\nu} + J_{\mu\nu} - K_{\mu\nu}^{\alpha} \;,\]


\[h_{\mu\nu} = \left( \mu | \hat{T} |\nu \right) + \left( \mu | \hat{V}_{nuc} |\nu \right) \;,\]

usually called the core Hamiltonian, is the sum of kinetic energy operator and nuclear attraction operator, and

\[J_{\mu\nu} = \sum_{\lambda\sigma} \left(\mu\nu|\lambda\sigma\right) \left(P_{\lambda\sigma}^{\alpha}+P_{\lambda\sigma}^{\beta}\right)\]


\[K_{\mu\nu}^{\alpha} = \sum_{\lambda\sigma} \left(\mu\lambda|\nu\sigma\right) P_{\lambda\sigma}^{\alpha}\]

define the Coulomb and exchange operators, respectively. In the equations above, \(\mathbf{P}\) labels the density matrix:

\[P_{\mu\nu}^{\alpha} = \sum_{i\in occ} C_{\mu i}^{\alpha} C_{i\nu}^{\alpha\dagger} \;,\]

where \(\mathbf{C}\) represents the MO coefficients:

\[|\psi^{\alpha}\rangle = \sum_{\mu} |\mu\rangle C_{\mu i}^{\alpha} \;.\]

Within the AO representation, the Hartree-Fock equations reduce to the Roothaan-Hall or Pople-Nesbet [1] equations:

\[\mathbf{F}^{\alpha} \mathbf{C}^{\alpha} = \boldsymbol{\varepsilon}^{\alpha} \mathbf{S} \mathbf{C}^{\alpha} \;,\]


\[S_{\mu\nu} = \left( \mu | \nu \right)\]

is the AO overlap matrix. Note that as the Coulomb and exchange operators depend on the electron density, solving the Hartree-Fock equation is a nonlinear procedure, which requires iteratively updating the MOs and the Fock operator until reaching self-consistency. Finally, with the converged MOs or density matrix, one may compute various ground-state properties. For example, the ground-state electronic energy is expressed as

\[E_{\rm HF} = \frac{1}{2} \left\{{\rm Tr}[\mathbf{h}(\mathbf{P}^{\alpha}+\mathbf{P}^{\beta})] + {\rm Tr}(\mathbf{F}^{\alpha}\mathbf{P}^{\alpha}) + {\rm Tr}(\mathbf{F}^{\beta}\mathbf{P}^{\beta}) \right\} \;.\] Periodic boundary conditions

PySCF also alows the user to perform SCF calculations for solids. With crystalline Gaussian-type AOs as the underlying single-partial basis (see Section 7.12.1), the molecular SCF code can be easily adapted to the cases where periodic boundary conditions (PBCs) are applied. Instead of solving only one set of Roothaan-Hall or Pople-Nesbet equtions for molecules, it is now necessary to solve them for each k point for solids:

\[\mathbf{F}(\mathbf{k}) \mathbf{C}(\mathbf{k}) = \boldsymbol{\varepsilon}(\mathbf{k}) \mathbf{S}(\mathbf{k}) \mathbf{C}(\mathbf{k}) \;,\]

where the Fock matrix is defined (within the restricted formalism) as

\[\mathbf{F}(\mathbf{k}) = \mathbf{T}(\mathbf{k}) + \mathbf{V}^{\rm PP}(\mathbf{k}) +\mathbf{J}(\mathbf{k}) - \frac{1}{2} \mathbf{K}(\mathbf{k}) + \mathbf{V}^{L+J}(\mathbf{k}) \;.\]

Here, \(\mathbf{V}^{\rm PP}\) denotes the pseudopotential contribution and \(\mathbf{V}^{L+J}\) deals with the divergence of local pseudopotential and Hartree potential (see below).

The one-electron overlap, kinetic energy, and local pseudopotential integrals are evaluated through numerical integrations on the real-space grid according to

\[S_{\mu\nu}(\mathbf{k}) = \int_\Omega d\mathbf{r} \phi_{\mu\mathbf{k}}^{*}(\mathbf{r}) \phi_{\nu\mathbf{k}}(\mathbf{r}) \;,\]
\[T_{\mu\nu}(\mathbf{k}) = -\frac{1}{2} \int_\Omega d\mathbf{r} \phi_{\mu\mathbf{k}}^{*}(\mathbf{r}) \boldsymbol{\nabla}_{\mathbf{r}}^2 \phi_{\nu\mathbf{k}}(\mathbf{r}) \;,\]


\[V_{\mu\nu}^{\rm L-PP}(\mathbf{k}) = \int_\Omega d\mathbf{r} \phi_{\mu\mathbf{k}}^{*}(\mathbf{r}) v^{\rm L-PP}(\mathbf{r}) \phi_{\nu\mathbf{k}}(\mathbf{r}) \;,\]

where \(\Omega\) labels the unit cell volume. The non-local part of the pseudopotential is computed in the reciprocal space:

\[V_{\mu\nu}^{\rm NL-PP}(\mathbf{k}) = \Omega \sum_{\mathbf{G},\mathbf{G}'} \phi_{\mu\mathbf{k}}^{*}(\mathbf{G}) v^{\rm NL-PP}(\mathbf{k}+\mathbf{G}, \mathbf{k}+\mathbf{G}') \phi_{\nu\mathbf{k}}(\mathbf{G}') \;,\]


\[v^{\rm NL-PP}(\mathbf{k}+\mathbf{G}, \mathbf{k}+\mathbf{G}') = \frac{1}{\Omega} \int d\mathbf{r} \int d\mathbf{r}' e^{-i(\mathbf{k}+\mathbf{G})\cdot\mathbf{r}} v^{\rm NL-PP}(\mathbf{r},\mathbf{r}') e^{ i(\mathbf{k}+\mathbf{G}^{'})\cdot\mathbf{r}'} \;.\]


The way that the pseudopotential integrals are computed differs in different density fitting schemes and for different pseudopotentials. Interested readers should refer to Section 7.12.2 and Section 7.12.3.

The Coulomb and exchange matrices are defined similarly as

\[J_{\mu\nu}(\mathbf{k}) = \int_{\Omega} d\mathbf{r} \phi_{\mu\mathbf{k}}^{*}(\mathbf{r}) v_{\rm H}(\mathbf{r}) \phi_{\nu\mathbf{k}}(\mathbf{r}) \;,\]


\[K_{\mu\nu}(\mathbf{k}) = \int_{\Omega} d\mathbf{r} \int d\mathbf{r}' \phi_{\mu\mathbf{k}}^{*}(\mathbf{r}) \frac{\rho(\mathbf{r}, \mathbf{r}')}{|\mathbf{r}-\mathbf{r}'|} \phi_{\nu\mathbf{k}}(\mathbf{r}') \;.\]

Here \(v_{\rm H}\) is the Hartree potential

\[v_{\rm H}(\mathbf{r}) = \frac{4\pi}{\Omega} \sum_{\mathbf{G}\neq \mathbf{0}} \frac{\rho(\mathbf{G})}{G^2} e^{i\mathbf{G}\cdot\mathbf{r}} \;,\]

and \(\rho(\mathbf{r}, \mathbf{r}')\) is the density matrix

\[\rho(\mathbf{r}, \mathbf{r}') = \sum_{\mathbf{k}} w_{\mathbf{k}} \sum_{\lambda\sigma} P_{\lambda\sigma}(\mathbf{k}) \phi_{\lambda\mathbf{k}}(\mathbf{r}) \phi_{\sigma\mathbf{k}}^{*}(\mathbf{r}') \;,\]

where \(w_{\mathbf{k}}\) represents the weight of each k point.

Note that the local part of the pseudopotential and the Hartree potential diverge at \(G=0\); however, their sum is not, which leads to the \(V^{\rm L+J}\) term (for charge neutral unit cell):

\[V_{\mu\nu}^{\rm L+J} (\mathbf{k}) = \frac{S_{\mu\nu}}{\Omega} \int d\mathbf{r} \left(v^{\rm L-PP}(\mathbf{r}) + \sum_{\alpha} \frac{Z_{\alpha}e^2}{r} \right) \;,\]

where \(Z_{\alpha}\) denotes the nuclear charge of the \(\alpha\)-th atom.


For details about how to compute the Coulomb (J) and exchange (K) integrals, see Section 7.12.2.

Finally, the total electronic energy differs from the molecular case only by a k-point summation:

\[E_{\rm HF} = \sum_{\mathbf{k}} w_{\mathbf{k}} E_{\rm HF}(\mathbf{k}) \;,\]


\[E_{\rm HF}(\mathbf{k}) = \frac{1}{2} \left\{ {\rm Tr}\left[\mathbf{h}(\mathbf{k}) (\mathbf{P}^{\alpha}(\mathbf{k})+\mathbf{P}^{\beta}(\mathbf{k}))\right] + {\rm Tr}\left[\mathbf{F}^{\alpha}(\mathbf{k}) \mathbf{P}^{\alpha}(\mathbf{k})\right] + {\rm Tr}\left[\mathbf{F}^{\beta}(\mathbf{k}) \mathbf{P}^{\beta}(\mathbf{k})\right] \right\} \;.\]

7.1.3. Initial guess

As the Roothaan-Hall and Pople-Nesbet equations are solved iteratively, an initial guess for the MOs or the density matrices must be supplied. Poor initial guess may cause slow convergence or even divergence of the procedure. Furthermore, when treating magnetic or open-shell systems, the initial guess must be carefully chosen in order to get the correct state.

There are several options available in PySCF for selecting the initial guess to solve the SCF problem. One can set the attribute mf.init_guess to the following values to generate the initial guess in different ways:

  • 'minao' (default)

    The initial guess density matrix is first generated based on the atomic natural orbital (ANO) basis [2][3][4][5][6][7], then projected onto the basis set used for the SCF calculation.

  • 'hcore'

    The core Hamiltonian is diagonalized to get the initial MOs.

  • 'atom'

    The initial guess density matrix is from the superposition of atomic HF density matrix. Commonly know as the ‘SAD’ method.

  • 'chk'

    Read the existing SCF results from the checkpoint file, then the density matrix is projected onto the basis set used for the new SCF calculation.

Alternatively, the user could manually set the initial guess density matrix for an SCF calculation by using the 'dm0' argument. For example, the followings script first computes the HF density matrix for \(\rm Cr^{6+}\) cation, which is then used as the initial guess for the HF calculation of \(\rm Cr\) atom.

# use cation to produce initial guess
mol = gto.Mole()
    symmetry = 'D2h',
    atom = [['Cr',(0, 0, 0)], ],
    basis = 'cc-pvdz',
    charge = 6,
    spin = 0,

mf = scf.RHF(mol)
dm1 = mf.make_rdm1()

mol.charge = 0
mol.spin = 6

mf = scf.RHF(mol)

More examples can be found in examples/scf/15-initial_guess.py.

7.1.4. Accelerating SCF convergence Direct Inversion in the Iterative Subspace (DIIS)

At convergence of an SCF calcuation, one should expect the density matrix commute with the Fock matrix:

\[\mathbf{SPF} - \mathbf{FPS} = \mathbf{0} \;.\]

Prior to convergence, it is possible to define an error vector as

\[\mathbf{e}_i \equiv \mathbf{S}\mathbf{P}_i\mathbf{F}_i - \mathbf{F}_i\mathbf{P}_i\mathbf{S} \;,\]

where \(\mathbf{F}_i\) is a linear combination of the Fock matrices in the previous SCF cycles:

\[\mathbf{F}_i = \sum_{k=i-L}^{i-1} c_k \mathbf{F}_k \;,\]

\(\mathbf{P}_i\) is obtained by diagonalizing \(\mathbf{F}_i\), and \(L\) is the size of the DIIS subspace, which can be modified by setting the mf.diis_space attribute (the default size is 8). The DIIS method [8][9] minimizes the square of the error vector with respect to the DIIS coefficients \(c_k\) under the constraint that \(\sum_k c_k = 1\). The Euler–Lagrange equation of such a constrained minimization problem reads:

\[\begin{split}\left( \begin{array}{cccc} \mathbf{e}_1\cdot\mathbf{e}_1 &\dots &\mathbf{e}_1\cdot\mathbf{e}_L &1 \\ \vdots &\ddots &\vdots &\vdots \\ \mathbf{e}_L\cdot\mathbf{e}_1 &\dots &\mathbf{e}_L\cdot\mathbf{e}_L &1 \\ 1 &\dots &1 &0 \end{array} \right) \left( \begin{array}{c} c_1 \\ \vdots \\ c_{L} \\ \lambda \end{array} \right) = \left( \begin{array}{c} 0 \\ \vdots \\ 0 \\ 1 \end{array} \right)\end{split}\]

PySCF also implements two other similar DIIS algorithms, namely, EDIIS [10] and ADIIS [11]. Interested readers should refer to the reference. An example of selecting different DIIS schemes can be found in examples/scf/24-tune_diis.py Co-iterative augmented hessian (CIAH) second order SCF solver [12]

7.1.5. Level shifting and smearing

7.1.6. Stability analysis

7.1.7. References


JA Pople and R_K Nesbet. Self-consistent orbitals for radicals. The Journal of Chemical Physics, 22(3):571–572, 1954.


Per-Olof Widmark, Per-Åke Malmqvist, and Björn O Roos. Density matrix averaged atomic natural orbital (ano) basis sets for correlated molecular wave functions. Theoretica chimica acta, 77(5):291–306, 1990.


Björn O Roos, Valera Veryazov, and Per-Olof Widmark. Relativistic atomic natural orbital type basis sets for the alkaline and alkaline-earth atoms applied to the ground-state potentials for the corresponding dimers. Theoretical Chemistry Accounts, 111(2-6):345–351, 2004.


Björn O Roos, Roland Lindh, Per-Åke Malmqvist, Valera Veryazov, and Per-Olof Widmark. Main group atoms and dimers studied with a new relativistic ano basis set. The Journal of Physical Chemistry A, 108(15):2851–2858, 2004.


Björn O Roos, Roland Lindh, Per-Åke Malmqvist, Valera Veryazov, and Per-Olof Widmark. New relativistic ano basis sets for transition metal atoms. The Journal of Physical Chemistry A, 109(29):6575–6579, 2005.


Björn O Roos, Roland Lindh, Per-Åke Malmqvist, Valera Veryazov, and Per-Olof Widmark. New relativistic ano basis sets for actinide atoms. Chemical physics letters, 409(4-6):295–299, 2005.


Björn O Roos, Roland Lindh, Per-Åke Malmqvist, Valera Veryazov, Per-Olof Widmark, and Antonio Carlos Borin. New relativistic atomic natural orbital basis sets for lanthanide atoms with applications to the ce diatom and luf3. The Journal of Physical Chemistry A, 112(45):11431–11435, 2008.


Péter Pulay. Convergence acceleration of iterative sequences. the case of scf iteration. Chemical Physics Letters, 73(2):393–398, 1980.


Peter Pulay. Improved scf convergence acceleration. Journal of Computational Chemistry, 3(4):556–560, 1982.


Konstantin N Kudin, Gustavo E Scuseria, and Eric Cances. A black-box self-consistent field convergence algorithm: one step closer. The Journal of chemical physics, 116(19):8255–8261, 2002.


Xiangqian Hu and Weitao Yang. Accelerating self-consistent field convergence with the augmented roothaan–hall energy function. The Journal of chemical physics, 132(5):054109, 2010.


Qiming Sun. Co-iterative augmented hessian method for orbital optimization. arXiv preprint arXiv:1610.08423, 2016.