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

create an instance of the RHF or UHF class
initiate the SCF optimization using a function call

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

dump orbitals (for instance into molden files)
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,

# Converge RHF
hf = RHF(lf, occ_model)
# the order of the arguments does not matter
hf_output = hf(kin, na, er, 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, na, er, 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. 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.3.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.3.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.4. 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

PyBEST’s internal format based on HDF5

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.4.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.4.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
# First, add gobasis as attribute
hf_output.gobasis = gobasis
# Ready to dump results
hf_output.to_file("water-scf.molden")

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

6.5. Restarting SCF calculations

6.5.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.

# 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,
    na,
    er,
    external,
    olp,
    orb_a,
    restart="pybest-results/checkpoint_scf.h5",
)

6.5.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,

# 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, na, er, 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.6. 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

#!/usr/bin/env python

import numpy as np

from pybest import context
from pybest.gbasis import get_gobasis
from pybest.gbasis.gobasis import (
    compute_eri,
    compute_kinetic,
    compute_nuclear,
    compute_nuclear_repulsion,
    compute_overlap,
)
from pybest.io.iodata import IOData
from pybest.linalg.dense import DenseLinalgFactory
from pybest.scf import AufbauOccModel
from pybest.wrappers.hf 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
gobasis = get_gobasis("cc-pvdz", fn_xyz)

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

# Compute integrals
olp = compute_overlap(gobasis)
kin = compute_kinetic(gobasis)
na = compute_nuclear(gobasis)
er = compute_eri(gobasis)
external = compute_nuclear_repulsion(gobasis)

# Create alpha orbitals
orb_a = lf.create_orbital()

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

# Converge RHF
hf = RHF(lf, occ_model)
# the order of the arguments does not matter
hf_output = hf(kin, na, er, external, olp, orb_a)

# Write SCF results to a molden file
# First, add gobasis as attribute
hf_output.gobasis = gobasis
# Ready to dump results
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,
    na,
    er,
    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, na, er, external, restart.olp, restart.orb_a)