Time-dependence of observables

In time-dependent Hartree–Fock, we are determining the time-dependence of the orbital rotation parameters (or electron transfer parameters). In spectroscopic applications, this time-dependence arises due to the system being subjected to an external time-dependent electromagnetic field as described by the Hamiltonian

\[ \hat{H}(t) = \hat{H}_0 + \hat{V}(t) \]

where, in response theory, the coupling term, \(\hat{V}(t)\), is treated as a perturbation and expanded in the frequency domain according to

\[ \hat{V}(t) = \sum_\omega \hat{V}^\omega_\alpha F^\omega_\alpha e^{-i \omega t} \]

The sum over \(\omega\) includes both positive and negative Fourier components of the field and we have

\[ \hat{V}^{-\omega}_\alpha = \Big[ \hat{V}^{\omega}_\alpha \Big]^\dagger; \qquad F^{-\omega}_\alpha =\Big[ F^{\omega}_\alpha \Big]^* \]

There is also an implied Einstein summation over the Cartesian molecular axes \(\alpha=\{x,y,z\}\).

Representing a monochromatic electric field

In the electric-dipole approximation, the coupling between the system and a monochromatic sinusoidal electric field is given by

\[ \hat{V}(t) = -\hat{\mu}_\alpha E^{\omega_0}_\alpha \sin(\omega_0 t) \]

In this case, the summation over \(\omega\) includes two terms namely \(\omega = \{\pm\omega_0\}\), and we have

\[ \hat{V}^{\omega_0}_\alpha = -\hat{\mu}_\alpha; \qquad F^{\omega_0}_\alpha = \frac{1}{2i} E^{\omega_0}_\alpha \]

Perturbation expansions

The time-dependent electron transfer amplitudes are expanded in orders of the perturbation

\[ \kappa_{ai}(t) = \kappa_{ai}^{(1)}(t) + \kappa_{ai}^{(2)}(t) + \ldots \]

where it has been used that the zeroth-order solution to the phase-isolated wave function is equal to the reference state, or, in other words, \(\kappa_{ai}^{(0)}(t) = 0\).

With use of Fourier expansions, we have

\[\begin{align*} \kappa_{ai}^{(1)}(t) & = \sum_\omega \kappa_{ai}^{(1)}(\omega) e^{-i\omega t} \\ \kappa_{ai}^{(2)}(t) & = \sum_{\omega_1, \omega_2} \kappa_{ai}^{(2)}(\omega_1, \omega_2) e^{-i(\omega_1 + \omega_2) t} \\ \end{align*}\]

et cetera for higher orders. The arguments are here used to separate the time-dependent functions from their Fourier amplitudes.

The terms in the generator of orbital rotations that are linear in the perturbation can be written

\[\begin{align*} \hat{\kappa}^{(1)}(t) & = \sum_{a}^\mathrm{virt} \sum_{i}^\mathrm{occ} \Big[ \kappa^{(1)}_{ai}(t) \, \hat{a}^{\dagger}_a \hat{a}_i + \big[\kappa^{(1)}_{ai}(t)\big]^\ast \, \hat{a}^{\dagger}_i \hat{a}_a \Big] \\ & = \sum_\omega \sum_{a}^\mathrm{virt} \sum_{i}^\mathrm{occ} \Big[ \kappa^{(1)}_{ai}(\omega) \, \hat{a}^{\dagger}_a \hat{a}_i + \big[\kappa^{(1)}_{ai}(-\omega)\big]^\ast \, \hat{a}^{\dagger}_i \hat{a}_a \Big] \, e^{-i\omega t} \\ & \equiv \sum_\omega \hat{\kappa}^{(1)}(\omega) \, e^{-i\omega t} \end{align*}\]

where, in the second term, we have made use of the fact that we are free to change sign of the frequency summation variable in the Fourier expansion.

Spin-adapted electron transfer operators

With use of spin-adapted electron transfer operators, we get

\[\begin{align*} \hat{\kappa}^{(1)}(\omega) & = \sum_{a}^\mathrm{virt} \sum_{i}^\mathrm{occ} \Big[ \kappa^{(1)}_{ai}(\omega) \, \hat{E}^{\dagger}_{ai} + \big[\kappa^{(1)}_{ai}(-\omega)\big]^\ast \, \hat{E}^{\dagger}_{ia} \Big] \\ \hat{E}^{\dagger}_{pq} & = \hat{a}^{\dagger}_{p\alpha} \hat{a}_{q\alpha} + \hat{a}^{\dagger}_{p\beta} \hat{a}_{q\beta} \end{align*}\]

where summations run over molecular orbitals.

Expectation values of observables

In time-dependent Hartree–Fock, the expectation value of an observable becomes

\[\begin{align*} \langle \hat{\Omega} \rangle(t) & = \langle \Psi(t) | \hat{\Omega} | \Psi(t) \rangle = \langle \bar{\Psi}(t) | \hat{\Omega} | \bar{\Psi}(t) \rangle = \langle 0 | e^{i\hat{\kappa}} \hat{\Omega} e^{-i\hat{\kappa}} | 0 \rangle \\ & = \langle 0 | \hat{\Omega} | 0 \rangle + i \langle 0 | \Big[ \hat{\kappa}(t), \hat{\Omega} \Big] | 0 \rangle - \frac{1}{2} \langle 0 | \Big[ \hat{\kappa}(t), \Big[ \hat{\kappa}(t), \hat{\Omega} \Big]\Big] | 0 \rangle + \cdots \end{align*}\]

where we have made use of the fact that the expectation value can equally well be determined with respect to the phase isolated wave function.

The contributions to the expectation value that are linear in the perturbation can be written

\[\begin{align*} \langle \hat{\Omega} \rangle^{(1)}(t) & = i \langle 0 | \Big[ \hat{\kappa}^{(1)}(t), \hat{\Omega} \Big] | 0 \rangle = i \sum_\omega \langle 0 | \Big[ \hat{\kappa}^{(1)}(\omega), \, \hat{\Omega} \Big] | 0 \rangle \, e^{-i\omega t} \equiv \sum_\omega \langle \hat{\Omega} \rangle^{(1)}(\omega) \, e^{-i\omega t} \end{align*}\]

The MO-representation of the \(\hat{\kappa}^{(1)}\)-operator can be determined from linear response theory [NRS18] and calculated with use of the VeloxChem program, and it is thereafter a trivial task to calculate the expectation value of its commutator with the observable \(\hat{\Omega}\).

The factor “\(i\)” multiplying the expectation value is typically suppressed in intermediate steps of the calculations as to return at a point of convenience.

Response functions

Response functions are defined from the order expansions of observables. The linear response functions are identified from the expression

\[ \langle \hat{\Omega} \rangle^{(1)}(\omega) = \langle \langle \hat{\Omega}; \hat{V}^\omega_\alpha \rangle \rangle_\omega \, F^\omega_\alpha \]

As an example, the tensor components of the electric-dipole polarizability can be determined as follows

\[ \alpha_{\alpha\beta}(-\omega;\omega) = - \langle \langle \hat{\mu}_\alpha; \hat{\mu}_\beta \rangle \rangle_\omega \]

Illustration

import numpy as np
import veloxchem as vlx

ethylene_xyz = """
6

C        0.67759997    0.00000000   0.00000000
C       -0.67759997    0.00000000   0.00000000
H        1.21655197    0.92414474   0.00000000
H        1.21655197   -0.92414474   0.00000000
H       -1.21655197   -0.92414474   0.00000000
H       -1.21655197    0.92414474   0.00000000
"""

molecule = vlx.Molecule.read_xyz_string(ethylene_xyz)
basis = vlx.MolecularBasis.read(molecule, "def2-svp", ostream=None)

nocc = molecule.number_of_alpha_electrons()
norb = basis.get_dimension_of_basis(molecule)
nvirt = norb - nocc

print("Number of occupied orbitals :", nocc)
print("Number of virtual orbitals  :", nvirt)
print("Number of molecular orbitals:", norb)

Number of occupied orbitals : 8
Number of virtual orbitals  : 40
Number of molecular orbitals: 48

molecule.show()

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

This system is subjected to an electric field along the \(x\)-direction with a angular frequency of \(\omega = 0.0656\) a.u. and an amplitude \(F^\omega_x\).

Within the electric-dipole approximation the perturbation operator is in this case given by

\[ \hat{V}^\omega = -\hat{\mu}_x F_x^\omega \]

and we will be concerned with the determination of the induced dipole moment along the C–C bond axis associated with the expectation value of the operator

\[ \hat{\Omega} = \hat{\mu}_x \]

First, we determine the Hartree–Fock reference state.

scf_drv = vlx.ScfRestrictedDriver()
scf_drv.ostream.mute()

scf_results = scf_drv.compute(molecule, basis)

Second, we get the MO representation of the operator \(\hat{\Omega}\).

dipole_drv = vlx.ElectricDipoleIntegralsDriver()
dipole_mats = dipole_drv.compute(molecule, basis)

mu_x_ao = -1.0 * dipole_mats.x_to_numpy()

C = scf_results["C_alpha"]

mu_x = np.einsum("ap, ab, bq -> pq", C, mu_x_ao, C)

Third, we determine the linear response parameters.

lrf_drv = vlx.LinearResponseSolver()
lrf_drv.ostream.mute()

lrf_drv.a_operator = "electric dipole"
lrf_drv.b_operator = "electric dipole"

lrf_drv.a_components = ["x"]
lrf_drv.b_components = ["x"]

lrf_drv.frequencies = [0.0656]

lrf_results = lrf_drv.compute(molecule, basis, scf_results)

In VeloxChem, the resulting linear response parameters are structured into 1D-arrays known as response vectors

\[\begin{split} \mathbf{N}^\alpha(\omega) = \begin{pmatrix} \mathbf{Z}^{\alpha; \omega} \\ \mathbf{Y}^{\alpha; \omega} \end{pmatrix} \end{split}\]

The response vectors depend on the angular frequency, \(\omega\), and they associated with separate Cartesian axes, \(\alpha\). The elements of the response vector relate to the wave function parameters

\[ \kappa_{ai}^{(1)}(\omega) = i Z_{ai}^{\alpha; \omega} F^\omega_\alpha ; \qquad \kappa_{ia}^{(1)}(-\omega) = i Y_{ia}^{\alpha; \omega} F^{-\omega}_\alpha \]

where “\(ai\)” and “\(ia\)” are thought of as MO-matrix indices for the \(\kappa\)-parameters but compound vector indices in the 1D-arrays. Here appears a factor of “\(i\)” that combined with the aforementioned one give an overall factor of -1 in the final result.

Spin-adapted electron transfer operators in VeloxChem

In VeloxChem, normalzed spin-adapted electron transfer operators are adopted

\[ \hat{E}_{pq} = \frac{1}{\sqrt{2}} \Big[ \hat{a}^{\dagger}_{p\alpha} \hat{a}_{q\alpha} + \hat{a}^{\dagger}_{p\beta} \hat{a}_{q\beta} \Big] \]

Among other things, with this choice the diagonal elements in the electronic Hessian will be close to orbital energy differences and not two times these values. Whatever the choice, final results will be identical as long as it is implemented in a consistent manner.

Response vectors are obtained from VeloxChem with use of the get_full_vector(index) method, where indices of “\(0\)” and “\(1\)” will instruct the method to return vectors \((\mathbf{Z} + \mathbf{Y})\) and \((\mathbf{Z} - \mathbf{Y})\), respectively.

Z = -0.5 * (
    lrf_results["solutions"][("x", 0.0656)].get_full_vector(0)
    + lrf_results["solutions"][("x", 0.0656)].get_full_vector(1)
) * np.sqrt(2)

Y = -0.5 * (
    lrf_results["solutions"][("x", 0.0656)].get_full_vector(0)
    - lrf_results["solutions"][("x", 0.0656)].get_full_vector(1)
) * np.sqrt(2)

A factor of -1 has here been introduced in recognition of the fact that the perturbation operator is equal to minus the electric dipole operator and the factor of \(\sqrt{2}\) is associated with the aforementioned normalization of electron transfer operators.

We scatter the response parameters into an MO-matrix.

def vector_to_matrix(Z, Y):

    zymat = np.zeros((norb, norb), dtype=float)

    idx = 0
    for i in range(nocc):
        for a in range(nocc, norb):
            zymat[a, i] = Z[idx]
            zymat[i, a] = Y[idx]
            idx += 1

    return zymat

zymat = vector_to_matrix(Z, Y)

Fourth, we determine the expectation value of the commutator between \(\hat{\kappa}^{(1)}(\omega)\) and \(\hat{\Omega}\) to determine the electric-dipole polarizability

\[ \alpha_{xx}(-\omega;\omega) = \langle \hat{\Omega} \rangle^{(1)}(\omega) / F^\omega_x \]

def commutator(A, B):

    return np.matmul(A, B) - np.matmul(B, A)

def expectation_value(A):  # closed shell state

    return 2 * np.einsum("ii ->", A[:nocc, :nocc])

alpha_xx = - expectation_value(commutator(zymat, mu_x))

We note that the factor of -1 comes from the left out factor of \(i^2\).

print(f"Polarizability: {alpha_xx:12.8f}")

Polarizability:  34.84380178

We compare our result with the reference value from VeloxChem.

print(f"Reference: {-lrf_results['response_functions'][('x', 'x', 0.0656)]:12.8f}")

Reference:  34.84380178

We note that the two results are in perfect agreement.