20.3. Response Multipole Moments

PyBEST provides pathways to calculate various single-electron properties, such as the dipole moment and quadrupole moment, by solving the response \(\Lambda\)-equations at pCCD and various post-pCCD levels of theory. If you use this module, please cite [pccd-dipole-moments-jctc-2024]. For response properties with the pCCD-LCC methods, please also cite [nowak2021].

The response property module can be called as

response_multipole = ResponseMultipole(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 = response_multipole(
            method_output,
            property_options={"moment": "dipole", "quadrupole"}
            )

The method_output is the output from the electronic structure calculation as IOData object, which contains the 1-RDM. 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 response 1-RDMs for the following post-Hartree-Fock methods:

• MP2 (relaxed and unrelaxed)

• pCCD

• LCC (LCCD, LCCSD, pCCD-LCCD, pCCD-LCCSD)

Note

For the LCC methods, the response 1-RDMs are only available when the \(\Lambda\)-equations are solved by setting the keyword argument lambda_equations=True in the LCC function call.

lcc=LCC(lf, occ_model)
lcc_result = lcc(kin, ne, eri, reference_output, lambda_equations=True)

where

kin:

the kinetic energy integrals.

ne:

the nucleus-electron attraction integrals.

eri:

the two-electron repulsion integrals.

20.3.1. Example calculation

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

Listing 20.3 data/examples/properties/water_response_multipole_pccd.py

from pybest import context
from pybest.cc import 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 ResponseMultipole
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)

# Call multipole moment property class for response density matrix-based calculations.
response_multipole = ResponseMultipole(lf, occ_model)

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


# Compute multipole moment.
# Call multipole moment property class for response density matrix-based calculations.
response_multipole = ResponseMultipole(lf, occ_model)

# 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 multipole 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,
# they can set this option to `False` in property_options.

hf_multipole = response_multipole(
    hf_output,
    property_options={
        "moment": ("dipole", "quadrupole"),
        "molecular_orbitals": False,
    },
)

pccd = RpCCD(lf, occ_model)
pccd_output = pccd(kin, ne, eri, hf_output)
pccd_multipole = response_multipole(
    pccd_output,
    property_options={
        "moment": ("dipole", "quadrupole"),
    },
)

# Note for users:

# For LCC methods, response Lambda equations are not solved by default.
# In order to access properties in this case, `lambda_equations=True` needs to be set specifically.

pccd_lccsd = RpCCDLCCSD(lf, occ_model)
pccd_lccsd_output = pccd_lccsd(
    kin, ne, eri, pccd_output, lambda_equations=True
)
pccd_lccsd_multipole = response_multipole(
    pccd_lccsd_output,
    property_options={
        "moment": ("dipole", "quadrupole"),
    },
)

Octupole moments can also be calculated by setting the moment key in the property_options dictionary to “octupole”. For example, the following code calculates the octupole moment at the HF and pCCD level of theory.

Listing 20.4 data/examples/properties/water_response_multipole_octupole.py

response_multipole = ResponseMultipole(lf, occ_model)

# Converge RHF
# ------------
hf = RHF(lf, occ_model)
hf_output = hf(kin, ne, eri, nuc, olp, orb_a)
hf_octupole = response_multipole(
    hf_output,
    property_options={
        "moment": "octupole",
        "molecular_orbitals": False,
    },
)

pccd = RpCCD(lf, occ_model)
pccd_output = pccd(kin, ne, eri, hf_output)
pccd_octupole = response_multipole(
    pccd_output,
    property_options={
        "moment": "octupole",
    },
)

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 get_com

# Calculate the center of mass (COM) coordinates
com_coordinates = get_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}
            )