Molecular Calculations with PySCF

The PySCF [1] interface decorates an existing restricted or unrestricted Kohn–Sham object. Molecular geometry, charge, spin, basis, grids, density fitting, occupations, and SCF controls remain PySCF concepts; CiderPress replaces the numerical XC evaluation and adds any model-specific energy terms.

NOTE: the decorated object sets ks.xc to a placeholder semilocal name ("PBE" for GGA-level models, "R2SCAN" for meta-GGA) so that PySCF routes the calculation through its semilocal machinery instead of treating it as a hybrid. The functional actually evaluated is the CIDER model, so ks.xc should not be read as the functional in use.

Choose the functional composition first

CIDER26XC models contain full exchange and correlation. Their PySCF initializer defaults are correct:

mf = make_cider_calc(dft.RKS(mol), "CIDER26XCCHEM")

CIDER23X and CIDER24X contain exchange. Specify the surrogate-hybrid composition explicitly:

mf = make_cider_calc(
    dft.RKS(mol),
    "CIDER23X_NL_MGGA_DTR",
    xmix=0.25,
    xkernel="GGA_X_PBE",
    ckernel="GGA_C_PBE",
)

The compositions associated with each packaged model are listed in Choosing a CIDER Functional.

Closed-shell workflow

A complete closed-shell example is:

 1#!/usr/bin/env python
 2"""Closed-shell molecular CIDER26XC calculation in PySCF."""
 3
 4import argparse
 5
 6from pyscf import dft, gto
 7
 8from ciderpress.pyscf.dft import make_cider_calc
 9
10MODELS = ("CIDER26XCCHEM", "CIDER26XCCHEMD4", "CIDER26XCSURFSCI")
11
12
13def main():
14    parser = argparse.ArgumentParser()
15    parser.add_argument("--model", choices=MODELS, default="CIDER26XCCHEM")
16    args = parser.parse_args()
17
18    mol = gto.M(
19        atom="""
20        O  0.000000  0.000000  0.117790
21        H  0.000000  0.755453 -0.471161
22        H  0.000000 -0.755453 -0.471161
23        """,
24        basis="def2-tzvp",
25        charge=0,
26        spin=0,
27    )
28
29    base = dft.RKS(mol)
30    base.grids.level = 3
31    mf = make_cider_calc(base, args.model)
32    mf = mf.density_fit(auxbasis="def2-universal-jfit")
33    mf.conv_tol = 1e-9
34    mf.max_cycle = 200
35    energy = mf.kernel()
36
37    if not mf.converged:
38        raise RuntimeError("CIDER SCF did not converge")
39
40    print(f"model = {args.model}")
41    print(f"total energy = {energy:.12f} Ha")
42    if hasattr(mf, "e_vdw_expected"):
43        print(f"SCF/base energy = {mf.e_tot_base:.12f} Ha")
44        print(f"expected dispersion = {mf.e_vdw_expected:.12f} Ha")
45        print(f"dispersion adjustment = {mf.e_vdw_delta:.12f} Ha")
46
47
48if __name__ == "__main__":
49    main()

The example uses the def2-tzvp basis and PySCF grid level 3. For a small end-to-end energetics workflow, examples/pyscf/compute_ae.py computes a molecular atomization energy with a chosen libxc or packaged CIDER functional.

Density fitting

Calling density_fit after make_cider_calc accelerates the Coulomb problem. CIDER26XC has a zero exact-exchange fraction. The auxiliary basis is supplied through the normal PySCF interface:

mf = make_cider_calc(dft.RKS(mol), "CIDER26XCCHEM")
mf = mf.density_fit(auxbasis="def2-universal-jfit")

Density fitting approximates the Coulomb contribution. The selected CIDER model and its D4 behavior remain unchanged.

Open-shell systems

Use UKS and set mol.spin to \(N_\alpha-N_\beta\). spin_square() returns \(\langle S^2\rangle\) and the corresponding multiplicity for the converged unrestricted solution.

For a difficult system, first converge a conventional functional with the same molecule, basis, grid, charge, and spin, then provide its density to the CIDER calculation. The following script uses a small O2 calculation to keep the restart example quick to run:

  1#!/usr/bin/env python
  2"""Demonstrate restart-ladder mechanics on a compact open-shell molecule.
  3
  4O2 keeps the example inexpensive.  The ladder is a template for a larger
  5open-shell calculation that needs a PBE warm start or conservative SCF
  6controls.
  7"""
  8
  9from pathlib import Path
 10
 11from pyscf import dft, gto
 12from pyscf.scf import diis as pyscf_diis
 13
 14from ciderpress.pyscf.dft import make_cider_calc
 15
 16LADDER = (
 17    {
 18        "label": "cdiis8",
 19        "DIIS": pyscf_diis.CDIIS,
 20        "diis_space": 8,
 21        "conv_tol": 1e-9,
 22        "level_shift": 0.0,
 23        "damp": 0.0,
 24    },
 25    {
 26        "label": "cdiis12",
 27        "DIIS": pyscf_diis.CDIIS,
 28        "diis_space": 12,
 29        "conv_tol": 1e-8,
 30        "level_shift": 0.0,
 31        "damp": 0.0,
 32    },
 33    {
 34        "label": "adiis",
 35        "DIIS": pyscf_diis.ADIIS,
 36        "diis_space": 12,
 37        "conv_tol": 1e-7,
 38        "level_shift": 0.0,
 39        "damp": 0.0,
 40    },
 41    {
 42        "label": "ediis",
 43        "DIIS": pyscf_diis.EDIIS,
 44        "diis_space": 12,
 45        "conv_tol": 1e-7,
 46        "level_shift": 0.0,
 47        "damp": 0.0,
 48    },
 49    {
 50        "label": "cdiis_shift02",
 51        "DIIS": pyscf_diis.CDIIS,
 52        "diis_space": 8,
 53        "conv_tol": 1e-7,
 54        "level_shift": 0.2,
 55        "damp": 0.0,
 56    },
 57    {
 58        "label": "cdiis_shift05",
 59        "DIIS": pyscf_diis.CDIIS,
 60        "diis_space": 8,
 61        "conv_tol": 1e-6,
 62        "level_shift": 0.5,
 63        "damp": 0.3,
 64    },
 65)
 66
 67
 68def configure(mf, rung):
 69    mf.DIIS = rung["DIIS"]
 70    mf.diis_space = rung["diis_space"]
 71    mf.conv_tol = rung["conv_tol"]
 72    mf.level_shift = rung["level_shift"]
 73    mf.damp = rung["damp"]
 74    mf.max_cycle = 500
 75
 76
 77def run_ladder(mf, warm_density):
 78    guesses = (("warm", warm_density), ("atom", mf.get_init_guess(key="atom")))
 79    for guess_name, density in guesses:
 80        for rung in LADDER:
 81            configure(mf, rung)
 82            mf.chkfile = str(Path(f"o2_cider_{rung['label']}_{guess_name}.chk"))
 83            energy = mf.kernel(dm0=density)
 84            if not mf.converged:
 85                continue
 86
 87            # A relaxed rung supplies a preconditioned density.  This example
 88            # follows it with the standard CDIIS control set.
 89            if rung["label"] != "cdiis8":
 90                tight_density = mf.make_rdm1()
 91                configure(mf, LADDER[0])
 92                mf.chkfile = "o2_cider_final.chk"
 93                energy = mf.kernel(dm0=tight_density)
 94                if not mf.converged:
 95                    continue
 96            return energy, rung["label"], guess_name
 97    raise RuntimeError("CIDER SCF did not converge with the suggested ladder")
 98
 99
100def main():
101    mol = gto.M(
102        atom="O 0 0 0; O 0 0 1.21",
103        basis="def2-svp",
104        charge=0,
105        spin=2,
106    )
107
108    baseline = dft.UKS(mol)
109    baseline.xc = "PBE"
110    baseline.grids.level = 3
111    baseline.conv_tol = 1e-9
112    baseline.max_cycle = 200
113    baseline.chkfile = "o2_pbe.chk"
114    baseline.kernel()
115    if not baseline.converged:
116        raise RuntimeError("Baseline PBE calculation did not converge")
117    warm_density = baseline.make_rdm1()
118
119    base = dft.UKS(mol)
120    base.grids.level = 3
121    mf = make_cider_calc(base, "CIDER26XCCHEM")
122    energy, rung, guess = run_ladder(mf, warm_density)
123
124    print(f"total energy = {energy:.12f} Ha")
125    print(f"initial density = {guess}")
126    print(f"successful preconditioning rung = {rung}")
127    print(f"<S^2>, multiplicity = {mf.spin_square()}")
128
129
130if __name__ == "__main__":
131    main()

The fallback ladder records each control set separately. A relaxed or level-shifted rung supplies a new starting density, followed in this example by a calculation with the listed standard CDIIS controls.

Checkpoints and interrupted calculations

Set a different chkfile for the baseline, each fallback rung, and the final CIDER calculation. A checkpoint can recover molecular orbitals and a density after interruption:

from pyscf.scf import chkfile

mol_from_chk, scf_data = chkfile.load_scf("cider.chk")
mo = scf_data["mo_coeff"]
occ = scf_data["mo_occ"]
dm0 = mf.make_rdm1(mo, occ)
energy = mf.kernel(dm0=dm0)

Recreate the selected CIDER model, then load the orbitals and occupations from the PySCF checkpoint. load_scf also returns the checkpoint’s Mole object; the stored orbital dimensions must match the new mean-field object.

D4-corrected energy

CIDER26XCCHEMD4 evaluates its electronic full-XC contribution during the SCF and D4 from the geometry afterward, so the dispersion term changes the total energy but not the density or the CIDER potential.

Because a user may already have attached a dispersion wrapper, CiderPress reconciles what is present with what the model expects rather than simply adding a term. After kernel(), every CIDER26XC calculation exposes:

mf.e_tot_base

Energy returned by the underlying SCF object, before the dispersion adjustment.

mf.e_vdw_present

Dispersion already contained in that energy, if a supported with_dftd3/with_dftd4 wrapper was attached.

mf.e_vdw_expected

Dispersion the selected model was trained with.

mf.e_vdw_delta

The correction e_vdw_expected - e_vdw_present.

The returned mf.e_tot is mf.e_tot_base + mf.e_vdw_delta, which contains the expected term exactly once. For CIDER26XCCHEM and CIDER26XCSURFSCI the expected dispersion is zero, so an attached wrapper is disabled and does not contribute to the returned energy.

For CIDER26XCCHEMD4, mf.nuc_grad_method() adds the analytical D4 derivative to the electronic CIDER and nuclear-repulsion terms. See Energies and Derivative Properties for the gradient interface.

CIDER24X

Install ciderpress[cider24] and use the same exchange-only composition as CIDER23X. These SDMX models use the density matrix and PyTorch-backed mapped evaluator. Keep the model on the device selected by its evaluator and avoid mixing CUDA and CPU PyTorch installations within one environment. A minimal SDMX calculation is shown in examples/pyscf/simple_sdmx.py.

Gradients and other methods

Analytical restricted and unrestricted nuclear gradients are available for the molecular NLDF path, including density-fitted calculations and the D4 correction selected by CIDER26XCCHEMD4. Use mf.nuc_grad_method() as shown in Energies and Derivative Properties. The gradient code raises NotImplementedError for models carrying SDMX features, so CIDER24X calculations give energies but not forces.

The CIDER mean-field object provides the energies and derivatives listed in Energies and Derivative Properties. Hessians, NMR, polarizability, coupled-cluster, multireference, and similar response or post-SCF methods are outside that interface. The documented analytical-gradient implementation covers the molecular NLDF path.

The PBE-seeded CDIIS, ADIIS, EDIIS, level-shift, and damping settings used by the restart example are listed in Handling SCF Convergence Issues.