19.6. Linear Response Transition Dipole Moment and Related Properties

The current version of PyBEST supports property calculations from pCCD-based (see Detailed input structure of pCCD) Linear Response methods (see The Linear Response Module). Thus, we will focus on a pCCD reference function below. The Transition Dipole Moment (TDM), in the absence of any external potential, can be obtained using only the \(\hat{A}\)-dependent part of the full linear response expression, given as:

\[\begin{split}\Gamma_{0 \rightarrow k}^{\hat{A}} &= \langle 0 | \hat{A} | k \rangle \\ &= \langle \Lambda |[\hat{A}, \, \hat{\tau}_k]|\mathrm{pCCD}\rangle \\ &\quad - \sum_{\nu} (-\hat{J} + \omega \boldsymbol{I})_{\mu\nu}^{-1} \langle \bar{\nu} | \hat{A} | \mathrm{pCCD} \rangle \langle \Lambda | [[\hat{H}_0, \hat{\tau}_\nu], \hat{\tau}_k] | \mathrm{pCCD} \rangle\end{split}\]

Using dipole operators as \(\hat{A}\), one can compute the transition dipole moments for selected excited states.

An intrinsic advantage of using pCCD as the reference wave function, unlike CCSD-based approaches, is that the Jacobian matrix becomes approximately symmetric. This symmetry allows the transition dipole moment to be approximated using only the right eigenvectors without significantly compromising accuracy.

This approximation enables reduced computational scaling to \(\mathcal{O}(o^2v^2)\).

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

19.6.1. Implemented Features

Molecular properties are determined from Linear Response calculations. Currently, PyBEST supports two different Linear Response flavors (see The Linear Response Module).

• LRpCCD: Linear Response based on a pCCD reference

• LRpCCDS: Linear Response based on a pCCD with single excitations

Each class is derived from its corresponding EOM class (see Quick Guide) and reuses most of the functionality, modifying only the construction of the Hamiltonian to match the true LR structure.

The corresponding LR property calculations support the following options when determining transition properties. They are collected as a dictionary and passed as the property_options argument during function call. Supported dictionary keys are:

operator_A

One-electron operator (e.g., dipole integrals)

operator_B

One-electron operator (e.g., dipole integrals)

coordinates

Used for origin-dependent properties (list containing x, y, z coordinates)

transition_dipole_moment

(bool) Whether to compute transition dipole moments

19.6.2. Performance and Applicability

These LR methods are especially useful for:

• Accurate and efficient description of singly and doubly excited states

• Cost-effective alternatives to EOM-CCSD for large systems

• Systems requiring particle-number-conserving excitation operators

19.6.3. Quick Guide: Transiton dipole moment for pCCD and pCCD+S

This is a short introduction to the Linear Response module. More information on the input and output structure can be found in the following sections. Similar to the previous modules, we will assume the following names for all PyBEST objects

lf

a LinalgFactory instance (see Preliminaries).

occ_model

an Aufbau occupation model of the AufbauOccModel class

kin

the kinetic energy integrals

ne

the nucleus-electron attraction integrals

eri

the two-electron repulsion integrals

The current version of PyBEST supports two different Linear Response flavors, both based on a pCCD reference function (The pCCD module). For the transition dipole moment calculation, you need to pass the integrals of the dipole operator (see Computing the multipole moment integrals) and the molecular coordinates of the origin. The code snippet below demonstrates how to generate or extract this data from previously generated PyBEST output pccd_output.

from pybest.gbasis import compute_dipole
from pybest.utility import get_com

# Compute dipole integrals for center of charge/mass
x, y, z = get_com(basis)
dipole = compute_dipole(basis, x=x, y=y, z=z)

# Return the origin from previously generated output
from pybest.wrappers.multipole import check_coord

# check_coord returns the origin of the dipole integrals
coord = check_coord(dipole, pccd_output)

The LRpCCD class allows you to determine the Linear Response for pair excitations only. A code snippet can be found below.

from pybest.properties import LRpCCDS

# For pCCD Linear Response (only electron pairs)

# Compute transition dipole moment
tm_dipole = LRpCCD(lf, occ_model)
results = tm_dipole(
    kin,
    ne,
    eri,
    jac_output,
    property_options={
        "operator_A": dipole,
        "operator_B": dipole,
        "coordinates": coord,
        "transition_dipole_moment": True,
    },
    printoptions={"nroot": 2},
)

To include single excitations in the Linear Response, you can use the corresponding LRpCCDS class.

from pybest.properties import LRpCCDS

# For pCCD+S Linear Response

# Compute transition dipole moment
tm_dipole = LRpCCDS(lf, occ_model)
results = tm_dipole(
    kin,
    ne,
    eri,
    jac_output,
    property_options={
        "operator_A": dipole,
        "operator_B": dipole,
        "coordinates": coord,
        "transition_dipole_moment": True,
    },
    printoptions={"nroot": 2},
)

The result includes excitation energies and electronic, nuclear, and total transition dipole moments, along with oscillator strength from the Jacobian matrix. The results are returned as an IOData container, while all results are saved to the pybest-results/checkpoint_LR_pCCD.h5 or pybest-results/checkpoint_LR_pCCDS.h5 files. Specifically, the IOData container contains the following attributes.

orb_a

Molecular orbitals from the pCCD reference

e_ref

Total reference energy (pCCD or pCCDS)

e_ee

Excitation energies (first entry is 0.0)

civ_ee

Eigenvectors (CI-like) of the LR Jacobian

t_dm

Transition dipole moments (if requested)

19.6.4. Keyword Arguments

The LR modules support the following keyword arguments:

19.6.4.1. General Options

nroot

(int) Number of excited states to compute (default: 1)

threshold

(float) Threshold for printing CI vector contributions (default: 0.1)

dump_cache

(bool) Cache intermediate matrices on disk (default: True for large systems)

19.6.4.2. Davidson Options

tolerance

(float) Convergence threshold for excitation energies (default: 1e-6)

tolerancev

(float) Threshold for CI eigenvectors (default: 1e-4)

maxiter

(int) Maximum number of iterations (default: 200)

nguessv

(int) Number of initial guess vectors (default: (nroot - 1)*4 + 1)

maxvectors

(int) Maximum number of vectors in Davidson space (default: (nroot - 1)*10)

todisk

(bool) Write intermediates to disk (default: False)

19.6.5. Example Usage

The example files for obtaining transition dipole moments using the linear response method are collected in the directory data/examples/properties.

19.6.5.1. Example for LR-pCCD transition dipole moments and oscillator strength

This is a basic example of how to calculate properties from linear response methods in PyBEST. This script performs an RHF calculation, followed by a pCCD and JacobianpCCD calculation, and an LR-pCCD calculation with a pCCD reference function on the water molecule using the cc-pVDZ basis set.

Listing 19.11 data/examples/properties/water_lr_tdm_pccd.py

from pybest import context
from pybest.ee_jacobian import JacobianpCCD
from pybest.gbasis import (
    compute_dipole,
    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 LRpCCD
from pybest.utility import get_com
from pybest.wrappers import RHF
from pybest.wrappers.multipole import check_coord

# Set up molecule, define basis set
fn_xyz = context.get_fn("test/water.xyz")
basis = get_gobasis("cc-pvdz", fn_xyz)

# Define occupation model, orbitals, and overlap
lf = DenseLinalgFactory(basis.nbasis)
occ_model = AufbauOccModel(basis)
orb_a = lf.create_orbital(basis.nbasis)
olp = compute_overlap(basis)


# Construct Hamiltonian
kin = compute_kinetic(basis)
ne = compute_nuclear(basis)
eri = compute_eri(basis)
nuc = compute_nuclear_repulsion(basis)
nuc = compute_nuclear_repulsion(basis)
# dipole moment
x, y, z = get_com(basis)
dipole = compute_dipole(basis, x=x, y=y, z=z)

# Do a Hartree-Fock calculation
hf = RHF(lf, occ_model)
hf_output = hf(kin, ne, eri, nuc, olp, orb_a)

# Do pCCD calculations
pccd = RpCCD(lf, occ_model)
pccd_output = pccd(kin, ne, eri, hf_output)
# check coordinates
coord = check_coord(dipole, pccd_output)

# Compute excitation energies using Jacobian matrix.

jac = JacobianpCCD(lf, occ_model)
jac_output = jac(kin, ne, eri, pccd_output, davidson=True)

# Compute transition dipole moment
tm_dipole = LRpCCD(lf, occ_model)
dipole_output = tm_dipole(
    kin,
    ne,
    eri,
    jac_output,
    property_options={
        "operator_A": dipole,
        "operator_B": dipole,
        "coordinates": coord,
        "transition_dipole_moment": True,
    },
    printoptions={"nroot": 2},
)

19.6.5.2. Example for LR-pCCD+S transition dipole moments and oscillator strength

This is a basic example of how to calculate properties from linear response methods in PyBEST. This script performs an RHF calculation, followed by a pCCD and JacobianpCCD+S calculation, and an LR-pCCD+S calculation with a pCCD reference function on the water molecule using the cc-pVDZ basis set.

Listing 19.12 data/examples/properties/water_lr_tdm_pccd_s.py

from pybest import context
from pybest.ee_jacobian import JacobianpCCDS
from pybest.gbasis import (
    compute_dipole,
    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 LRpCCDS
from pybest.utility import get_com
from pybest.wrappers import RHF
from pybest.wrappers.multipole import check_coord

# Set up molecule, define basis set
fn_xyz = context.get_fn("test/water.xyz")
basis = get_gobasis("cc-pvdz", fn_xyz)

# Define occupation model, orbitals, and overlap
lf = DenseLinalgFactory(basis.nbasis)
occ_model = AufbauOccModel(basis)
orb_a = lf.create_orbital(basis.nbasis)
olp = compute_overlap(basis)


# Construct Hamiltonian
kin = compute_kinetic(basis)
ne = compute_nuclear(basis)
eri = compute_eri(basis)
nuc = compute_nuclear_repulsion(basis)
# dipole moment
x, y, z = get_com(basis)
dipole = compute_dipole(basis, x=x, y=y, z=z)

# Do a Hartree-Fock calculation
hf = RHF(lf, occ_model)
hf_output = hf(kin, ne, eri, nuc, olp, orb_a)

# Do pCCD calculations
pccd = RpCCD(lf, occ_model)
pccd_output = pccd(kin, ne, eri, hf_output)
# check coordinates
coord = check_coord(dipole, pccd_output)

# Compute excitation energies using Jacobian matrix.

jac = JacobianpCCDS(lf, occ_model)
jac_output = jac(kin, ne, eri, pccd_output, davidson=True)

# Compute transition dipole moment
tm_dipole = LRpCCDS(lf, occ_model)
dipole_output = tm_dipole(
    kin,
    ne,
    eri,
    jac_output,
    property_options={
        "operator_A": dipole,
        "operator_B": dipole,
        "coordinates": coord,
        "transition_dipole_moment": True,
    },
    printoptions={"nroot": 7},
)