Quickstart

Coulomb matrix computation

This example demonstrates the complete workflow: building a molecule with PySCF, converting bases to libintx structures, constructing the DF metric (P|Q), and computing the Coulomb matrix J.

import numpy as np
from pyscf import gto, scf, df
import libintx
from libintx.tools import pyscf_to_libintx_basis

# Build a PySCF molecule
mol = gto.Mole()
mol.atom = '''
O  0.000000  0.000000  0.000000
H  0.000000 -0.757160  0.586260
H  0.000000  0.757160  0.586260
'''
mol.basis = 'def2-svp'
# Spherical harmonic GTOs are not yet supported:
mol.cart = True
mol.build()

# Build an auxiliary basis and its two-center metric (P|Q)
auxmol = df.addons.make_auxmol(mol)  # set up the default auxiliary basis
V = auxmol.intor('int2c2e', hermi=1)  # 2-center 2-electron integrals matrix

# Convert PySCF bases to libintx Basis
basis = pyscf_to_libintx_basis(mol)
aux_basis = pyscf_to_libintx_basis(auxmol)

# Create the J-engine and compute J
jeng = libintx.Jengine(basis, aux_basis, V, vinv_type='cholesky', vinv_device='cuda')

D = np.eye(mol.nao)

J = jeng.compute_j(D, scale=0.5)

# Check the results
mf = scf.RHF(mol).density_fit()
J_ref = mf.get_j(dm=D)

print(f'|J_ref - J| = {np.linalg.norm(J_ref-J)}')

Schwarz Screening

Schwarz screening can dramatically reduce compute time, especially for large systems. You can enable screening using per-shell two-center two-electron \((P|P)\) and four-center two-electron \((ij|ij)\) estimates.

Computation:

\((P|P)\) can be computed from the diagonal AO blocks of (P|Q) aggregated per auxiliary shell. It should be a vector of length nsh_aux, where nsh_aux is the number of shells in the auxiliary basis (len(aux_basis)).
\((ij|ij)\) can be initialized via PySCF’s RHF helper. It should be a matrix of shape (nsh, nsh), where nsh is the number of shells in the primary basis (len(basis)).

Note

Some PySCF bases are general-contracted; convert them to segmented-contraction variants before computing screening quantities.

An example code:

import numpy as np
from pyscf import scf
import libintx

# Compute screening quantities
pp = np.sqrt(np.array([np.max(np.abs(V[i:j,i:j])) for i,j in aux_basis.aos_by_shell()]))
ijij = scf.RHF(mol).init_direct_scf().q_cond

# Enable screening
jeng.set_screening(pp, ijij, cutoff=1e-12)

Solver Options

To use CUDA-based solver, pass vinv_device=’cuda’. Optionally adjust vinv_mem_fraction (default 0.5) to limit the maximum fraction of GPU memory for the solver.

jeng = libintx.Jengine(
    basis, aux_basis, V,
    vinv_device='cuda',
    vinv_mem_fraction=0.5
)

If memory limits are exceeded, the Cholesky/LU factors will be stored on the host and loaded to the GPU when needed. Consider distributed Cholesky solver (cholesky_mp) for very large molecules (requires MPI).