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}
\]
And the property gradient
\[\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*}\]
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*}\]
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()
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(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
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.36669362
Associated orbital energy differences: 20.76493335 0.70960170
Maximum element of B-block of E2 : 0.14823088
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
\[\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()
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]}
\]
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.