Python bindings

import ferric gives you the same engine as the ferric CLI, as plain function calls that return result objects holding floats and numpy arrays. To install it, see Installation. The wheel is enough; you do not need a clone of the repository to run anything on this page.

If you already know PySCF, For PySCF users maps the calls you know onto ferric's and lists where the two behave differently.

Every snippet below uses an inline geometry. Run them in order, since later ones reuse water, o2, bs, bs_dz and aux from earlier ones. The snippets through "Properties and charges" were run when this page was written (2026-09-23), and any output shown is what they printed. The one-line calls in the reference tables were not run.

Molecules and basis sets

import ferric

water = ferric.Molecule.from_xyz_string("""3
water
O   0.000000   0.000000   0.117790
H   0.000000   0.755453  -0.471161
H   0.000000  -0.755453  -0.471161
""")
bs = ferric.BasisSet.bundled("sto-3g")

The string is standard XYZ: an atom count on the first line, a comment line, then one symbol x y z line per atom. The count must be the first line, so start the string with """3. Starting it with """ and a newline gives an empty first line, and parsing fails with bad atom count. Molecule.from_xyz(path) reads the same format from a file. A leading @ on a symbol (@H) makes a ghost atom, which carries basis functions but no nucleus or electrons.

Charge and spin belong to the molecule, not to the SCF call:

o2 = ferric.Molecule.from_xyz_string("""2
O2 triplet
O 0.0 0.0 0.0
O 0.0 0.0 1.208
""", charge=0, multiplicity=3)

multiplicity is 2S+1, so a triplet is 3. (PySCF's spin is 2S, so the same triplet there is spin=2.) Both constructors default to charge=0, multiplicity=1. run_uhf and run_rohf take no spin argument; they read it from the molecule. The geometry-changing drivers (run_frequencies, run_saddle, run_irc) also accept a multiplicity= keyword.

Units

print(water.coords()[0])        # (0.0, 0.0, 0.11779)            Ångström
print(water.coords_bohr()[0])   # (0.0, 0.0, 0.22259084021251865) Bohr

The Sharp bits page has the full units table, including the QM/MM accessors.

Bundled basis sets

BasisSet.bundled(name) loads a basis compiled into the library. Names are case-insensitive. An unknown name raises ValueError. These 25 sets are available (cc-pvdz-rifit is also accepted, as an alias of cc-pvdz-ri):

The list comes from the bundled() match in crates/ferric-core/src/basis.rs. The source records two coverage gaps: cc-pvtz lacks K, and cc-pvtz-rifit lacks K and Ca. An RI run on a missing element errors; it is not silently patched. The Python API has no loader for basis files on disk; the Rust library parses BSE-JSON and Gaussian-94 files (see the Rust API).

Most drivers take a BasisSet object. The drivers that move atoms (run_optimize, run_frequencies, run_saddle, run_irc, run_qmmm, run_optimize_qmmm) take the basis name as a string instead, because they rebuild the basis at every geometry.

Ground state

rhf = ferric.run_rhf(water, bs)
print(rhf)
print(f"RHF: {rhf.energy:.10f} Ha, converged={rhf.converged}")

RHF Energy: -74.9631468000 Ha (converged: true, 8 iterations)
RHF: -74.9631468000 Ha, converged=True

Open shell, using the triplet o2 built above:

bs_dz = ferric.BasisSet.bundled("cc-pvdz")
uhf  = ferric.run_uhf(o2, bs_dz)
rohf = ferric.run_rohf(o2, bs_dz)
print(f"UHF  {uhf.energy:.10f}  converged={uhf.converged}")
print(f"ROHF {rohf.energy:.10f}  converged={rohf.converged}")

UHF  -149.6276689907  converged=True
ROHF -149.6079865946  converged=True

PySCF 2.12 on the same geometry and basis gives −149.6276689907 (UHF) and −149.6079865946 (ROHF). run_rohf returns a UhfResult: it has α and β densities and orbital energies, which coincide for the spatial orbitals.

Kohn–Sham DFT:

dft = ferric.run_dft(water, bs_dz, functional="b3lyp")
print(f"B3LYP {dft.total_energy:.10f}")

run_dft is closed-shell only, and its default functional is LDA. It uses density fitting for Coulomb by default (def2-universal-jkfit); run_rhf does not. Pass df_j_aux="exact" for conventional four-centre Coulomb when you compare against an exact-Coulomb code. dispersion="d3bj" adds Grimme D3(BJ); the result then carries e_scf, e_dispersion and their sum in total_energy. with_gradient=True also returns the analytic nuclear gradient from dft.gradient(). There is no open-shell run_dft. Unrestricted Kohn–Sham is reachable from Python only inside other drivers: run_frequencies(reference="uhf", xc=...), run_u_gw(xc=...), run_qmmm(method="uks") and run_cdft(functional=...).

run_rhf also takes implicit solvent (solvent=78.4 or solvent="water", IEF-PCM), point charges and a uniform field. See its docstring (help(ferric.run_rhf)) for the full SCF knob set, which matches the CLI [scf] section.

Check .converged, and know what it means

run_rhf, run_uhf, run_rohf and run_qmmm return a result whether or not the SCF converged. A non-converged result is not an error; it is a result with converged = False, and its energy is a plausible, wrong number. run_dft, run_ksdft and the correlated drivers (MP2, CC, RPA, GW, TDDFT) behave differently: they raise if their reference SCF does not converge.

converged = True means the SCF reached a stationary point. It does not mean the lowest one. The O2 triplet above with sto-3g instead of cc-pvdz shows this. run_uhf with the default MINAO guess converges to −147.63397 Ha, which is a saddle of the orbital Hessian 1.33 mHa above the UHF minimum at −147.635296 Ha (PySCF's default guess lands on the same saddle). Without stability_descent nothing in the result flags it. run_uhf(o2, bs, stability_descent=True) (CLI: [scf] stability_descent = true on kind = "uhf") follows the downhill eigenvector and reaches the minimum. guess="hcore" converges (converged=True) to a much higher stationary point, −147.3789 Ha, which lies above ferric's own ROHF (−147.63219 Ha). A UHF energy above the ROHF energy for the same molecule cannot be a ground state, so comparing the two is a cheap check for open-shell work.

Correlation

Correlated drivers run their own reference SCF internally, so they take the molecule and basis rather than an SCF result. They also take an explicit auxiliary (RI) basis. There is no automatic choice.

aux = ferric.BasisSet.bundled("cc-pvdz-ri")

rhf = ferric.run_rhf(water, bs_dz)
mp2 = ferric.run_rimp2(water, bs_dz, aux)
cc  = ferric.run_ccsd_t(water, bs_dz, aux)

print(f"RHF      {rhf.energy:.10f}")
print(f"RI-MP2   {mp2.total_energy:.10f}  (corr {mp2.mp2_corr:.10f})")
print(f"CCSD(T)  {rhf.energy + cc.correlation_energy + cc.t_correction:.10f}  total")
print(f"         {cc.correlation_energy:.10f}  CCSD correlation")
print(f"         {cc.t_correction:.10f}  (T)")

closed-shell CCSD converged in 10 iterations. E_corr = -0.2135061893
RHF      -76.0267679974
RI-MP2   -76.2308014541  (corr -0.2040334567)
CCSD(T)  -76.2433412449  total
         -0.2135061893  CCSD correlation
         -0.0030670582  (T)

The first line is progress output that the CCSD solver prints to stdout.

CcResult holds only correlation_energy and t_correction (which is None for run_ccd and run_ccsd). It carries no reference energy, so the total above adds run_rhf(...).energy by hand. The MP2-family results all carry a total_energy, and most also carry the reference energy.

Other members of the family use the same call shape:

att   = ferric.run_attenuated_rimp2(water, bs_dz, aux, omega=0.420)  # Å⁻¹
scs   = ferric.run_scs_mp2(water, bs_dz, aux)
sos   = ferric.run_laplace_sos_mp2(water, bs_dz, aux)
oo    = ferric.run_oo_rimp2(water, bs_dz, aux)
mp3   = ferric.run_mp3(water, bs_dz, aux)

frozen_core= is accepted by every correlated driver. run_terfc_rimp2 has no CLI method.kind of its own (the CLI reaches it through att-rimp2 with att_operator = "terfc"); the reference table below marks the CLI kind of every driver. run_rimp2, run_drpa and run_linlccd compute their method exactly unless local="amplitude-threshold" and eps= are passed (see Exact and local correlation).

Response and excited states

rpa   = ferric.run_pdep_rpa(water, bs_dz, aux)            # RPA correlation
gw    = ferric.run_gw(water, bs_dz, aux)                  # G0W0@HF by default
tddft = ferric.run_tddft(water, bs_dz, aux, n_roots=3, method="tda")

print(rpa.total_energy, rpa.eigensolver_converged)
print(gw.mo_indices, gw.eps_qp)          # quasiparticle energies, Hartree
print(tddft.excitation_energies)         # Hartree

run_gw runs method="g0w0" on an HF reference unless you pass xc= (for example xc="pbe"); by default it corrects HOMO−2 through LUMO+2. Check gw.outer_converged and gw.qp_converged before using the numbers. run_u_gw is the open-shell version.

run_tddft is closed-shell only. With functional=... it includes the (ia|f_xc|jb) exchange-correlation kernel; meta-GGA, VV10 and range-separated functionals are refused because their kernel is not built. With no functional, it is CIS (method="tda") or TDHF (method="casida"). Both methods match PySCF TDA/TDDFT to 1e-3 eV on the systems listed in Capabilities and validation.

Properties and charges

The property functions take the molecule, the basis and a converged RhfResult or DftResult. They work on closed-shell results.

import numpy as np

rhf = ferric.run_rhf(water, bs)                      # water / STO-3G

q_lowdin = ferric.lowdin_charges(water, bs, rhf)
q_resp   = ferric.resp_charges(water, bs, rhf)
print(np.round(q_lowdin, 4), np.round(q_resp, 4))

# ESP 3 Bohr above each atom. Points are an (N, 3) array in Bohr.
pts = np.array(water.coords_bohr()) + np.array([0.0, 0.0, 3.0])
print(np.round(ferric.esp_at_points(water, bs, rhf, pts), 6))

[-0.2525  0.1263  0.1263] [-0.6176  0.3088  0.3088]
[-0.064684 -0.067154 -0.067154]

The charge family is mulliken_charges, lowdin_charges, hirshfeld_charges, chelpg_charges and resp_charges, all returning one charge per atom in units of e. hirshfeld_charges builds its proatoms from free-atom SCF densities in the molecule's basis, as the CLI does; proatom="slater" selects the single-exponential Slater proatom instead, which is 0.23–0.72 e away on H2O, CO and CH3OH. resp_charges is a single-stage restrained fit, not the multi-stage, multi-conformer RESP procedure. esp_at_atoms gives the potential at each nucleus. hirshfeld_polarizability returns per-atom 3×3 polarizability tensors (Bohr³) and needs an RI basis.

For NPZ export of ML-ready features (MO coefficients, PDEP eigenvectors, ESP, charges, polarizabilities, C6 coefficients) in one run, use the CLI's [rpa] export_npz section; see Input file (TOML).

What comes back

Results are Python objects with plain attributes for scalars and methods for arrays:

D = rhf.density()             # numpy.ndarray, (n_bf, n_bf), AO basis
e = rhf.orbital_energies()    # numpy.ndarray, ascending, Hartree
C = rhf.mo_coefficients()     # numpy.ndarray, (n_bf, n_mo)

Matrices and tensors come back as numpy.ndarray. Per-atom lists (charges, ESP values) come back as Python lists; wrap them in np.asarray if you need arrays. run_drpa, run_linlccd and tune_omega return plain dicts, and run_drpa_scan returns a list of them.

AO-basis matrices follow libint2's basis-function conventions, which are not PySCF's. MEASURED on CO/cc-pVDZ: total and orbital energies agree with PySCF to 3e-12 and 8e-9 Ha, and the two AO density matrices differ element by element by up to 1.5. Compare invariant quantities, not raw AO matrices.

Memory, threads and MPI

Most drivers accept memory_budget_gb (GiB). It sets the same per-allocation limits as the CLI's [memory] budget_gb: an allocation that does not fit is spilled to disk, recomputed on demand (the DFT grid AO cache) or refused with an error naming it (for example run_rimp2 and run_ccsd). Unlike the CLI, Python installs no shared ledger, so each check compares its own allocation with the whole budget rather than with what other live allocations have left. Two checks instead subtract the process's current resident memory (RSS) first and allow 90% of the remainder: the KS-DFT decision to store or recompute the grid AO cache, and the UKS Newton/TRAH f_xc kernel's second grid cache. Because RSS includes everything already resident, those two see less than the full budget. It is not a cap on total process memory; see Sharp bits.

mp2 = ferric.run_rimp2(water, bs_dz, aux, memory_budget_gb=8.0)

import ferric pins OpenBLAS to one thread unless OPENBLAS_NUM_THREADS is already set. ferric parallelizes with rayon instead. run_rhf, run_uhf, run_rohf, run_dft, run_rimp2, run_ccsd, run_ccsd_t, run_gw and run_tddft release the GIL, so independent jobs submitted from a ThreadPoolExecutor run in parallel. The other drivers hold it. For throughput across many molecules, prefer many single-threaded processes.

Do not run a Python script under mpirun. The bindings expose no rank or world-size accessor, so every rank runs the whole script. Distributed-memory runs go through the CLI built from source with MPI; see Installation.

Full reference

The module registers 61 public functions and 39 classes. That count excludes _cli_main, the entry point behind the ferric console command. It also exports two constants: DEFAULT_TEMPERATURE_K (298.15) and BOLTZMANN_HARTREE_PER_K. The list below was taken from the registration block of #[pymodule] fn ferric in crates/ferric-python/src/lib.rs, and each purpose line is condensed from that item's doc comment, or from its code where it has none. help(ferric.<name>) shows the full docstring and signature.

"CLI" gives the matching method.kind, task or TOML section, or "—" when the capability is Python-only. How well each one is validated is on Capabilities and validation.

Molecules and basis sets

SCF and DFT

Constrained DFT

See Constrained DFT for a worked example.

Geometry, vibrations and reaction paths

QM/MM

See QM/MM for a worked example.

MP2 family

Exact and local correlation

run_rimp2, run_drpa and run_linlccd name a method and compute it exactly by default. The local approximation is asked for with keywords that mirror the CLI's [local] section, under the same rules, raised as ValueError before any SCF:

• local= is None (exact; "none" is the same) or "amplitude-threshold"; any other value is an error.
• eps= is required with local="amplitude-threshold": the threshold is part of the model and has no default. eps=0 reproduces the exact method. eps= without local= is an error, not ignored.
• compute_reference=True (local only) also computes the exact method and reports it: local["e_corr_canonical_ri"] (MP2), e_corr_plasmon_canonical (dRPA), local["e_corr_exact"] (LinLCCD). Off, the key is present and None. It is a full exact calculation, so it removes any cost saving.
• run_rimp2(integral_direct=True, ...) selects the integral-direct local MP2, with its locality maps aux_radius, virt_radius (Bohr), ao_tail, schwarz_skip, batch_merge, gate_cal and virt_schwarz_kappa (the CLI names and defaults). They are errors without integral_direct=True. kappa is an error on the local MP2.

Every result says which model it is: RiMp2Result.local and the dicts' "local" key are None for the exact method, else a dict with scheme, eps, keep_fraction and integral_direct plus solver counters. All local paths are closed shell.

Coupled cluster

All use RI integrals from the auxbasis argument and a closed-shell RHF reference.

RPA, GW and excited states

Properties and charges

Each takes (mol, basis_set, result) with a converged closed-shell RhfResult or DftResult; esp_at_points also takes the points and hirshfeld_polarizability an RI basis.

The [rpa] compute_* keys write into the CLI's NPZ export.

Integrals and orbitals (prototyping)

These take omega/r0 in raw Bohr units, unlike the run_* drivers.

Conformer ensembles