# 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$

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

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

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]}$

In order, we have introduced:

The Hermitian electronic Hessian

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

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

The (unknown) response vector

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

with

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

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

with

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

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

and let us adopt the electric-dipole approximation

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

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

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

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

and

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

and, in this case, the block dimension equals

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

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

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

$\begin{split} \mathbf{E}^{[2]} = \begin{pmatrix} \mathbf{A} & \mathbf{B} \\ \mathbf{B}^* & \mathbf{A}^* \\ \end{pmatrix} \end{split}$
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


• 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.36669362
Associated orbital energy differences:   20.76493335    0.70960170

Maximum element of B-block of E2 :   0.21128138


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)


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

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

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

mu_z_ao = -1.0 * dipole_mats.z_to_numpy()

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

"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]}$
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]}$
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.