6.1. Preliminaries

The SCF module in PyBEST contains various solvers to solve the Schrödinger equation self-constistently. The same solvers are used in restricted and unrestricted Hartree-Fock calculations.

To perform an SCF calculation, PyBEST requires

  • In case of molecular calculation, the molecular geometry, basis set, and molecular Hamiltonian need to be present.

  • For model Hamiltonians, the corresponding Hamiltonian matrix elements have to be available.

See also Defining Basis Sets, Molecular Geometries, and Hamiltonians for more details.

Note

Before setting up an SCF calculation, make sure that you have followed the instructions in Defining Basis Sets, Molecular Geometries, and Hamiltonians.

In PyBEST, all electronic structure methods feature a similar structure. The individual steps of an SCF calculation are

  1. create an instance of the RHF or UHF class

  2. initiate the SCF optimization using a function call

After the SCF is converged, the final data can be processed to

  1. dump orbitals (for instance into molden files)

  2. pass output to post-HF methods

Throughout this documentation, all electronic structure methods will label all objects according to the PyBEST name convention. See also Input/Output Operations: the IOData Container for a list of used acronyms in PyBEST.

Note that all variables listed below have to be initialized and/or evaluated beforehand (see also Defining Basis Sets, Molecular Geometries, and Hamiltonians). Below, we assume that the following quantities have been defined prior entering the SCF module

lf:

the LinalgFactory used to handle all tensor contractions

occ_model:

the occupation model used

kin:

the kinetic energy integrals

ne:

the nucleus-electron attraction integrals

eri:

the two-electron repulsion integrals

external:

some external potential, for instance the nuclear repulsion term

olp:

the overlap integrals of the AO basis

orb_a:

the molecular orbitals for alpha spin (or for restricted orbitals)

orb_b:

the molecular orbitals for beta spin (if present)

6.2. The Hartree-Fock wrapper

PyBEST provides a wrapper that allows the user to easily define and perform Hartree-Fock calculations. It automatically generates some initial guess orbitals, chooses some DIIS flavor, and dumps and returns all Hartree-Fock output data required for restarts and post-processing.

6.2.1. Restricted Hartree-Fock calculations

The code snippet below shows how to perform a restricted Hartree-Fock calculation in PyBEST,

hf = RHF(lf, occ_model)

# the order of the arguments does not matter
hf_output = hf(kin, ne, eri, external, olp, orb_a)

Note

All arguments (all integrals and orbitals) can be passed in any order. The RHF class determines their type automatically.

By default, PyBEST performs a core Hamiltonian guess, where the one-electron Hamiltonian is diagonalized. Specifically, the Hamiltonian contains only all one-electron terms that have been passed to the function call hf(), while all two-electron interactions are neglected. Due to the lack of the electron-electron interaction terms, the eigenfunctions of the one-electron Hamiltonian can be obtained by simply diagonalizing the core Hamiltonian. These orbitals can then be used as guess orbitals for any SCF procedure.

Note

The Hartree-Fock method overwrites the initial guess orbitals orb_a with the corresponding canonical HF orbitals. It also updates the occupation numbers and the orbital energies.

The RHF class returns an instance of the IOData container class, which stores all important Hartree-Fock results as attributes. This container can be passed as argument to other (post-Hartree-Fock or post-processing) methods to provide for instance (guess) orbitals, reference energies, etc. Specifically, the hf_output output container possesses the following attributes by default

hf_output.converged:

if converged, returns True (boolean)

hf_output.olp:

the overlap integrals

hf_output.dm_1:

the Hartree-Fock 1-particle reduced density matrix (RDM) in the AO basis

hf_output.orb_a:

the canonical Hartree-Fock orbitals (stored as a true/deep copy)

hf_output.e_tot:

the total Hartree-Fock energy

hf_output.e_ref:

the energy of the Hartree-Fock reference determinant (equivalent to e_tot)

hf_output.e_kin:

the kinetic energy

hf_output.e_ne:

the nuclear attraction energy

hf_output.e_hartree:

the Coulomb energy

hf_output.e_x_hf:

the Hartree-Fock exchange energy

hf_output.e_core:

the energy associated with some external potential

By default, the results of a Hartree-Fock calculation are stored to disk using PyBEST’s internal format and can be used for restarts. The default filename is checkpoint_scf.h5.

Note

All checkpoint files are stored in PyBEST’s results directory called pybest-results. This directory can be changed globally to any directory specified by the user using the filemanager,

from pybest import filemanager

# "some_other_directory" will be created related to the PyBEST run script
filemanager.result_dir = "some_other_directory"

6.2.2. Unrestricted Hartree-Fock calculations

To perform unrestricted Hartree-Fock calculation in PyBEST, the above code snippet has to be modified only slightly,

# Converge UHF
hf = UHF(lf, occ_model)
hf_output = hf(kin, ne, eri, external, olp, orb_a, orb_b)

Thus, an instance of the UHF class has to be created and the beta orbitals have to be passed as an additional argument.

Note

PyBEST assumes that the first orbitals passed as arguments correspond to alpha electrons, while the second set of orbitals is associated with beta electrons.

The corresponding IOData container (return value) also contains the beta orbitals and the 1-particle RDM of the alpha and beta electrons as attribute,

hf_output.dm_1_a:

the Hartree-Fock 1-particle RDM in the AO basis for alpha electrons

hf_output.dm_1_b:

the Hartree-Fock 1-particle RDM in the AO basis for beta electrons

hf_output.orb_b:

the canonical Hartree-Fock orbitals (stored as a true/deep copy) for beta electrons

6.3. Hartree-Fock calculations in the presence of point charges and static embedding potentials

Restricted Hartree-Fock calculations can be done in PyBEST in the presence of user-specified point charges and static embedding potentials obtained from other programs.

6.3.1. Restricted Hartree-Fock with point charges

If you use this part of the code, please cite [chakraborty2023].

The compute_external_charges() and compute_point_charges() functions calculate the electrostatic interactions between the external point charges and the nuclei and electrons, respectively. These terms can be computed in the following way

from pybest.gbasis.gobasis import compute_nuclear_pc, compute_point_charges
# compute the interaction between point charges and nuclei
external_pc = compute_nuclear_pc(gobasis, charges)
# compute the interaction between point charges and electrons
# (a ``TwoIndex`` object)
pc = compute_point_charges(gobasis, charges)

The charges argument is an instance of ExternalCharges (see Loading point charges and embedding potentials from another file for more details).

Then, Hartree-Fock calculations in the presence of external point charges can be performed in the following way,

hf = RHF(lf, occ_model)
# The order of the arguments does not matter
hf_output= hf(kin, eri, external, external_pc, pc, olp, orb_a)

6.3.2. Restricted Hartree-Fock with static embedding

If you use this part of the code, please cite [chakraborty2023].

The matrix elements of the static embedding potential (\(\nu^{emb}(r)\)) represented in the atomic orbital basis set are calculated as,

\[\nu_{ij}=\bra{\phi_{i}}\nu^{emb}(r)\ket{\phi_{j}} \approx \sum_{k}w_{k}\nu^{emb}(r_{k})\phi_{i}(r_{k})\phi_{j}(r_{k})\]

where \(i,j\) labels are for atomic orbitals, and the embedding potential is approximated on a grid with grid points \(r_{k}\) and corresponding weights \(w_{k}\) [gomes2008embedding]. The compute_static_embedding() function determines the effect of the embedding potential. The following code-snippet is an example of how to use this function,

from pybest.gbasis.gobasis import compute_static_embedding
# a ``TwoIndex`` object
emb = compute_static_embedding(gobasis, embedding_potential)

and the computed terms can be added to construct the effective Hamiltonian in Hartree-Fock calculations as follows

hf = RHF(lf, occ_model)
# The order of the arguments does not matter
hf_output_emb = hf(kin, ne, eri, external, emb, olp, orb_a)

6.4. Steering the SCF optimization

The convergence behaviour of the SCF algorithm can be steered using keyword arguments. The procedure is similar for both unrestricted and restricted Hartree-Fock calculations.

6.4.1. Setting the convergence threshold

By default the convergence threshold is set to 1e-8. To modify this value, adjust the RHF or UHF function call by adding the threshold keyword argument as follows,

# set a tighter convergence threshold in RHF
hf_output = hf(kin, na, er, external, olp, orb_a, threshold=1e-10)

The maximum number of SCF iterations steps (default 128) can be adjusted using the maxiter keyword,

# set a tighter convergence threshold and increase number of iterations in RHF
hf_output = hf(kin, na, er, external, olp, orb_a, threshold=1e-10, maxiter=200)

6.4.2. Choosing the DIIS solver

PyBEST supports four different SCF algorithms:

  • A simple SCF solver without any convergence acceleration: plain or None

  • The traditional SCF solver employing Pulay mixing: cdiis (default)

  • An SCF solver using energy-based DIIS: ediis

  • An SCF solver combing Pulay mixing and energy-based DIIS: ediis2

If plain or None is selected, an ordinary SCF solver is chosen that does not feature any convergence acceleration techniques. Instead, in each iteration, the Fock matrix is built and then diagonalized. By default, the cdiis options is selected, which corresponds to the (traditional) commutator-based direct inversion of the iterative subspace (CDIIS) algorithm, which is also known as Pulay mixing [pulay1980]. This SCF flavour is typically very efficient, albeit sensitive to the initial guess. ediis invokes the energy direct inversion of the iterative subspace (EDIIS) method [kudin2002]. This DIIS flavour works well for the first few iterations, becomes, however, numerically unstable close to convergence. It is the preferred SCF algorithm in cases the initial guess is poor. Finally, the ediis2 method combines the CDIIS and EDIIS algorithms [kudin2002]. This hybrid method aims at joining the benefits of both approaches.

The DIIS flavor can be changed using the keyword diis in the function call (similar for both restricted and unrestricted orbitals),

# select plain solver
hf_output = hf(kin, na, er, external, olp, orb_a, diis='plain') # or diis=None
# select CDIIS (default)
hf_output = hf(kin, na, er, external, olp, orb_a, diis='cdiis')
# select EDIIS
hf_output = hf(kin, na, er, external, olp, orb_a, diis='ediis')
# select EDIIS2
hf_output = hf(kin, na, er, external, olp, orb_a, diis='ediis2')

The behavior of the DIIS algorithm can be adjusted using additional keywords. For instance the convergence threshold and number of iterations can be controlled as described above (Setting the convergence threshold). In addition, the number of DIIS vectors can be adjusted (default value is 6). To do so, use the keyword nvector,

# select CDIIS (default) with at most 8 vectors
hf_output = hf(kin, na, er, external, olp, orb_a, nvector=8)

6.5. Dumping SCF results

At each SCF iteration step, the current results (orbitals, densities, etc.) are dumped to disk and stored for later restarts, visualization, or post-processing. If converged, those results are updated one last time and stored in the checkpoint file checkpoint_scf.h5 of the pybest-results directory (by default).

After SCF convergence, you can also dump all Hartree-Fock results to user-specified file and file format. Currently, PyBEST supports the following file formats

  1. PyBEST’s internal format based on HDF5

  2. the Molden format: Generating molden files

The internal format has the advantage that as much information as provided by the user can be stored in full precision. For internal storage, PyBEST uses the IOData container (see Input/Output Operations: the IOData Container for more details).

For visualization purposes, the Molden format is preferred. Furthermore, the Molden format can be read by other quantum chemistry codes. Unfortunately, Molden only supports basis functions up to including G functions.

Since the return value of the RHF and UHF function call is an IOData container, dumping the corresponding (RHF or UHF) results to different file formats is particularly easy. To dump results to disk the to_file() method of the IOData container has to be used, where the file ending specifies the file format.

6.5.1. PyBEST’s internal format

To dump into PyBEST’s internal file format, use the extension .h5

# Dump data to file using PyBEST's internal format
hf_output.to_file('water-scf.h5')

6.5.2. The Molden format

To dump Molden files, we have to update the IOData container to include also the gobasis information (by default, it is omitted from the container)

# Write SCF results to a molden file
hf_output.to_file("water-scf.molden")

# Restart RHF
hf = RHF(lf, occ_model)

See also Filling the IOData container on how to modify an IOData container on the fly.

6.6. Restarting SCF calculations

6.6.1. Restarting from a checkpoint file

One straightforward way to restart a Hartree-Fock calculation from a PyBEST checkpoint file is to use the restart keyword, which specifies the file name of the previously dumped SCF solution. PyBEST will automatically read in the new orbitals and perform a projection.

# results are stored in pybest-results directory
hf_output = hf(
    kin,
    ne,
    eri,
    external,
    olp,
    orb_a,
    restart="pybest-results/checkpoint_scf.h5",
)

# Read in orbitals and modify them
restart = IOData.from_file("pybest-results/checkpoint_scf.h5")

6.6.2. Restarting from perturbed orbitals

In case you want to perturb the orbitals prior restart, like swapping their order (Swapping orbitals) or 2x2 rotations (Rotating orbitals), you have to read in the corresponding checkpoint file form disk manually, perform the orbital manipulation you desire, and then pass the modified orbitals to the RHF or UHF function call,

# Swap HOMO-LUMO orbitals (Python indexing, starts with 0)
swaps = np.array([[4, 5]])
# Only swap AO/MO coefficients
restart.orb_a.swap_orbitals(swaps, skip_occs=True, skip_energies=True)

# Do RHF calculation using those orbitals
hf = RHF(lf, occ_model)
hf_output = hf(kin, ne, eri, external, restart.olp, restart.orb_a)

Note

In order to guarantee that the RHF or UHF solvers swap the orbitals in question, skip_occs and skip_energies have to be set to True. Otherwise, orbital swapping has no effect at all.

6.7. Example Python scripts

The following, we list a basic example script for a restricted Hartree-Fock calculation for the water molecule.

Listing 6.1 data/examples/hf/rhf_water_dense.py
import numpy as np

from pybest import context
from pybest.gbasis import (
    compute_eri,
    compute_kinetic,
    compute_nuclear,
    compute_nuclear_repulsion,
    compute_overlap,
    get_gobasis,
)
from pybest.io import IOData
from pybest.linalg import DenseLinalgFactory
from pybest.occ_model import AufbauOccModel
from pybest.wrappers import RHF

# Hartree-Fock calculation
# ------------------------

# Load the coordinates from file.
# Use the XYZ file from PyBEST's test data directory.
fn_xyz = context.get_fn("test/water.xyz")

# Create a Gaussian basis set
factory = get_gobasis("cc-pvdz", fn_xyz)

# Create a linalg factory
lf = DenseLinalgFactory(factory.nbasis)

# Compute integrals
olp = compute_overlap(factory)
kin = compute_kinetic(factory)
ne = compute_nuclear(factory)
eri = compute_eri(factory)
external = compute_nuclear_repulsion(factory)

# Create alpha orbitals
orb_a = lf.create_orbital()

# Decide how to occupy the orbitals (5 alpha electrons)
occ_model = AufbauOccModel(factory)

# Converge RHF
hf = RHF(lf, occ_model)

# the order of the arguments does not matter
hf_output = hf(kin, ne, eri, external, olp, orb_a)

# Write SCF results to a molden file
hf_output.to_file("water-scf.molden")

# Restart RHF
hf = RHF(lf, occ_model)
# we still need to provide some inital guess orbitals
# results are stored in pybest-results directory
hf_output = hf(
    kin,
    ne,
    eri,
    external,
    olp,
    orb_a,
    restart="pybest-results/checkpoint_scf.h5",
)

# Read in orbitals and modify them
restart = IOData.from_file("pybest-results/checkpoint_scf.h5")

# Swap HOMO-LUMO orbitals (Python indexing, starts with 0)
swaps = np.array([[4, 5]])
# Only swap AO/MO coefficients
restart.orb_a.swap_orbitals(swaps, skip_occs=True, skip_energies=True)

# Do RHF calculation using those orbitals
hf = RHF(lf, occ_model)
hf_output = hf(kin, ne, eri, external, restart.olp, restart.orb_a)