Seminumerical exchange (SGX)#

Modules: pyscf.sgx

The SGX module implements seminumerical computation of the exchange matrix (and also the Coulomb matrix).

Introduction#

Direct computation of the Hartree-Fock exchange matrix in the atomic orbital basis scales poorly with system size. To achieve better scaling, one three-dimensional integral in the 6-dimensional two-electron integrals can be computed analytically, while the other can be evaluated on a real-space grid, as proposed by Friesner [88]. The PySCF implementation resembles the chain-of-spheres (COSX) algorithm of Neese et al. [89], but uses slightly different default grids and a modified P-junction screening algorithm. Overlap fitting is used to reduce aliasing errors [90]. Before screening of negligible integrals, SGX scales as \(O(N^3)\) with system size, as opposed to the \(O(N^4)\) scaling of analytical exchange [89]. This makes its performance favorable for large basis sets and large systems. In principle, both SGX and analytical exchange scale asymptotically as \(O(N^2)\) if negligible integrals are screened out, and as \(O(N)\) if the system has a gap and negligible density matrix components are screened. However, the linear scaling is difficult to achieve in practice and only seen for some very large systems.

Usage and Example#

Any scf.hf.SCF object mf can be converted to an equivalent object that computes the Coulomb and exchange matrices with SGX instead of analytical integration by calling sgx.sgx_fit(mf).

#!/usr/bin/env python

'''
This example shows how to use pseudo spectral integrals in SCF calculation.
'''

from pyscf import gto
from pyscf import scf
from pyscf import sgx
mol = gto.M(
    atom='''O    0.   0.       0.
            H    0.   -0.757   0.587
            H    0.   0.757    0.587
    ''',
    basis = 'ccpvdz',
)
# Direct K matrix for comparison
mf = scf.RHF(mol)
mf.kernel()

# Using SGX for J-matrix and K-matrix, without using P-junction screening
mf = sgx.sgx_fit(scf.RHF(mol), pjs=False)
mf.kernel()

# Using RI for Coulomb matrix while K-matrix is constructed with SGX method
mf.with_df.dfj = True
mf.kernel()

# Turn on P-junction screening to accelerate large calculations
# (uses algorithm similar to COSX)
mf.with_df.pjs = True
mf.kernel()

# direct_scf_sgx turns on direct SCF for SGX
# JK matrix is rebuilt from scratch every rebuild_nsteps steps
mf.direct_scf_sgx = True
# If grids_level_i == grids_level_j, no grid switch occurs
mf.with_df.grids_level_i = 1
mf.kernel()

# If dfj is off at runtime, it is turned on and a user warning is issued
# because SGX-J cannot be used with P-junction screening.
mf.with_df.dfj = False
mf.kernel()

# Use direct J-matrix evaluation (slow, primarily for testing purposes)
mf = sgx.sgx_fit(scf.RHF(mol), pjs=False)
mf.with_df.direct_j = True
mf.with_df.dfj = False
mf.kernel()

In this case, the error of DF-J+SGX-K compared to analytical exchange is about 0.03 mEh. The line

mf.with_df.dfj = True

specifies to compute the Coulomb matrix using density fitting (DF-J) while using SGX for the exchange matrix.

Passing pjs=True to sgx_fit sets dfj=True and also turns on P-junction screening (optk=True), which means that the three-center integrals are screened by the density matrix. If a three-center integral gets contracted only with negligibly small density matrix elements, it is not computed. Because the J-matrix cannot be screened in this manner, optk=True only activates P-junction screening when used with dfj=True. If dfj=False, optk is ignored.

When optk=True, the P-junction screening threshold is determined by two parameters: mf.with_df.sgx_tol_energy and mf.with_df.sgx_tol_potential. These parameters correspond to approximate upper bounds on the error of the exchange energy and K-matrix due to P-junction screening, respectively. By default, sgx_tol_energy=sgx_tol_potential="auto", in which case sgx_tol_energy is set to mf.direct_scf_tol and sgx_tol_potential=sgx_tol_energy**0.5. For more aggressive screening that is still numerically stable and sufficiently precise for most cases, we recommend sgx_tol_energy=mf.conv_tol or 0.1*mf.conv_tol and sgx_tol_potential="auto". The first thing to try if you have numerical issues is setting mf.rebuild_nsteps=1. This turns off incremental building of the K-matrix, which is faster but can accumulate errors between SCF steps.

Direct evaluation of the J matrix can be used by setting mf.with_df.direct_j = True, but this is much slower than SGX-J or DF-J and only intended for testing purposes.

Adjustable Parameters#

Calling the sgx_fit function on an scf.hf.SCF object returns an equivalent SGXHF object with a with_df attribute that handles SGX integration. To use a non-default auxiliary basis (for dfj=True), auxbasis can be specified in the call to sgx_fit. The other input option to sgx_fit is pjs. If pjs=False (default), dfj=optk=False, otherwise dfj=optk=True. In addition, there are various adjustable parameters for SGX that can be set after initialization.

with_df attributes:

  • grids_level_i: The grid level to use for initial SCF iterations.

  • grids_level_f: The grid level to use for final SCF iterations.

  • dfj: If True, density fitting is used for the J-matrix. If False, SGX is used.

  • optk: If True, P-junction screening is used. Threshold determined by sgx_tol_energy and sgx_tol_potential. Only applies if dfj=True.

  • sgx_tol_energy and sgx_tol_potential: Sets the approximate upper bound on the energy and K-matrix errors due to density matrix screening. Ignored if optk=False or dfj=False. See the Usage and Example section above for instructions on how to set these, as well as the docstrings for more details.

  • use_opt_grids: Whether to use optimized grids for SGX based on COSX in ORCA [91]. Default is True.

  • fit_ovlp: Whether to numerically fit the overlap matrix to improve numerical precision. Default is True.

  • grids_thrd: When the values of all atomic orbitals, multiplied by the qudrature weight, are less than this threshold for a given grid point, it is removed from the SGX numerical grid.

  • grids_switch_thrd: The threshold for the magnitude of the change in the density matrix that triggers the switch from the initial grid specified by grids_level_i to the final grid specified by grids_level_f.

  • bound_aglo: Determines how to estimate the maximum value of the three-center SGX integrals. Options are "ovlp", "sample", and "sample_pos" (default). "ovlp" assumes the maximum value of each integral is roughly equal to the maximum overlap between orbitals in a shell (with a value of 1 for shells on the same atom to account for orthogonality). "sample" constructs a coarse grid for each atom pair and computes the maximum value of the integrals over these grids. "sample_pos" does the same but also contructs an approximate estimate for the integral inside each grid block (position-dependent estimate), which can produce a small speedup for global exchange and a larger speedup for short-range exchange. Usually the default value "sample_pos" should be fine.

  • direct_j: If True, direct evaluation of the J-matrix is used (slow, for testing only).

SGXHF attributes:

  • auxbasis: The auxiliary basis for J-matrix density fitting. Used if dfj=True.

  • rebuild_nsteps: By default, when mf.direct_scf=True, the K-matrix is constructed incrementally from the change in the density matrix between SCF steps (as is the case for other exchange algorithms in PySCF). The SGX J and K matrices are then rebuilt from scratch every rebuild_nsteps steps (default 5). Set to 1 to rebuild from scratch at every step, and set to greater than the max number of cycles to turn off rebuilding the JK matrix (Warning: This can cause creeping numerical error). The SGX matrix is built from scratch when the numerical grid changes, regardless of the value of rebuild_nsteps, because the same density matrix will have different energies for different grid sizes.