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)

(1)

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}

(2)

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]^*

(3)

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

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

(6)

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*}

(7)

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*}

(8)

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.

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*}

(10)

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*}

(11)

The MO-representation of the $\hat{\kappa}^{(1)}$ -operator can be determined from linear response theory Norman et al. (2018) 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

(12)

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

(13)

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

Loading...

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

(14)

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

(15)

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

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

(16)

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

(17)

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.

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

(19)

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.