`vp.scf.run`

Restricted or unrestricted Hartree-Fock for a molecule in a Gaussian basis set, with PySCF's names, shapes and units.

Part of the distribution vipyr-chem; it runs from a checkout of the repository. Planned on PyPI as pip install "vipyr[chem]".

The molecule is data: atomic numbers, positions in angstrom, a basis set, the charge and the spin. The run builds the one- and two-electron integrals, starts from the core-Hamiltonian guess, and iterates the self-consistent field with DIIS until both the energy change and the orbital gradient are below their tolerances. The defaults are water in STO-3G.

The results use PySCF's conventions: e_tot is the total energy in hartree, including the nuclear repulsion e_nuc; mo_energy, mo_coeff and dm have the shapes pyscf.scf.RHF and UHF give them.

Signature

vp.scf.run(*, atoms=[8, 1, 1], coords=..., basis="sto-3g", charge=0, spin=0, method="rhf", jk="auto", e_tol=1.0e-12, grad_tol=1.0e-9, max_iter=100, diis_size=8, diis_start=1)

Parameters

atoms : ndarray[int64], default [8, 1, 1]: Atomic numbers, one per atom. Range: each element >= 1, <= 18. 1 (H) to 18 (Ar), the elements the basis tables cover; given with coords.
coords : ndarray[float64], default the program's 3×3 array: Positions [natm, 3], angstrom. Units: angstrom. Range: shape (len(atoms), 3). One row of three per atom; given with atoms.
basis : str, default "sto-3g": The basis set. Range: one of sto-3g, 6-31g, 6-31g*, cc-pvdz, cc-pvtz. Spelled as PySCF reads a basis name: any case, -, _ and spaces ignored, * read as s.
charge : int, default 0: Net charge, in units of e. Units: units of e. Range: >= -18, <= 18.
spin : int, default 0: Alpha minus beta electrons (2S). Range: >= 0. 2S, the number of unpaired electrons: 0 for RHF; with UHF its parity must match the electron count, or the program stops with a message.
method : str, default "rhf": RHF or UHF. Range: one of rhf, uhf. Case does not matter.
jk : str, default "auto": How J and K are built. Range: one of auto, incore, direct. Case does not matter.
e_tol : float, default 1.0e-12: The SCF converges once the energy changes by less than this (hartree). Units: hartree. Range: >= 0. 0 asks for an exact repeat of the energy, which the run may not reach before max_iter.
grad_tol : float, default 1.0e-9: And the largest element of the orbital gradient (the commutator of the Fock and density matrices) is below this. Range: >= 0.
max_iter : int, default 100: The SCF iteration limit. Range: >= 1. A run that reaches it returns converged = False, and the program's test that the run converged fails in check().
diis_size : int, default 8: How many Fock matrices the DIIS extrapolation mixes. Range: >= 1.
diis_start : int, default 1: Iterations DIIS leaves out (1: the guess's, as PySCF). Range: >= 0.

atoms, coords are given together or not at all.
coords has shape (len(atoms), 3).

Returns

A frozen SCFResult. Iterating it gives the fields below in order; it also has params, tests, log and check().

e_tot : float, hartree: Total energy, hartree.
e_nuc : float, hartree: Nuclear repulsion, hartree.
converged : bool: Whether both convergence tests passed.
iterations : int: SCF iterations run.
mo_energy : ndarray[float64], hartree: Orbital energies, hartree: [n] RHF, [2, n] UHF (rank 1 when method is rhf, rank 2 when method is uhf).
mo_coeff : ndarray[float64]: Orbital coefficients, AO by MO: [n, n] RHF, [2, n, n] UHF (rank 2 when method is rhf, rank 3 when method is uhf).
dm : ndarray[float64]: The AO density: total for RHF, [2, n, n] alpha and beta for UHF (rank 2 when method is rhf, rank 3 when method is uhf).

Raises

vipyr.ParameterError: A ValueError: an argument outside its range, choices, rank or shape, or one of a group given without the rest; its .param names the argument.
TypeError: An argument of the wrong type, or an unknown keyword.
vipyr.VivBuildError: An ImportError: no C compiler, or the program's C does not build.
RuntimeError: The program stops with its own message (a combination it rejects, or no convergence where the program treats that as an error).

Notes

Hartree-Fock writes the wavefunction as one Slater determinant and makes the energy stationary over the orbitals. In a basis of $n$ atomic orbitals with overlap matrix $S$ , restricted Hartree-Fock solves the Roothaan-Hall equations

F (D) C = S C ε, D = 2 C_{occ} C_{occ}^{T},

where the Fock matrix $F = h + J (D) - \frac{1}{2} K (D)$ depends on the density $D$ it determines, so the equations are solved by iteration. Unrestricted Hartree-Fock keeps one density per spin (the Pople-Nesbet equations).

Each iteration builds $J$ and $K$ from the electron-repulsion integrals (McMurchie-Davidson, lib/hf), in core or direct as jk says, then diagonalizes $F$ in an orthonormal basis. Canonical orthogonalization drops overlap eigenvalues at or below $1 0^{- 8}$ , so nmo can be smaller than nao for a near-dependent basis. Pulay's DIIS extrapolates the Fock matrix from the last diis_size iterations, minimizing the norm of the commutator $F D S - S D F$ .

The run has converged when the energy changes by less than e_tol and the largest element of $F D S - S D F$ is below grad_tol. With $N$ basis functions the two-electron integrals take $O (N^{4})$ time per build (in core, $O (N^{4})$ memory as well).

Limits: atoms H to Ar (the basis tables) and the five basis sets listed; no point-group symmetry and no density fitting. lib/hf has a level shift and damping, but this function does not expose them.

References

Q. Sun et al., PySCF: the Python-based simulations of chemistry framework, WIREs Computational Molecular Science 8, e1340 (2018). The upstream code this program is checked against. doi:10.1002/wcms.1340
C. C. J. Roothaan, New developments in molecular orbital theory, Reviews of Modern Physics 23, 69 (1951). doi:10.1103/RevModPhys.23.69
J. A. Pople and R. K. Nesbet, Self-consistent orbitals for radicals, Journal of Chemical Physics 22, 571 (1954). doi:10.1063/1.1740120
P. Pulay, Convergence acceleration of iterative sequences. The case of SCF iteration, Chemical Physics Letters 73, 393 (1980). doi:10.1016/0009-2614(80)80396-4
L. E. McMurchie and E. R. Davidson, One- and two-electron integrals over Cartesian Gaussian functions, Journal of Computational Physics 26, 218 (1978). doi:10.1016/0021-9991(78)90092-X

Examples

import vipyr as vp
# water in STO-3G; atoms= (atomic numbers) and coords= (angstrom) for another molecule
r = vp.scf.run(basis="sto-3g", method="rhf")
print(r.e_tot, r.converged, r.mo_energy[:3])
assert r.check()

-74.96306312972935 True [-20.24196697  -1.26816105  -0.61738544]

An open-shell molecule: the OH radical (one unpaired electron) by UHF.

import numpy as np
import vipyr as vp
r = vp.scf.run(atoms=[8, 1], coords=np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 0.97]]),
               basis="sto-3g", method="uhf", spin=1)
print(round(r.e_tot, 8), r.converged, r.mo_energy.shape)

-74.36266919 True (2, 6)

Each example ran as written when these pages were generated, and the output shown is what it printed.

Checked against

Reference: PySCF 2.5.0, run live by study/ports/pyscf/test_vivscf_vs_pyscf.py (VIVSCF_BACKEND=export)
Cases: study/ports/pyscf/port.toml, the 11 cases
Energy tolerance: 10⁻¹⁰ hartree, |E - E_PySCF|
Mo energy tolerance: 10⁻⁸ hartree, each orbital energy

Credits

Checked against PySCF 2.5.0 (Apache-2.0), which runs as an external oracle; no PySCF code is copied. The Gaussian basis sets come from the Basis Set Exchange.

Source

Manifest: study/ports/pyscf/vivscf.export.toml
Viv: study/ports/pyscf/vivscf.viv
The defs it reaches and the calls between them: map of the code
Watch it run, line by line: /explore/walk/pyscf-port-water

vp.scf.run