PyFock

About The Project

PyFock is a pure Python quantum chemistry package that enables efficient Kohn-Sham density functional theory (DFT) calculations for molecular systems. Unlike traditional quantum chemistry codes written in Fortran or C/C++, PyFock is written entirely in Python—including all performance-critical molecular integral evaluations—while achieving computational efficiency comparable to established codes like PySCF and Psi4.

What Makes PyFock Different?

• 100% Pure Python: All code, including computationally intensive molecular integrals, is written in Python
• High Performance: Achieves efficiency comparable to C/C++ backends through Numba JIT compilation, NumPy, NumExpr, SciPy, and CuPy
• GPU Acceleration: Leverages CUDA via Numba and CuPy for up to 14× speedup on large systems
• Easy Installation: Simple pip install on all major operating systems (Linux, macOS, Windows)
• Accessible: Designed for education, prototyping, and community development
• Near-Quadratic Scaling: ~O(N²·⁰⁵) scaling through density fitting with Cauchy-Schwarz screening
• Gaussian-Type Orbitals: Employs GTOs as basis functions for molecular calculations
• Efficient Parallelization: Multi-core CPU support and multi-GPU acceleration capabilities

Performance Highlights

• Numerical accuracy: Consistent with PySCF (< 10⁻⁷ Ha)
• Parallel efficiency: Comparable to state-of-the-art C++ backends on multicore CPUs
• GPU speedup: Up to 14× faster than 4-core CPU execution for large systems
• Scaling: Near-quadratic ~O(N²·⁰⁵) for electron repulsion integrals (Coulomb term)
• XC evaluation: Sub-quadratic scaling ~O(N¹·²⁵⁻¹·⁵) for exchange-correlation contributions

Key Features

• ✅ Pure Python Implementation: Including molecular integral evaluations (overlap, kinetic, nuclear attraction, electron repulsion integrals)
• ✅ Density Fitting: Efficient density fitting approximation with Cauchy-Schwarz screening
• ✅ GPU Acceleration: Full GPU support for integral evaluation, XC term, and matrix operations. use_gpu=True runs the XC term in single precision until the relative energy change drops below 5e-7 and in double precision from there on (dynamic_precision, on by default since it is ~1.5x faster for the same converged energy; automatically disabled for meta-GGAs, whose tau single precision resolves too poorly)
• ✅ Multiple Integration Schemes:
  • Classical Taketa-Huzinaga-O-ohata scheme
  • Rys quadrature method (roots 1–10) for efficient ERI evaluation
  • Obara-Saika method for ERI evaluation
• ✅ XC Functionals: Support for LDA, GGA and meta-GGA functionals natively and optionally via LibXC integration, plus the Skala neural functional
• ✅ DIIS Convergence: Direct inversion of iterative subspace for SCF acceleration
• ✅ Parallel Execution: Multi-core CPU and multi-GPU support via Numba and Joblib
• ✅ Modular Design: Standalone integral modules for benchmarking and embedding
• ✅ Web-based GUI: Interactive interface for visualization and input generation
• ✅ Cartesian and Spherical Basis: Support for both CAO and SAO representations
• ✅ Effective Core Potentials: Support for evaluation of ECP integrals
• ✅ Analytical gradients & forces: Fast analytical nuclear gradients for density-fitted DFT — one-electron (overlap/kinetic/nuclear), DF Coulomb (3c2e + 2c2e), and XC for LDA, GGA and meta-GGA (native or LibXC) — matching PySCF forces and faster
• ✅ ASE Calculator: Optional ASE interface (geometry optimization and the wider ASE ecosystem), using analytical forces by default
• ✅ XC integration grids: Treutler-Ahlrichs radial and Lebedev angular grids with region-wise angular pruning and Becke partitioning (levels 0-9), built by a parallel Numba kernel or on the GPU (~6x faster, same grid); energies agree with PySCF to the SCF convergence threshold in our benchmarks; numgrid grids remain available as an alternative scheme
• ✅ SANO initial guess: superposition of atomic natural-orbital densities (ANO-RCC-MB natural orbitals projected onto the calculation basis, no atomic SCF needed) as the default SCF starting point; it roughly halves the number of SCF iterations relative to the core-Hamiltonian guess and prunes the XC grid accurately
• ✅ Cross-Platform: Works on Linux, macOS, and Windows

Installation

Basic Installation

PyFock can be easily installed via pip:

pip install pyfock

Installing from GitHub (Latest Development Version)

To get the latest development version directly from GitHub:

pip install git+https://github.com/manassharma07/pyfock.git

Or clone the repository and install locally:

git clone https://github.com/manassharma07/pyfock.git
cd pyfock
pip install -e .

Installing LibXC (Optional Dependency)

PyFock can use LibXC for exchange-correlation functionals not available natively in PyFock. The installation method depends on your system:

Using Conda (Recommended - Easiest Method)

conda install -c conda-forge pylibxc -y

On Ubuntu/Debian

sudo apt-get install libxc-dev
pip install pylibxc2

On macOS

brew install libxc
pip install pylibxc2

Note: The conda method is recommended as it works reliably across all platforms.

Optional Dependencies

None of these are needed to import or run PyFock. Each is available as a pip extra:

Extras combine as usual:

pip install "pyfock[ase,skala]"

and work the same when installing from a clone:

pip install -e ".[ase,skala]"

GPU acceleration is kept separate, because the right wheel depends on your CUDA version:

pip install cupy-cuda11x  # Replace 11x with your CUDA version

Two notes:

• ASE is only required when you actually use PyFockCalculator; PyFock imports and runs without it.
• pyfock[skala] does not install the skala package. PyFock reads the published TorchScript checkpoint directly with torch.jit.load, which avoids that package's own dependencies on PySCF (no Windows wheels) and e3nn. dftd3 is included because Skala is parametrised together with a DFT-D3 correction, and it is what reproduces Skala's published numbers. See Skala: the neural exchange-correlation functional.

Quick Start

Here's a minimal example to get you started with PyFock:

from pyfock import Basis, Mol, DFT

# Define molecule from XYZ file
mol = Mol(coordfile='h2o.xyz')

# Set up basis sets
basis = Basis(mol, {'all': Basis.load(mol=mol, basis_name='def2-SVP')})
auxbasis = Basis(mol, {'all': Basis.load(mol=mol, basis_name='def2-universal-jfit')})

# Create DFT calculation object
dftObj = DFT(mol, basis, auxbasis, xc='PBE')

# Set calculation parameters
dftObj.conv_crit = 1e-7
dftObj.max_itr = 20
dftObj.ncores = 4

# Run SCF calculation
energy, dmat = dftObj.scf()
print(f"Total Energy: {energy} Ha")

xc also accepts LibXC IDs instead of a name — xc=[101, 130] is the same PBE as xc='PBE'.

Usage

Computing Molecular Integrals

PyFock provides standalone access to all molecular integrals:

from pyfock import Integrals, Basis, Mol

mol = Mol(coordfile='h2o.xyz')
basis = Basis(mol, {'all': Basis.load(mol=mol, basis_name='def2-SVP')})

# One-electron integrals
S_ovlp = Integrals.overlap_mat_symm(basis)
V_kin = Integrals.kin_mat_symm(basis)
V_nuc = Integrals.nuc_mat_symm(basis, mol)

# Two-electron integrals (classical scheme)
ERI_slow = Integrals.conv_4c2e_symm(basis)

# Two-electron integrals (Rys quadrature - faster)
ERI_fast = Integrals.rys_4c2e_symm(basis)

# Three-center integrals for density fitting
ERI_3c2e = Integrals.rys_3c2e_symm(basis, auxbasis)

# Two-center integrals
ERI_2c2e = Integrals.rys_2c2e_symm(basis)

GPU-Accelerated Integrals

# GPU versions (returns CuPy arrays in device memory)
S_ovlp_gpu = Integrals.overlap_mat_symm_cupy(basis)
V_kin_gpu = Integrals.kin_mat_symm_cupy(basis)
V_nuc_gpu = Integrals.nuc_mat_symm_cupy(basis, mol)
ERI_3c2e_gpu = Integrals.rys_3c2e_symm_cupy(basis, auxbasis)

Converting Between Cartesian and Spherical Basis

# Convert from Cartesian to Spherical atomic orbitals
V_kin_CAO = Integrals.kin_mat_symm(basis)
c2sph_mat = basis.cart2sph_basis()
V_kin_SAO = np.dot(c2sph_mat, np.dot(V_kin_CAO, c2sph_mat.T))

Subset Evaluation

# Evaluate integrals for a subset of basis functions
S_ovlp_subset = Integrals.overlap_mat_symm(basis, slice=[0, 5, 0, 5])
# slice = [row_start, row_end, col_start, col_end]

Skala: the neural exchange-correlation functional

Skala is a machine-learned exchange-correlation functional from Microsoft Research AI for Science that reaches hybrid-like accuracy at semilocal cost. Pass its name as the xc argument and nothing else changes:

from pyfock import Basis, DFT, Grids, Mol

mol      = Mol(coordfile='H2O.xyz')
basis    = Basis(mol, {'all': Basis.load(mol=mol, basis_name='def2-SVP')})
auxbasis = Basis(mol, {'all': Basis.load(mol=mol, basis_name='def2-universal-jfit')})

dftObj = DFT(mol, basis, auxbasis, xc='skala-1.1', grids=Grids(mol, level=3))
dftObj.conv_crit = 1e-8
energy, dmat = dftObj.scf()

Available names are skala-1.1 (recommended), skala-1.1-rev1, skala-1.1-rev0 and skala-1.0. See examples/ex44_Skala_neural_functional.py for a runnable version, with and without the dispersion correction.

Getting the model

pip install "pyfock[skala]" is all that is needed (see Optional Dependencies) — notably not the skala package itself, which is why this works on Windows despite Skala's own packaging being Linux/macOS only.

The 2.4 MB checkpoint is downloaded from Hugging Face on first use and cached afterwards. Its SHA-256 is verified against the digests published by Microsoft before it is loaded, because TorchScript deserialization executes code from the file. To run without network access — or without huggingface_hub at all — download skala-1.1-rev1.fun from huggingface.co/microsoft/skala-1.1 and set:

export SKALA_LOCAL_MODEL_PATH=/path/to/skala-1.1-rev1.fun

Note that this path bypasses hash verification, so only point it at a file you trust.

What is different about it

Every other functional in PyFock is pointwise: the energy density at a grid point depends only on the density at that point. Skala is not — its non-local layers aggregate over the points of each atomic grid, so it consumes the whole grid at once and returns the total XC energy as a single number, and the potential comes from automatic differentiation of that number. PyFock therefore evaluates it in three passes (density over the whole grid → one model call → potential), rather than in the single fused blocked loop used for LDA/GGA/meta-GGA. The derivatives the model returns are exactly the vrho/vsigma/vtau intermediates the meta-GGA path already contracts with the AO values, so the potential assembly itself is unchanged.

Consequences worth knowing:

• Grids must be the native 'treutler' ones (the default). Skala needs both the Becke-partitioned weights and the raw single-atom weights.
• Do not density-prune the grid. The model's non-local features are integrals over each atom's full grid.
• The GPU runs the model, the CPU runs the rest. Pass skala_gpu=True and the neural functional is evaluated on the device while the integrals, the Coulomb term and the assembly of Vxc stay on the CPU — see Skala on the GPU. PyFock's full GPU SCF path (use_gpu=True) does not support Skala and raises a clear error; it would need a CuPy↔Torch bridge that is not wired up yet.
• Closed-shell (restricted) only, like the rest of PyFock's DFT.
• Analytical nuclear gradients are supported, so Skala works for geometry optimization like any other functional — see Analytical Forces & Geometry Optimization. Unlike the semilocal path, the Skala gradient includes the full grid response (Becke weight derivatives and the grid-translation term); its features are integrals over each atomic grid, so a frozen grid would cost ~1e-2 Ha/Bohr rather than the ~1e-4 it costs a meta-GGA.
• Dispersion is off by default — pass dispersion=True to include it, see below.
• The model carries about 1e-9 Ha of its own numerical noise. Presenting the network with differently shaped batches (a different max_points_per_chunk) shifts the energy at that level. Within one calculation the chunking is fixed, so the shift is systematic rather than random and SCF convergence to 1e-8 is unaffected — but do not expect two runs with different chunk sizes to agree bit for bit.

Dispersion

Skala 1.1 is parametrised together with a DFT-D3 correction — the checkpoint declares it, and dftObj.skala.d3_settings() returns 'b3lyp5', meaning D3(BJ) damping with B3LYP5 parameters and no three-body term. Leaving it out is fine for comparisons against Skala's own reference energies (which exclude it) but wrong for anything where dispersion matters, such as non-covalent interactions or conformer ranking.

D3 depends only on the atomic numbers, the coordinates and those damping parameters — not on the density. It never enters the Kohn-Sham matrix and cannot change the SCF, so PyFock evaluates it once after convergence and adds it to the total energy:

pip install dftd3

dftObj = DFT(mol, basis, auxbasis, xc='skala-1.1', dispersion=True)
energy, dmat = dftObj.scf()      # dispersion-corrected total energy
print(dftObj.Edisp)              # the correction on its own, in Hartree

dispersion=True uses whatever the functional declares. For any other functional, name the parametrisation yourself — dispersion='pbe', dispersion='b3lyp' — and tune the damping through dftObj.dispersion_version (default 'd3bj') and dftObj.dispersion_atm (default False). pyfock.Dispersion also exposes d3_energy and d3_energy_and_gradient directly, the latter giving dE_disp/dR for forces. See examples/ex43_D3_dispersion_correction.py.

The ASE calculator uses the same backend by default:

PyFockCalculator(functional='PBE', dispersion=True, dispersion_kwargs={'xc': 'pbe'})

torch-dftd (pip install pyfock[dispersion-gpu]) remains available there for GPU runs, where it evaluates the correction on the device. Select it with backend='torch-dftd', or implicitly by asking for a non-CPU device:

dispersion_kwargs={'xc': 'pbe', 'backend': 'torch-dftd', 'device': 'cuda'}

Note that Skala's own skala.dispersion module routes D3 through PySCF (pyscf.dispersion.dftd3 or dftd3.pyscf), which is why PyFock does not use it. PyFock calls dftd3.interface from the same simple-dftd3 package instead — identical numbers, no PySCF, and Windows wheels are available.

Validating against the published reference energies

Skala ships the total energies behind its benchmark report as benchmark/reference/measurements.json (a Git-LFS file, ~7.8 MB): 1221 converged runs over skala-1.1, r2scan, m06-2x and b3lyp5, each at def2-SVP/TZVP/QZVP, for ~30 molecules from GMTKN55 and a conformer benchmark. The geometries are pulled from grimme-lab/GMTKN55 at a pinned commit.

One detail is easy to get wrong: the Skala energies there include the D3 correction, because the runner builds the method through SkalaKS, whose constructor defaults to with_dftd3=True. The plain PySCF functionals in the same set do not. Comparing a bare PyFock Skala energy against them leaves a residual exactly equal to the dispersion energy.

benchmarks_tests/validate_skala_reference.py downloads both, runs PyFock under the matching protocol (spherical orbitals via dftObj.sao = True, density fitting with def2-universal-jkfit, grid level 3, conv_crit = 5e-6, dispersion=True for Skala only) and tabulates the differences:

python validate_skala_reference.py --molecules H2O --basis def2-tzvp

It runs r2SCAN alongside Skala as a control: PyFock's r2SCAN is independently validated, so whatever it shows against the same reference measures the grid, basis and density-fitting baseline rather than anything to do with Skala. Over H2, H2O, H2O2 and H3N at def2-SVP and def2-TZVP:

PyFock reproduces the published Skala energies to within the reference data's own run-to-run spread. tests/test_skala.py keeps one of these values (H2 at def2-SVP) as a permanent regression check.

Cost

Skala is several times more expensive per SCF iteration than a semilocal functional. Measured with benchmarks_tests/benchmark_skala.py (def2-SVP, level-3 grid, 4 CPU cores, density fitting):

One XC evaluation on an identical density and grid, which isolates the functional from any difference in iteration count:

Two things are worth reading out of this. First, the three-pass restructuring that non-locality forces on PyFock is not where the time goes: the extra density and potential passes cost 0.07 s of H2O's 0.79 s, and the neural network is essentially the whole overhead. Second, the relative cost falls as the system grows (21.7x to 10.8x against PBE), because the model's cost tracks the number of grid points while the semilocal functionals also carry AO work that grows with the basis.

Skala on the GPU

Skala is a PyTorch model, so it can be evaluated on the GPU without PyFock's own CuPy machinery. Pass skala_gpu=True and the model runs on the device while everything else — the integrals, the Coulomb term and the assembly of Vxc from the model's cotangents — stays on the CPU:

dftObj = DFT(mol, basis, auxbasis, xc='skala-1.1', dispersion=True, skala_gpu=True)

It works the same way through the ASE calculator, and forces come along with it:

atoms.calc = PyFockCalculator(functional='skala-1.1', skala_gpu=True, dispersion=True)

All this needs is a CUDA build of PyTorch (pip install torch --index-url https://download.pytorch.org/whl/cu128 for current cards); PyFock downloads the CUDA checkpoint, which is a separate file from the CPU one. use_gpu=True, the flag for PyFock's full GPU SCF path, is a different thing and still raises a clear error for Skala.

benchmarks_tests/benchmark_skala_gpu.py measures the speedup and the CPU/GPU agreement, for single points and for geometry optimizations.

Initial Guess for the SCF

The SCF starts by default from the SANO guess (superposition of atomic natural-orbital densities): the spherically averaged density of every free atom, written in the ANO-RCC-MB minimal basis whose contracted functions are the atomic natural orbitals of Roos and co-workers, projected onto the calculation basis and renormalized shell by shell. It is the PyFock analogue of PySCF's minao guess, costs a fraction of a second and needs no atomic SCF. The core-Hamiltonian guess remains available:

dft_obj = DFT(mol, basis, auxbasis, xc='PBE', dmat_guess_method='sano')   # default
dft_obj = DFT(mol, basis, auxbasis, xc='PBE', dmat_guess_method='core')   # core-Hamiltonian guess

# The guess density itself (CAO representation) and its references
from pyfock import Guess
dmat_guess, info = Guess.sano_dmat(mol, basis)
print(Guess.sano_citation_text(info['nuclear_charges']))

The SCF output lists the references to cite for the guess. Compared with the core guess it saves 12-63 % of the SCF iterations, 41 % over a 22-case benchmark (e.g. decane/def2-SVP/PBE 18 -> 8, cholesterol 27 -> 10), and, because PyFock prunes the XC grid with the starting density, it also removes a grid-pruning error of the too-compact core-Hamiltonian density (3e-4 Ha for decane, more for heavy atoms). See docs/sano_guess.md for the method, the benchmark table and the tests.

XC Integration Grids

The exchange-correlation term is integrated on Treutler-Ahlrichs radial grids combined with Lebedev angular grids (tables from numgrid), pruned region by region and joined by Becke partitioning with Treutler's atomic-size adjustment. gridsLevel runs from 0 (coarsest) to 9 (finest), default 3. In our benchmarks the resulting SCF energies agree with PySCF's at the same level to the convergence threshold, and the grid is built in about the time PySCF needs (decane: 0.11 s vs 0.14 s).

dftObj = DFT(mol, basis, auxbasis, xc='PBE', gridsLevel=3)            # default: grids_scheme='treutler'
dftObj.grids_options = {'pruning': None, 'size_adjustment': 'becke'} # variants of the default scheme
dftObj.grids_options = {'points_per_element': {'C': (75, 302)}}      # (n_rad, n_ang) overrides per element
dftObj.grids_scheme = 'numgrid'; dftObj.grids_preset = 'compact'     # numgrid grids ('dense': the old PyFock grids)
dftObj.use_pyscf_grids = True                                        # or let PySCF build the grid (PySCF must be installed)

A GPU calculation (use_gpu=True) also builds the grid on the GPU: the points, their atom indices and their box order are identical to the CPU build and the weights agree to ~1e-13, at about 6x the speed from roughly 30 atoms on (decane 0.13 s -> 0.02 s, taxol 3.7 s -> 0.68 s). It falls back to the CPU with a warning when CuPy or a CUDA device is missing; grids_options={'use_gpu': False} turns it off. The density pruning that follows the build runs on the GPU too, with the same AO kernel the XC term uses: it keeps exactly the same points as the CPU pruning at about 4.5x the speed (cholesterol 2.5 s -> 0.55 s).

Standalone grids come from Grids(mol, level=3) (coords, weights and the atom index of every point, grouped into 1.2 Bohr boxes; use_gpu=True builds them on the GPU). See docs/xc_grids.md for the construction, its validation and a study of what the numgrid grids can and cannot do.

Density-Fitting Coulomb Algorithms and Memory Budget

The Coulomb term is evaluated with density fitting and Schwarz screening. Three algorithms are available through DFT.DF_algo:

• 11 (default): shell-blocked Rys evaluation with block-sparse storage (CPU/GPU). It honours a memory budget for the stored integrals:

• 12: multipole-accelerated density fitting (CPU, pure functionals). The near-field three-center integrals are evaluated and stored as in 11; well-separated pairs of a basis-function-product distribution and an auxiliary function are handled through multipole expansions (exact finite moments of the primitive pair products and of the auxiliary functions, box-level expansions truncated at lmax). Energies agree with 11 to ~1e-8 Hartree. Both the speed-up and the memory saving grow with the extent of the molecule, since the far field is what they come from (def2-SVP, 4 cores, three-center build plus all per-iteration Coulomb work):

  System	far field	time vs 11	stored integrals + moments
  Caffeine (14 bohr across)	4.0%	0.9x	0.18 -> 0.19 GB
  Cholesterol (35 bohr)	32.4%	1.5x	1.61 -> 1.28 GB
  Icosane C20H42 (48 bohr)	52.6%	1.4x	0.58 -> 0.37 GB
  Tetracontane C40H82 (95 bohr)	74.7%	2.7x	2.41 -> 0.80 GB

  The crossover is around 25 to 30 bohr of molecular extent; below it a group of distributions is mostly not expanded at all, because it is only expanded when it has enough far-field auxiliary functions to pay for the expansion. {'low_memory': True} drops the pre-translated box-centred moments and re-derives them each iteration, which takes icosane from 0.37 to 0.28 GB for about 20% slower Coulomb iterations and identical energies. Parameters are set through dft_obj.multipole_options, e.g. {'precision': 1e-10, 'lmax': 12, 'box_size': 2.5, 'separation': 4.0}; see docs/df_algo12_multipoles.md.

• 10: the previous default, per-function Rys evaluation with sparse triangular storage of the screened three-center integrals. It gives the same energies (identical to ~1e-10 Hartree in the CPU comparisons). It remains selectable for benchmarking.

dft_obj.max_memory_ints3c2e = 2.0   # GB; None (default) = store everything, 0 = recompute every SCF iteration

With a budget smaller than the significant integrals, the most expensive shell-pair blocks are kept in memory and the cheaper ones are re-evaluated in every SCF cycle. With use_gpu=True, algorithm 11 evaluates and contracts shell blocks on the GPU using CuPy and Numba-CUDA. The memory budget limits cached device values. Direct pairs use an additional temporary buffer (up to 1 GB by default); the plan summary reports its size. Cached blocks remain on the device even if keep_ints3c2e_in_gpu=False; use max_memory_ints3c2e to control their storage. The GPU integral builder supports orbital and auxiliary shells through i (l=6). See GPU design and validation for validation coverage, benchmark commands, memory accounting, and limitations of the surrounding driver.

Analytical Forces & Geometry Optimization

After a converged DFT calculation, analytical nuclear gradients (and forces) are available directly via DFT_Grad (density fitting; LDA/GGA/meta-GGA and Skala; CPU):

from pyfock import DFT_Grad

# dftObj must already be converged (dftObj.scf() called)
grad = DFT_Grad(dftObj)
result = grad.calculate()
forces = result["forces"]      # (natoms, 3) in Ha/Bohr
gradient = result["gradient"]  # = -forces

For geometry optimization, use the ASE calculator (requires ase). It uses the analytical forces by default and falls back to finite differences only for configurations the analytical gradients do not yet cover (e.g. HF, no DF):

from ase import Atoms
from ase.optimize import BFGS
from pyfock import PyFockCalculator

water = Atoms("OHH", positions=[[0, 0, 0.119], [0, 0.763, -0.477], [0, -0.763, -0.477]])
water.calc = PyFockCalculator(functional="PBE", basis="def2-SVP",
                              auxbasis="def2-universal-jfit", ncores=4)
BFGS(water).run(fmax=0.02)

By default each step runs in a fresh subprocess, which keeps a crashing or non-converging step from taking the optimizer down and leaves a full PyFock output on disk for every step. The price is that every step repeats the cold start — the imports, loading the Skala checkpoint and warming up TorchScript, and a CUDA context when skala_gpu is on — which for small molecules costs more than the step itself. (PyFock's own Numba kernels are compiled with cache=True, so they reload from disk and are a small part of it.) run_in_process=True runs the steps here instead and keeps all of that alive between them:

PyFockCalculator(functional="skala-1.1", run_in_process=True)

It is the faster choice for Skala by a wide margin; keep the subprocess for long unattended runs. BLAS reads its thread count from the environment when numpy is first imported, so with run_in_process=True set OMP_NUM_THREADS at the top of your script rather than relying on ncores alone. See examples/ex45_ASE_geometry_optimization_with_Skala.py.

Grid response

By default the XC gradient treats the quadrature grid as fixed — its dependence on the nuclear positions is not differentiated. Passing grid_response=True adds the two terms that removes: the grid points of an atom translating with it, and the Becke weights depending on every nucleus. The forces then become exactly translationally invariant:

DFT_Grad(dftObj, grid_response=True).calculate()
PyFockCalculator(functional="PBE", basis="def2-SVP", grid_response=True)   # same flag in ASE

It costs one extra pass over the partitioning and needs the native ('treutler') grids. For semilocal functionals it is a small correction and off by default, so existing results are unchanged. For Skala it is mandatory and on by default — its features are integrals over each atomic grid, so a frozen grid would be wrong by ~1e-2 Ha/Bohr, the size of the forces themselves.

Generating Visualization Files

from pyfock import Utils

# Generate cube files for molecular orbitals and density
Utils.write_density_cube(dftObj, filename='benzene_density.cube')

Graphical User Interface

PyFock includes a web-based GUI for interactive calculations and visualization:

🌐 Try it online: https://pyfock-gui.bragitoff.com

GUI Features

• Interactive 3D Visualization: View molecules and molecular orbitals using Py3Dmol
• Easy Configuration: Select basis sets, functionals, and calculation parameters
• Automatic Cube File Generation: HOMO, LUMO, and density visualizations
• Input Script Generator: Export Python code for local execution
• PySCF Validation: Built-in comparison with PySCF for accuracy verification
• Molecule Library: Pre-loaded common molecules or custom XYZ input

Running GUI Locally

The GUI source code is available on GitHub and can be run locally:

git clone https://github.com/manassharma07/PyFock-GUI.git
cd PyFock-GUI
pip install -r requirements.txt
streamlit run app.py

Tutorials

Interactive Jupyter Notebooks

🚀 Coming Soon: Interactive tutorials on Kaggle and Google Colab

• Kaggle Notebook: Introduction to PyFock
• Kaggle Notebook: Advanced Features and GPU accelerated computations
• Google Colab Notebook: Introduction to PyFock
• Kaggle Notebook: Benchmarking PyFock against PySCF
• Kaggle Notebook: PyFock GPU Dynamic Precision
• Kaggle Notebook: PySCF vs. PyFock GPU Benchmark

Documentation

📚 Full Documentation: https://pyfock-docs.bragitoff.com

• Getting Started Guide
• API Reference
• Examples and Tutorials

Roadmap

• Density Fitting with Cauchy-Schwarz screening
• GPU acceleration for integrals and XC evaluation
• DIIS convergence acceleration
• Web-based GUI
• Rys quadrature (roots 1–10)
• Analytical nuclear gradients & forces (density fitting; LDA/GGA/meta-GGA; CPU)
• ASE calculator & geometry optimization
• Analytical gradients on GPU and for non-DF / ECP calculations
• Electron dynamics & Excited state calculations (RT-TDDFT)
• Periodic boundary conditions
• Hybrid functionals with exact exchange (native B3LYP/PBE0 and LibXC hybrids, RI-K via DF_algo=11; CPU)
• Skala neural exchange-correlation functional, with analytical gradients (CPU; model on GPU via skala_gpu=True)
• Grid-response XC gradients (Becke weight derivatives, Grids.becke_weight_gradient)
• Multi-GPU parallelization
• Basis set optimization tools

See the open issues for a full list of proposed features and known issues.

Contributing

Contributions are what make the open-source community such an amazing place to learn, inspire, and create. Any contributions you make are greatly appreciated.

If you have a suggestion that would make PyFock better, please fork the repo and create a pull request. You can also simply open an issue with the tag "enhancement".

1. Fork the Project
2. Create your Feature Branch (git checkout -b feature/AmazingFeature)
3. Commit your Changes (git commit -m 'Add some AmazingFeature')
4. Push to the Branch (git push origin feature/AmazingFeature)
5. Open a Pull Request

Don't forget to give the project a star! ⭐ Thanks!

License

Distributed under the MIT License. See LICENSE file for more information.

Citation

If you use PyFock in your research, please cite:

@article{sharma2026pyfock,
  title        = {PyFock: A Just-In-Time Compiled Gaussian Basis DFT Python Code for CPU and GPU Architectures},
  author       = {Sharma, Manas and Sierka, Marek},
  journal = {The Journal of Physical Chemistry A},
  year = {2026},
  month = {08},
  issn = {1089-5639},
  doi = {10.1021/acs.jpca.6c03727},
  url = {https://doi.org/10.1021/acs.jpca.6c03727},
  eprint = {https://pubs.acs.org/jpcafh/article-pdf/doi/10.1021/acs.jpca.6c03727/67108190/acs.jpca.6c03727.pdf},
}

Journal Article Link

PyPI Package: https://pypi.org/project/pyfock/

Contact

Manas Sharma

• Email: manas.sharma@uni-jena.de
• Website: manas.bragitoff.com
• LinkedIn: linkedin.com/in/manassharma07
• Project Homepage: https://pyfock.bragitoff.com
• Project Link: https://github.com/manassharma07/pyfock

Built With