16.4. Computing the quadrupole moment

Similar to dipoles, PyBEST provides a method-independent wrapper to calculate the total (nuclear + electric) quadrupole moment of some electronic wave functions. The code snippet below shows how to determine the total quadrupole moment for some electronic structure method whose results are stored in the IOData container method_output. The quadrupole moment integrals are stored in the quadrupole instance as a tuple (see Computing the multipole moment integrals for more details).

quadrupole = compute_quadrupole(gobasis, x=0, y=0, z=0)

The above function returns each component of the electric quadrupole moment integrals in the order {olp, \(\mu_x\), \(\mu_y\), \(\mu_z\), \(\mu_{xx}\), \(\mu_{xy}\), \(\mu_{xz}\), \(\mu_{yy}\), \(\mu_{yz}\), \(\mu_{zz}\)} . As of now, all quadrupole moment calculations are done strictly taking the origin of the cartesian coordinate system (0,0,0) as the only reference point.

Note

In order to calculate the quadrupole moment, the 1-RDM has to be available and stored in the IOData container that is passed to the method.

In general, the compute_quadrupole() function allows you to calculate the quadrupole moment for a 1-RDM expressed in the atomic or molecular orbital basis. This functionality can be steered using the molecular_orbitals keyword argument. In total, 3 keyword arguments can be set

molecular_orbitals

(boolean) if True, the 1-RDM is transformed back to the atomic orbital basis (default 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 a factor of 2.0

name

(str) if provided, the name is printed in the header of the std output.

16.4.1. 1-RDMs expressed in the atomic orbital basis

PyBEST stores the following 1-RDMs in the atomic orbital basis:

  • The RHF/UHF 1-RDM

The code snippet below shows how to determine the quadrupole moment from an RHF wave function,

# Quadrupole moment integrals of the atomic basis set wrt [0.0, 0.0, 0.0]
# -----------------------------------------------------------------------
quadrupole = compute_quadrupole(gobasis, x=0, y=0, z=0)

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

# Compute quadrupole moment (nuclear and electronic parts)
# --------------------------------------------------------
quadrupole_moment = compute_quadrupole_moment(quadrupole, hf_output)

16.4.2. 1-RDMs expressed in the molecular orbital basis

All 1-RDMs determined in post-Hartree-Fock methods are stored in the molecular orbital basis. Thus, they need to be transformed back to the atomic orbital basis prior to the calculation of the multipole integrals. This back-transformation is performed automatically if the keyword argument molecular_orbitals is set to True.

The current version of PyBEST supports 1-RDMs for

  • MP2 (relaxed and unrelaxed)

  • OOpCCD

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

The code snippet below shows how to determine the quadrupole moment from an OOpCCD wave function assuming the OOpCCD results are stored in oopccd_output,

# Quadrupole moment integrals of the atomic basis set wrt [0.0, 0.0, 0.0]
# -----------------------------------------------------------------------
quadrupole = compute_quadrupole(gobasis, x=0, y=0, z=0)

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

# Converge pCCD
# -------------
oopccd = ROOpCCD(lf, occ_model)
oopccd_output = oopccd(kin, ne, eri, hf_output)

# Compute quadrupole moment, activate keyword for MOs
# ---------------------------------------------------
quadrupole_oopccd = compute_quadrupole_moment(
    quadrupole, oopccd_output, molecular_orbitals=True
)

The main difference is that we have to set the keyword argument molecular_orbitals to True.

16.4.3. Example Python scripts

16.4.3.1. The quadrupole moment of the water molecule determined from a restricted Hartree-Fock calculation

This is a basic example of how to calculate the quadrupole moment in PyBEST. This script performs a restricted Hartree-Fock calculation for the water molecule using the cc-pVDZ basis set.

Listing 16.7 data/examples/multipole/quadrupole_water_scf.py
from pybest import context
from pybest.gbasis import (
    compute_eri,
    compute_kinetic,
    compute_nuclear,
    compute_nuclear_repulsion,
    compute_overlap,
    compute_quadrupole,
    get_gobasis,
)
from pybest.linalg import DenseLinalgFactory
from pybest.occ_model import AufbauOccModel
from pybest.wrappers import RHF, compute_quadrupole_moment

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

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

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

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

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

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

# Decide how to occupy the orbitals (default Aufbau occupation)
# -------------------------------------------------------------
occ_model = AufbauOccModel(gobasis)

# Quadrupole moment integrals of the atomic basis set wrt [0.0, 0.0, 0.0]
# -----------------------------------------------------------------------
quadrupole = compute_quadrupole(gobasis, x=0, y=0, z=0)

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

# Compute quadrupole moment (nuclear and electronic parts)
# --------------------------------------------------------
quadrupole_moment = compute_quadrupole_moment(quadrupole, hf_output)

16.4.3.2. The quadrupole moment of the water molecule determined from a restricted pCCD calculation

This is a basic example of how to calculate the quadrupole moment in PyBEST with the RpCCD method. This script performs a restricted Hartree-Fock calculation followed by an orbital-optimized pCCD calculation for the water molecule using the cc-pVDZ basis set.

Listing 16.8 data/examples/multipole/quadrupole_water_pccd.py
from pybest import context
from pybest.gbasis import (
    compute_eri,
    compute_kinetic,
    compute_nuclear,
    compute_nuclear_repulsion,
    compute_overlap,
    compute_quadrupole,
    get_gobasis,
)
from pybest.geminals import ROOpCCD
from pybest.linalg import DenseLinalgFactory
from pybest.occ_model import AufbauOccModel
from pybest.wrappers import RHF, compute_quadrupole_moment

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

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

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

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

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

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

# Decide how to occupy the orbitals (default Aufbau occupation)
# -------------------------------------------------------------
occ_model = AufbauOccModel(gobasis)


# Quadrupole moment integrals of the atomic basis set wrt [0.0, 0.0, 0.0]
# -----------------------------------------------------------------------
quadrupole = compute_quadrupole(gobasis, x=0, y=0, z=0)

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

# Converge pCCD
# -------------
oopccd = ROOpCCD(lf, occ_model)
oopccd_output = oopccd(kin, ne, eri, hf_output)

# Compute quadrupole moment, activate keyword for MOs
# ---------------------------------------------------
quadrupole_oopccd = compute_quadrupole_moment(
    quadrupole, oopccd_output, molecular_orbitals=True
)

16.4.3.3. The quadrupole moment of the water molecule from RpCCD-LCCSD calculation

The example code given below shows how to determine the quadrupole moment from a pCCD-LCCSD wave function. The script details performing a restricted Hatree-Fock calculation followed by an orbital-optimized pCCD calculation and then an LCCSD calculation on top of that. Notice that the \(\Lambda\)-amplitudes are optimized by setting the keyword argument lambda_equations=True in the LCC function call.

Listing 16.9 data/examples/multipole/quadrupole_water_pccd_lccsd.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,
    compute_quadrupole,
    get_gobasis,
)
from pybest.geminals import ROOpCCD
from pybest.linalg import DenseLinalgFactory
from pybest.occ_model import AufbauOccModel
from pybest.wrappers import RHF, compute_quadrupole_moment

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

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

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

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

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

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

# Decide how to occupy the orbitals (default Aufbau occupation)
# -------------------------------------------------------------
occ_model = AufbauOccModel(gobasis)

# Quadrupole moment integrals of the atomic basis set wrt [0.0, 0.0, 0.0]
# -----------------------------------------------------------------------
quadrupole = compute_quadrupole(gobasis, x=0, y=0, z=0)

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

# Converge pCCD
# -------------
oopccd = ROOpCCD(lf, occ_model)
oopccd_output = oopccd(kin, ne, eri, hf_output)

# LCCSD on top of pCCD
# --------------------
lccsd = RpCCDLCCSD(lf, occ_model)
lccsd_output = lccsd(kin, ne, eri, oopccd_output, lambda_equations=True)

# Compute quadrupole moment from pCCD-LCCSD, activate keyword for MOs
# -------------------------------------------------------------------
quadrupole_lccsd = compute_quadrupole_moment(
    quadrupole, lccsd_output, molecular_orbitals=True
)