14.3. Computing the electric dipole moment

PyBEST provides a method-independent wrapper to calculate the electric dipole moment of some electronic wave function. The code snippet below shows how to determine the electric 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 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 mo keyword argument. In total, 3 keyword arguments can be set

mo

(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.

14.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(gobasis)

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

14.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 ``mo`` is set to True.

The current version of PyBEST supports 1-RDMs for

  • MP2 (relaxed and unrelaxed)

  • OOpCCD

The code snippet below shows how to determine the dipole moment from an OOpCCD wave function,

# Compute dipole moment, activate keyword for MOs
# -----------------------------------------------
dipole_moment = compute_dipole_moment(dipole, oopccd_output, mo=True)

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

14.3.3. Example Python scripts

14.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 14.3 data/examples/multipole/dipole_water_scf.py
from pybest import context
from pybest.gbasis import get_gobasis
from pybest.gbasis.gobasis import (
    compute_dipole,
    compute_eri,
    compute_kinetic,
    compute_nuclear,
    compute_nuclear_repulsion,
    compute_overlap,
)
from pybest.linalg.dense import DenseLinalgFactory
from pybest.scf import AufbauOccModel
from pybest.wrappers.hf import RHF
from pybest.wrappers.multipole import compute_dipole_moment
from pybest.utils import get_com


# 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)
na = compute_nuclear(gobasis)
er = 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 (5 alpha electrons)
# -----------------------------------------------------
occ_model = AufbauOccModel(5)

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

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

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

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

14.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. 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 14.4 data/examples/multipole/dipole_water_pccd.py
from pybest import context
from pybest.gbasis import get_gobasis
from pybest.gbasis.gobasis import (
    compute_dipole,
    compute_eri,
    compute_kinetic,
    compute_nuclear,
    compute_nuclear_repulsion,
    compute_overlap,
)
from pybest.geminals import ROOpCCD
from pybest.linalg.dense import DenseLinalgFactory
from pybest.scf import AufbauOccModel
from pybest.wrappers.hf import RHF
from pybest.wrappers.multipole import compute_dipole_moment
from pybest.utils import get_com


# 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)
na = compute_nuclear(gobasis)
er = 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 (5 alpha electrons)
# -----------------------------------------------------
occ_model = AufbauOccModel(5)

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

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

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

# Converge pCCD
# -------------
oopccd = ROOpCCD(lf, occ_model)
oopccd_output = oopccd(kin, na, er, hf_output)

# Compute dipole moment, activate keyword for MOs
# -----------------------------------------------
dipole_moment = compute_dipole_moment(dipole, oopccd_output, mo=True)