9. The Symmetry Adapted Perturbation Theory module
Calculations of interaction energies in quantum chemistry are particularly challenging, primarily because of the strong cancellation between the total energies of the monomers and the dimer, which leads to a loss of accuracy. The main philosophy of SAPT is first, to perform separate calculations for unperturbed monomers, in order to obtain wave functions, density matrices, or orbitals, which are then used to calculate the interaction energy directly, as we describe below.
9.1. Supported features
The SAPT module implements the SAPT0 un-coupled flavour for spin-restricted HF orbitals. Both DenseLinalgFactory and CholeskyLinalgFactory are supported.
All implemented expressions are based on the MO formulation with the \(S^2\) approximation for the exchange terms.
9.2. Theory background
In Symmetry Adapted Perturbation Theory, it is assumed that the whole system consists of two subsystems that interact weakly, i.e. without covalent bonds involved. In particular, the Hamiltonian of such a system can be written as
where the first two terms denote the Hamiltonians of the monomers A and B, respectively, and \(V\) is the interaction between them: the electrostatic repulsion between the electrons and nuclei of the interacting species, and the attraction of the electrons of A with the nuclei of B and vice versa. Such a partitioning of the Hamiltonian defines the perturbation series in powers of \(V\), with a zeroth-order wave function of the form of \(\Psi_{\rm A} \Psi_{\rm B}\), where \(\Psi_{\rm A}\) and \(\Psi_{\rm B}\) are eigenstates of A and B, respectively. Since the permutational symmetry of the zeroth-order wave function is lower than the symmetry of the full Hamiltonian, a symmetry forcing procedure is also needed. In the simplest formulation (SRS - symmetrized Rayleigh-Schrödinger perturbation series), the antisymmetrizer is used to enforce the correct symmetry. For more details, see reviews [jeziorski1994] and [patkowski2020].
PyBEST uses the so-called SAPT0 approximation, where the wave functions of the monomers used for perturbation theory are Slater determinants, and in eq. (9.1) Fock operators are used instead of \(H_{\rm A}\) and \(H_{\rm B}\).
The perturbation series in SAPT in the present formulation is truncated at the second-order and the interaction energy is defined as the sum of the following contributions:
These contributions can be attributed to physically meaningful effects that can be later exploited by the user. The first-order term \(E^{(1)}_{\rm elst}\) is the energy of the electrostatic interaction of the subsystems. The so-called exchange term \(E^{(1)}_{\rm exch}\) is responsible for the leading repulsion between the interacting species due to Pauli repulsion. The induction term \(E^{(2)}_{\rm ind}\) represents the leading interaction between the mutually polarized molecules A and B. Finally, the dispersion energy can be viewed as the interaction of instantaneous dipole moments caused by the correlated motion of the electrons around the molecules. In PyBEST, you can access all these terms as well as the interaction energy itself.
9.3. Quick Guide: running SAPT0(RHF)
The SAPT module provides a sapt_utils wrapper that allows for computing
DCBS counter-poised corrected RHF.
The list of configurable parameters is as follows:
- basis:
String representation of desired basis set
- fn_geo:
Path to dimer geometry in .xyz format, the geometry of the monomers is assumed to be located at {dimer}a.xyz and {dimer}b.xyz
- occupations:
A tuple of fragments occupation given in the (monA, monB, dimer) order
- fourindex_type:
A string specifying the type of the linalg factory; Dense or Cholesky available
- er_threshold:
A float specifying the accuracy of CD-ERI if fourindex_type is Cholesky
- solver:
A string specifying the type of the SCF solver, default to “ediis2”
dimer, monA, monB = sapt_utils.prepare_cp_hf(
basis="cc-pvdz",
fn_geo="dimer_xyz_file_path",
occupations=(9, 5, 14),
fourindex_type="dense",
solver="ediis2"
)
prepare_cp_hf returns IOData containers suitable to be used in the SAPT0 solvers
Note
Only expert users should modify the matrix elements stored in the containers, e.g. modifying one-electron potentials of the monomers, to customize SAPT0 behaviour.
Invoking SAPT0(RHF) is then straightforward
sapt0 = SAPT0(monA, monB)
sapt0(monA, monB, dimer)
9.4. Example Python scripts
Some complete examples can be found in the directory
data/examples/sapt.
9.4.1. SAPT0(RHF) using dense ERI
This is a basic example on how to perform a SAPT0(RHF) calculation in PyBEST using DenseLinalgFactory.
from pybest import context
from pybest.sapt import sapt_utils, SAPT0
###############################################################################
## Perform DCBS HF calculations ###############################################
###############################################################################
# get the XYZ file from PyBEST's test data directory
dimer_xyz = context.get_fn("test/water_dimer.xyz")
# monA is 'test/water_dimera.xyz'
# monB is 'test/water_dimerb.xyz'
basis_set = "cc-pvdz"
# occupation model specific kwargs (e.g., charge); not required here
occupation = ({}, {}, {})
# call wrapper routine that returns fragments IOData objects
dimer, monA, monB = sapt_utils.prepare_cp_hf(
"cc-pvdz",
dimer_xyz,
occupation,
fourindex_type="dense",
)
###############################################################################
## Build SAPT0 solver and run computation #####################################
###############################################################################
# construct solver
sapt0_solver = SAPT0(monA, monB)
# call solver routine
sapt0_solver(monA, monB, dimer)
# dictionary with corrections values
corrections = sapt0_solver.result
# access value of each correction
e10_elst = corrections["E^{(10)}_{elst}"]
e10_exch = corrections["E^{(10)}_{exch}(S^2)"]
e20_ind = corrections["E^{(20)}_{ind},unc"]
e20_exch_ind = corrections["E^{(20)}_{exch-ind}(S^2),unc"]
e20_disp = corrections["E^{(20)}_{disp},unc"]
e20_exch_disp = corrections["E^{(20)}_{exch-disp}(S^2),unc"]
9.4.2. SAPT0(RHF) using Cholesky ERI
This is a basic example on how to perform a SAPT0(RHF) calculation in PyBEST using CholeskyLinalgFactory.
from pybest import context
from pybest.sapt import sapt_utils, SAPT0
###############################################################################
## Perform DCBS HF calculations ###############################################
###############################################################################
# get the XYZ file from PyBEST's test data directory
dimer_xyz = context.get_fn("test/water_dimer.xyz")
# monA is 'test/water_dimera.xyz'
# monB is 'test/water_dimerb.xyz'
basis_set = "cc-pvdz"
# occupation model specific kwargs (e.g., charge); not required here
occupation = ({}, {}, {})
# call wrapper routine that returns fragments IOData objects
dimer, monA, monB = sapt_utils.prepare_cp_hf(
"cc-pvdz",
dimer_xyz,
occupation,
# use cd-eri and set threshold to 1e-6
fourindex_type="cholesky",
er_threshold=1e-6,
)
###############################################################################
## Build SAPT0 solver and run computation #####################################
###############################################################################
# construct solver
sapt0_solver = SAPT0(monA, monB)
# call solver routine
sapt0_solver(monA, monB, dimer)
# dictionary with corrections values
corrections = sapt0_solver.result
# access value of each correction
e10_elst = corrections["E^{(10)}_{elst}"]
e10_exch = corrections["E^{(10)}_{exch}(S^2)"]
e20_ind = corrections["E^{(20)}_{ind},unc"]
e20_exch_ind = corrections["E^{(20)}_{exch-ind}(S^2),unc"]
e20_disp = corrections["E^{(20)}_{disp},unc"]
e20_exch_disp = corrections["E^{(20)}_{exch-disp}(S^2),unc"]