Water Hartree-Fock energies compared with PySCF

Run restricted Hartree-Fock on water in five Gaussian basis sets with vp.scf.run, the Viv port of PySCF's SCF, and compare every total energy and orbital energy with PySCF 2.5.0's own numbers for the same molecule. The energy falls toward the Hartree-Fock limit as the basis grows. At every size the total energy agrees with PySCF to a few times 10⁻¹³ hartree, and the orbital energies to about 2 × 10⁻⁹, the size the orbital-gradient tolerance allows.

The PySCF numbers come from study/ports/pyscf/oracle/, written by running PySCF 2.5.0 itself (oracle/make_oracle.sh), so this page does not need PySCF installed. The geometry is PySCF's water example, in angstrom.

import json
import numpy as np
import vipyr as vp
from vipyr.figure import figure

bases = ["sto-3g", "6-31g", "6-31g*", "cc-pvdz", "cc-pvtz"]
oracle_file = {
    "sto-3g": "water_sto3g",
    "6-31g": "water_631g",
    "6-31g*": "water1_631gs",
    "cc-pvdz": "water_ccpvdz",
    "cc-pvtz": "water_ccpvtz",
}
pyscf = {
    b: json.load(open("study/ports/pyscf/oracle/%s.json" % oracle_file[b]))
    for b in bases
}
print(pyscf["sto-3g"]["atom"])

O 0 0 0; H 0 -0.757 0.587; H 0 0.757 0.587

Total energy in five basis sets

Each call compiles the program once (cached), runs the SCF to an energy change below 10⁻¹² and an orbital gradient below 10⁻⁹, and returns the energy in hartree.

runs = {b: vp.scf.run(basis=b) for b in bases}
print(
    "%-8s %5s %22s %22s %9s"
    % ("basis", "AOs", "vipyr E (hartree)", "PySCF E (hartree)", "|diff|")
)
for b in bases:
    e, ref = runs[b].e_tot, float(pyscf[b]["e_tot"])
    print(
        "%-8s %5d %22.12f %22.12f %9.1e"
        % (b, pyscf[b]["nao"], e, ref, abs(e - ref))
    )
print("all converged:", all(r.converged for r in runs.values()))

fig = figure(
    "Water RHF energy against basis size", "basis functions", "E (hartree)"
)
nao = [pyscf[b]["nao"] for b in bases]
fig.plot(nao, [runs[b].e_tot for b in bases], "vipyr")
fig.plot(nao, [float(pyscf[b]["e_tot"]) for b in bases], "PySCF 2.5.0")
fig.style("vipyr", "s1 circle")
fig.style("PySCF 2.5.0", "s2 dash ring")

basis      AOs      vipyr E (hartree)      PySCF E (hartree)    |diff|
sto-3g       7       -74.963063129729       -74.963063129729   1.6e-13
6-31g       13       -75.983948498105       -75.983948498106   1.4e-13
6-31g*      18       -76.009099108597       -76.009099108597   2.3e-13
cc-pvdz     24       -76.026765673120       -76.026765673120   4.3e-14
cc-pvtz     58       -76.057114083120       -76.057114083119   2.3e-13
all converged: True

Orbital energies

The energy is one number; the orbital energies are all of them. The worst difference from PySCF over every occupied and virtual orbital, per basis. vipyr stops when the orbital gradient is below 10⁻⁹, and PySCF's reference run stopped below 1e-8. An orbital energy carries an error about that size. The total energy is stationary at the solution, so its error is of the order of the gradient squared, and rounding is what limits the agreement instead.

worst = []
for b in bases:
    ref = np.array([float(x) for x in pyscf[b]["mo_energy"]])
    d = np.max(np.abs(np.asarray(runs[b].mo_energy) - ref))
    worst.append(d)
    print("%-8s %3d orbitals   worst |diff| %.1e hartree" % (b, len(ref), d))

fig = figure(
    "Largest orbital-energy difference from PySCF",
    "basis functions",
    "hartree",
)
fig.plot(nao, worst, "max |vipyr - PySCF|")
fig.plot(
    nao,
    [abs(runs[b].e_tot - float(pyscf[b]["e_tot"])) for b in bases],
    "|E difference|",
)
fig.style("max |vipyr - PySCF|", "s1 circle")
fig.style("|E difference|", "s2 square")
fig.yscale("log")

sto-3g     7 orbitals   worst |diff| 1.9e-09 hartree
6-31g     13 orbitals   worst |diff| 1.1e-10 hartree
6-31g*    18 orbitals   worst |diff| 1.7e-09 hartree
cc-pvdz   24 orbitals   worst |diff| 2.3e-09 hartree
cc-pvtz   58 orbitals   worst |diff| 2.0e-09 hartree

check() on each run

check() runs the program's own tests on the result.

for b in bases:
    rep = runs[b].check()
    print("%-8s passed=%s  (%d tests)" % (b, rep.passed, len(runs[b].tests)))

sto-3g   passed=True  (3 tests)
6-31g    passed=True  (3 tests)
6-31g*   passed=True  (3 tests)
cc-pvdz  passed=True  (3 tests)
cc-pvtz  passed=True  (3 tests)

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

gallery/science/plot_water_hartree_fock.py, run as written when these pages were built. The output and the figures above are from that run.

The Viv program it calls is study/ports/pyscf/vivscf.viv: its atlas, the defs and the calls between them, or open it in the workbench.