20.4. Expectation-value Multipole Moments

If you use this module, please cite [chakraborty2025expectation].

PyBEST provides a computationally cheaper alternative to calculate the multipole moments, such as the dipole moment and quadrupole moment, by using the expectation-value approach. This method exploits an alternate excitation operator \(\hat{S}\) defined as

\[e^{\hat{S}} |\Phi_0\rangle = \frac{e^{\hat{T}^{\dagger}} e^{\hat{T}} |\Phi_0\rangle} {\langle \Phi_0 | e^{\hat{T}^{\dagger}} e^{\hat{T}} | \Phi_0 \rangle}\]
\[\langle \hat{O} \rangle_{\mathrm{XCC}} = \langle \Phi_0 | e^{S^\dagger} e^{-T} \hat{O} e^{T} e^{-S^\dagger} | \Phi_0 \rangle\]

Based on the order-by-order expansion of \(\hat{S}\), and using appropriate truncation scheme, an approximated form of the 1-RDM can be obtained directly from the cluster amplitudes without the need for solving the \(\Lambda\)-equations. The different blocks of the 1-RDM, then, can be expressed as

\[ \begin{align}\begin{aligned}\begin{split}\gamma_i^a = t_i^a + \sum_{jb} t_j^b \left(2 t_{ij}^{ab} - t_{ij}^{ba}\right)\qquad \text{[Occ--Virt]} \\\end{split}\\\begin{split}\gamma_k^i = \delta_i^k - \sum_{jab}t_{ij}^{ab}\left(2 t_{kj}^{ab} - t_{kj}^{ba}\right) \qquad \text{[Occ--Occ]} \\\end{split}\\\begin{split}\gamma_c^a =\sum_{ijb} t_{ij}^{ab} \left(2 t_{ij}^{cb} - t_{ij}^{bc}\right) \qquad \text{[Virt--Virt]}\\\end{split}\end{aligned}\end{align} \]

The expectation-value property module can be called as

expectation_multipole = ExpectationMultipole(lf, occ_model)

where

lf:

A LinalgFactory instance (see Preliminaries).

occ_model:

An Aufbau occupation model of the AufbauOccModel class.

Thereafter, the available multipole moments can be calculated as

multipole = expectation_multipole(
            method_output,
            property_options={"moment": "dipole", "quadrupole"}
            )

The method_output is the output from the electronic structure calculation as IOData object. The available options for the property_options dictionary are listed below. The moment key is mandatory, while the other keys are optional and have default values.

moment:

(str, list, tuple) the type of multipole moment to calculate, “dipole” or “quadrupole” or both as a tuple or list of strings.

molecular_orbitals:

(boolean) if True, the 1-RDM is transformed back to the atomic orbital basis (default True). This option is important for post-Hartree-Fock methods where the 1-RDM is stored in the molecular orbital basis. For the Hartree-Fock method, the 1-RDM is already stored in the atomic orbital basis, so this option should be set to False.

scale:

(float) the scaling factor for the 1-RDM (default 2.0). All internal 1-RDMs in PyBEST are stored for one set of orbitals. Thus, even in the restricted case, only one spin-component is saved. The total 1-RDM (spin-free) is obtained by multiplying by the scaling factor.

origin:

(list, tuple) the coordinates of the reference point for the multipole moment calculation (in a.u.). If not provided, the center of nuclear charge (COC) is used as the reference point.

The current version of PyBEST supports expectation-value 1-RDMs for the following post-Hartree-Fock methods:

  • pCCD

  • LCC (pCCD-LCCD, pCCD-LCCSD)

  • fpCC (fpCCD, fpCCSD)

  • CCD

  • CCSD

20.4.1. Example calculation

Here we show an example of calculating the dipole and quadrupole moments for the water molecule at the pCCD and various post-pCCD levels of theory.

Listing 20.5 data/examples/properties/water_expectation_multipole.py
from pybest import context
from pybest.cc import RCCSD, RfpCCSD, RpCCDLCCSD
from pybest.gbasis import (
    compute_eri,
    compute_kinetic,
    compute_nuclear,
    compute_nuclear_repulsion,
    compute_overlap,
    get_gobasis,
)
from pybest.geminals import RpCCD
from pybest.linalg import DenseLinalgFactory
from pybest.occ_model import AufbauOccModel
from pybest.properties import ExpectationMultipole
from pybest.wrappers import RHF

# Define coordinate file
fn_xyz = context.get_fn("test/water.xyz")

# Create a Gaussian basis set
basis = get_gobasis("cc-pvdz", fn_xyz)

# Create a linalg factory
lf = DenseLinalgFactory(basis.nbasis)

# Compute integrals in the atom-centered Gaussian basis set
kin = compute_kinetic(basis)
ne = compute_nuclear(basis)
eri = compute_eri(basis)
nuc = compute_nuclear_repulsion(basis)

# Compute overlap matrix of atom-centered basis set
olp = compute_overlap(basis)

# Create orbitals that store the AO/MO coefficients
orb_a = lf.create_orbital()
occ_model = AufbauOccModel(basis)

# Converge RHF
hf = RHF(lf, occ_model)
hf_output = hf(kin, ne, eri, nuc, olp, orb_a)

# Call multipole moment property class for expectation-value density matrix-based calculations.
expectation_multipole = ExpectationMultipole(lf, occ_model)

# Compute expectation-value multipole moment.
# Note for users:
# 1. The "moment" option can be set to either "dipole", "quadrupole", or both as a tuple
# or list of strings.
# 2. The "origin" option specifies the origin for the dipole moment calculation.
# If not provided, the center of nuclear charge of the chemical system will be used as
# the default origin.
# 3. Post-HF density matrices are saved in the MO basis, so the "molecular_orbitals"
# option is set to `True` by default. If the user wants to use AO basis density matrices for some reasons,
# they can set this option to `False` in property_options.

# Compute expectation-value multipole moment using pCCD wave function.
pccd = RpCCD(lf, occ_model)
pccd_output = pccd(kin, ne, eri, hf_output)
pccd_multipole = expectation_multipole(
    pccd_output,
    property_options={
        "moment": ("dipole", "quadrupole"),
    },
)
# Compute expectation-value multipole moment using CCSD wave function.
ccsd = RCCSD(lf, occ_model)
ccsd_output = ccsd(kin, ne, eri, hf_output)
ccsd_multipole = expectation_multipole(
    ccsd_output,
    property_options={
        "moment": ("dipole", "quadrupole"),
    },
)
# Compute expectation-value multipole moment using fpCCSD wave function.
fpccsd = RfpCCSD(lf, occ_model)
fpccsd_output = fpccsd(kin, ne, eri, pccd_output)
fpccsd_multipole = expectation_multipole(
    fpccsd_output,
    property_options={
        "moment": ("dipole", "quadrupole"),
    },
)
# Compute expectation-value multipole moment using pCCD-LCCSD wave function.
pccd_lccsd = RpCCDLCCSD(lf, occ_model)
pccd_lccsd_output = pccd_lccsd(kin, ne, eri, pccd_output)
pccd_lccsd_multipole = expectation_multipole(
    pccd_lccsd_output,
    property_options={
        "moment": ("dipole", "quadrupole"),
    },
)

Note

1. The same procedure can be followed for the optimized pCCD natural orbital basis. In that case, the following changes should be made in the above code:

from pybest.geminals import ROOpCCD

# Converge pCCD in the optimized natural orbital basis
pccd_opt = ROOpCCD(lf, occ_model)
pccd_output = pccd_opt(kin, ne, eri, hf_output)

2. The origin can be changed by providing the coordinates (in a.u.) of the reference point in the property_options dictionary. For example, to use the center of mass (COM) as the reference point, the following code can be used:

from pybest.utils import calculate_com

# Calculate the center of mass (COM) coordinates
com_coordinates = calculate_com(molecule)

# Calculate the multipole moments with COM as the reference point
multipole = response_multipole(
            method_output,
            property_options={"moment": "dipole", "quadrupole", "origin": com_coordinates}
            )