16.3. Computing the dipole moment¶

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

dipole_moment = compute_dipole_moment(dipole, method_output)

The above function returns each component of the electric dipole moment in the order {x,y,z}.

Note

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

In general, the compute_dipole() function allows you to calculate the dipole 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.3.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 dipole moment from an RHF wave function,

# Get center of mass
# ------------------
x, y, z = get_com(factory)

# electric dipole moment integrals of the atomic basis set wrt COM
# ----------------------------------------------------------------
dipole = compute_dipole(factory, x=x, y=y, z=z)

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

# Compute dipole moment
# ---------------------
dipole_moment = compute_dipole_moment(dipole, hf_output)

Note

The get_com() function calculates the center of charge (COM, in our terminology) of the molecular structure. Coordinates (in bohr) of any other reference point can be provided by the user by setting the values of x, y, and z directly, while calling the compute_dipole() function.

dipole = compute_dipole(gobasis, x=x_ref, y=y_ref, z=z_ref)

The nuclear dipole moments are also calculated and the total dipole moments are printed upon calling the compute_dipole() function.

16.3.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 dipole moment from an OOpCCD wave function assuming the OOpCCD results are stored in oopccd_output,

# -----------------------------------------------
dipole_oopccd = compute_dipole_moment(
    dipole, oopccd_output, molecular_orbitals=True
)

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

16.3.3. Example Python scripts¶

16.3.3.1. The dipole moment of the water molecule determined from a restricted Hartree-Fock calculation¶

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

Listing 16.4 data/examples/multipole/dipole_water_scf.py¶

from pybest import context
from pybest.gbasis import (
    compute_dipole,
    compute_eri,
    compute_kinetic,
    compute_nuclear,
    compute_nuclear_repulsion,
    compute_overlap,
    get_gobasis,
)
from pybest.linalg import DenseLinalgFactory
from pybest.occ_model import AufbauOccModel
from pybest.utils import get_com
from pybest.wrappers import RHF, compute_dipole_moment

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

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

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

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

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

# 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(factory)

# Get center of mass
# ------------------
x, y, z = get_com(factory)

# electric dipole moment integrals of the atomic basis set wrt COM
# ----------------------------------------------------------------
dipole = compute_dipole(factory, x=x, y=y, z=z)

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

# Compute dipole moment
# ---------------------
dipole_moment = compute_dipole_moment(dipole, hf_output)

16.3.3.2. The dipole moment of the water molecule determined from a restricted pCCD calculation¶

This is a basic example of how to calculate the dipole 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.5 data/examples/multipole/dipole_water_pccd.py¶

from pybest import context
from pybest.cc import RpCCDLCCSD
from pybest.gbasis import (
    compute_dipole,
    compute_eri,
    compute_kinetic,
    compute_nuclear,
    compute_nuclear_repulsion,
    compute_overlap,
    get_gobasis,
)
from pybest.geminals import ROOpCCD
from pybest.linalg import DenseLinalgFactory
from pybest.occ_model import AufbauOccModel
from pybest.utils import get_com
from pybest.wrappers import RHF, compute_dipole_moment

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

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

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

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

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

# 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(factory)

# Get center of mass
# ------------------
x, y, z = get_com(factory)

# electric dipole moment integrals of the atomic basis set wrt COM
# ----------------------------------------------------------------
dipole = compute_dipole(factory, x=x, y=y, z=z)

# 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 dipole moment, activate keyword for MOs
# -----------------------------------------------
dipole_oopccd = compute_dipole_moment(
    dipole, oopccd_output, molecular_orbitals=True
)

lccsd = RpCCDLCCSD(lf, occ_model)
lccsd_output = lccsd(kin, ne, eri, oopccd_output, lambda_equations=True)

dipole_lccsd = compute_dipole_moment(
    dipole, lccsd_output, molecular_orbitals=True
)

16.3.3.3. The dipole moment of the water molecule from RpCCD-LCCSD calculation¶

The example code given below shows how to determine the dipole 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 \(\Lambda\)-amplitudes are optimized by setting the keyword argument lambda_equations=True in the LCC function call.

Listing 16.6 data/examples/multipole/dipole_water_pccd_lccsd.py¶

from pybest import context
from pybest.cc import RpCCDLCCSD
from pybest.gbasis import (
    compute_dipole,
    compute_eri,
    compute_kinetic,
    compute_nuclear,
    compute_nuclear_repulsion,
    compute_overlap,
    get_gobasis,
)
from pybest.geminals import ROOpCCD
from pybest.linalg import DenseLinalgFactory
from pybest.occ_model import AufbauOccModel
from pybest.utils import get_com
from pybest.wrappers import RHF, compute_dipole_moment

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

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

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

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

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

# 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(factory)

# Get center of mass
# ------------------
x, y, z = get_com(factory)

# electric dipole moment integrals of the atomic basis set wrt COM
# ----------------------------------------------------------------
dipole = compute_dipole(factory, x=x, y=y, z=z)

# 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 dipole moment from pCCD-LCCSD, activate keyword for MOs
# ---------------------------------------------------------------
dipole_lccsd = compute_dipole_moment(
    dipole, lccsd_output, molecular_orbitals=True
)