Energy determination of parameters

import numpy as np
import veloxchem as vlx

Order expansion of parameters¶

The wave function parameters are expanded in orders of the external field

\kappa_{ai}(F) = \kappa_{ai}^{(1)} + \kappa_{ai}^{(2)} + \ldots

(1)

such that $\kappa_{ai}= 0$ in the limit of $F = 0$ .

Conversely,

\kappa_{ai}^{(1)} = \left. \frac{ d \kappa_{ai} }{ d F_{\alpha} } \right|_{F=0} F_{\alpha}

(2)

Linear response equation¶

The first-order response can be determined from the field-derivative of the variational condition

\begin{align*} 0 & = \Big[ \frac{ d }{ d F_{\alpha} } \frac{\partial E(\kappa, F)}{\partial \kappa_{bj}^*} \Big]_{F=0} \\ & = \Big[ \frac{\partial^2 E(\kappa, F)}{\partial \kappa_{bj}^* \partial F_\alpha} + \sum_{ai} \Big( \frac{ \partial^2 E(\kappa, F) }{ \partial \kappa_{bj}^* \partial \kappa_{ai} } \frac{d \kappa_{ai}}{d F_\alpha} + \frac{ \partial^2 E(\kappa, F) }{ \partial \kappa_{bj}^* \partial \kappa_{ai}^* } \frac{d \kappa_{ai}^*}{d F_\alpha} \Big) \Big]_{F=0} \end{align*}

(3)

We here view the $ai$ -pair of unoccupied/occupied orbital indices as a compound index $n$ (and similarly for $bj$ and $m$ ), inviting us to introduce a compact matrix notation to arrive at the linear response equation.

\mathbf{E}^{[2]} \mathbf{N} = - \mathbf{B}^{[1]}

(4)

In order, we have introduced:

The Hermitian electronic Hessian

\mathbf{E}^{[2]} = \begin{pmatrix} \mathbf{A} & \mathbf{B} \\ \mathbf{B}^* & \mathbf{A}^* \\ \end{pmatrix}

(5)

with blocks

A_{mn} = \left. \frac{ \partial^2 E(\kappa, F) }{ \partial \kappa_{m}^* \partial \kappa_{n} } \right|_{F = 0} ; \qquad B_{mn} = \left. \frac{ \partial^2 E(\kappa, F) }{ \partial \kappa_{m}^* \partial \kappa_{n}^* } \right|_{F = 0}

(6)

The (unknown) response vector

\mathbf{N} = \begin{pmatrix} \mathbf{Z} \\ \mathbf{Z}^* \end{pmatrix}

(7)

with

Z_n = \left. \frac{d \kappa_n}{d F_\alpha} \right|_{F=0}

(8)

And the property gradient

\mathbf{B}^{[1]} = \begin{pmatrix} \mathbf{g} \\ \mathbf{g}^* \\ \end{pmatrix}

(9)

with

g_m = \left. \frac{\partial^2 E(\kappa, F)}{\partial \kappa_{m}^* \partial F_\alpha} \right|_{F=0}

(10)

Baker–Campbell–Hausdorff expansion¶

We have

\begin{align*} E(\kappa,F) & = \langle 0 | e^{i\hat{\kappa}} \hat{H} e^{-i\hat{\kappa}} | 0 \rangle \\ & = \langle 0 | \hat{H} | 0 \rangle + i \langle 0 | [\hat{\kappa}, \hat{H}] | 0 \rangle - \frac{1}{2} \langle 0 | [\hat{\kappa}, [\hat{\kappa}, \hat{H}]] | 0 \rangle + \cdots \end{align*}

(11)

and let us adopt the electric-dipole approximation

\hat{H} = \hat{H}_0 - \hat{\mu}_\alpha F_\alpha

(12)

It is then straightforward to evaluate that electronic Hessian and the property gradient. We get

\begin{align*} A_{mn} & = - \langle 0 | [\hat{a}^\dagger_j \hat{a}_b, [\hat{a}^\dagger_a \hat{a}_i, \hat{H}_0]] | 0 \rangle \\ & = \langle 0_j^b | \hat{H}_0 | 0_i^a \rangle - \langle 0 | \hat{H}_0 | 0 \rangle \delta_{ab} \delta_{ij} \\ B_{mn} & = - \langle 0 | [\hat{a}^\dagger_j \hat{a}_b, [\hat{a}^\dagger_i \hat{a}_a, \hat{H}_0]] | 0 \rangle \\ & = - \langle 0_{ij}^{ab} | \hat{H}_0 | 0 \rangle ( 1 - \delta_{ab} ) ( 1 - \delta_{ij} ) \end{align*}

(13)

and

\begin{align*} g_m & = i \langle 0 | [\hat{a}^\dagger_j \hat{a}_b, -\hat{\mu}_\alpha ] | 0 \rangle \\ & = - i \langle 0_j^b | \hat{\mu}_\alpha | 0 \rangle \\ & = -i \langle \psi_j | \hat{\mu}_\alpha | \psi_b \rangle \end{align*}

(14)

Spin-adapted excitation operators¶

For closed-shell reference states in VeloxChem, we are using normalized spin-adapted excitation operators that will result in

\begin{align*} A_{mn} &= (\varepsilon_a - \varepsilon_i) \delta_{ab} \delta_{ij} + 2 (jb|ia) - (ji|ba) \\ B_{mn} &= - 2 (jb|ia) + (ja|ib) \end{align*}

(15)

and

\begin{align*} g_m & = -i \sqrt{2} ( \phi_j | \hat{\mu}_\alpha | \phi_b ) \end{align*}

(16)

and, in this case, the block dimension equals

n_\mathrm{dim} = n_\mathrm{occ} \times n_\mathrm{virt}

(17)

Illustration¶

Setup¶

h2o_xyz = """3

O    0.000000000000        0.000000000000        0.000000000000
H    0.000000000000        0.740848095288        0.582094932012
H    0.000000000000       -0.740848095288        0.582094932012
"""

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

molecule = vlx.Molecule.read_xyz_string(h2o_xyz)
basis = vlx.MolecularBasis.read(molecule, "6-31g", ostream=None)

scf_results = scf_drv.compute(molecule, basis)

epsilon = scf_results["E_alpha"]  # orbital energies

erimo_drv = vlx.MOIntegralsDriver()  # ERI blocks in the MO basis

ovov = erimo_drv.compute_in_memory(molecule, basis, scf_drv.mol_orbs, "chem_ovov")
oovv = erimo_drv.compute_in_memory(molecule, basis, scf_drv.mol_orbs, "chem_oovv")

print("(ov|ov):", ovov.shape)
print("(oo|vv):", oovv.shape)

(ov|ov): (5, 8, 5, 8)
(oo|vv): (5, 5, 8, 8)

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

print("Number of occupied orbitals:", nocc)
print("Number of virtual orbitals :", nvirt)
print("Block dimension            :", ndim)

Number of occupied orbitals: 5
Number of virtual orbitals : 8
Block dimension            : 40

Electronic Hessian¶

We determine the electronic Hessian

\mathbf{E}^{[2]} = \begin{pmatrix} \mathbf{A} & \mathbf{B} \\ \mathbf{B}^* & \mathbf{A}^* \\ \end{pmatrix}

(18)

def delta(i, j):

    if i == j:
        delta = 1
    else:
        delta = 0

    return delta

E2 = np.zeros((2 * ndim, 2 * ndim))  # real valued

A = np.zeros((ndim, ndim))
B = np.zeros((ndim, ndim))

m = 0
for j in range(nocc):
    for b in range(nvirt):
        n = 0
        for i in range(nocc):
            for a in range(nvirt):
                A[m, n] = (
                    (epsilon[a + nocc] - epsilon[i]) * delta(a, b) * delta(i, j)
                    + 2 * ovov[j, b, i, a]
                    - oovv[j, i, b, a]
                )
                B[m, n] = -2 * ovov[j, b, i, a] + ovov[j, a, i, b]
                n += 1
        m += 1

E2[:ndim, :ndim] = A
E2[:ndim, ndim:] = B
E2[ndim:, :ndim] = B
E2[ndim:, ndim:] = A

Two observations can be made:

matrix block $A$ is diagonal dominant with diagonal elements approximately equal to orbital energy differences (in particular for core electron transitions)
elements of $B$ are comparatively small

core, homo, lumo = 0, 4, 5  # Python indexing of orbitals

print(f"{'1s-LUMO':>48s}{'HOMO-LUMO':>17s}")
print(" " * 41 + "=" * 25)
print(
    "Diagonal elements of A-block of E2   :"
    + f"{A[0,0]:14.8f}"
    + f"{A[(nocc-1)*nvirt,(nocc-1)*nvirt]:14.8f}"
)
print(
    "Associated orbital energy differences:"
    + f"{epsilon[lumo] - epsilon[core]:14.8f}"
    + f"{epsilon[lumo] - epsilon[homo]:14.8f}\n"
)

print(f"Maximum element of B-block of E2 : {np.max(B):12.8f}")

                                         1s-LUMO        HOMO-LUMO
                                         =========================
Diagonal elements of A-block of E2   :   20.34024945    0.36669368
Associated orbital energy differences:   20.76493334    0.70960177

Maximum element of B-block of E2 :   0.21128143

We ensure that our result equals the electronic Hessian obtained with the get_e2() method of the LinearResponseEigenSolver class.

lres_drv = vlx.LinearResponseEigenSolver()
lres_drv.ostream.mute()

E2_ref = lres_drv.get_e2(molecule, basis, scf_results)

np.testing.assert_allclose(E2, E2_ref, atol=1.0e-6)

Property gradient¶

We determine the property gradient

\mathbf{B}^{[1]} = \begin{pmatrix} \mathbf{g} \\ \mathbf{g}^* \\ \end{pmatrix}

(19)

for the operator $\hat{\mu}_z$ .

dipole_mats = vlx.compute_electric_dipole_integrals(molecule, basis)

mu_z_ao = dipole_mats[2]

C = scf_results["C_alpha"]

mu_z = np.einsum("ap, ab, bq -> pq", C, mu_z_ao, C)

B1 = np.zeros(2 * ndim, dtype=complex)  # imaginary valued

g = []
for j in range(nocc):
    for b in range(nocc, norb):
        g.append(1j * np.sqrt(2) * mu_z[j, b])
g = np.array(g)

B1[:ndim] = g
B1[ndim:] = g.conj()

We ensure that our result equals the reference property gradient obtained with the get_prop_grad() method of the LinearResponseSolver class.

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

B1_ref = lrs_drv.get_prop_grad(
    "electric_dipole", "z", molecule, basis, scf_drv.scf_tensors
)[0]

np.testing.assert_allclose(B1.imag, B1_ref, atol=1.0e-12)

Response vector¶

We determine the response vector

\mathbf{N} = - \Big[ \mathbf{E}^{[2]} \Big]^{-1}\mathbf{B}^{[1]}

(20)

N = -np.matmul(np.linalg.inv(E2), B1)

We ensure that our result equals the reference response vector obtained with the get_full_vector() method.

lrs_drv.a_operator = "electric dipole"
lrs_drv.b_operator = "electric dipole"

lrs_drv.a_components = ["z"]
lrs_drv.b_components = ["z"]

lrs_drv.frequencies = [0.0]

lrs_results = lrs_drv.compute(molecule, basis, scf_results)

N_ref = np.zeros(2 * ndim)

N_ref[:ndim] = -lrs_results["solutions"][("z", 0.0)].get_full_vector(1)
N_ref[ndim:] = lrs_results["solutions"][("z", 0.0)].get_full_vector(1)

np.testing.assert_allclose(N.imag, N_ref, atol=1.0e-6)

We note that the imaginary “ $i$ ” is left out in the response vector provided by VeloxChem.

Response function¶

\langle \langle \hat{A}; \hat{B} \rangle \rangle = - \Big[\mathbf{A}^{[1]}\Big]^\dagger \Big[ \mathbf{E}^{[2]} \Big]^{-1} \mathbf{B}^{[1]}

(21)

A1 = B1  # property gradient of the observable

print(f"Response function: {np.dot(A1.conj(), N).real:10.6f}")

Response function:  -4.298489

print(f"Reference value: {lrs_results['response_functions'][('z', 'z', 0.0)]:10.6f}")

Reference value:  -4.298489

We note that our result is in perfect agreement with the reference value obtained from the VeloxChem program.