PySCF: Python Chemistry Simulation Framework

• PySCF is a Python-based electronic structure framework featuring a modular design that leverages pure Python and optimized C kernels for efficient quantum simulations.
• It supports a wide range of methods from mean-field to post-mean-field and extends to differentiable quantum chemistry, GPU, and IPU acceleration for high-performance computing.
• The platform enables flexible workflows across molecular, periodic, and custom Hamiltonians, facilitating novel method development and practical computational insights.

PySCF, the Python-based Simulations of Chemistry Framework, is a general-purpose electronic structure platform designed from the ground up to emphasize code simplicity, both to aid new method development, as well as for flexibility in computational workflow. It supports first-principles simulations of molecules and solids, finite size systems, extended systems with periodic boundary conditions, low dimensional periodic systems, and custom Hamiltonians, using mean-field and post-mean-field methods with standard Gaussian basis functions. Over the past decade it has developed into a widely used open-source platform for electronic structure theory and quantum chemical method development, while retaining a Python/C hybrid architecture in which most of the logic is written in pure Python and computationally critical paths are implemented with heavily optimized C routines (Sun et al., 2017, Sun et al., 2020, Sun et al., 14 Mar 2026).

1. Design philosophy and software structure

PySCF’s stated design philosophy is “Simplicity > Generality > Efficiency.” Almost everything is implemented in Python, with only the most time-critical routines in hand-tuned C, notably via the Libcint Gaussian-integral library and custom C kernels for tensor contractions. The package is Python-first and modular: nearly all of PySCF is written in pure Python; sub-operations such as integral generation, Fock builds, orbital transforms, and density fitting are exposed as small functional APIs following NumPy/SciPy conventions; and functions have minimal side-effects, so calls can be made anywhere, in any order, without corrupting a monolithic program state (Sun et al., 2017, Sun et al., 2020).

The internal organization is based on loosely coupled, small modules. In the current architecture, the main packages include gto and pbc/gto for basis-set definitions, atomic structures, and Gaussian integrals via Libcint; df and pbc/df for density fitting; scf, dft, pbc/scf, and pbc/dft for Hartree–Fock and DFT in molecular and periodic settings; tdscf for time-dependent HF/DFT; mp for MP2; gw for RPA and GW; cc for coupled cluster; and fci for full configuration interaction. These modules share common data structures such as molecule or crystal objects, integrals, density matrices, and orbital coefficients, so that methods can be combined with only a few line changes while still presenting a uniform Pythonic interface (Zhang et al., 2022).

A central architectural consequence is that the Python layer orchestrates high-performance kernels rather than replacing them. NumPy is used for array operations, SciPy and MPI4Py optionally for linear algebra and parallelism, and the Python driver layer spends much less than one percent of runtime in Python overhead for large calculations. This organization underlies a recurring pattern in the PySCF literature: method development occurs at the level of readable Python code, whereas integral evaluation, tensor contractions, and related kernels remain delegated to optimized backends (Sun et al., 2020).

2. Methodological scope and mathematical foundations

PySCF supports a broad hierarchy of electronic-structure methods. In the molecular and periodic Gaussian-AO framework, the core capabilities include Hartree–Fock and generalized Kohn–Sham DFT; MP2 and MP3; coupled-cluster theory up through CCSDT and CCSDTQ in small bases, as well as the perturbative triples correction; multireference approaches including CASSCF, state-averaged MCSCF, NEVPT2, DSRG-MRPT2, selected CI, DMRG interfaces, and MC-PDFT; excited-state formalisms including TDDFT, ppRPA, ADC, GW, Bethe–Salpeter, and EOM-CCSD; AFQMC; continuum solvation; QM/MM with periodic Ewald embedding; and ab initio molecular dynamics. Since 2020, notable additions have included periodic plane-wave HF/DFT/MP2/CCSD infrastructure, local CCSD(T) gradients via automatic differentiation, a full AFQMC implementation with CISD-based trial states, new CCSDT and CCSDTQ energy routines with shared-memory parallelization, extended ADC modules, and GPU-accelerated DFT and TDDFT via GPU4PySCF (Sun et al., 2020, Sun et al., 14 Mar 2026).

At the mean-field level, the second-quantized molecular Hamiltonian is written as

$$\hat H = \sum_{pq} h_{pq}\,a_p^\dagger a_q + \tfrac12\sum_{pqrs}(pq|rs)\,a_p^\dagger a_q^\dagger a_s a_r,$$

with one-electron integrals $h_{pq}$ and two-electron integrals $(pq|rs)$. The closed-shell Hartree–Fock equations take the Roothaan–Hall form

$F\,C = S\,C\,\varepsilon,$

where $F$ is the Fock matrix, $S$ the AO overlap matrix, $C$ the MO coefficient matrix, and $\varepsilon$ the orbital energies. In Kohn–Sham DFT, the total energy functional is

$$E[\rho] = T_s[\rho] + \int V_{ne}(r)\rho(r)\,dr + \tfrac12\iint \frac{\rho(r_1)\rho(r_2)}{r_{12}}\,dr_1dr_2 + E_{xc}[\rho,\nabla\rho,\ldots],$$

with exchange–correlation terms supplied through Libxc and XCFun (Sun et al., 2017, Sun et al., 2020).

Post-mean-field formalisms follow the same pattern of exposed equations and reusable intermediates. The MP2 correlation energy is given in PySCF’s overviews in the usual occupied–virtual denominator form, while CCSD solves projected amplitude equations for singles and doubles and yields the CCSD energy from the connected similarity-transformed Hamiltonian. EOM-CCSD then diagonalizes the effective Hamiltonian in the relevant excitation manifold. The documentation and review papers emphasize that these are implemented in a way that preserves textbook transparency at the Python level while retaining optimized kernels and asynchronous I/O for large runs (Sun et al., 2017, Sun et al., 2020, Sun et al., 14 Mar 2026).

3. Molecular, periodic, and programmable workflows

A distinctive feature of PySCF is that the same framework supports molecules, solids, low-dimensional periodic systems, and custom Hamiltonians. For periodic boundary conditions, PySCF employs a Bloch-summed Gaussian basis,

$$\phi_{\mathbf k,\chi}(\mathbf r) = \sum_{\mathbf T} e^{\,i\mathbf k\cdot\mathbf T}\,\chi(\mathbf r - \mathbf T),$$

with explicit $h_{pq}$0-point sampling for SCF and correlated methods. The same core algorithms are reused across molecular and periodic calculations, differing primarily in the handling of translational symmetry, Coulomb sums, and density-fitting representations (Sun et al., 2017, Sun et al., 2020).

Density fitting is a unifying theme across these workflows. In molecular calculations, four-center integrals may be replaced by three-index tensors; in solids, auxiliary representations may be Gaussian auxiliary AOs, plane-wave functions on a grid, or mixed Gaussian plus plane-wave forms. The ten-year overview further notes three complementary density-fitting formalisms for periodic boundary conditions—FFTDF, GDF, and ISDF—as well as range-separated variants such as RSGDF and RSJK for hybrid DFT and post-HF at $h_{pq}$1 and modest $h_{pq}$2-meshes (Sun et al., 2020, Sun et al., 14 Mar 2026).

Programmability is built directly into the interface. Custom Hamiltonians can be introduced by overwriting mf.get_hcore, mf.get_ovlp, and mf._eri; advanced options such as density fitting, exact two-component relativity, and second-order SCF can be layered through decorator-style calls; and major steps in the SCF procedure are exposed as Python functions that can be overridden or decorated at runtime. The framework also supports MPI parallelism through Python daemons in MPI4PySCF, so that the Python API remains identical whether running in serial OpenMP mode or distributed-memory MPI mode. This combination of monkey-patching, plugin-style composition, and shared abstractions is one of the main reasons PySCF has been used as a development environment rather than only as an end-user code (Sun et al., 2017, Sun et al., 2020).

A common misconception is that a Python-centered interface implies a rigid or pedagogical code. The published design documents point instead to a system in which Mole or Cell objects, mean-field drivers, post-SCF solvers, and low-level helper modules remain light-weight and composable, permitting both routine workflows and highly nonstandard problem formulations without altering the core library (Sun et al., 2020, Sun et al., 14 Mar 2026).

4. Differentiable quantum chemistry with PySCFAD

The differentiable extension of PySCF, PySCFAD, introduces automatic differentiation into the framework using JAX. Its implementation strategy consists of replacing NumPy and SciPy calls by their JAX analogues, treating expensive C-kernels such as ERIs via Libcint and XC derivatives via Libxc as primitive AD operations by registering their analytic derivatives, and using implicit differentiation for iterative solvers such as SCF and coupled-cluster amplitude equations rather than unrolling all iterations. In this construction, the restricted Hartree–Fock energy functional in an AO basis is

$h_{pq}$3

and for a fixed-point condition $h_{pq}$4 the derivative with respect to an external parameter $h_{pq}$5 is expressed as

$h_{pq}$6

The paper explicitly relates this reverse-mode treatment to the conventional Z-vector approach (Zhang et al., 2022).

The reported applications span both molecules and solids. The paper demonstrates orbital-optimized RPA, response properties such as static Raman and IR intensities at CCSD, EOM-EE-CCSD excitation energies, derivative couplings at CIS or FCI, and the HF stress tensor in a Si primitive cell. For BH at CCSD/aug-cc-pVDZ, Raman frequencies, activities, and depolarization ratios agree with reference to better than $h_{pq}$7; for H$h_{pq}$8 on 6-31G$h_{pq}$9, the lowest four singlet EOM-CCSD states match conventional EOM-CCSD to within machine precision; and for the Si primitive cell with GTH-SZV pseudopotential, AD stress-tensor results agree with finite differences to better than $(pq|rs)$0 over a range of lattice parameters (Zhang et al., 2022).

The published performance characterization is deliberately conservative. Objective evaluations and gradient computations in PySCFAD plus JAX are typically $(pq|rs)$1–$(pq|rs)$2 slower than pure PySCF, mainly due to disabled ERI permutation symmetry, while overall memory use remains comparable if implicit differentiation is used because there is no need to store all SCF or CC iterations. The same body of work identifies several open issues: the JAX–NumPy dichotomy complicates control flow and dynamic constructs; ERI eight-fold permutation symmetries are not yet enforced in AD; out-of-core workflows require new data-history tracking; higher-order AD for ERI exponents and contraction coefficients beyond first derivatives is missing; and fully differentiable eigensolvers for degenerate cases remain an open challenge. A related optimization study shows that AD gradients versus hand-coded CPHF derivatives can agree to approximately $(pq|rs)$3 a.u. with developer effort of approximately 20 lines of JAX code, illustrating the intended trade-off between prototyping speed and low-level optimization (Zhang et al., 2022, Pu et al., 7 Jun 2025).

5. Hardware acceleration and performance engineering

Recent PySCF development has placed strong emphasis on modern hardware, algorithmic optimization, and Python tooling. GPU4PySCF mirrors the PySCF namespace and exposes .to_gpu() and .to_cpu() helper methods that recursively convert objects and large arrays between NumPy and CuPy forms. In GPU4PySCF v1.0, GPU-accelerated functionality includes direct and density-fitting Hartree–Fock and Kohn–Sham DFT, geometry optimization, vibrational frequency analysis, implicit solvent models, nonlocal VV10 correlation, three-center integrals for density fitting, initial GPU-powered Fock builds for MP2 and CCSD, CHELPG charges, GIAO NMR shielding constants, Fukui functions, and electrostatic potentials. Under the hood, the code re-implements key methods such as kernel(), nuc_grad_method(), and Hessian(), dispatching contractions to cuTENSOR or custom CUDA kernels and evaluating XC contributions on device with an experimental CUDA-enabled libXC (Wu et al., 2024).

The mathematical core of GPU acceleration is described in two forms. In the density-fitting Fock build,

$(pq|rs)$4

while grid-based XC evaluation is written as

$(pq|rs)$5

A lower-level account focuses on GPU Rys-quadrature ERIs for contracted basis sets up to $(pq|rs)$6 functions, integral-direct HF Fock builds, and GPU nuclear gradients via double contraction. In that implementation, the ERI kernels use shell-pair screening, batching, and eight-fold permutational symmetry, and benchmark calculations show a significant speedup of two orders of magnitude with respect to the multi-threaded CPU Hartree–Fock code of PySCF, with performance comparable to GAMESS and QUICK on a single NVIDIA A100 GPU (Li et al., 2024).

Application-level benchmark data are similarly explicit. One GPU4PySCF paper reports that, when performing DFT calculations on modern GPU platforms, GPU4PySCF delivers 30 times speedup over a 32-core CPU node, resulting in approximately 90% cost savings for most DFT tasks; a broader benchmark set reports SCF speedup of $(pq|rs)$7 on average, gradient speedups of $(pq|rs)$8–$(pq|rs)$9, analytic Hessian speedups of $F\,C = S\,C\,\varepsilon,$0–$F\,C = S\,C\,\varepsilon,$1, and DFT plus PCM speedups of $F\,C = S\,C\,\varepsilon,$2–$F\,C = S\,C\,\varepsilon,$3 for SCF, $F\,C = S\,C\,\varepsilon,$4–$F\,C = S\,C\,\varepsilon,$5 for gradients, and $F\,C = S\,C\,\varepsilon,$6–$F\,C = S\,C\,\varepsilon,$7 for Hessians (Wu et al., 2024).

Performance engineering within PySCF-based simulations is not limited to GPUs. A 2025 guide discusses improved initial guesses, second-order SCF via the newton solver, multigrid integration, density fitting approximations, Numba JIT for custom kernels, and automatic differentiation via PySCFAD. For RKS(PBE0) on an NVIDIA A100 with density fitting, the reported timings scale from $F\,C = S\,C\,\varepsilon,$8 s on CPU versus $F\,C = S\,C\,\varepsilon,$9 s on GPU for $F$0 AOs to $F$1 s on CPU versus $F$2 s on GPU for $F$3 AOs, with observed scaling $F$4 and $F$5. The same paper frames these results as a practical guide for enhancing the use of PySCF’s capabilities in quantum chemistry research (Pu et al., 7 Jun 2025).

A distinct accelerator line re-implements the SCF loop in JAX for Graphcore IPUs. PySCF$F$6 kept Libcint and Libxc but ported main Python and NumPy routines to JAX, compiled a fully fused SCF kernel through XLA, and exploited IPU MIMD parallelism for two-electron integrals. On a representative 46-AO molecule, the reported speed-up versus CPU float64 is $F$7 for the total SCF iteration and for SCF to convergence. This implementation was used to generate QM1B with one billion training examples in approximately $F$8 IPU-hours, corresponding to about five days on 320 IPUs, whereas the estimated CPU time was approximately two years for comparably sized datasets (Mathiasen et al., 2023).

6. Extensions, ecosystem, and frontier applications

PySCF’s architecture has enabled a substantial ecosystem of companion projects. The ten-year overview identifies a primary pyscf repository for the stable modules most users need, a pyscf-forge repository for experimental features, and companion repositories such as GPU4PySCF, PySCFAD, and MPI4PySCF. Earlier reviews also document native interfaces or independent external projects involving geomeTRIC and pyberny, libdftd3, Block, CheMPS2, Dice, SHCI, NECI, CPPE, OpenFermion, Qiskit-Aqua, OpenQEMIST, QMCPACK, QWalk, pyQMC, HANDE-QMC, QSoME, pDMET, PyDMFET, LSQC, and PySCF-NAO. This suggests an ecosystem in which PySCF functions as both a standalone electronic-structure code and a backend for embedding, stochastic methods, quantum computing, and machine-learning workflows (Sun et al., 2020, Sun et al., 14 Mar 2026).

One specialized extension is PyFLOSIC, an open-source Python implementation of the Fermi–Löwdin orbital self-interaction correction based on PySCF. It is designed as a lightweight layer on top of PySCF, re-using almost every foundational routine and adding only the pieces necessary to carry out FLO-SIC. The package wraps the ordinary PySCF DFT workflow so that, at each SCF iteration, it constructs Fermi orbitals from the current density, orthonormalizes them to yield FLOs, evaluates the orbital-by-orbital self-interaction correction, and adds the corresponding SIC potential to the PySCF Fock matrix. Because it piggy-backs on PySCF’s handling of basis sets, integrals, numerical grids, and Libxc or XCFun functionals, it immediately gains support for a wide range of basis and quadrature choices. The reported basis-set convergence study on SO$F$9 and SF$S$0 showed that FLO-SIC requires very large basis sets and very large unpruned grids to converge atomization energies to approximately $S$1 kcal/mol, and the paper explicitly notes that older DFO and DFO+ basis sets were barely double-$S$2 in quality (Schwalbe et al., 2019).

Another example uses PySCF in strong magnetic fields through GQCP-generated GIAO integrals. In that workflow, GQCP builds complex one-electron, overlap, and two-electron integrals that include the magnetic-field Hamiltonian and London or GIAO phase factors, and these are then monkey-patched into PySCF by overriding get_hcore, get_ovlp, and mf._eri. Because PySCF’s SCF and CC stacks already handle complex arrays, the same infrastructure can then run RHF and CCSD with reduced ERI symmetry. The reported cost increase is roughly a factor of $S$3–$S$4 in disk or memory and about $S$5–$S$6 in CPU to build and store the ERI tensor, a roughly $S$7 wall-time increase for the SCF stage, and in a 10-orbital H$S$8 example about 30 s at zero field versus about 120 s at $S$9 a.u. on a single 16-core node (Zalialiutdinov et al., 9 Mar 2025).

The PySCF core itself has also expanded into higher-order coupled-cluster territory. A 2026 implementation adds open-source RCCSDT, RCCSDT(Q), RCCSDTQ, and UCCSDT within PySCF, combining compact triangular storage of the highest-order amplitudes with the multithreaded tensor contraction backend pytblis. The shared-memory implementation delivers near-ideal thread scaling up to 90 cores, while MPI-based distributed-memory algorithms extend RCCSDT, RCCSDT(Q), and RCCSDTQ to approximately 3,000 CPU cores on 32 nodes. The reported practical reach is CCSDT(Q) calculations with approximately 100 correlated electrons in 450 orbitals and CCSDTQ calculations with approximately 50 correlated electrons in 115 orbitals (Jin et al., 1 Jul 2026).

Taken together, these developments show a consistent pattern rather than a single methodological direction. PySCF remains a framework in the strict software sense: a common set of abstractions, data structures, and interfaces on top of which differentiable programming, GPU kernels, IPU-oriented data generation, self-interaction corrections, magnetic-field Hamiltonians, and high-order coupled-cluster solvers can be attached with limited disruption to the core design. The published limitations are correspondingly specific rather than conceptual—out-of-core differentiation, fully differentiable eigensolvers in degenerate cases, broader post-HF GPU support, multi-GPU algorithms, and further integration of experimental modules remain active areas of work (Zhang et al., 2022, Wu et al., 2024, Pu et al., 7 Jun 2025, Sun et al., 14 Mar 2026).