API Reference

This page contains the API reference for the prop_sapt package.

Dimer Class

Definition of a core Dimer (molecule-like) class.

class prop_sapt.molecule.Dimer(geometry, reference='RHF', memory=None, **kwargs)

Bases: sinfinity

Core dimer (molecule-like) class.

Initialises and stores all necessary values for ‘propSAPT’ calculations. Performs SCF calculations for the monomers.

The fragment structure of the geometry (atoms, charges and multiplicities of monomers A and B and of any extra fragments) is in fragments, a prop_sapt.fragments.DimerFragments.

The settings the dimer is built with are kept: psi4_options, the Psi4 options changed from their defaults, and fitting_bases, the fitting basis sets as built (basis.blend()) by option, “DF_BASIS_SCF” and “DF_BASIS_SAPT”. Changing the options afterwards does not change the dimer.

Memory comes from Psi4’s setting (psi4.set_memory), as for any Psi4 calculation.

Parameters:

• geometry (str) – Dimer geometry in Psi4 format, with monomers A and B as its first two fragments; any further fragments are ghost atoms.

• reference (str, optional) – “RHF” or “RKS”. Defaults to “RHF”.

• memory (float | None, optional) – Deprecated, set psi4.set_memory instead. Memory of the JK object in GB; None, the default, follows Psi4’s memory.

Keyword Arguments:

• functional (str) – Functional of an RKS reference. Defaults to the default_functional option.

• grac_A (float) – GRAC shift of monomer A (Psi4’s dft_grac_shift). Defaults to no shift.

• grac_B (float) – GRAC shift of monomer B. Defaults to no shift.

clear_theta_amplitudes()

Free the theta amplitudes of the exchange density matrices, see get_theta_amplitudes. The exchange-dispersion ones are as large as the dispersion amplitudes. calc_densities frees them when it is done; call this after calling calc_density_matrix for both monomers yourself.

evaluate_density(dm: ndarray | Sequence[ndarray], points: ndarray, basisset: BasisSet | None = None, chunk_size: int | None = None) → ndarray

Evaluate the density of AO density matrix dm at points. No spin factor is applied: for the total density of an RHF monomer pass 2 * Da.

Parameters:

• dm (np.ndarray | Sequence[np.ndarray]) – AO density matrix, shape (nbf, nbf), or several of them, as a sequence or an array of shape (M, nbf, nbf).

• points (np.ndarray) – Points in bohr, shape (…, 3).

• basisset (psi4.core.BasisSet | None, optional) – Basis set of dm. Defaults to the dimer basis set.

• chunk_size (int | None, optional) – Number of points evaluated at once. Defaults to the point_eval_chunk_size option.

Returns:

Density at the points, shape (…) for a single matrix

and (M, …) for several.

Return type:

np.ndarray

evaluate_esp(dm: ndarray | Sequence[ndarray], points: ndarray, nuclei: bool = False, basisset: BasisSet | None = None) → ndarray

Evaluate the electrostatic potential of AO density matrix dm at points, from analytic potential integrals. Electrons contribute negatively.

Parameters:

• dm (np.ndarray | Sequence[np.ndarray]) – AO density matrix, shape (nbf, nbf), or several of them, as a sequence or an array of shape (M, nbf, nbf).

• points (np.ndarray) – Points in bohr, shape (…, 3).

• nuclei (bool, optional) – Add the potential of the nuclei of the basis set’s molecule, use it for total (not interaction-induced) densities. Defaults to False.

• basisset (psi4.core.BasisSet | None, optional) – Basis set of dm. Defaults to the dimer basis set.

Returns:

Electrostatic potential at the points, shape (…) for a single

matrix and (M, …) for several.

Return type:

np.ndarray

evaluate_orbital(coeffs: ndarray, points: ndarray, basisset: BasisSet | None = None, chunk_size: int | None = None) → ndarray

Evaluate orbitals given by MO coefficients coeffs at points. Orbitals are the columns of coeffs, as in Psi4’s Ca().

Parameters:

• coeffs (np.ndarray) – MO coefficients, shape (nbf,) for one orbital or (nbf, k) for k orbitals.

• points (np.ndarray) – Points in bohr, shape (…, 3).

• basisset (psi4.core.BasisSet | None, optional) – Basis set of coeffs. Defaults to the dimer basis set.

• chunk_size (int | None, optional) – Number of points evaluated at once. Defaults to the point_eval_chunk_size option.

Returns:

Orbital values at the points, shape (…) for coeffs of shape

(nbf,) and (k, …) for coeffs of shape (nbf, k).

Return type:

np.ndarray

fragment_operators(weights: str | Mapping[str, ndarray] = 'becke', grid: QuadratureGrid | None = None, dipole: bool = True, chunk_size: int | None = None, partition_of_unity: bool = True, **grid_kwargs) → FragmentOperators

Compute the fragment population operators of monomers A and B in the dimer basis set, see prop_sapt.partition.FragmentOperators.

Parameters:

• weights (str | Mapping[str, np.ndarray], optional) – Fragment weights: “becke” for the Becke cells of each monomer’s atoms, “hirshfeld” for rho_X / (rho_A + rho_B) from the unperturbed monomer densities, or the weight of each fragment at the points of grid. Defaults to “becke”.

• grid (QuadratureGrid | None, optional) – Quadrature grid, required for weights given as arrays. Defaults to molecular_grid of the dimer.

• dipole (bool, optional) – Also compute the first-moment matrices, needed for dipoles. Defaults to True.

• chunk_size (int | None, optional) – Number of points evaluated at once. Defaults to the point_eval_chunk_size option.

• partition_of_unity (bool, optional) – Whether the fragments of weights together cover the whole system, i.e. their weights sum to 1 at every point. If True, this is checked, and every population is checked to add up to the electrons of the density, Tr(D S), see FragmentOperators.from_grid. Pass False for only some of the fragments. Defaults to True.

• **grid_kwargs – Passed to prop_sapt.partition.molecular_grid when no grid is given, e.g. radial_points or spherical_points.

Returns:

Operators of fragments “A” and “B” for named weights,

or of the fragments of weights.

Return type:

FragmentOperators

get_cpscf_ra() → ndarray

Grab CPSCF response amplitudes for monomer A or calculate them if they haven’t been calculated already.

Returns:

Response amplitudes for monomer A.

Return type:

np.ndarray

get_cpscf_sb() → ndarray

Grab CPSCF response amplitudes for monomer B or calculate them if they haven’t been calculated already.

Returns:

Response amplitudes for monomer B.

Return type:

np.ndarray

get_psi4_basisset() → BasisSet

Get ‘psi4.core.BasisSet’ object of the dimer: the one it was built with, whatever the BASIS option is now.

Returns:

Psi4 basis set object of the dimer.

Return type:

psi4.core.BasisSet

get_psi4_molecule() → Molecule

Get ‘psi4.core.Molecule’ object of the dimer.

Returns:

Psi4 molecule object of the dimer.

Return type:

psi4.core.Molecule

get_psi4_wavefunction() → Wavefunction

Get ‘psi4.core.Wavefunction’ object of the dimer, in the basis set the dimer was built with, whatever the BASIS option is now.

Returns:

Psi4 wavefunction object of the dimer.

Return type:

psi4.core.Wavefunction

get_theta_amplitudes(key: str, build: Callable[[Dimer], ndarray | tuple[ndarray, ...]]) → ndarray | tuple[ndarray, ...]

Grab the theta amplitudes key of an exchange density matrix, or calculate them with build if they haven’t been calculated already.

The density matrix of each monomer needs the amplitudes of both, so they are kept for the other monomer’s calc_density_matrix call. They are read-only. Free them with clear_theta_amplitudes.

Parameters:

• key (str) – Density component, e.g. “exch-ind_S2”.

• build (Callable[[Dimer], np.ndarray | tuple[np.ndarray, ...]]) – Calculates the amplitudes for this dimer.

Returns:

The amplitudes, as returned by build.

Return type:

np.ndarray | tuple[np.ndarray, …]

make_cube(matrix: ndarray | Sequence[ndarray], obj_type: str | Iterable[str] = 'density', **kwargs) → Cube | list[Cube]

Create a Cube object with volumetric data of the given matrix calculated on a grid.

Parameters:

• matrix (np.ndarray | Sequence[np.ndarray]) – Either density matrix, orbital coefficients vector, or list of such objects.

• obj_type (str | Iterable[str], optional) – Specifies the matrix argument type(s), either “density”, “orbital”, “esp” or a list with such strings. Defaults to “density”.

• **kwargs – Passed to prop_sapt.cubes.make_cube, e.g. grid_step or esp_nuclei. basisset defaults to the dimer’s basis set, see get_psi4_basisset.

Returns:

Resulting cube(s).

Return type:

Cube | list[Cube]

save_cube(matrix: ndarray | Sequence[ndarray], obj_type: str | Iterable[str] = 'density', filename: str | PathLike | Iterable[str | PathLike] = 'density.cube', **kwargs)

Evaluate matrix values on a grid and save to a .cube file.

Parameters:

• matrix (np.ndarray | Sequence[np.ndarray]) – Either density matrix, orbital coefficients vector, or list of such objects.

• obj_type (str | Iterable[str], optional) – Specifies the matrix argument type(s), either “density”, “orbital”, “esp” or a list with such strings. Defaults to “density”.

• filename (str | os.PathLike | Iterable[str | os.PathLike], optional) – File path(s) to save resulting cube(s), one per matrix. Defaults to “density.cube”.

Raises:

ValueError – If the number of filenames differs from the number of matrices.

store_theta_amplitudes(key: str, amplitudes: ndarray | tuple[ndarray, ...]) → None

Keep the theta amplitudes key, made read-only, for get_theta_amplitudes.

Parameters:

• key (str) – Density component, e.g. “exch-ind_S2”.

• amplitudes (np.ndarray | tuple[np.ndarray, ...]) – The amplitudes.

Properties

Calculations of SAPT one-electron properties.

prop_sapt.properties.calc_induced_dipole(dimer: Dimer, response: bool | str = True) → DataFrame

Calculate interaction-induced dipole moment along X, Y and Z axes.

The CPSCF equations of the three axes are solved together, a batch per monomer at a time, which costs much less than solving them axis by axis.

Parameters:

• dimer (Dimer) – A dimer system.

• response (bool | str, optional) – Induction amplitudes, True (coupled), False (uncoupled) or “both”; see perform_property_contractions. Defaults to True.

Returns:

DataFrame with results.

Return type:

pd.DataFrame

prop_sapt.properties.calc_property(dimer: Dimer, prop: str | ndarray, response: bool | str = True, **kwargs) → DataFrame

Calculation of interaction-induced property.

The property is the field derivative of the energy calc_sapt_energy gives for the same response.

Known issue: the S^inf exchange-dispersion term x2_exch-disp is not yet the field derivative of EXCH-DISP2. For the dipole of He-H2 in aug-cc-pVTZ it is 7.2e-5 a.u. off the finite-field derivative (about 18%), and the total x_induced carries that error. Every term of x_induced_S2 agrees with its finite-field derivative. See “Known limitations” in the usage guide.

Parameters:

• dimer (Dimer) – A dimer system.

• prop (str | np.ndarray) – Property to calculate. Either name of already implemented property, e.g. dipole moment (prop=”dipole”) or an array with matrix representation, in AO basis, of the operator corresponding to a property to calculate.

• response (bool | str, optional) – Induction amplitudes of the second-order induction terms. True: the coupled (CPSCF) ones, x2_ind,r, x2_exch-ind,r_S2 and x2_exch-ind,r (plus the uncoupled x2_ind). False: the uncoupled ones, x2_ind, x2_exch-ind_S2 and x2_exch-ind, with no CPSCF equations solved for them. “both”: both sets. The totals take the coupled terms unless response is False. The monomers relax in the field in every mode, so the other terms do not depend on it. Defaults to True.

Keyword Arguments:

results (str) – Path of the CSV file the results are saved to. Defaults to “results.csv”.

Returns:

DataFrame with results.

Return type:

pd.DataFrame

prop_sapt.properties.perform_property_contractions(dimer: Dimer, property_matrix: ndarray, response: bool | str = True) → DataFrame

Calculates property contributions. Returns pandas.DataFrame with values of separate contributions and the summed-up total.

The property is the field derivative of the energy of calc_sapt_energy with the same response. The monomers relax in the field in every mode (the first-order, dispersion and exchange-dispersion terms do not depend on it); response chooses the induction amplitudes of the second-order induction terms:

• True: the coupled (CPSCF) ones, x2_ind,r, x2_exch-ind,r_S2 and x2_exch-ind,r. The uncoupled x2_ind is reported too.

• False: the uncoupled ones, x2_ind, x2_exch-ind_S2 and x2_exch-ind. No CPSCF equations are solved for the induction amplitudes or their derivatives.

• “both”: both sets.

The totals x_induced and x_induced_S2 take the coupled terms unless response is False, as calc_sapt_energy does.

Parameters:

• dimer (Dimer) – A dimer system.

• property_matrix (np.ndarray) – A matrix representation of the operator corresponding to a property to calculate, in AO basis.

• response (bool | str, optional) – True, False or “both”, see above. Defaults to True.

Returns:

DataFrame with results.

Return type:

pd.DataFrame

Densities

prop_sapt.densities.BATCHED_COMPONENTS = (('exch', <function _exch_density_builder>, 'First-order Exchange'), ('ind', <function _ind_density_builder>, 'Second-order Induction'), ('disp', <function _disp_density_builder>, 'Second-order Dispersion'), ('exch-ind_S2', <function exch_ind_s2_density_builder>, 'Second-order Exchange-Induction (S^2)'), ('exch-ind', <function exch_ind_density_builder>, 'Second-order Exchange-Induction (S^inf)'), ('exch-disp_S2', <function exch_disp_s2_density_builder>, 'Second-order Exchange-Dispersion (S^2)'))

Density components whose CPSCF equations calc_density_matrix solves in batches: key, builder (see prop_sapt.cpscf_batching) and the name in the output.

prop_sapt.densities.DENSITY_COMPONENTS = ('pol', 'exch', 'ind', 'exch-ind_S2', 'exch-ind', 'disp', 'exch-disp_S2', 'Ind', 'Disp', 'total')

Keys of the density components returned by calc_density_matrix, in order; calc_densities returns them as “<key>_A”, “<key>_B” and, in the AO basis, “<key>”.

prop_sapt.densities.MATRIX_FORMATS = ('npz', 'npy')

“npz” for the self-describing bundle (see prop_sapt.artifacts), “npy” for bare files (deprecated).

Type:

Formats of the matrices saved by calc_densities

prop_sapt.densities.calc_densities(dimer: Dimer, orbital_basis='AO', save_cubes=False, cubes_dir='.', save_matrices=False, matrices_dir='.', matrix_format: str | None = None, metadata: Mapping | None = None, compress_matrices: bool | None = None) → dict[str, ndarray]

Calculate interaction-induced changes in the density matrix of dimer (and both monomers).

The matrices are for one spin (the Da convention): the total density of the closed-shell monomers is twice their sum. Saved matrices are always in the AO basis (transformed from the MO basis if needed); saved cubes hold the total density.

Parameters:

• dimer (Dimer) – A dimer system.

• orbital_basis (str, optional) – Select orbital basis of the returned density matrix, either “AO” or “MO”. Defaults to “AO”.

• save_cubes (bool, optional) – Whether to save density cubes. Defaults to False.

• cubes_dir (str, optional) – Directory to save density cubes. Defaults to “.”.

• save_matrices (bool, optional) – Whether to save density matrices. Defaults to False.

• matrices_dir (str, optional) – Directory to save density matrices. Defaults to “.”.

• matrix_format (str | None, optional) – “npz” saves the matrices, the monomer wavefunctions and their metadata into one self-describing bundle, densities.npz, read with prop_sapt.load_densities. “npy” saves bare delta_dm_<key>.npy and wfn_mon_<A/B>.npy files; it is deprecated. Defaults to “npz”.

• metadata (Mapping | None, optional) – JSON-serialisable entries stored in the bundle’s metadata under “extra”. Defaults to None.

• compress_matrices (bool | None, optional) – Deprecated, use matrix_format: True is “npz” and False is “npy”.

Returns:

Density matrices keyed “<component>_A”, “<component>_B”

and, in the AO basis, “<component>” for their sum, with the components of DENSITY_COMPONENTS. In the MO basis there is no sum: each monomer’s matrices are in its own MO basis, and the two are not orthogonal.

Return type:

dict[str, np.ndarray]

prop_sapt.densities.calc_density_matrix(dimer: Dimer, monomer: str, orbital_basis='AO') → dict[str, ndarray]

Calculate interaction-induced change in the density matrix of the ‘monomer’.

The theta amplitudes of the second-order exchange components are shared by the two monomers, so they are kept on dimer for the other monomer’s call; dimer.clear_theta_amplitudes() frees them (calc_densities does this itself).

The CPSCF equations of the components are solved together, with one Dimer.cpscf_batch call per monomer at a time: the first monomer’s call solves its components and the theta amplitudes, then the exchange-induction components that need them (16 equations of the two calls in 4 batches).

Parameters:

• dimer (Dimer) – A dimer system.

• monomer (str) – Select a monomer, either “A” or “B”.

• orbital_basis (str, optional) – Select orbital basis of the returned density matrix, either “AO” or “MO”. Defaults to “AO”.

Returns:

Density matrix components, keyed by DENSITY_COMPONENTS:

”pol” and “exch”: first-order polarisation and exchange corrections; “ind”, “exch-ind_S2” (S^2) and “exch-ind” (S^inf): second-order induction and exchange-induction; “disp” and “exch-disp_S2” (S^2): second-order dispersion and exchange-dispersion; “Ind” and “Disp”: induction and dispersion totals with S^2 exchange; “total”: sum of “pol”, “exch”, “Ind” and “Disp”. Every component is a symmetric matrix.

Return type:

dict[str, np.ndarray]

Warns:

InvariantWarning – If a component is not symmetric within the check_symmetry_atol option, see prop_sapt.checks.check_symmetric.

prop_sapt.densities.density_mo_to_ao(dimer: Dimer, monomer: str, density_matrix: ndarray) → ndarray

Transform density matrix from MO to AO.

Parameters:

• dimer (Dimer) – A dimer system.

• monomer (str) – Select a monomer, either “A” or “B”.

• density_matrix (np.ndarray) – Density matrix to transform.

Returns:

Transformed density matrix.

Return type:

np.ndarray

Density Artifacts

Self-describing density artifacts: interaction-induced density matrices in one .npz bundle that records what they are, and a loader that checks them before use.

calc_densities(…, save_matrices=True) writes the bundle densities.npz; load_densities reads it back. The bundle holds

• metadata: a JSON document (a 0-d string array) with the conventions, the method, the basis set, the Psi4 and propSAPT options, the versions, and the geometry with its fragments (prop_sapt.fragments.DimerFragments);

• delta_dm_<key>: the density matrices, e.g. delta_dm_pol_A;

• wfn_mon_A and wfn_mon_B: the monomer wavefunctions, Psi4’s Wavefunction.to_file() data with its arrays stored as entries wfn_mon_A/... and everything else in a JSON manifest.

Nothing in the bundle is pickled: it is written and read with allow_pickle=False.

Conventions (CONVENTIONS): the matrices are in the AO basis of the dimer basis set, in Psi4 AO order, for one spin (the Da convention of the matrices calc_densities returns). The total density of a closed-shell system is twice the stored matrix, see DensitySet.total; the .cube files written by calc_densities hold the total density.

Classes

DensitySet

Density matrices and metadata read from a bundle, with the monomer wavefunctions.

ArtifactError

Raised for a malformed bundle and for a bundle that does not match the expected geometry or basis set.

Functions

save_densities

Write density matrices of a dimer into a bundle.

load_densities

Read a bundle, checking its conventions and, optionally, its geometry and basis set.

artifact_metadata

The metadata of a dimer recorded in the bundle.

exception prop_sapt.artifacts.ArtifactError

Bases: ValueError

Raised when a bundle is malformed, or does not match what it is checked against.

prop_sapt.artifacts.CONVENTIONS = {'ao_order': 'psi4', 'orbital_basis': 'AO', 'spin': 'alpha', 'total_density_factor': 2}

Conventions of the stored matrices, recorded in and checked on every bundle.

prop_sapt.artifacts.DENSITIES_FILENAME = 'densities.npz'

Filename of the bundle written by calc_densities.

prop_sapt.artifacts.DENSITIES_FORMAT = 'propsapt-densities'

Value of the format metadata entry.

class prop_sapt.artifacts.DensitySet(matrices: dict[str, ndarray], metadata: dict, fragments: DimerFragments, path: Path)

Bases: object

Density matrices read from a bundle by load_densities, with their metadata.

check(dimer=None, geometry: DimerFragments | str | PathLike | None = None, basisset: BasisSet | None = None, atol: float = 1e-06)

Check that the bundle belongs to a dimer, a geometry and/or a basis set.

Parameters:

• dimer (Dimer | None, optional) – Check the geometry, fragments and basis set against this dimer. Defaults to None.

• geometry (DimerFragments | str | os.PathLike | None, optional) – Check the geometry and fragments against these fragments, a geometry string or a geometry file (see prop_sapt.fragments.parse_fragments). Defaults to None.

• basisset (psi4.core.BasisSet | None, optional) – Check the basis set against this one. Defaults to None.

• atol (float, optional) – Largest coordinate difference allowed, in bohr. Defaults to 1e-6.

Raises:

ArtifactError – Naming the first difference found.

property components: tuple[str, ...]

The density components in the bundle, e.g. “pol”, “exch”, …

fragments: DimerFragments

Geometry and fragments of the dimer.

matrices: dict[str, ndarray]

AO basis, Psi4 AO order, one spin, as calc_densities returns them. See total for the total density.

Type:

Density matrices by key (e.g. “pol_A”, “pol_B”, “pol”)

metadata: dict

The metadata of the bundle, see prop_sapt.artifacts.

path: Path

Path of the bundle.

total(key: str) → ndarray

Total (both spins) density matrix of key, twice the stored matrix.

Parameters:

key (str) – Matrix key, e.g. “pol_A”.

Returns:

The total density matrix.

Return type:

np.ndarray

wavefunction(monomer: str) → Wavefunction

Monomer wavefunction rebuilt with Psi4’s Wavefunction.from_file.

Psi4 builds the orbital basis set by its name, so a basis set given by per-element overrides or read from a basis set directory has to be available to Psi4 in the current session; wavefunction_data needs no basis set.

Parameters:

monomer (str) – “A” or “B”.

Raises:

RuntimeError – If Psi4 cannot rebuild the wavefunction, e.g. its basis set.

Returns:

The monomer wavefunction.

Return type:

psi4.core.Wavefunction

wavefunction_data(monomer: str) → dict

Monomer wavefunction data as written by Psi4’s Wavefunction.to_file(), read from the bundle without building any Psi4 object.

Parameters:

monomer (str) – “A” or “B”.

Returns:

e.g. data["matrix"]["Da"] or data["vector"]["epsilon_a"].

Return type:

dict

prop_sapt.artifacts.FORMAT_VERSION = 1

Version of the bundle layout and metadata, raised on incompatible changes.

prop_sapt.artifacts.artifact_metadata(dimer) → dict

Metadata of a dimer recorded in a bundle: format, versions, conventions, method, basis set, options, geometry and fragments. The basis sets and the Psi4 options are those the dimer was built with, whatever the options are now.

Parameters:

dimer (Dimer) – The dimer.

Returns:

JSON-serialisable metadata, without the entries describing the stored

matrices, which save_densities adds.

Return type:

dict

prop_sapt.artifacts.load_densities(path: str | PathLike, dimer=None, *, geometry: DimerFragments | str | PathLike | None = None, basisset: BasisSet | None = None, atol: float = 1e-06) → DensitySet

Read a density bundle written by calc_densities or save_densities.

The bundle is always checked for its format and conventions and for the shapes of its matrices. Give dimer, geometry and/or basisset to also check that it belongs to them, see DensitySet.check.

Parameters:

• path (str | os.PathLike) – The bundle, or a directory holding densities.npz.

• dimer (Dimer | None, optional) – Check geometry, fragments and basis set against this dimer. Defaults to None.

• geometry (DimerFragments | str | os.PathLike | None, optional) – Check geometry and fragments against these fragments, a geometry string or a geometry file. Defaults to None.

• basisset (psi4.core.BasisSet | None, optional) – Check the basis set against this one. Defaults to None.

• atol (float, optional) – Largest coordinate difference allowed, in bohr. Defaults to 1e-6.

Raises:

ArtifactError – If the bundle is malformed, uses unknown conventions, or does not match dimer, geometry or basisset.

Returns:

Matrices, metadata and fragments of the bundle.

Return type:

DensitySet

prop_sapt.artifacts.save_densities(path: str | PathLike, matrices: Mapping[str, ndarray], dimer, metadata: Mapping | None = None, wavefunctions: bool = True) → Path

Write density matrices of a dimer into a bundle.

Parameters:

• path (str | os.PathLike) – Path of the bundle; “.npz” is appended if missing.

• matrices (Mapping[str, np.ndarray]) – AO density matrices of the dimer basis set by key, one spin (the Da convention). Keys ending in “_A” or “_B” are monomer components, see DensitySet.components.

• dimer (Dimer) – The dimer the matrices belong to.

• metadata (Mapping | None, optional) – Additional top-level metadata entries; they cannot replace the entries written here. Defaults to None.

• wavefunctions (bool, optional) – Also store the monomer wavefunctions. Defaults to True.

Raises:

ValueError – If a matrix is not (nbf, nbf), or metadata repeats an entry.

Warns:

InvariantWarning – If a matrix is not symmetric within the check_symmetry_atol option, see prop_sapt.checks.check_symmetric; it is saved anyway.

Returns:

Path of the written bundle.

Return type:

Path

Invariant Checks

Invariant checks: cheap checks of properties that propSAPT results must have.

A failed check warns with InvariantWarning, attributed to the first caller outside propSAPT, and writes the same message into the Psi4 output file; it never stops a calculation. Turn the warnings into errors with warnings.simplefilter("error", InvariantWarning). Each check takes its tolerance from a propSAPT option, and setting the option to None turns the check off.

Classes

InvariantWarning

Warned when a result violates an invariant by more than its tolerance.

Functions

check_symmetric

Density matrices are symmetric, max|D - D^T| within the check_symmetry_atol option. calc_density_matrix checks every component, save_densities every matrix it saves.

check_partition_of_unity

Fragment weights sum to 1 wherever the basis functions reach, within the check_partition_atol option. FragmentOperators.from_grid checks the weights unless given partition_of_unity=False.

check_charge_closure

Fragment populations add up to the total, within the check_closure_atol option. FragmentOperators.population and charge check them against Tr(D S) when the fragment weights are a partition of unity.

exception prop_sapt.checks.InvariantWarning

Bases: RuntimeWarning

Warned when a result violates an invariant by more than its tolerance.

prop_sapt.checks.check_charge_closure(populations: Mapping[str, float | ndarray], total: float | ndarray, atol: float | None = None, where: str = '') → float | ndarray | None

Check that the populations of fragments add up to the total, sum_X N_X = N, and warn for the densities whose populations do not.

For a density matrix D and FragmentOperators the populations are Tr(D W_X) and the total is the electron count Tr(D S); a closure error then means that the quadrature grid or the fragment weights are not accurate enough for that density. Fragment charges and the total charge, e.g. integrated from cube files, work the same way.

Parameters:

• populations (Mapping[str, float | np.ndarray]) – Population of each fragment, a float or an array with one value per density.

• total (float | np.ndarray) – Total population, with the shape of each population.

• atol (float | None, optional) – Largest |sum_X N_X - N| allowed, in electrons. Defaults to the check_closure_atol option; the check is off if that is None.

• where (str, optional) – What the populations belong to, for the message. Defaults to “”.

Raises:

ValueError – If the shapes do not match, or the tolerance is not a number >= 0.

Returns:

The closure error sum_X N_X - N, a float or an

array with one value per density; None if the check is off.

Return type:

float | np.ndarray | None

prop_sapt.checks.check_partition_of_unity(points: ndarray, weights: ndarray, fragment_weights: Mapping[str, ndarray], basisset: BasisSet, atol: float | None = None, where: str = '', chunk_size: int | None = None) → float | None

Check that fragment weights sum to 1 wherever the basis functions reach, and warn if they do not.

The defect is the largest element of

Delta[mu,nu] = sum_q w_q (sum_X w_X(r_q) - 1) phi_mu(r_q) phi_nu(r_q),

the difference between the sum of the fragments’ weighted-overlap matrices and the overlap matrix on the same grid. Points that no basis function reaches do not count, so weights that stop summing to 1 far from the molecule (LISA weights do, tens of bohr away) pass. With nonnegative quadrature weights the defect is at most max|sum_X w_X - 1| times the largest diagonal element of the overlap matrix on the grid (about 1 for normalized basis functions), so weights within the tolerance at every point pass without evaluating any basis function.

Parameters:

• points (np.ndarray) – Points of the quadrature grid in bohr, shape (…, 3).

• weights (np.ndarray) – Quadrature weights, with the shape of points without its last axis.

• fragment_weights (Mapping[str, np.ndarray]) – Weight w_X of each fragment at the points, with the same shape as weights.

• basisset (psi4.core.BasisSet) – Basis set of the functions phi.

• atol (float | None, optional) – Largest defect allowed. Defaults to the check_partition_atol option; the check is off if that is None.

• where (str, optional) – What the weights belong to, for the message. Defaults to “”.

• chunk_size (int | None, optional) – Number of points evaluated at once. Defaults to the point_eval_chunk_size option.

Raises:

ValueError – If the shapes do not match, or the tolerance is not a number >= 0.

Returns:

The defect if it is above the tolerance, else None.

Return type:

float | None

prop_sapt.checks.check_symmetric(matrices: ndarray | Mapping[str, ndarray], atol: float | None = None, where: str = '') → dict[str, float]

Check that density matrices are symmetric, max|D - D^T| <= atol, and warn about those that are not.

The density matrix of a real wavefunction is symmetric. An antisymmetric part does not enter the density or the trace with a symmetric operator, but it corrupts eigendecompositions (numpy.linalg.eigh reads one triangle) and formats that store one triangle, such as fchk files.

Parameters:

• matrices (np.ndarray | Mapping[str, np.ndarray]) – A square matrix, or square matrices by name.

• atol (float | None, optional) – Largest element of |D - D^T| allowed. Defaults to the check_symmetry_atol option; the check is off if that is None.

• where (str, optional) – What the matrices belong to, for the message, e.g. “monomer A”. Defaults to “”.

Raises:

ValueError – If a matrix is not square, or the tolerance is not a number >= 0.

Returns:

max|D - D^T| of each matrix above the tolerance, by name

(“matrix” for a single matrix); empty if all pass or the check is off.

Return type:

dict[str, float]

Fragments

Fragment structure of a dimer geometry: the atoms of monomers A and B and of any extra fragments, with their charges and multiplicities.

Fragments follow the Dimer convention: the first fragment of the geometry (the atoms before the first --) is monomer A, the second is monomer B, and every further fragment, e.g. midbond functions, only adds basis functions to both monomer calculations.

Coordinates are in bohr, in the frame the Dimer uses: the input orientation, moved to the center of mass unless the geometry says no_com.

Classes

Fragment

Atoms, charge and multiplicity of one fragment.

DimerFragments

Atoms of a dimer geometry (elements, ghost flags, coordinates) and their fragments.

Functions

parse_fragments

Fragment structure of a Psi4 geometry string or a .p4m file, without running any calculation.

class prop_sapt.fragments.DimerFragments(symbols: tuple[str, ...], atomic_numbers: tuple[int, ...], real: tuple[bool, ...], coordinates: ndarray, fragments: tuple[Fragment, ...])

Bases: object

Atoms of a dimer geometry and the fragments they belong to.

Build it with parse_fragments from a geometry, or with from_molecule from a Psi4 molecule; Dimer.fragments holds the one of a dimer.

property A: Fragment

Monomer A, the first fragment.

property B: Fragment

Monomer B, the second fragment.

atomic_numbers: tuple[int, ...]

Atomic number of each atom, ghost atoms included.

check_same_system(other: DimerFragments, atol: float = 1e-06)

Check that other describes the same atoms, fragments and coordinates.

Parameters:

• other (DimerFragments) – Fragments to compare with.

• atol (float, optional) – Largest coordinate difference allowed, in bohr. Defaults to 1e-6.

Raises:

ValueError – Naming the first difference found, values of self first.

coordinates: ndarray

Coordinates in bohr, shape (natom, 3).

property extra: tuple[Fragment, ...]

The fragments after the monomers, which carry basis functions only.

fragment_of(atom: int) → str

Label of the fragment holding atom.

Parameters:

atom (int) – Atom index in the dimer geometry, counted from 0.

Returns:

“A”, “B”, or “X1”, “X2”, … for the extra fragments.

Return type:

str

fragments: tuple[Fragment, ...]

Monomer A, monomer B, then the extra fragments, in the order of the geometry.

classmethod from_dict(data: Mapping) → DimerFragments

Build the fragments from a dictionary written by DimerFragments.to_dict.

Parameters:

data (Mapping) – Dictionary from to_dict.

Raises:

ValueError – If an entry is missing or the entries are inconsistent.

Returns:

The fragments.

Return type:

DimerFragments

classmethod from_molecule(molecule: Molecule) → DimerFragments

Read the fragments of a Psi4 molecule, after updating its geometry.

For the frame of a Dimer, fix the orientation and the center of mass of a molecule not yet updated, as parse_fragments does.

Parameters:

molecule (psi4.core.Molecule) – Molecule with at least two fragments.

Raises:

ValueError – If the molecule has fewer than two fragments.

Returns:

The fragments.

Return type:

DimerFragments

property natom: int

Number of atoms, ghost atoms included.

real: tuple[bool, ...]

False for ghost atoms (Gh(X) or @X), which carry basis functions only.

symbols: tuple[str, ...]

Element symbol of each atom, e.g. “He”; ghost atoms keep their element.

to_dict() → dict

Return the fragments as a plain, JSON-serialisable dictionary.

Returns:

Read back with DimerFragments.from_dict; coordinates are nested lists

in bohr, which JSON stores without loss.

Return type:

dict

class prop_sapt.fragments.Fragment(label: str, atoms: tuple[int, ...], charge: int, multiplicity: int)

Bases: object

One fragment of a dimer geometry.

atoms: tuple[int, ...]

Indices of the fragment’s atoms in the dimer geometry, counted from 0.

charge: int

Charge of the fragment.

label: str

“A” and “B” for the monomers, “X1”, “X2”, … for the extra fragments.

multiplicity: int

Spin multiplicity of the fragment.

prop_sapt.fragments.parse_fragments(geometry: str | PathLike) → DimerFragments

Read the fragment structure of a dimer geometry without running any calculation.

The coordinates are those a Dimer built from the same geometry uses.

Parameters:

geometry (str | os.PathLike) – Geometry in Psi4 format, or the path of a file holding one (e.g. a .p4m file). A string without line breaks that names an existing file is read as a path.

Raises:

ValueError – If the geometry has fewer than two fragments.

Returns:

The fragments.

Return type:

DimerFragments

SAPT Driver

SAPT Driver Module This module provides a function to calculate the SAPT energy components for a dimer system.

prop_sapt.sapt_driver.calc_sapt_energy(dimer: Dimer, **kwargs) → Series

Calculate SAPT energy for a given dimer.

This function computes the SAPT energy components for a dimer system: SAPT0 for an RHF reference, SAPT(DFT) for an RKS one. The dispersion is the uncoupled one in both cases. With RHF that is the SAPT0 dispersion. With RKS it is not the coupled dispersion that Psi4’s SAPT(DFT) uses by default, so SAPT(DFT) dispersion and total energies differ from Psi4’s. Coupled dispersion is not implemented. The results are returned as a pandas Series and, unless results=False, saved to a CSV file, “sapt_results.csv” by default.

The SAPT energy components are calculated as follows:

• ELST1: First-order electrostatic energy

• EXCH1: First-order exchange energy

• EXCH1(S^2): First-order exchange energy with S^2 correction

• IND2,R: Second-order induction energy (response)

• EXCH-IND2,R: Exchange-induction energy (response)

• DISP2: Second-order dispersion energy

• EXCH-DISP2: Exchange-dispersion energy

Parameters:

dimer (Dimer) – A dimer system for which to calculate the SAPT energy.

Keyword Arguments:

• results (str | bool) – Path of the CSV file the results are saved to. Defaults to “sapt_results.csv”, which True also gives; False (or None) saves nothing.

• response (bool | str) – Induction amplitudes: True, the coupled (CPSCF) ones; False, the uncoupled ones; “both”. Defaults to True.

Returns:

A pandas Series containing the SAPT energy components.

Return type:

pd.Series

prop_sapt.sapt_driver.print_sapt_summary(results: Series, **kwargs)

Prints a summary of the SAPT energy components.

Parameters:

results (pd.Series) – A pandas Series containing the SAPT energy components.

Finite-Field Properties

prop_sapt.finite_field.calc_ff_dipole(geometry: str, method: str, cp_correction: bool = False, field_strength: float = 0.001, formula: str = 'two_point') → DataFrame

Calculate the finite field dipole moment.

Parameters:

• geometry (str) – Geometry of the molecule in Psi4 format.

• method (str) – Electronic structure method to use (e.g., “HF”, “DFT”).

• prop (str | np.ndarray) – Property to calculate (e.g., “dipole”).

• cp_correction (bool, optional) – Whether to apply counterpoise correction. Defaults to False.

• field_strength (float, optional) – Strength of the external electric field in atomic units. Defaults to 0.001.

• formula (str, optional) – Finite difference formula to use (“two_point” or “four_point”). Defaults to “two_point”.

Returns:

DataFrame containing the finite field dipole moment results.

Return type:

pd.DataFrame

prop_sapt.finite_field.calc_ff_sapt_dipole(geometry: str, reference: str = 'RHF', field_strength: float = 0.001, formula: str = 'two_point', sapt_kwargs: Dict | None = None, **kwargs) → DataFrame

Calculate the finite field SAPT dipole moment.

Parameters:

• geometry (str) – Geometry of the dimer in Psi4 format.

• reference (str, optional) – Reference wavefunction type. Defaults to “RHF”.

• field_strength (float, optional) – Strength of the external electric field. Defaults to 0.001.

• formula (str, optional) – Finite difference formula to use (“two_point” or “four_point”). Defaults to “two_point”.

• sapt_kwargs (Optional[Dict], optional) – Additional arguments for SAPT calculations. Defaults to None.

Returns:

DataFrame containing the results of the finite field SAPT calculation.

Return type:

pd.DataFrame

prop_sapt.finite_field.calculate_ff_dipole_along_axis(geometry: str, axis: str, method: str, cp_correction: bool = False, field_strength: float = 0.001, formula: str = 'two_point') → Series

prop_sapt.finite_field.calculate_ff_sapt_dipole_along_axis(geometry: str, axis: str, reference: str = 'RHF', field_strength: float = 0.001, formula: str = 'two_point', sapt_kwargs: Dict | None = None, **kwargs) → Series

Calculate the finite field SAPT dipole moment along a specified axis.

Parameters:

• geometry (str) – Geometry of the dimer in Psi4 format.

• axis (str) – Axis along which to calculate the dipole moment (e.g., “X”, “Y”, “Z”).

• reference (str, optional) – Reference wavefunction type. Defaults to “RHF”.

• field_strength (float, optional) – Strength of the external electric field. Defaults to 0.001.

• formula (str, optional) – Finite difference formula to use (“two_point” or “four_point”). Defaults to “two_point”.

• sapt_kwargs (Optional[Dict], optional) – Additional arguments for SAPT calculations. Defaults to None.

Returns:

Series containing the results of the finite field SAPT calculation.

Return type:

pd.Series

prop_sapt.finite_field.finite_field_property(geometry: str, method: str, prop: str | ndarray, cp_correction: bool = False, field_strength: float = 0.001, formula: str = 'two_point', **kwargs) → DataFrame

Calculate the interaction-induced property using finite field method.

The energies of method (anything psi4.energy runs, e.g. “scf”, a functional, “mp2” or “ccsd”) are computed in a uniform field along each axis (Psi4’s PERTURB_H) and differentiated, so the result includes the orbital relaxation.

Parameters:

• geometry (str) – Geometry of the molecule in Psi4 format.

• method (str) – Electronic structure method, as named for psi4.energy.

• prop (str | np.ndarray) – Property to calculate (e.g., “dipole”).

• cp_correction (bool, optional) – Differentiate the counterpoise-corrected interaction energy (bsse_type=”cp”), giving the interaction-induced property of a dimer. Without it the total energy is differentiated; Psi4’s field acts on the electrons only, so that gives the electronic part of the dipole, which depends on the origin. Defaults to False.

• field_strength (float, optional) – Strength of the external electric field in atomic units. Defaults to 0.001.

• formula (str, optional) – Finite difference formula to use (“two_point” or “four_point”). Defaults to “two_point”.

Keyword Arguments:

results (str) – Path of the CSV file the results are saved to. Defaults to “results-ff-<method>.csv”.

Raises:

• ValueError – If the property is not implemented.

• ValueError – If the calculation fails.

Returns:

DataFrame containing the results of the finite field property calculation.

Return type:

pd.DataFrame

prop_sapt.finite_field.finite_field_sapt(geometry: str, prop: str | ndarray, reference: str = 'RHF', field_strength: float = 0.001, formula: str = 'two_point', sapt_kwargs: Dict | None = None, **kwargs) → DataFrame

Calculate the interaction-induced property using finite field SAPT.

Parameters:

• geometry (str) – Geometry of the dimer in Psi4 format.

• prop (str | np.ndarray) – Property to calculate (e.g., “dipole”).

• field_strength (float, optional) – Strength of the external electric field. Defaults to 0.001.

• formula (str, optional) – Finite difference formula to use (“two_point” or “four_point”). Defaults to “two_point”.

• sapt_kwargs (Optional[Dict], optional) – Additional arguments for SAPT calculations. Defaults to None.

Raises:

• ValueError – If the property is not implemented.

• ValueError – If the SAPT calculation fails.

Returns:

DataFrame containing the results of the finite field SAPT calculation.

Return type:

pd.DataFrame

Runner

Configurable runner: one configuration file, one geometry, the selected propSAPT tasks.

Used by the prop-sapt run command (see prop_sapt.cli), and from Python:

>>> from prop_sapt import run
>>> results = run("config.yaml", "geos/geo1.p4m")

See prop_sapt.config for the configuration format.

prop_sapt.runner.add_basis_dirs(basis_dirs) → None

Prepends directories to the PSIPATH environment variable, where Psi4 looks for <basis name>.gbs files, skipping those already on it.

Parameters:

basis_dirs (Iterable[str]) – Directories holding basis set files.

prop_sapt.runner.as_run_config(config: RunConfig | str | PathLike | Mapping) → RunConfig

Returns config as a RunConfig.

Parameters:

config (RunConfig | str | os.PathLike | Mapping) – A resolved configuration, the path of a configuration file, or a raw configuration mapping (relative paths in it are resolved against the current working directory).

Returns:

The resolved configuration.

Return type:

RunConfig

prop_sapt.runner.build_dimer(config: RunConfig, geometry: str) → Dimer

Creates the Dimer for geometry as described by config.

Psi4 has to be set up first, see setup_psi4.

Parameters:

• config (RunConfig) – Run configuration.

• geometry (str) – Dimer geometry in Psi4 format.

Returns:

The dimer, with the monomer SCF calculations done.

Return type:

Dimer

prop_sapt.runner.run(config: RunConfig | str | PathLike | Mapping, geometry_file: str | PathLike) → dict

Runs the tasks of a configuration for one dimer geometry.

Every file is written into the configured output directory, and a record of the run (configuration, Psi4 options, versions, timings, the geometry’s fragments) into propsapt_run.json there.

Parameters:

• config (RunConfig | str | os.PathLike | Mapping) – Run configuration, see as_run_config.

• geometry_file (str | os.PathLike) – Path to the dimer geometry in Psi4 format (e.g. a .p4m file).

Returns:

The result of each task that ran, keyed by task name.

Return type:

dict

prop_sapt.runner.setup_psi4(config: RunConfig, output_file: str) → list[str]

Sets up Psi4 for a calculation described by config.

Resets the Psi4 and propSAPT options first, so the calculation depends on config alone. Then sets memory, threads and the output file, the basis set search path, the Psi4 options, the per-element basis overrides and the propSAPT options.

The directories in config.basis_dirs are prepended to the PSIPATH environment variable, which Psi4 searches for <basis name>.gbs files.

Parameters:

• config (RunConfig) – Run configuration.

• output_file (str) – Path of the Psi4 output file.

Returns:

The basis options that were overridden per element.

Return type:

list[str]

Run Configuration

Run configuration for the propSAPT runner (prop-sapt run).

A run is described by a small YAML or JSON file. Every key is optional except tasks:

reference: RHF                # "RHF" or "RKS"
functional: pbe0              # RKS only; defaults to the `default_functional` option
grac: {A: 0.0851, B: 0.0999}  # RKS only, GRAC shifts in Hartree

basis: aug-cc-pvtz
df_basis: aug-cc-pvtz         # stem: DF_BASIS_SCF = <stem>-jkfit, DF_BASIS_SAPT = <stem>-ri
df_basis_scf: null            # full basis names, taking precedence over the stem
df_basis_sapt: null
basis_overrides: basis_overrides.json   # path, or inline {element: {basis_key: basis}}
basis_dirs: [basis]           # extra directories searched for <basis name>.gbs files

memory: 2                     # in GB, or a string such as "2 GB" or "500 MB"
threads: 2

psi4_options: {scf_type: direct, e_convergence: 1.0e-10}
propsapt_options: {r_convergence: 1.0e-10}

output:
  dir: "."                    # "{geometry}" expands to the geometry file stem
  psi4_output: null           # defaults to output_hf.dat / output_propsapt_dft_<functional>.dat

tasks:                        # true for defaults, a mapping to change them
  monomer_references: {cubes: true, fchk: true}
  sapt_energy: {results: sapt_results.csv, response: true}
  property: {prop: dipole, results: property_results.csv, response: true}
  densities: {orbital_basis: AO, save_cubes: true, save_matrices: true,
              matrix_format: npz}

Values are taken from the built-in defaults, then the file, then the command-line flags. A –df-basis flag sets both fitting bases from its stem, replacing df_basis_scf and df_basis_sapt given in the file.

DF_BASIS_SAPT is always set: from df_basis_sapt, else <df_basis>-ri, else <basis>-ri. DF_BASIS_SCF is set from df_basis_scf, else <df_basis>-jkfit; with neither given it is left to Psi4, which picks <basis>-jkfit where that exists and a universal JKFIT basis otherwise (e.g. for He with aug-cc-pvtz). The densities task saves its matrices into one self-describing bundle, densities.npz (matrix_format: npz, read with prop_sapt.load_densities), or as bare .npy files (matrix_format: npy, deprecated). The older compress_matrices: true/false setting still works, as npz/npy, with a deprecation warning.

Relative paths inside the file (basis_overrides, basis_dirs) are resolved against the directory of the file; the output directory is resolved against the current working directory.

Basis sets missing from Psi4’s library can be added as Psi4-format files named <basis name>.gbs in one of the basis_dirs, e.g. written from ccRepo or the Basis Set Exchange with prop-sapt basis fetch (see prop_sapt.basis_fetch).

exception prop_sapt.config.ConfigError

Bases: ValueError

Raised when a run configuration cannot be read or is invalid.

class prop_sapt.config.RunConfig(reference: str = 'RHF', functional: str | None = None, grac_A: float | None = None, grac_B: float | None = None, basis: str = 'aug-cc-pvtz', df_basis_scf: str | None = None, df_basis_sapt: str = 'aug-cc-pvtz-ri', basis_overrides: dict = <factory>, basis_dirs: tuple = (), memory_gb: float = 2.0, threads: int = 1, psi4_options: dict = <factory>, propsapt_options: dict = <factory>, output_dir_template: str = '.', psi4_output: str = 'output_hf.dat', tasks: dict = <factory>, source: str | None = None)

Bases: object

Resolved and validated run configuration. Build it with resolve_config.

applied_basis_overrides() → dict

Returns the per-element basis overrides that apply to this run, those of the basis options it sets (see apply_basis_overrides).

basis: str = 'aug-cc-pvtz'

basis_dirs: tuple = ()

basis_overrides: dict

df_basis_sapt: str = 'aug-cc-pvtz-ri'

df_basis_scf: str | None = None

effective_psi4_options() → dict

Returns every Psi4 option the runner sets: the basis keys, save_jk and psi4_options.

functional: str | None = None

grac_A: float | None = None

grac_B: float | None = None

memory_gb: float = 2.0

property method_tag: str

Short method label used in the default filenames, “hf” or “dft_<functional>”.

output_dir(geometry_stem: str) → str

Returns the output directory for a geometry, with “{geometry}” expanded.

output_dir_template: str = '.'

propsapt_options: dict

psi4_options: dict

psi4_output: str = 'output_hf.dat'

reference: str = 'RHF'

source: str | None = None

tasks: dict

threads: int = 1

to_dict() → dict

Returns the configuration as a plain, JSON-serialisable dictionary.

prop_sapt.config.load_config(path: str | PathLike) → dict

Reads a run configuration file.

Files ending in “.json” are read as JSON, anything else as YAML.

Parameters:

path (str | os.PathLike) – Path to the configuration file.

Raises:

ConfigError – If the file cannot be read or parsed, or does not hold a mapping.

Returns:

The raw configuration, to be passed to resolve_config.

Return type:

dict

prop_sapt.config.resolve_config(raw: Mapping, base_dir: str | PathLike | None = None, overrides: Mapping | None = None, source: str | None = None) → RunConfig

Validates a raw configuration and fills in the defaults.

Parameters:

• raw (Mapping) – Raw configuration, e.g. from load_config.

• base_dir (str | os.PathLike | None, optional) – Directory that relative paths in raw are resolved against, normally the directory of the config file. Defaults to the current working directory.

• overrides (Mapping | None, optional) – Values that take precedence over raw, keyed by the names in OVERRIDE_KEYS. None values are ignored. Defaults to None.

• source (str | None, optional) – Path of the config file, kept for provenance.

Raises:

ConfigError – If the configuration is invalid.

Returns:

The resolved configuration.

Return type:

RunConfig

Basis Set Overrides

Per-element basis set overrides.

Some elements lack the auxiliary basis sets that pair with the orbital basis in Psi4 (e.g. aug-cc-pvtz has no -jkfit / -ri set for Zn). For those elements another basis is substituted, leaving every other element on the default. Overrides are given element-first, as a mapping or as a JSON file holding one:

{
    "Zn": {"DF_BASIS_SCF": "def2-tzvpd-ri", "DF_BASIS_SAPT": "def2-tzvpd-ri"}
}

where each key under an element is a Psi4 basis option (“BASIS”, “DF_BASIS_SCF”, “DF_BASIS_SAPT”, …). An override is applied only for a basis option that the calculation sets, so one file can serve several kinds of calculation: entries for options the calculation does not use (e.g. “DF_BASIS_MP2” in a propSAPT run) are skipped.

prop_sapt.basis_overrides.apply_basis_overrides(options: Mapping, overrides: Mapping[str, Mapping[str, str]]) → list[str]

Registers per-element basis overrides with Psi4.

Call after psi4.set_options(options). For each basis option that is both overridden and set in options, a custom basis is registered with psi4.basis_helper: every element gets options[basis_option] and the overridden elements get their substitutes. psi4.basis_helper then points the basis option at the custom basis.

Parameters:

• options (Mapping) – Psi4 options already passed to psi4.set_options.

• overrides (Mapping[str, Mapping[str, str]]) – Overrides as returned by load_basis_overrides.

Returns:

The basis options that were overridden.

Return type:

list[str]

prop_sapt.basis_overrides.load_basis_overrides(spec: str | PathLike | Mapping | None, base_dir: str | PathLike | None = None) → dict[str, dict[str, str]]

Reads per-element basis overrides and groups them by basis option.

Parameters:

• spec (str | os.PathLike | Mapping | None) – Path to a JSON file, or a mapping {element: {basis_option: basis_name}}. None means no overrides.

• base_dir (str | os.PathLike | None, optional) – Directory that a relative path is resolved against. Defaults to the current working directory.

Raises:

• OSError – If the file cannot be read.

• ValueError – If the file is not valid JSON, or a basis set name cannot be assigned to an element by Psi4 (e.g. a “bse:” name).

• TypeError – If the overrides are not a mapping of mappings of strings.

Returns:

Overrides as {BASIS_OPTION: {element: basis_name}}, with

basis options in upper case.

Return type:

dict[str, dict[str, str]]

Fetching Basis Sets

Fetch basis sets missing from Psi4’s library into Psi4-format files (prop-sapt basis fetch).

The files are written once, e.g. on a workstation or login node, into a directory that runs then list in basis_dirs (see prop_sapt.config), so calculations read fixed inputs and need no network or extra packages. Two sources are available, each through an optional package that is installed separately and only where this command is used:

ccrepo ccRepo, http://www.grant-hill.group.shef.ac.uk/ccrepo/, through ccrepo_utility

(its dependencies are not installed with it, so they are listed as well): pip install git+https://github.com/Sheffield-Theoretical-Chemistry/ccrepo_utility colorlog tqdm requests

bse Basis Set Exchange, https://www.basissetexchange.org, through basis_set_exchange:

pip install basis_set_exchange

Each file is named after its basis set, by default “<source>-<basis set>”, e.g. “ccrepo-aug-cc-pVTZ_MP2FIT” in ccrepo-aug-cc-pvtz_mp2fit.gbs. A file must not take the name of a basis set in Psi4’s library: Psi4 reads the first file of a name it finds, so a file holding a few elements would hide the library basis set for every other element.

exception prop_sapt.basis_fetch.BasisFetchError

Bases: RuntimeError

Raised when a basis set cannot be fetched or written.

prop_sapt.basis_fetch.fetch_basis_file(source: str, name: str, elements: Sequence[str], directory: str | Path = 'basis', as_name: str | None = None, force: bool = False) → tuple[Path, str]

Fetches a basis set and writes it into a Psi4-format file.

Parameters:

• source (str) – Source of the basis set, a key of SOURCES (“ccrepo” or “bse”).

• name (str) – Name of the basis set at the source, e.g. “aug-cc-pVTZ_MP2FIT”.

• elements (Sequence[str]) – Element symbols to fetch the basis set for.

• directory (str | Path, optional) – Directory to write the file into. Defaults to “basis”.

• as_name (str | None, optional) – Name of the basis set in propSAPT/Psi4, which also names the file. Defaults to “<source>-<name>”.

• force (bool, optional) – Whether to overwrite an existing file. Defaults to False.

Raises:

BasisFetchError – If the source or its package is unavailable, the basis set is missing for an element, the name would hide a Psi4 library basis set, or the file exists.

Returns:

Path of the written file and the basis set name to use in a run

configuration.

Return type:

tuple[Path, str]

Command-Line Interface

Command-line interface of propSAPT:

prop-sapt run CONFIG GEOMETRY [GEOMETRY ...] [--memory 32 --threads 8 ...]
prop-sapt basis fetch SOURCE NAME --elements He,Li [--dir basis]

run runs the tasks of a configuration file (see prop_sapt.config) for each geometry in turn; flags override the values in the file. basis fetch writes a basis set missing from Psi4’s library into a file for basis_dirs (see prop_sapt.basis_fetch). Exit status: 0 on success, 1 if a calculation failed, 2 for an invalid configuration, command line or basis set request.

prop_sapt.cli.build_parser() → ArgumentParser

Builds the prop-sapt argument parser.

prop_sapt.cli.main(argv: list[str] | None = None) → int

Entry point of the prop-sapt command.

Parameters:

argv (list[str] | None, optional) – Command-line arguments without the program name. Defaults to sys.argv[1:].

Returns:

Exit status.

Return type:

int

propSAPT Options

Options module for the propSAPT package.

Package-wide default values for keyword arguments that are otherwise hardcoded in propSAPT routines. An explicitly passed keyword argument always takes precedence over the value stored here. Use prop_sapt.set_propsapt_options and prop_sapt.get_propsapt_options to modify and inspect the current values.

Available options:

default_functional (str): DFT functional used for the monomer calculations

with the RKS reference, if functional is not given to Dimer. Defaults to “pbe0”.

r_convergence (float): Residual convergence threshold of the CPSCF solver.

Defaults to 1.0e-10.

cpscf_maxiter (int): Maximum number of CPSCF iterations.

Defaults to 50.

cube_grid_step (float): Step size (in bohr) of the grid built by make_cube,

if no grid is given. Defaults to 0.2.

cube_grid_overage (float): Extent (in bohr) of the grid built by make_cube

beyond the outermost atoms, if no grid is given. Defaults to 4.0.

cube_iso_sum_level (float): Fraction of the density (or of the squared orbital)

enclosed by the isosurface, used for the isovalues reported in cube files. Defaults to 0.85.

point_eval_chunk_size (int): Number of points for which basis function values are

computed at once by the prop_sapt.evaluate functions and make_cube; the memory used is about point_eval_chunk_size * nbf floats. Defaults to 2000.

quadrature_radial_points (int): Number of radial points per atom of the molecular

quadrature grid built by prop_sapt.partition.molecular_grid. Defaults to 75.

quadrature_spherical_points (int): Number of Lebedev points per radial shell of that

grid. Defaults to 590.

quadrature_bs_radius_alpha (float): Scale of the Bragg-Slater radii setting the radial

extent of that grid, Psi4’s DFT_BS_RADIUS_ALPHA. Psi4’s default of 1.0 is too short for products of diffuse basis functions. Defaults to 2.0.

check_symmetry_atol (float | None): Largest element of |D - D^T| allowed for a

density matrix by prop_sapt.checks.check_symmetric, which calc_density_matrix and save_densities run; a larger one warns. None turns the check off. Defaults to 1.0e-10.

check_partition_atol (float | None): Largest defect of fragment weights allowed by

prop_sapt.checks.check_partition_of_unity, which FragmentOperators.from_grid runs: the largest element of sum_q w_q (sum_X w_X - 1) phi_mu phi_nu. A larger one warns. None turns the check off. Defaults to 1.0e-8.

check_closure_atol (float | None): Largest |sum_X Tr(D W_X) - Tr(D S)| in

electrons allowed by prop_sapt.checks.check_charge_closure, which FragmentOperators.population and charge run; a larger one warns. None turns the check off. Defaults to 1.0e-6.

class prop_sapt.options.PropsaptOptions

Bases: object

Container for the package-wide propSAPT options.

Provides read-only, dict-like access (options[key], key in options, iteration over keys, len). Options should be modified only through prop_sapt.set_propsapt_options. The module-level PROPSAPT_OPTIONS instance is the one used throughout the package.

__contains__(key: str) → bool

Checks whether key is a known option.

__getitem__(key: str)

Returns the current value of the option key.

__iter__() → Iterator[str]

Iterates over the option keys.

__len__() → int

Returns the number of options.

get_options_dict() → dict

Returns a copy of the current options.

Modifying the returned dictionary does not change the propSAPT options.

Returns:

Dictionary containing option_key: option_value pairs.

Return type:

dict

prop_sapt.options.check_if_valid_option(option_key: str) → None

Checks whether option_key is a valid propSAPT option key.

Parameters:

option_key (str) – Option key to check (case-insensitive).

Raises:

• TypeError – Raised when a non-string option_key is given.

• KeyError – Raised when invalid option_key is given.

prop_sapt.options.get_propsapt_options() → dict

Returns the current package-wide propSAPT options.

Modifying the returned dictionary does not change the propSAPT options, use prop_sapt.set_propsapt_options instead.

Returns:

Dictionary containing option_key: option_value pairs.

Return type:

dict

prop_sapt.options.reset_propsapt_options(option_keys: Collection[str] | None = None) → None

Resets the package-wide propSAPT options to the default values.

Parameters:

option_keys (Collection[str] | None, optional) – Option keys to be set to default values (e.g. a list, tuple, set or dict), or a single key as a string. If None, all options are being set to default values. Defaults to None.

Raises:

• TypeError – Raised when option_keys is not a collection (e.g. a generator) or a non-string option key is given.

• KeyError – Raised when an invalid option key is given.

• ValueError – Raised when option_keys is empty.

prop_sapt.options.set_propsapt_options(options: dict) → None

Sets package-wide propSAPT options.

The new values are used as defaults by all subsequent propSAPT calls, unless overridden by an explicitly passed keyword argument. Option keys are case-insensitive. Options not present in options keep their current values. See the prop_sapt.options module docstring for the available options.

Example

>>> prop_sapt.set_propsapt_options({"r_convergence": 1.0e-8, "cpscf_maxiter": 100})

Parameters:

options (dict) – Dictionary containing option_key: option_value pairs.

Raises:

• TypeError – Raised when a non-string option_key is found in options.

• KeyError – Raised when invalid option_key is found in options.

Utility Functions

class prop_sapt.utils.CalcTimer(name: str)

Bases: object

Class to time calculations.

__enter__()

Starts the timer.

Returns:

Returns the CalcTimer instance to allow for use in a with statement.

Return type:

CalcTimer

__exit__(*args)

Stops the timer. And prints out the time taken for the operation into the output file.

__init__(name: str)

Initializes timer with a name.

Parameters:

name (str) – Name of the timed operation.

prop_sapt.utils.changed_psi4_options() → dict | None

The Psi4 options changed from their defaults, None if Psi4 cannot list them.

prop_sapt.utils.energy_printer(name: str, value: float, output: str = 'psi4', **kwargs)

Prints out the energy value in mH and kcal/mol.

Parameters:

• name (str) – Name of the energy term.

• value (float) – Value of the energy term in Hartree.

• output (str, optional) – Output method. Defaults to “psi4” - Psi4’s output file.

Raises:

• ValueError – If an unknown output type is provided.

• ValueError – If the output filename is not provided when output is ‘file’.

prop_sapt.utils.package_version(name: str) → str

Returns the installed version of a package, or “unknown”.

Parameters:

name (str) – Distribution name of the package, e.g. “prop_sapt”.

Returns:

The version, or “unknown” if the package is not installed.

Return type:

str

prop_sapt.utils.prepare_path(file_path: str) → str

Create directories (and subdirectories) from file_path if they don’t exist.

Parameters:

file_path (str) – Path to check or create.

Returns:

The file_path given.

Return type:

str

prop_sapt.utils.trace_memory_peak(func: Callable)

Decorator. Trace peak memory usage of ‘func’.

Parameters:

func (Callable) – A function to trace peak memory usage.

Point Evaluation

Evaluation of densities, orbitals and the electrostatic potential at arbitrary points.

All functions take points in bohr, as an array of shape (…, 3), and the Psi4 basis set in which the matrices are expressed (Psi4 AO order, basisset.nbf() functions). The values are returned with the shape of points without its last axis, preceded by one axis over the matrices or orbitals when several are given.

Basis function values are computed for blocks of chunk_size points at a time, which bounds the memory to chunk_size * nbf values; the default is the point_eval_chunk_size option.

Functions

evaluate_basis_functions

Values of all basis functions at the points.

evaluate_density

Density of one or more AO density matrices at the points.

evaluate_orbital

Values of one or more orbitals (columns of MO coefficients) at the points.

evaluate_esp

Electrostatic potential of one or more AO density matrices at the points.

prop_sapt.evaluate.evaluate_basis_functions(points: ndarray, basisset: BasisSet, chunk_size: int | None = None) → ndarray

Evaluate all basis functions of basisset at points.

Parameters:

• points (np.ndarray) – Points in bohr, shape (…, 3).

• basisset (psi4.core.BasisSet) – Basis set to evaluate.

• chunk_size (int | None, optional) – Number of points evaluated at once. Defaults to the point_eval_chunk_size option.

Returns:

Basis function values, shape (…, nbf), in Psi4 AO order.

Return type:

np.ndarray

prop_sapt.evaluate.evaluate_density(dm: ndarray | Sequence[ndarray], points: ndarray, basisset: BasisSet, chunk_size: int | None = None) → ndarray

Evaluate the density of AO density matrix dm at points, rho(r) = sum_{mu,nu} phi_mu(r) dm[mu,nu] phi_nu(r).

No spin factor is applied: for the total density of an RHF wavefunction pass 2 * Da.

Parameters:

• dm (np.ndarray | Sequence[np.ndarray]) – AO density matrix, shape (nbf, nbf), or several of them, as a sequence or an array of shape (M, nbf, nbf).

• points (np.ndarray) – Points in bohr, shape (…, 3).

• basisset (psi4.core.BasisSet) – Basis set of dm.

• chunk_size (int | None, optional) – Number of points evaluated at once. Defaults to the point_eval_chunk_size option.

Returns:

Density at the points, shape (…) for a single matrix

and (M, …) for several.

Return type:

np.ndarray

prop_sapt.evaluate.evaluate_esp(dm: ndarray | Sequence[ndarray], points: ndarray, basisset: BasisSet, nuclei: bool = False) → ndarray

Evaluate the electrostatic potential of AO density matrix dm at points, from analytic potential integrals. Electrons contribute negatively.

As for evaluate_density, no spin factor is applied to dm.

Parameters:

• dm (np.ndarray | Sequence[np.ndarray]) – AO density matrix, shape (nbf, nbf), or several of them, as a sequence or an array of shape (M, nbf, nbf).

• points (np.ndarray) – Points in bohr, shape (…, 3).

• basisset (psi4.core.BasisSet) – Basis set of dm.

• nuclei (bool, optional) – Add the potential of the nuclei of basisset.molecule(), use it for total (not interaction-induced) densities. Ghost atoms carry no charge, and a nucleus lying exactly on a point is skipped at that point. Defaults to False.

Returns:

Electrostatic potential at the points, shape (…) for a single matrix

and (M, …) for several.

Return type:

np.ndarray

prop_sapt.evaluate.evaluate_orbital(coeffs: ndarray, points: ndarray, basisset: BasisSet, chunk_size: int | None = None) → ndarray

Evaluate orbitals given by MO coefficients coeffs at points, psi_i(r) = sum_mu phi_mu(r) coeffs[mu,i].

Orbitals are the columns of coeffs, as in Psi4’s Ca(): pass C[:, i] for orbital i, or C[:, i:j] for several.

Parameters:

• coeffs (np.ndarray) – MO coefficients, shape (nbf,) for one orbital or (nbf, k) for k orbitals.

• points (np.ndarray) – Points in bohr, shape (…, 3).

• basisset (psi4.core.BasisSet) – Basis set of coeffs.

• chunk_size (int | None, optional) – Number of points evaluated at once. Defaults to the point_eval_chunk_size option.

Returns:

Orbital values at the points, shape (…) for coeffs of shape (nbf,)

and (k, …) for coeffs of shape (nbf, k).

Return type:

np.ndarray

Fragment Populations

Fragment population operators: weighted-overlap matrices of a partition of space.

Given a molecular quadrature grid (points r_q, weights w_q) and fragment weights w_X(r_q) at its points, the weighted-overlap matrix of fragment X is

W_X[mu,nu] = sum_q w_q w_X(r_q) phi_mu(r_q) phi_nu(r_q),

computed once per geometry. The number of electrons any AO density matrix D puts on fragment X is then a trace, Tr(D W_X), with no further grid work. The first-moment matrices (r_q times the same weights) give fragment dipoles in the same way.

Points are in bohr. Matrices are in the AO basis of the given Psi4 basis set, in Psi4 AO order, and no spin factor is applied, as in prop_sapt.evaluate: the density matrices of calc_densities follow the Da convention, pass 2 * D for the total (RHF) density.

Fragment weights come from any partitioning: the Becke partition of Psi4’s molecular grid (becke_fragment_weights), Hirshfeld-like weights from the monomer densities (hirshfeld_fragment_weights), or weights computed elsewhere, e.g. LISA or MBIS weights from horton_part on its own grid, passed as arrays with that grid’s points and weights.

Classes

QuadratureGrid

Points, weights and parent atoms of a molecular quadrature grid.

FragmentOperators

Weighted-overlap and first-moment matrices of the fragments, with populations, charges and dipoles of density matrices.

Functions

molecular_grid

Psi4’s Becke-partitioned molecular quadrature grid of a basis set’s molecule.

becke_fragment_weights

Becke fragment weights on a QuadratureGrid, from the atoms of each fragment.

hirshfeld_fragment_weights

Hirshfeld-like fragment weights from fragment densities, e.g. of the monomers.

class prop_sapt.partition.FragmentOperators(overlaps: Mapping[str, ndarray], first_moments: Mapping[str, ndarray] | None = None, ao_overlap: ndarray | None = None)

Bases: object

Weighted-overlap matrices W_X and first-moment matrices R_X of a set of fragments,

W_X[mu,nu] = integral(w_X(r) phi_mu(r) phi_nu(r) dr), R_X[k,mu,nu] = integral(w_X(r) r_k phi_mu(r) phi_nu(r) dr),

with r about the coordinate origin. Build them with from_grid; populations, charges and dipoles of density matrices are then traces.

When the fragments cover the whole system, i.e. their weights are a partition of unity (the default of from_grid), the weights are checked to sum to 1, and every population and charge call checks that the populations of all fragments add up to the electrons Tr(dm S), see closure_error and prop_sapt.checks.

__init__(overlaps: Mapping[str, ndarray], first_moments: Mapping[str, ndarray] | None = None, ao_overlap: ndarray | None = None)

Parameters:

• overlaps (Mapping[str, np.ndarray]) – W_X of each fragment, shape (nbf, nbf).

• first_moments (Mapping[str, np.ndarray] | None, optional) – R_X of each fragment, shape (3, nbf, nbf). Without them dipole is not available. Defaults to None.

• ao_overlap (np.ndarray | None, optional) – The analytic AO overlap matrix S, shape (nbf, nbf). Giving it declares that the fragment weights are a partition of unity, i.e. the fragments cover the whole system: population and charge then check that the populations of all fragments add up to Tr(dm S), and closure_error is available. Defaults to None.

property ao_overlap: ndarray | None

The analytic AO overlap matrix S, read-only, if the fragment weights are a partition of unity; None otherwise.

charge(fragment: str, dm: ndarray | Iterable[ndarray]) → float | ndarray

Electronic charge of AO density matrix dm on fragment, -Tr(dm W_X). For an interaction-induced density this is the charge moved onto the fragment. No spin factor is applied: for the total density of an RHF matrix pass 2 * Da.

Parameters:

• fragment (str) – Name of the fragment.

• dm (np.ndarray | Iterable[np.ndarray]) – AO density matrix, shape (nbf, nbf), or several of them, as a sequence or an array of shape (M, nbf, nbf).

Returns:

Charge, a float for a single matrix

and shape (M,) for several.

Return type:

float | np.ndarray

Warns:

InvariantWarning – As population.

closure_error(dm: ndarray | Iterable[ndarray]) → float | ndarray

Electrons of AO density matrix dm that the fragments miss or count twice, sum_X Tr(dm W_X) - Tr(dm S): the quadrature error of the grid and the defect of the fragment weights, for this density. No spin factor is applied to dm.

Parameters:

dm (np.ndarray | Iterable[np.ndarray]) – AO density matrix, shape (nbf, nbf), or several of them, as a sequence or an array of shape (M, nbf, nbf).

Raises:

ValueError – If the operators have no ao_overlap, i.e. the fragment weights are not a partition of unity.

Returns:

Closure error in electrons, a float for a single matrix

and shape (M,) for several.

Return type:

float | np.ndarray

dipole(fragment: str, dm: ndarray | Iterable[ndarray], origin: ndarray | None = None) → ndarray

Electronic dipole moment of AO density matrix dm on fragment, -integral(w_X(r) (r - origin) rho(r) dr) = -Tr(dm R_X) + origin Tr(dm W_X), in atomic units. Electrons contribute negatively, as in Psi4’s ao_dipole, so at the coordinate origin the dipoles of fragments that cover the whole system sum to the electronic dipole Tr(dm @ ao_dipole). No spin factor is applied to dm.

Parameters:

• fragment (str) – Name of the fragment.

• dm (np.ndarray | Iterable[np.ndarray]) – AO density matrix, shape (nbf, nbf), or several of them, as a sequence or an array of shape (M, nbf, nbf).

• origin (np.ndarray | None, optional) – Origin of the dipole in bohr, shape (3,). Defaults to the coordinate origin.

Returns:

Dipole, shape (3,) for a single matrix and (M, 3) for several.

Return type:

np.ndarray

first_moment(fragment: str) → ndarray

First-moment matrices R_X of fragment, shape (3, nbf, nbf), read-only, with r about the coordinate origin and no charge sign.

property fragments: tuple[str, ...]

Names of the fragments.

classmethod from_grid(points: ndarray, weights: ndarray, fragment_weights: Mapping[str, ndarray], basisset: BasisSet, dipole: bool = True, chunk_size: int | None = None, partition_of_unity: bool = True) → FragmentOperators

Compute the fragment operators by quadrature, in one pass over the points.

With partition_of_unity (the default) the fragment weights are first checked to sum to 1 wherever the basis functions reach (prop_sapt.checks.check_partition_of_unity), and the operators get the analytic overlap matrix, so that every population is checked to close (closure_error). A failed check warns, see prop_sapt.checks.

Parameters:

• points (np.ndarray) – Points of the quadrature grid in bohr, shape (…, 3).

• weights (np.ndarray) – Quadrature weights, with the shape of points without its last axis.

• fragment_weights (Mapping[str, np.ndarray]) – Weight w_X of each fragment at the points, with the same shape as weights.

• basisset (psi4.core.BasisSet) – Basis set of the operators.

• dipole (bool, optional) – Also compute the first-moment matrices, needed for dipole. Defaults to True.

• chunk_size (int | None, optional) – Number of points evaluated at once. Defaults to the point_eval_chunk_size option.

• partition_of_unity (bool, optional) – Whether the fragments together cover the whole system, i.e. fragment_weights sum to 1 at every point. If True, this is checked here, and the operators get the analytic overlap matrix S, so that every population and charge call checks that the populations of all fragments add up to Tr(dm S). Pass False for weights that do not sum to 1, e.g. only some of the fragments or atoms, or a grid covering part of space: nothing is checked then, and closure_error is not available. Defaults to True.

Returns:

Operators of the fragments in the basis of basisset.

Return type:

FragmentOperators

property nbf: int

Number of basis functions of the operators.

overlap(fragment: str) → ndarray

Weighted-overlap matrix W_X of fragment, shape (nbf, nbf), read-only. The sum over fragments that cover the whole system is the overlap matrix, up to quadrature error.

population(fragment: str, dm: ndarray | Iterable[ndarray]) → float | ndarray

Number of electrons of AO density matrix dm on fragment, Tr(dm W_X). No spin factor is applied: for the total density of an RHF matrix pass 2 * Da.

Parameters:

• fragment (str) – Name of the fragment.

• dm (np.ndarray | Iterable[np.ndarray]) – AO density matrix, shape (nbf, nbf), or several of them, as a sequence or an array of shape (M, nbf, nbf).

Returns:

Population, a float for a single matrix

and shape (M,) for several.

Return type:

float | np.ndarray

Warns:

InvariantWarning – For a partition of unity, if the populations of all fragments miss Tr(dm S) by more than the check_closure_atol option, see closure_error.

class prop_sapt.partition.QuadratureGrid(points: ndarray, weights: ndarray, atom_index: ndarray)

Bases: object

Molecular quadrature grid: the integral of f over space is sum(weights * f(points)).

atom_index: ndarray

Index of the atom whose Becke cell holds each point, shape (n,), in the order of the atoms of the grid’s molecule.

points: ndarray

Points in bohr, shape (n, 3).

weights: ndarray

Quadrature weights, shape (n,), with the Becke partition of space between the atoms included.

prop_sapt.partition.becke_fragment_weights(grid: QuadratureGrid, fragment_atoms: Mapping[str, Iterable[int]]) → dict[str, ndarray]

Becke fragment weights: 1 at the points in the Becke cells of the fragment’s atoms, 0 elsewhere. With the grid’s weights they integrate over the fragment’s cells.

Every atom that carries grid points must belong to one fragment, ghost atoms included, so that the weights sum to 1 at every point.

Parameters:

• grid (QuadratureGrid) – Grid from molecular_grid.

• fragment_atoms (Mapping[str, Iterable[int]]) – Atom indices of each fragment, e.g. {“A”: range(0, 3), “B”: range(3, 6)}.

Returns:

Weights of each fragment, shape (n,).

Return type:

dict[str, np.ndarray]

prop_sapt.partition.hirshfeld_fragment_weights(points: ndarray, densities: Mapping[str, ndarray], basisset: BasisSet, chunk_size: int | None = None) → dict[str, ndarray]

Hirshfeld-like fragment weights, w_X(r) = rho_X(r) / sum_Y rho_Y(r), from fragment densities, e.g. the unperturbed monomer densities of a dimer.

Negative density values (rounding) count as 0, and where the total density is below 1e-30 the weights are split equally between the fragments.

Parameters:

• points (np.ndarray) – Points in bohr, shape (…, 3).

• densities (Mapping[str, np.ndarray]) – AO density matrix of each fragment, shape (nbf, nbf), e.g. {“A”: 2 * Da_A, “B”: 2 * Da_B}.

• basisset (psi4.core.BasisSet) – Basis set of the density matrices.

• chunk_size (int | None, optional) – Number of points evaluated at once. Defaults to the point_eval_chunk_size option.

Returns:

Weights of each fragment, with the shape of points

without its last axis.

Return type:

dict[str, np.ndarray]

prop_sapt.partition.molecular_grid(basisset: BasisSet, radial_points: int | None = None, spherical_points: int | None = None, nuclear_scheme: str = 'TREUTLER', bs_radius_alpha: float | None = None) → QuadratureGrid

Build Psi4’s molecular quadrature grid around the atoms of basisset.molecule().

Every atom, ghost atoms included, gets a Treutler radial grid times a Lebedev spherical grid, unpruned, and space is divided between the atoms by nuclear_scheme. The grid options are set here and do not follow the global Psi4 DFT options.

Parameters:

• basisset (psi4.core.BasisSet) – Basis set whose molecule carries the grid.

• radial_points (int | None, optional) – Radial points per atom. Defaults to the quadrature_radial_points option.

• spherical_points (int | None, optional) – Lebedev points per radial shell, one of Psi4’s Lebedev orders. Defaults to the quadrature_spherical_points option.

• nuclear_scheme (str, optional) – Partition of space between the atoms, Psi4’s DFT_NUCLEAR_SCHEME (“TREUTLER”, “BECKE”, “STRATMANN”, …). Defaults to “TREUTLER”.

• bs_radius_alpha (float | None, optional) – Scale of the Bragg-Slater radii setting the radial extent, Psi4’s DFT_BS_RADIUS_ALPHA. Defaults to the quadrature_bs_radius_alpha option.

Returns:

Points, weights and parent atoms of the grid.

Return type:

QuadratureGrid

Cube Files

Cube files related things.

Cube file reference: https://paulbourke.net/dataformats/cube/

class prop_sapt.cubes.Cube(*, origin: ndarray, atoms: list[tuple[int, float, ndarray]], x_vector: ndarray, y_vector: ndarray, z_vector: ndarray, volumetric_data: ndarray, comment1: str = '', comment2: str = '', n_atoms: int | None = None, n_x: int | None = None, n_y: int | None = None, n_z: int | None = None)

Bases: object

Volumetric data on a regular grid, as stored in a .cube file.

All arguments are keyword-only. Coordinates and grid vectors are in bohr. Arrays are converted to float NumPy arrays, without a copy when they already are.

Keyword Arguments:

• origin (np.ndarray) – Origin of the grid, shape (3,).

• atoms (list) – List of tuples describing molecule. Should look like: (atomic_number, charge, xyz_coordinates).

• x_vector (np.ndarray) – Grid step in x direction, shape (3,).

• y_vector (np.ndarray) – Grid step in y direction, shape (3,).

• z_vector (np.ndarray) – Grid step in z direction, shape (3,).

• volumetric_data (np.ndarray) – Values at the grid points, shape (n_x, n_y, n_z).

• comment1 (str, optional) – First comment line of .cube file. Defaults to “”.

• comment2 (str, optional) – Second comment line of .cube file. Defaults to “”.

• n_atoms (int | None, optional) – Number of atoms. Defaults to len(atoms).

• n_x (int | None, optional) – Number of grid points in x direction. Defaults to volumetric_data.shape[0].

• n_y (int | None, optional) – Number of grid points in y direction. Defaults to volumetric_data.shape[1].

• n_z (int | None, optional) – Number of grid points in z direction. Defaults to volumetric_data.shape[2].

Cubes on the same grid can be added and subtracted, and a cube multiplied or divided by a scalar: rho - rho_pol, 0.5 * rho, -rho. Grids are checked as in add_cubes. Results take the grid and atoms of the left cube, share no arrays with the operands and have empty comments.

Raises:

ValueError – If an array has the wrong shape, or a given count does not match atoms or volumetric_data.

atoms: list[tuple[int, float, ndarray]]

comment1: str = ''

comment2: str = ''

copy() → Cube

Create a copy of the Cube object.

Returns:

Copied cube.

Return type:

Cube

classmethod from_file(filename: str | PathLike) → Cube

Read cube data from filename.

Parameters:

filename (str | os.PathLike) – Path to file.

Returns:

Loaded cube.

Return type:

Cube

integrate(weight: Cube | ndarray | None = None) → float | ndarray

Integrate the volumetric data over the grid: the sum over voxels times voxel_volume, optionally with a weight at every voxel.

For a density cube, cube.integrate() is the number of electrons in the box and cube.integrate(w_A) the electrons that the weight w_A assigns to fragment A. The L2 norm of the data is np.sqrt(cube.integrate(cube.volumetric_data)). Integrals on a cube carry the error of its voxel grid; for fragment populations from density matrices see prop_sapt.partition.FragmentOperators.

Parameters:

weight (Cube | np.ndarray | None, optional) – Weight at every voxel, either a Cube on the same grid or an array whose last three axes are (n_x, n_y, n_z). Leading axes of an array stack several weights, which are integrated at once. Defaults to None, no weight.

Raises:

ValueError – If weight is a cube on another grid or an array of another shape.

Returns:

The integral, or for stacked weights an array of their

leading shape.

Return type:

float | np.ndarray

moved_charge() → float

Calculate the total charge moved by the density described by the cube.

Returns:

Total charge moved.

Return type:

float

n_atoms: int | None = None

n_x: int | None = None

n_y: int | None = None

n_z: int | None = None

origin: ndarray

save(filename: str | PathLike)

Save a Cube object into filename.

Parameters:

filename (str | os.PathLike) – File path.

volumetric_data: ndarray

property voxel_volume: float

Volume of one voxel in bohr^3, the triple product of the grid vectors, |(x_vector x y_vector) . z_vector|.

x_vector: ndarray

y_vector: ndarray

z_vector: ndarray

prop_sapt.cubes.add_cubes(cubes: list[Cube] | tuple[Cube, ...]) → Cube

Calculate a sum of volumetric data of multiple cubes, cube_1 + cube_2 + ....

Cube grids have to be the same for this operation. Data about molecule geometry is taken from the first cube; the result shares no arrays with the operands and has empty comments (a single cube is copied as it is).

Parameters:

cubes (list[Cube] | tuple[Cube, ...]) – List or tuple of cubes to add together.

Raises:

ValueError – If no cubes are given, or the cube grids differ.

Returns:

Resulting Cube.

Return type:

Cube

prop_sapt.cubes.calculate_isocontour(volumetric_data: ndarray | Cube, threshold: float | None = None, obj_type: str = 'density') → tuple[float, float]

Calculate isocontour values for a given threshold, assumed as the density fraction.

Parameters:

• volumetric_data (np.ndarray | Cube) – Data for isovalue calculations.

• threshold (float | None, optional) – Fraction of the density to be inside of the isosurface described by isovalues. Defaults to the cube_iso_sum_level option.

• obj_type (str, optional) – Specifies the type of volumetric data, either “density” or “orbital”. Defaults to “density”.

Returns:

Isovalues tuple.

Return type:

tuple[float, float]

prop_sapt.cubes.cosine_similarity(cube_1: Cube, cube_2: Cube) → float

Calculate cosine similarity between two cubes on the same grid.

The cosine similarity is calculated as:

cosine_similarity = (A . B) / (||A|| ||B||)

utilizing the L2 norm.

Parameters:

• cube_1 (Cube) – First cube.

• cube_2 (Cube) – Second cube.

Raises:

ValueError – If the cube grids differ.

Returns:

Cosine similarity between the two cubes.

Return type:

float

prop_sapt.cubes.make_cube(mol: Molecule, matrix: ndarray | Sequence[ndarray], obj_type: str | Iterable[str] = 'density', **kwargs) → Cube | list[Cube]

Create a Cube object with volumetric data of the given matrix calculated on a grid for a psi4 Molecule object. If iterables with matrices and types are given as matrix and obj_type returns a list with a Cube for each element of the list.

Parameters:

• mol (psi4.core.Molecule) – A corresponding molecule.

• matrix (np.ndarray | Sequence[np.ndarray]) – Either density matrix, orbital coefficients vector (shape (nbf,) or one column (nbf, 1) of Psi4’s Ca()), or list of such objects.

• obj_type (str | Iterable[str], optional) – Specifies the matrix argument type(s), must be either “density”, “orbital”, “esp” or a list with such strings. Defaults to “density”. For “esp” the matrix is a density matrix and the cube holds the electrostatic potential it generates, from analytic potential integrals (electrons negative).

• **kwargs – basisset (psi4.core.BasisSet) is the basis set the matrices are in. Without it, one is built for mol from the current BASIS option; pass it for matrices of a Dimer (or use Dimer.make_cube, which does). grid (dict) is a grid from prepare_grid, reused as given; otherwise one is built around mol with grid_step and grid_overage (bohr), which default to the cube_grid_step and cube_grid_overage options. esp_nuclei (bool) adds the potential of the nuclear charges of mol to every “esp” cube, use it for total (not interaction-induced) densities; a nucleus lying exactly on a grid point is skipped at that point. Defaults to False. chunk_size (int) is the number of grid points for which basis function values are computed at once. Defaults to the point_eval_chunk_size option. iso_sum_level (float) is the fraction of the density (or squared orbital) enclosed by the isovalues written into the cube comments. Defaults to the cube_iso_sum_level option. Values at the grid points come from prop_sapt.evaluate.

Returns:

Resulting cube(s).

Return type:

Cube | list[Cube]

prop_sapt.cubes.prepare_grid(geometry: ndarray, grid_step: float | tuple, grid_overage: float | tuple) → dict

Prepares a simple scalar grid in 3D space.

Parameters:

• geometry (np.ndarray) – Molecule geometry. Should have shape (N_atom, 3).

• grid_step (float | tuple) – Step size of the grid (in bohr). If float then the same step is applied in all directions. If tuple, should have length 3, specifies step size for x, y and z directions separately.

• grid_overage (float | tuple) – Overage of the grid (in bohr). If float then the same overage in all directions. If tuple, should have length 3, specifies grid overage for x, y and z directions separately.

Returns:

Grid info.

Return type:

dict

prop_sapt.cubes.read_cube_file(filename: str | PathLike) → Cube

Reads the data from .cube file filename.

The volumetric data may be laid out with any number of values per line. Cubes of several orbitals (a negative number of atoms), with several values per voxel, or in angstrom (negative numbers of points) are not supported.

Parameters:

filename (str | os.PathLike) – Path to file.

Raises:

• NotImplementedError – For the cube variants that are not supported.

• ValueError – If the file holds other than n_x * n_y * n_z values.

Returns:

Loaded cube.

Return type:

Cube

prop_sapt.cubes.save_cube_file(cube: Cube, filename: str | PathLike)

Saves a Cube object into the cube file: filename.

Parameters:

• cube (Cube) – Cube object to be saved.

• filename (str | os.PathLike) – File path.

Raises:

ValueError – If the counts of cube do not match its data, e.g. after volumetric_data was replaced by an array of another shape.

prop_sapt.cubes.subtract_cubes(cube_1: Cube, cube_2: Cube) → Cube

Calculate a difference between volumetric data of two cubes, cube_1 - cube_2.

Cube grids have to be the same for this operation. Data about molecule geometry is taken from cube_1; the result shares no arrays with the operands and has empty comments.

Parameters:

• cube_1 (Cube) – Cube serving as minuend.

• cube_2 (Cube) – Cube serving as subtrahend.

Raises:

ValueError – If the cube grids differ.

Returns:

Resulting Cube.

Return type:

Cube

prop_sapt.cubes.translate_cube(cube: Cube, vector: ndarray, inplace: bool = False) → Cube

Translates a given Cube in 3D space by a given vector.

Parameters:

• cube (Cube) – A Cube object to be translated.

• vector (np.ndarray) – An array of length 3, a translation vector.

• inplace (bool) – Make an inplace modification if True, else returns a copy. Defaults to False.

Returns:

Resulting Cube object.

Return type:

Cube

Forces

Hellmann-Feynman forces on nuclei from one-particle electron density matrices.

Functions

calc_force_at_nuc

Total force on a nucleus: electronic contribution via AO electric field integrals plus (optionally) nuclear-nuclear repulsion. Accepts a single density matrix or a sequence of matrices, computing the integrals only once.

calc_force_nuc_nuc

Purely nuclear-nuclear repulsion force on a given nucleus (no electron density needed).

prop_sapt.forces.calc_force_at_nuc(mol: Molecule, rho: ndarray | Sequence[ndarray], atidx: int, add_nuclear: bool = True, **kwargs) → ndarray | list[ndarray]

Computes the Hellmann-Feynman force on nucleus atidx from one or more electron density matrices.

The electronic contribution is evaluated analytically via the electric field integrals at the nuclear position: F_i = Z_A * Tr[rho * <mu|(r_i - R_Ai)/|r - R_A|^3|nu>]. The nuclear-nuclear repulsion is added by default.

Accepts a single density matrix or a sequence of density matrices. When a sequence is provided, the electric field integrals are computed only once and reused for all matrices.

Parameters:

• mol (psi4.core.Molecule) – Psi4 molecule defining the geometry and nuclear charges.

• rho (np.ndarray | Sequence[np.ndarray]) – One-particle density matrix in the AO basis, shape (nbf, nbf), or a sequence of such matrices.

• atidx (int) – Zero-based index of the nucleus on which the force is computed.

• add_nuclear (bool, optional) – If True, adds the nuclear-nuclear repulsion force. Set to False to obtain the electronic contribution only. Defaults to True.

• **kwargs – Accepts basisset (psi4.core.BasisSet) to reuse a pre-built basis set instead of constructing one from the molecule.

Raises:

• ValueError – If rho is an empty sequence.

• TypeError – If mol is not a psi4.core.Molecule.

Returns:

Force vector of shape (3,) in atomic units (Hartree/Bohr)

when a single density matrix is given; a list of such vectors when a sequence is given.

Return type:

np.ndarray | list[np.ndarray]

prop_sapt.forces.calc_force_nuc_nuc(mol: Molecule, atidx: int) → ndarray

Computes the nuclear-nuclear repulsion force on nucleus atidx from all other nuclei.

Uses Coulomb’s law: F_a = sum_{b != a} Z_a * Z_b * (R_a - R_b) / |R_a - R_b|^3.

Parameters:

• mol (psi4.core.Molecule) – Psi4 molecule defining the geometry and nuclear charges.

• atidx (int) – Zero-based index of the nucleus on which the force is computed.

Returns:

Nuclear-nuclear force vector of shape (3,) in atomic units (Hartree/Bohr).

Return type:

np.ndarray