README.md

PyIntPak

PyIntPak is a readable, non-recursive Python program for computing one-, two-, three-, and four-electron integrals over Gaussian basis functions and Gaussian geminals.¹²³⁴⁵ It is the standalone Python successor to IntPak, with NumPy and SciPy replacing the original Fortran and SLATEC numerical path.

PyIntPak can be run with the same style of molecule and basis arguments as QuAcK:

python3 PyIntPak.py -x H2 -b 6-31g

This command reads mol/H2.xyz and basis/6-31g, then writes the integral files and a calculation report to int/. QuAcK is not required.

Current status: the one- and two-electron paths are suitable for producing QuAcK input. Three- and four-electron tensors are implemented and tested, but grow as N^6 and N^8; start with small systems when enabling them.

Installation

PyIntPak requires Python 3.10 or newer, NumPy, and SciPy.

git clone https://github.com/pfloos/PyIntPak.git
cd PyIntPak
python3 -m venv .venv
source .venv/bin/activate
python -m pip install .

On Windows PowerShell, activate the environment with:

.venv\Scripts\Activate.ps1

Installing the package provides a pyintpak command. Either form can be used from the repository root:

python3 PyIntPak.py -x H2 -b 6-31g
pyintpak -x H2 -b 6-31g

Run python3 PyIntPak.py -h or pyintpak -h to see every option.

Common examples

Use another bundled molecule and basis set:

python3 PyIntPak.py -x N2 -b cc-pvdz

Use your own XYZ geometry and NWChem/Basis Set Exchange format basis file:

python3 PyIntPak.py -x /path/to/water.xyz -b /path/to/aug-cc-pvdz

XYZ coordinates are interpreted as angstrom by default. Use --bohr when the coordinates are already in bohr:

python3 PyIntPak.py -x /path/to/molecule.xyz -b cc-pvdz --bohr

Specify charge, spin multiplicity, and frozen-core electrons:

python3 PyIntPak.py -x O2 -b cc-pvdz -c -1 -m 2 -fc 2

Generate Cartesian rather than real-spherical integrals, and also write the formatted two-electron file:

python3 PyIntPak.py -x He -b 6-31g --cartesian --formatted_2e

Real-spherical functions are the default. Their normalization, phases, and ordering match PySCF and QuAcK through J shells.

Selecting integrals

The file input/int controls the geminal exponent and enabled integral families. Comments and blank lines are ignored; the five data lines are:

1.0
F F F
T T F F
T F F
F F F

They mean, in order:

1. Slater-geminal exponent zeta;
2. overlap, kinetic, and nuclear-attraction switches;
3. ERI, f12, Yukawa, and error-function Coulomb switches;
4. three-electron Types 1, 2, and 3;
5. four-electron Types 1, 2, and 3.

For example, the third line above enables ordinary ERIs and Slater-geminal two-electron integrals:

T T F F

Every switch is honored in both legacy and QuAcK-style -x/-b runs. Dipole matrices are generated as a QuAcK auxiliary only when at least one one-electron integral is enabled. All formatted and binary two-electron output uses physicist notation.

Pass --input PATH to use a different option directory. If no input/int is present during a QuAcK-style run, PyIntPak defaults to overlap, kinetic, nuclear-attraction, dipole, and electron-repulsion integrals. Disabling files needed by a subsequent QuAcK calculation is allowed, but naturally produces an incomplete QuAcK input set.

PyIntPak does not delete files already present in the output directory. When changing selections, use a fresh output directory or consult IntPak.out to distinguish datasets produced by the current run from older files.

Operator definitions

PyIntPak writes descriptive operator names to both the terminal and int/IntPak.out. The notation is

• Cij = 1/r_ij for a Coulomb factor;
• fij = exp(-zeta r_ij) for a Slater factor represented by the six-Gaussian expansion.

Selection	Descriptive name	Operator
2e ERI	Electron repulsion	C12
2e F12	Slater geminal	f12
2e Yukawa	Screened Coulomb	C12 f12
2e Erf	Long-range Coulomb	erf(zeta r12) C12
3e Type 1	Coulomb-geminal chain	C12 f23
3e Type 2	Cyclic	C12 f13 f23
3e Type 3	Pure-geminal chain	f13 f23
4e Type 1	Unlinked-geminal chain	C12 f14 f23
4e Type 2	Trident	C12 f13 f14
4e Type 3	Linked-geminal chain	C12 f13 f34

The three- and four-electron names follow Fig. 1 of Barca and Loos.⁵

Output files

The default output directory is int/. A normal QuAcK-style run can produce:

File	Contents
IntPak.out	System, geometry, settings, operator definitions, counts, filenames, and timings
nBas.dat	Number of atomic-orbital basis functions
ENuc.dat	Nuclear-repulsion energy
Ov.dat	Overlap matrix
Kin.dat	Kinetic-energy matrix
Nuc.dat	Nuclear-attraction matrix
x.dat, y.dat, z.dat	Cartesian dipole matrices, when one-electron integrals are enabled
ERI.bin	Complete native-double ERI tensor in QuAcK physicist-index order
ERI.dat	Optional formatted ERIs, enabled with --formatted_2e
F12.dat, Yuk.dat, Erf.dat	Selected two-electron geminal operators
3eInt_TypeN.dat	Selected three-electron dataset
4eInt_TypeN.dat	Selected four-electron dataset

Choose another integral directory or report path with:

python3 PyIntPak.py -x H2 -b 6-31g \
  --output /path/to/integrals \
  --report /path/to/PyIntPak.out

PyIntPak checkpoints the report while it runs. If a long calculation is interrupted, the report retains the last completed dataset and a running status rather than disappearing.

Using the files with QuAcK

When selected in input/int, PyIntPak writes the metadata, one-electron matrices, dipoles, and ERI.bin expected by QuAcK. To write directly into a QuAcK checkout:

python3 PyIntPak.py \
  --working_dir /path/to/QuAcK \
  -x N2 -b cc-pvdz

Then run QuAcK normally:

/path/to/QuAcK/bin/QuAcK /path/to/QuAcK

PyIntPak also creates QuAcK-compatible input/molecule and input/eta_opt.dat. The --mmap_2e argument is accepted for command-line compatibility; PyIntPak always streams the ERI tensor through a memory map.

Screening, symmetry, and debugging

QuAcK-style runs enable exact permutation symmetry and Cauchy-Schwarz screening at 1e-15 by default. Change or disable screening with:

python3 PyIntPak.py -x N2 -b cc-pvdz --screening-threshold 1e-12
python3 PyIntPak.py -x N2 -b cc-pvdz --no-screening

Disable or explicitly enable symmetry with --no-symmetry or --symmetry. For exhaustive debugging:

python3 PyIntPak.py -x H2 -b 6-31g --no-symmetry --no-screening

The --limit N option stops each enabled multi-electron dataset after N contracted integrals. A limited ERI tensor is deliberately named ERI.partial.bin and cannot be read as a complete QuAcK tensor. An explicit scratch output directory is therefore required:

python3 PyIntPak.py \
  --input /path/to/debug-options \
  --output /tmp/pyintpak-check \
  -x N2 -b cc-pvdz --formatted_2e --limit 10

Basis and molecule libraries

The repository includes the QuAcK molecule library and the 6-31g and cc-pvdz basis files. Both -x and -b accept either a library name or a path. Add XYZ files under mol/ and NWChem-format basis files under basis/.

Additional basis sets can be installed from the Basis Set Exchange:

python -m pip install basis-set-exchange
cd basis
python install.py

Python API

The command-line program is the supported end-user interface. The integral kernels and contraction drivers can also be imported for research and testing:

from intpak import electron_repulsion, overlap

The public exports are listed in intpak/__init__.py. Their arguments use the small data classes in intpak/models.py; consult the tests for compact working examples.

Validation

Run the complete suite from the repository root:

python -m unittest discover -s tests -v

The tests separate historical IntPak compatibility from independent mathematical validation. References include closed-form values, dense Gaussian integration, direct quadrature, finite differences, exact symmetry identities, and black-box command-line runs. Continuous integration tests Python 3.10, 3.12, and 3.13.

See VALIDATION.md for scientific coverage and numerical tolerances, tests/README.md for test-maintenance guidance, and PORTING.md for implementation status.

Numerical implementation

• Integral recurrence dependency graphs are evaluated iteratively; Python functions do not recursively call themselves.
• The generalized Boys sequence uses SciPy's incomplete-gamma function and an equivalent confluent-hypergeometric representation near zero.
• NumPy provides vectorized Boys evaluation and memory-mapped ERI tensors.
• Cartesian recurrences are transformed at the contracted-shell boundary using cached real solid-harmonic matrices.
• Output records are streamed rather than retaining entire three- and four-electron tensors in memory.

Known limitations

• The four-electron implementation is derived from the published equations because the corresponding original IntPak computation was incomplete. It is validated independently by direct quadrature, center derivatives, and symmetry tests.
• Effective-core-potential integrals are not implemented. Basis files with ECP sections are rejected explicitly.
• Multi-electron datasets become extremely large. Start with a small molecule, a small basis, and --limit when enabling three- or four-electron operators.

The original Fortran implementation remains available in the IntPak repository.

License

PyIntPak is distributed under the GNU General Public License v3.0 or later. See LICENSE.

Contributor: Pierre-Francois Loos

Footnotes

1. P. M. W. Gill, Adv. Quantum Chem. 25, 141 (1994). ↩

2. R. Ahlrichs, Phys. Chem. Chem. Phys. 8, 3072 (2006). ↩

3. G. M. J. Barca, P. F. Loos and P. M. W. Gill, J. Chem. Theory Comput. 12, 1735 (2016). ↩

4. G. M. J. Barca and P. F. Loos, Adv. Quantum Chem. 76, 147 (2018). ↩

5. G. M. J. Barca and P. F. Loos, J. Chem. Phys. 147, 024103 (2017). ↩ ↩²