-
Hartree Fock Theory (notebook)
-
Second Quantization (notebook)
-
Many-Body Perturbation Theory (notebook)
-
Configuration Interaction with Single Excitations (notebook)
-
Static Response Properties (notebook)
Static Molecular Response Properties¶
In this project, we explore the response of a molecular system to a static external electric field, which is related to molecular properties such as the electric dipole moment, polarizability, and hyperpolarizability. Consider the following expression where we have expanded the energy of a molecular system as a Taylor series in applied static external electric field ($\epsilon$) as
$$\begin{align} E = E_0 - \mu_\alpha \epsilon_\alpha- \frac{1}{2} \alpha_{\alpha\beta} \epsilon_\alpha \epsilon_\beta - \frac{1}{6} \beta_{\alpha\beta\gamma} \epsilon_\alpha \epsilon_\beta\epsilon_\gamma + ... \end{align}$$
Here, $E_0$ is the unperturbed energy of the system, $\mu_\alpha$ represents cartesian components of the electric dipole moment, $\alpha_{\alpha\beta}$, $\xi_{\alpha\beta}$ is an element of the the electric dipole polarizability tensor, and $\beta_{\alpha\beta\gamma}$ is an element of the first hyperpolarizability tensor. Note that here and throughout this tutorial, repeated labels imply summation. Note also that this expression could be generalized to account for static external magnetic fields as well, in which case terms involving the magnetizability, NMR shielding tensor, etc. would appear. From this energy expression, we can see that the static response properties can be evaluated as derivatives of the energy with respect to the external field, i.e.,
$$\begin{align} \mu_\alpha &= -\left .\frac{d E}{d \epsilon_\alpha}\right |_{\epsilon = 0}\\ \alpha_{\alpha\beta} &= -\left .\frac{d^2 E}{d \epsilon_\alpha d \epsilon_\beta}\right |_{\epsilon = 0}\\ \beta_{\alpha\beta\gamma} &= -\left .\frac{d^3 E}{d \epsilon_\alpha d \epsilon_\beta d \epsilon_\gamma}\right |_{\epsilon = 0} \end{align}$$
Below, we will derive expressions and implement Python code to evaluate the static electric dipole polarizability for a molecular system described at the Hartree-Fock level of theory. All of the required molecular integrals will be obtained from the Psi4 package, and we will generate the required equations and Python expressions with the help of the p$^\dagger$q package.
The Perturbed Hartree-Fock Wave Function¶
At the Hartree-Fock level of theory, the $N$-electron wave function is approximated as an antisymmetrized product of $N$ one-electron functions called molecular orbitals (MOs), which is called a Slater determinant that we will label as $|\Phi\rangle$. When considering the response of this wave function to an external perturbation, it is convenient to define a perturbed wave function, which can be achieved via the introduction of orbital rotation parametes, $\kappa_{ia}$, as
$$\begin{align} |\Phi({\kappa})\rangle = \text{exp}\left (\kappa_{ia} ( \hat{a}^\dagger_a \hat{a}_i - \hat{a}^\dagger_i \hat{a}_a ) \right ) | \Phi\rangle \end{align}$$
The Hamiltonian corresponding to the system interacting with an external electric field is
$$\begin{align} \hat{H} = \hat{H}_0 - \hat{\mu}_\alpha \epsilon_\alpha \end{align}$$
where $\hat{H}_0$ represents the unperturbed Hamiltonian, and $\hat{\mu}_\alpha$ is a component of the electric dipole operator. To obtain expressions of the static response properties, we proceed by taking derivatives of the energy expression involving this Hamiltonian and the perturbed Hartree-Fock wave function,
$$\begin{align} E(\kappa, \epsilon) = \langle \Phi(\kappa)|\hat{H}_0 - \hat{\mu}_\alpha \epsilon_\alpha | \Phi(\kappa)\rangle \end{align}$$
Note that the energy depends explicitly (through the Hamiltonian) and implicitly (through the wave function parameters) on the field, and we should account for this fact via the chain rule when we differentiate, i.e.,
$$\begin{align} \frac{d}{d\epsilon_\alpha} = \frac{\partial }{\partial \epsilon_\alpha} + \frac{\partial \kappa_{ia}}{\partial \epsilon_\alpha}\frac{\partial }{\partial \kappa_{ia}} \end{align}$$
At this point, it is also useful to go ahead and consider the derivative of the perturbed wave function with respect to the external perturbation.
$$\begin{align} \left . \frac{d}{d \epsilon_\alpha} |\Phi({\kappa})\rangle\right |_{\epsilon=0} &= \left .\frac{\partial}{\partial \epsilon_a}|\Phi({\kappa})\rangle\right |_{\epsilon=0} + \frac{\partial |\Phi({\kappa})\rangle }{\partial \kappa_{ia}} \left . \frac{\partial \kappa_{ia}}{\partial \epsilon_\alpha} \right |_{\epsilon=0}\\ &= \left . ( \hat{a}^\dagger_a \hat{a}_i - \hat{a}^\dagger_i \hat{a}_a )|\Phi({\kappa_0})\rangle\frac{\partial \kappa_{ia}}{\partial \epsilon_\alpha} \right |_{\epsilon=0} \end{align}$$ and $$\begin{align} \left . \frac{d}{d \epsilon_\alpha} \langle \Phi({\kappa})| \right |_{\epsilon=0}&= \left . \left ( \frac{d}{d \epsilon_\alpha} |\Phi({\kappa})\rangle \right )^\dagger \right |_{\epsilon=0}\\ &= -\langle \Phi(\kappa_0)| ( \hat{a}^\dagger_a \hat{a}_i - \hat{a}^\dagger_i \hat{a}_a )\left . \frac{\partial \kappa_{ia}}{\partial \epsilon_\alpha}\right |_{\epsilon=0} \end{align}$$
Here, the partial derivative of the Hartree-Fock wave function with respect to the field parameter vanishes because there is no explicit field dependence in this term, and the notation $\kappa_0$ refers to the orbital rotation parameters determined at zero field, which are the optimal Hartree-Fock parameters.
Expressions for the Hartree-Fock Response Properties¶
Now, the dipole moment is given by
$$\begin{align} \mu_\alpha &= -\left .\frac{d E(\kappa, \epsilon)}{d \epsilon_\alpha}\right |_{\epsilon = 0} \\ &= -\left . \left (\frac{\partial }{\partial \epsilon_\alpha} + \frac{\partial \kappa_{ia}}{\partial \epsilon_\alpha}\frac{\partial }{\partial \kappa_{ia}} \right ) \langle \Phi(\kappa)|\hat{H}_0 - \hat{\mu}_\alpha \epsilon_\alpha | \Phi(\kappa)\rangle\right |_{\epsilon = 0} \\ %&= -\left .\frac{\partial}{\partial \epsilon_\alpha} \langle \Phi(\kappa)|\hat{H}_0| \Phi(\kappa)\rangle\right |_{\epsilon = 0} +\left .\frac{\partial}{\partial \epsilon_\alpha} \langle \Phi(\kappa)| \hat{\mu}_\alpha \epsilon_\alpha | \Phi(\kappa)\rangle\right |_{\epsilon = 0} \\ &= \langle \Phi(\kappa_0)|( \hat{a}^\dagger_a \hat{a}_i - \hat{a}^\dagger_i \hat{a}_a )\hat{H}_0| \Phi(\kappa_0)\rangle\left . \frac{\partial \kappa_{ia}}{\partial \epsilon_\alpha}\right |_{\epsilon=0} - \langle \Phi(\kappa_0)|\hat{H}_0 ( \hat{a}^\dagger_a \hat{a}_i - \hat{a}^\dagger_i \hat{a}_a )|\Phi(\kappa_0)\rangle\left . \frac{\partial \kappa_{ia}}{\partial \epsilon_\alpha}\right |_{\epsilon=0} + \langle \Phi(\kappa_0)|\hat{\mu}|\Phi(\kappa_0)\rangle \\ &= -\langle \Phi(\kappa_0)|[\hat{H}_0,( \hat{a}^\dagger_a \hat{a}_i - \hat{a}^\dagger_i \hat{a}_a )]| \Phi(\kappa_0)\rangle\left . \frac{\partial \kappa_{ia}}{\partial \epsilon_\alpha}\right |_{\epsilon=0} + \langle \Phi(\kappa_0)|\hat{\mu}_\alpha|\Phi(\kappa_0)\rangle \\ &= \langle \Phi(\kappa_0)|\hat{\mu}_\alpha|\Phi(\kappa_0)\rangle \end{align}$$ where, in the last line, we have recognized that the commutator expression corresponds to the orbital gradient, which is zero for the optimal Hartree-Fock parameters. So, we have found that the electric dipole moment is simply given by the expectation value of the dipole operator with respect to the Hartree-Fock wave function (with the optimal field-free parameters).
Now, let us proceed to the case of the dipole polarizability, which is related to the second derivative of the energy with respect to the external field. We have
$$\begin{align} \alpha_{\alpha\beta} &= -\left .\frac{d^2 E(\kappa, \epsilon)}{d \epsilon_\alpha d \epsilon_\beta}\right |_{\epsilon = 0} \\ &= \left .\left (\frac{\partial }{\partial \epsilon_\beta} + \frac{\partial \kappa_{ia}}{\partial \epsilon_\beta}\frac{\partial }{\partial \kappa_{ia}} \right ) \langle \Phi(\kappa)|\hat{\mu}_\alpha|\Phi(\kappa)\rangle \right |_{\epsilon=0}\\ &= \langle \Phi(\kappa_0)|[\hat{\mu}_\alpha,(\hat{a}^\dagger_a \hat{a}_i - \hat{a}^\dagger_i \hat{a}_a)]|\Phi(\kappa_0)\rangle \left .\frac{\partial \kappa_{ia}}{\partial \epsilon_\beta} \right |_{\epsilon=0} \\ &= 2 \mu_{ia}^{\alpha} \left .\frac{\partial \kappa_{ia}}{\partial \epsilon_\beta} \right |_{\epsilon=0} \end{align}$$
where
$$\begin{align} \mu_{ia}^{\alpha} = \langle \Phi(\kappa_0)|\hat{\mu}_\alpha \hat{a}^\dagger_a \hat{a}_i | \Phi(\kappa_0)\rangle = \langle \Phi(\kappa_0)|\hat{a}^\dagger_i \hat{a}_a\hat{\mu}_\alpha | \Phi(\kappa_0)\rangle \end{align}$$
Oh, no! In order to determine the electric dipole polarizability, we must first determine the first-order response of the wave function parameters to the perturbation, $\left .\frac{\partial \kappa_{ia}}{\partial \epsilon_\beta} \right |_{\epsilon=0}$! We can find an equation to determine this quantity from the orbital gradient for the perturbed Hamiltonian, which should be zero at zero field, i.e.,
$$\begin{align} 0 = \left . \langle \Phi(\kappa_0)|[\hat{H}_0 - \mu_\alpha \epsilon_\alpha,( \hat{a}^\dagger_a \hat{a}_i - \hat{a}^\dagger_i \hat{a}_a )]| \Phi(\kappa_0)\rangle \right |_{\epsilon=0} \end{align}$$
We differentiate this expression with respect to the external field to obtain
$$\begin{align} 0 &= \left . \frac{d}{d \epsilon_\beta} \langle \Phi(\kappa)|[\hat{H}_0 - \mu_\alpha \epsilon_\alpha,( \hat{a}^\dagger_a \hat{a}_i - \hat{a}^\dagger_i \hat{a}_a )]| \Phi(\kappa)\rangle \right |_{\epsilon=0} \\ &= \left . \left ( \frac{\partial}{\partial \epsilon_\beta} + \frac{\partial \kappa_{jb}}{\partial \epsilon_\beta} \frac{\partial }{\partial \kappa_{jb}}\right ) \langle \Phi(\kappa)|[\hat{H}_0 - \mu_\alpha \epsilon_\alpha,( \hat{a}^\dagger_a \hat{a}_i - \hat{a}^\dagger_i \hat{a}_a )]| \Phi(\kappa)\rangle \right |_{\epsilon=0} \\ &= \langle \Phi(\kappa_0)|[[\hat{H}_0,( \hat{a}^\dagger_a \hat{a}_i - \hat{a}^\dagger_i \hat{a}_a )],(\hat{a}^\dagger_b \hat{a}_j - \hat{a}^\dagger_j \hat{a}_b)]| \Phi(\kappa_0)\rangle \left . \frac{\partial \kappa_{jb}}{\partial \epsilon_\beta} \right |_{\epsilon = 0} - \langle \Phi(\kappa_0)|[\mu_\alpha \epsilon_\alpha,( \hat{a}^\dagger_a \hat{a}_i - \hat{a}^\dagger_i \hat{a}_a )]| \Phi(\kappa_0)\rangle \end{align}$$
or
$$\begin{align} \langle \Phi(\kappa_0)|[[\hat{H}_0,( \hat{a}^\dagger_a \hat{a}_i - \hat{a}^\dagger_i \hat{a}_a )],(\hat{a}^\dagger_b \hat{a}_j - \hat{a}^\dagger_j \hat{a}_b)]| \Phi(\kappa_0)\rangle \left . \frac{\partial \kappa_{jb}}{\partial \epsilon_\beta} \right |_{\epsilon = 0} = 2 \mu_{ia}^\alpha \end{align}$$
So, in order to determine the first-order response of the wave function, we must solve a set of linear equations involving the orbital hessian (the double commutator expression) and the dipole integrals. Defining
$$\begin{align} h_{jb, ia} = \langle \Phi(\kappa_0)|[[\hat{H}_0,( \hat{a}^\dagger_a \hat{a}_i - \hat{a}^\dagger_i \hat{a}_a )],(\hat{a}^\dagger_b \hat{a}_j - \hat{a}^\dagger_j \hat{a}_b)]| \Phi(\kappa_0)\rangle \end{align}$$
and assuming we can store and invert this tensor, we could determine the response vector as
$$\begin{align} \left . \frac{\partial \kappa_{jb}}{\partial \epsilon_\beta} \right |_{\epsilon = 0} = 2 h_{jb, ia} ^{-1} \mu_{ia}^\alpha \end{align}$$
Here, we have expressed all quantities in a basis of spin orbitals. If we explicitly account for spin, then we have two response equations for $\alpha$- and $\beta$-spin:
$$\begin{align} \left . \frac{\partial \kappa_{j_\sigma b_\sigma}}{\partial \epsilon_\beta} \right |_{\epsilon = 0} = 2 h^{-1}_{j_\sigma b_\sigma, i_\tau a_\tau} \mu_{i_\tau a_\tau}^\alpha \end{align}$$
We can use the p$^\dagger$q package to evaluate the different spin blocks of $h_{jb,ia}$ using the following code.
import pdaggerq
from pdaggerq.parser import contracted_strings_to_tensor_terms
pq = pdaggerq.pq_helper('fermi')
pq.add_double_commutator( 1.0, ['f'], ['a*(a)', 'a(i)'], ['a*(b)', 'a(j)'])
pq.add_double_commutator(-1.0, ['f'], ['a*(i)', 'a(a)'], ['a*(b)', 'a(j)'])
pq.add_double_commutator(-1.0, ['f'], ['a*(a)', 'a(i)'], ['a*(j)', 'a(b)'])
pq.add_double_commutator( 1.0, ['f'], ['a*(i)', 'a(a)'], ['a*(j)', 'a(b)'])
pq.add_double_commutator( 1.0, ['v'], ['a*(a)', 'a(i)'], ['a*(b)', 'a(j)'])
pq.add_double_commutator(-1.0, ['v'], ['a*(i)', 'a(a)'], ['a*(b)', 'a(j)'])
pq.add_double_commutator(-1.0, ['v'], ['a*(a)', 'a(i)'], ['a*(j)', 'a(b)'])
pq.add_double_commutator( 1.0, ['v'], ['a*(i)', 'a(a)'], ['a*(j)', 'a(b)'])
pq.simplify()
for s1 in ['a', 'b']:
for s2 in ['a', 'b']:
spin_labels = {
'j': s1,
'b': s1,
'i': s2,
'a': s2,
}
terms = pq.strings(spin_labels = spin_labels)
terms = contracted_strings_to_tensor_terms(terms)
print('')
print('# %s-%s spin block' % (s1, s2))
print('')
for term in terms:
print('#', term)
print("%s" % (term.einsum_string(update_val='h_'+s1+s2, output_variables=('j', 'b', 'i', 'a'))))
print()
# a-a spin block
# -1.00 d(a,b)*f_aa(j,i)
h_aa += -1.00 * einsum('ab,ji->jbia', kd_aa[va, va], f_aa[oa, oa])
# 1.00 d(i,j)*f_aa(a,b)
h_aa += 1.00 * einsum('ij,ab->jbia', kd_aa[oa, oa], f_aa[va, va])
# -1.00 d(b,a)*f_aa(i,j)
h_aa += -1.00 * einsum('ba,ij->jbia', kd_aa[va, va], f_aa[oa, oa])
# 1.00 d(j,i)*f_aa(b,a)
h_aa += 1.00 * einsum('ji,ba->jbia', kd_aa[oa, oa], f_aa[va, va])
# -1.00 <j,i||a,b>_aaaa
h_aa += -1.00 * einsum('jiab->jbia', g_aaaa[oa, oa, va, va])
# 1.00 <j,a||b,i>_aaaa
h_aa += 1.00 * einsum('jabi->jbia', g_aaaa[oa, va, va, oa])
# 1.00 <i,b||a,j>_aaaa
h_aa += 1.00 * einsum('ibaj->jbia', g_aaaa[oa, va, va, oa])
# -1.00 <a,b||j,i>_aaaa
h_aa += -1.00 * einsum('abji->jbia', g_aaaa[va, va, oa, oa])
# a-b spin block
# 1.00 <j,i||b,a>_abab
h_ab += 1.00 * einsum('jiba->jbia', g_abab[oa, ob, va, vb])
# 1.00 <j,a||b,i>_abab
h_ab += 1.00 * einsum('jabi->jbia', g_abab[oa, vb, va, ob])
# 1.00 <b,i||j,a>_abab
h_ab += 1.00 * einsum('bija->jbia', g_abab[va, ob, oa, vb])
# 1.00 <b,a||j,i>_abab
h_ab += 1.00 * einsum('baji->jbia', g_abab[va, vb, oa, ob])
# b-a spin block
# 1.00 <i,j||a,b>_abab
h_ba += 1.00 * einsum('ijab->jbia', g_abab[oa, ob, va, vb])
# 1.00 <a,j||i,b>_abab
h_ba += 1.00 * einsum('ajib->jbia', g_abab[va, ob, oa, vb])
# 1.00 <i,b||a,j>_abab
h_ba += 1.00 * einsum('ibaj->jbia', g_abab[oa, vb, va, ob])
# 1.00 <a,b||i,j>_abab
h_ba += 1.00 * einsum('abij->jbia', g_abab[va, vb, oa, ob])
# b-b spin block
# -1.00 d(a,b)*f_bb(j,i)
h_bb += -1.00 * einsum('ab,ji->jbia', kd_bb[vb, vb], f_bb[ob, ob])
# 1.00 d(i,j)*f_bb(a,b)
h_bb += 1.00 * einsum('ij,ab->jbia', kd_bb[ob, ob], f_bb[vb, vb])
# -1.00 d(b,a)*f_bb(i,j)
h_bb += -1.00 * einsum('ba,ij->jbia', kd_bb[vb, vb], f_bb[ob, ob])
# 1.00 d(j,i)*f_bb(b,a)
h_bb += 1.00 * einsum('ji,ba->jbia', kd_bb[ob, ob], f_bb[vb, vb])
# -1.00 <j,i||a,b>_bbbb
h_bb += -1.00 * einsum('jiab->jbia', g_bbbb[ob, ob, vb, vb])
# 1.00 <j,a||b,i>_bbbb
h_bb += 1.00 * einsum('jabi->jbia', g_bbbb[ob, vb, vb, ob])
# 1.00 <i,b||a,j>_bbbb
h_bb += 1.00 * einsum('ibaj->jbia', g_bbbb[ob, vb, vb, ob])
# -1.00 <a,b||j,i>_bbbb
h_bb += -1.00 * einsum('abji->jbia', g_bbbb[vb, vb, ob, ob])
In the code snippet, 'f' and 'v' refer to the Fock operator ($\hat{f}$) and fluctuation potential ($\hat{v}$), respectively, and $\hat{H} = \hat{f} + \hat{v}$. In the output, 'd(a,b)' is a Kronecker delta function, 'f(a,b)' is an element of the Fock matrix, and '<j,a||b,i>' refers to an antisymmetrized electron repulsion integral (ERI) in physicists' notation.
At this point, we have all of the ingredients we need to compute the electric dipole polarizability. The following code runs a Hartree-Fock calculation using Psi4 so that we can extract all of the required integrals.
import psi4
import numpy as np
from numpy import einsum
# set molecule
mol = psi4.geometry("""
0 1
O 0.00000 0.00000 0.11572
H 0.00000 0.74879 -0.46288
H 0.00000 -0.74879 -0.46288
symmetry c1
""")
psi4.set_options({'basis': 'cc-pvdz',
'scf_type': 'pk',
'e_convergence': 1e-12,
'd_convergence': 1e-12})
psi4.core.be_quiet()
# compute the Hartree-Fock energy and wavefunction
scf_e, wfn = psi4.energy('SCF', return_wfn=True)
Note that we asked the energy routine to return a wavefunction. This object contains all of the important information from the Hartree-Fock computation, including the number of electrons, the number of orbitals, the AO/MO transformation matrices, and more! The code below grabs everything we will need to evaluate the polarizability.
# Grab data from wavfunction
# all orbitals
Ca = wfn.Ca()
Cb = wfn.Cb()
# use Psi4's MintsHelper to generate ERIs and dipole integrals
mints = psi4.core.MintsHelper(wfn.basisset())
# build the integrals in chemists' notation
g_aaaa = np.asarray(mints.mo_eri(Ca, Ca, Ca, Ca))
g_aabb = np.asarray(mints.mo_eri(Ca, Ca, Cb, Cb))
g_bbbb = np.asarray(mints.mo_eri(Cb, Cb, Cb, Cb))
# antisymmetrized integrals in physicists' notation
g_aaaa = g_aaaa.transpose(0, 2, 1, 3) - g_aaaa.transpose(0, 2, 3, 1)
g_bbbb = g_bbbb.transpose(0, 2, 1, 3) - g_bbbb.transpose(0, 2, 3, 1)
g_abab = g_aabb.transpose(0, 2, 1, 3)
noa = wfn.nalpha()
nob = wfn.nbeta()
oa = slice(None, noa)
ob = slice(None, nob)
va = slice(noa, None)
vb = slice(nob, None)
nmo = wfn.nmo()
kd_aa = np.eye(nmo)
kd_bb = np.eye(nmo)
# orbital energies
f_aa = np.diag(wfn.epsilon_a())
f_bb = np.diag(wfn.epsilon_b())
# dipole integrals
mu = mints.ao_dipole()
mu_x = np.asarray(mu[0])
mu_y = np.asarray(mu[1])
mu_z = np.asarray(mu[2])
# transform the dipole integrals to the MO basis
Ca = np.asarray(Ca)
Cb = np.asarray(Cb)
mua_x = Ca.T @ mu_x @ Ca
mua_y = Ca.T @ mu_y @ Ca
mua_z = Ca.T @ mu_z @ Ca
mub_x = Cb.T @ mu_x @ Cb
mub_y = Cb.T @ mu_y @ Cb
mub_z = Cb.T @ mu_z @ Cb
Now, we can evaluate the matrix representation of the double commutator using the code generated by p$^\dagger$q. Note that we will need to modify the autogenerated code slightly to initialize each spin block of $h_{jb,ia}$.
# a-a spin block
# -1.00 d(a,b)*f_aa(j,i)
h_aa = -1.00 * einsum('ab,ji->jbia', kd_aa[va, va], f_aa[oa, oa])
# 1.00 d(i,j)*f_aa(a,b)
h_aa += 1.00 * einsum('ij,ab->jbia', kd_aa[oa, oa], f_aa[va, va])
# -1.00 d(b,a)*f_aa(i,j)
h_aa += -1.00 * einsum('ba,ij->jbia', kd_aa[va, va], f_aa[oa, oa])
# 1.00 d(j,i)*f_aa(b,a)
h_aa += 1.00 * einsum('ji,ba->jbia', kd_aa[oa, oa], f_aa[va, va])
# -1.00 <j,i||a,b>_aaaa
h_aa += -1.00 * einsum('jiab->jbia', g_aaaa[oa, oa, va, va])
# 1.00 <j,a||b,i>_aaaa
h_aa += 1.00 * einsum('jabi->jbia', g_aaaa[oa, va, va, oa])
# 1.00 <i,b||a,j>_aaaa
h_aa += 1.00 * einsum('ibaj->jbia', g_aaaa[oa, va, va, oa])
# -1.00 <a,b||j,i>_aaaa
h_aa += -1.00 * einsum('abji->jbia', g_aaaa[va, va, oa, oa])
# a-b spin block
# 1.00 <j,i||b,a>_abab
h_ab = 1.00 * einsum('jiba->jbia', g_abab[oa, ob, va, vb])
# 1.00 <j,a||b,i>_abab
h_ab += 1.00 * einsum('jabi->jbia', g_abab[oa, vb, va, ob])
# 1.00 <b,i||j,a>_abab
h_ab += 1.00 * einsum('bija->jbia', g_abab[va, ob, oa, vb])
# 1.00 <b,a||j,i>_abab
h_ab += 1.00 * einsum('baji->jbia', g_abab[va, vb, oa, ob])
# b-a spin block
# 1.00 <i,j||a,b>_abab
h_ba = 1.00 * einsum('ijab->jbia', g_abab[oa, ob, va, vb])
# 1.00 <a,j||i,b>_abab
h_ba += 1.00 * einsum('ajib->jbia', g_abab[va, ob, oa, vb])
# 1.00 <i,b||a,j>_abab
h_ba += 1.00 * einsum('ibaj->jbia', g_abab[oa, vb, va, ob])
# 1.00 <a,b||i,j>_abab
h_ba += 1.00 * einsum('abij->jbia', g_abab[va, vb, oa, ob])
# b-b spin block
# -1.00 d(a,b)*f_bb(j,i)
h_bb = -1.00 * einsum('ab,ji->jbia', kd_bb[vb, vb], f_bb[ob, ob])
# 1.00 d(i,j)*f_bb(a,b)
h_bb += 1.00 * einsum('ij,ab->jbia', kd_bb[ob, ob], f_bb[vb, vb])
# -1.00 d(b,a)*f_bb(i,j)
h_bb += -1.00 * einsum('ba,ij->jbia', kd_bb[vb, vb], f_bb[ob, ob])
# 1.00 d(j,i)*f_bb(b,a)
h_bb += 1.00 * einsum('ji,ba->jbia', kd_bb[ob, ob], f_bb[vb, vb])
# -1.00 <j,i||a,b>_bbbb
h_bb += -1.00 * einsum('jiab->jbia', g_bbbb[ob, ob, vb, vb])
# 1.00 <j,a||b,i>_bbbb
h_bb += 1.00 * einsum('jabi->jbia', g_bbbb[ob, vb, vb, ob])
# 1.00 <i,b||a,j>_bbbb
h_bb += 1.00 * einsum('ibaj->jbia', g_bbbb[ob, vb, vb, ob])
# -1.00 <a,b||j,i>_bbbb
h_bb += -1.00 * einsum('abji->jbia', g_bbbb[vb, vb, ob, ob])
With each spin block defined, we are almost ready to invert the double commutator matrix, but we should note that the spin blocks are coupled, i.e.,
$$\begin{align} h = \begin{pmatrix} h_{j_\alpha b_\alpha, i_\alpha a_\alpha} & h_{j_\alpha b_\alpha, i_\beta a_\beta} \\ h_{j_\beta b_\beta, i_\alpha a_\alpha} & h_{j_\beta b_\beta, i_\beta a_\beta} \\ \end{pmatrix} \end{align}$$
As such, we should invert the super matrix, rather than any of the individual blocks. We could then unpack the spin blocks of the inverted matrix, but it will be easier to work with the inverted super matrix.
# number of alpha- and beta-spin virtual orbitals
nva = nmo - noa
nvb = nmo - nob
# reshape tensors
h_aa = h_aa.reshape(noa*nva, noa*nva)
h_ab = h_ab.reshape(noa*nva, nob*nvb)
h_ba = h_ba.reshape(nob*nvb, noa*nva)
h_bb = h_bb.reshape(nob*nvb, nob*nvb)
# pack into super matrix
h = np.block([[h_aa, h_ab], [h_ba, h_bb]])
# invert super matrix
hinv = np.linalg.inv(h)
Given the inverted double commutator matrix and the dipole integrals in the MO basis, we can finally evaluate the response vectors and the elements of the dipole polarizability tensor. Recall, the response vectors are given by
$$\begin{align} \left . \frac{\partial \kappa_{j_\sigma b_\sigma}}{\partial \epsilon_\beta} \right |_{\epsilon = 0} = 2 h^{-1}_{j_\sigma b_\sigma, i_\tau a_\tau} \mu_{i_\tau a_\tau}^\alpha \end{align}$$
and the dipole polarizability is
$$\begin{align} \alpha_{\alpha\beta} = 2 \mu_{i_\sigma a_\sigma}^{\alpha} \left .\frac{\partial \kappa_{i_\sigma a_\sigma}}{\partial \epsilon_\beta} \right |_{\epsilon=0} \end{align}$$
where we have added spin labels in the polarizability expression. The following code evaluates the response vectors and the elements of the polarizability tensor. Note that matrix-vector products and dot products are most easily evaluated if we construct vectors that contain the relevant elements the alpha-spin and beta-spin dipole integral matrices.
# combine alpha- and beta-spin dipole integrals into a single vector the length of hinv
mu_x_vec = np.hstack((mua_x[oa, va].flatten(), mub_x[ob, vb].flatten()))
mu_y_vec = np.hstack((mua_y[oa, va].flatten(), mub_y[ob, vb].flatten()))
mu_z_vec = np.hstack((mua_z[oa, va].flatten(), mub_z[ob, vb].flatten()))
# response vectors
kappa_x = 2 * einsum('pq,q->p', hinv, mu_x_vec)
kappa_y = 2 * einsum('pq,q->p', hinv, mu_y_vec)
kappa_z = 2 * einsum('pq,q->p', hinv, mu_z_vec)
print('')
print(' ==> Static Dipole Polarizability <==')
print('')
alpha_xx = 2 * einsum('p,p->', mu_x_vec, kappa_x)
alpha_xy = 2 * einsum('p,p->', mu_x_vec, kappa_y)
alpha_xz = 2 * einsum('p,p->', mu_x_vec, kappa_z)
alpha_yy = 2 * einsum('p,p->', mu_y_vec, kappa_y)
alpha_yz = 2 * einsum('p,p->', mu_y_vec, kappa_z)
alpha_zz = 2 * einsum('p,p->', mu_z_vec, kappa_z)
print(' alpha(xx) = %8.4f' % (alpha_xx))
print(' alpha(xy) = %8.4f' % (alpha_xy))
print(' alpha(xz) = %8.4f' % (alpha_xz))
print(' alpha(yy) = %8.4f' % (alpha_yy))
print(' alpha(yz) = %8.4f' % (alpha_yz))
print(' alpha(zz) = %8.4f' % (alpha_zz))
print('')
==> Static Dipole Polarizability <==
alpha(xx) = 3.0444
alpha(xy) = -0.0000
alpha(xz) = -0.0000
alpha(yy) = 6.6932
alpha(yz) = -0.0000
alpha(zz) = 4.9785
Compare these values to the computed Hartree-Fock/cc-pVDZ values provided on the Computational Chemistry Comparison and Benchmark DataBase. How did we do?? Be sure to double check the units.
Iterative Solution of the Response Equations¶
Inverting the Hessian matrix is computationally demanding, requiring $\mathcal{O}(o^3v^3)$ floating-point operators, where $o$ and $v$ represent the number of occupied and unoccupied orbitals, respectively. Since the cost of solving for the Hartree-Fock energy scales only with the fourth power of system size, we should probably avoid the use of any operations in the corresponding response theory that exhibit worse scaling. Here, we solve
$$\begin{align} \left . h_{j_\sigma b_\sigma, i_\tau a_\tau} \frac{\partial \kappa_{j_\sigma b_\sigma}}{\partial \epsilon_\beta} \right |_{\epsilon = 0} = 2 \mu_{i_\tau a_\tau}^\alpha \end{align}$$
iteratively for the response vector, using the Jacobi method. Given an initial guess for the response vector, we obtain an updated guess from
$$\begin{align} \left . \frac{\partial \kappa_{i_\sigma a_\sigma}}{\partial \epsilon_\beta} \right |_{\epsilon = 0} = h^{-1}_{i_\sigma a_\sigma, i_\sigma a_\sigma} \left ( 2 \mu_{i_\tau a_\tau}^\alpha - h^\prime_{j_\sigma b_\sigma, i_\tau a_\tau} \left . \frac{\partial \kappa_{j_\sigma b_\sigma}}{\partial \epsilon_\beta} \right |_{\epsilon = 0} \right ) \end{align}$$
where the prime denotes that the sum excludes the diagonal elements of the hessian. In the code below, note that we have dampened the iterations because they do not converge otherwise.
kappa_x = np.zeros_like(mu_x_vec)
kappa_y = np.zeros_like(mu_y_vec)
kappa_z = np.zeros_like(mu_z_vec)
old_kappa_x = np.zeros_like(mu_x_vec)
old_kappa_y = np.zeros_like(mu_y_vec)
old_kappa_z = np.zeros_like(mu_z_vec)
# strip diagonals out of the hessian
h_nodiag = h.copy()
for ia in range (0, len(h)):
h_nodiag[ia, ia] = 0.0
h_diag = h.diagonal()
damp = 0.5
print('')
print(' ==> First-Order Response Equations <==')
print('')
print(' iter |dk|')
for i in range (0, 100):
new_kappa_x = 1.0 / h_diag * ( mu_x_vec - np.einsum('ij,j->i', h_nodiag, kappa_x) )
new_kappa_y = 1.0 / h_diag * ( mu_y_vec - np.einsum('ij,j->i', h_nodiag, kappa_y) )
new_kappa_z = 1.0 / h_diag * ( mu_z_vec - np.einsum('ij,j->i', h_nodiag, kappa_z) )
kappa_x = damp * old_kappa_x + (1.0 - damp) * new_kappa_x
kappa_y = damp * old_kappa_y + (1.0 - damp) * new_kappa_y
kappa_z = damp * old_kappa_z + (1.0 - damp) * new_kappa_z
conv = np.linalg.norm(kappa_x - old_kappa_x)
conv += np.linalg.norm(kappa_y - old_kappa_y)
conv += np.linalg.norm(kappa_z - old_kappa_z)
old_kappa_x = kappa_x.copy()
old_kappa_y = kappa_y.copy()
old_kappa_z = kappa_z.copy()
print(' %4i %12.8f' % (i, conv))
if conv < 1e-8 :
break
if i == 99 :
print('')
print(' warning! response equations did not converge!')
print('')
else:
print('')
print(' response equations converged!')
print('')
alpha_xx = 4 * np.einsum('i,i->', mu_x_vec, kappa_x)
alpha_xy = 4 * np.einsum('i,i->', mu_x_vec, kappa_y)
alpha_xz = 4 * np.einsum('i,i->', mu_x_vec, kappa_z)
alpha_yy = 4 * np.einsum('i,i->', mu_y_vec, kappa_y)
alpha_yz = 4 * np.einsum('i,i->', mu_y_vec, kappa_z)
alpha_zz = 4 * np.einsum('i,i->', mu_z_vec, kappa_z)
print(' alpha(xx) = %8.4f' % (alpha_xx))
print(' alpha(xy) = %8.4f' % (alpha_xy))
print(' alpha(xz) = %8.4f' % (alpha_xz))
print(' alpha(yy) = %8.4f' % (alpha_yy))
print(' alpha(yz) = %8.4f' % (alpha_yz))
print(' alpha(zz) = %8.4f' % (alpha_zz))
print('')
==> First-Order Response Equations <==
iter |dk|
0 1.42781193
1 0.59140965
2 0.28247428
3 0.14133161
4 0.07329898
5 0.03896359
6 0.02111298
7 0.01161663
8 0.00647116
9 0.00364119
10 0.00206562
11 0.00117962
12 0.00067731
13 0.00039063
14 0.00022612
15 0.00013129
16 0.00007643
17 0.00004458
18 0.00002605
19 0.00001525
20 0.00000894
21 0.00000524
22 0.00000308
23 0.00000181
24 0.00000106
25 0.00000063
26 0.00000037
27 0.00000022
28 0.00000013
29 0.00000008
30 0.00000004
31 0.00000003
32 0.00000002
33 0.00000001
response equations converged!
alpha(xx) = 3.0444
alpha(xy) = -0.0000
alpha(xz) = -0.0000
alpha(yy) = 6.6932
alpha(yz) = -0.0000
alpha(zz) = 4.9785
We get the same answer as when explicitly inverting the hessian! Great! Now, we also not that it may not always be desirable to even store the hessian matrix. Fortunately, the hessian matrix itself is not really necessary. All we need is the action of the hessian on the response vector,
$$\begin{align} \sigma_{i_\tau a_\tau} = h_{j_\sigma b_\sigma, i_\tau a_\tau} \left . \frac{\partial \kappa_{j_\sigma b_\sigma}}{\partial \epsilon_\beta} \right |_{\epsilon = 0} \end{align}$$
Recall that
$$\begin{align} h_{jb, ia} = \langle \Phi(\kappa_0)|[[\hat{H}_0,( \hat{a}^\dagger_a \hat{a}_i - \hat{a}^\dagger_i \hat{a}_a )],(\hat{a}^\dagger_b \hat{a}_j - \hat{a}^\dagger_j \hat{a}_b)]| \Phi(\kappa_0)\rangle \end{align}$$
so we can write
$$\begin{align} \sigma_{i_\tau a_\tau} = \langle \Phi(\kappa_0)|[[\hat{H}_0,( \hat{a}^\dagger_a \hat{a}_i - \hat{a}^\dagger_i \hat{a}_a )],(\hat{a}^\dagger_b \hat{a}_j - \hat{a}^\dagger_j \hat{a}_b)]| \Phi(\kappa_0)\rangle \left . \frac{\partial \kappa_{j_\sigma b_\sigma}}{\partial \epsilon_\beta} \right |_{\epsilon = 0} \end{align}$$
We can generate equations for this matrix vector product using the unitary coupled-cluster (CC) functionality in the p$^\dagger$q. It turns out that the singles amplitudes in unitary CC theory have a similar structure as the response vector, i.e., $\hat{t}_1 = t_{bj} (\hat{a}^\dagger_b \hat{a}_j - \hat{a}^\dagger_j \hat{a}_b)$. If we define
$$\begin{align} t_{bj} = \left .\frac{\partial \kappa_{j_\sigma b_\sigma}}{\partial \epsilon_\beta} \right |_{\epsilon = 0} \end{align}$$
then we can write
$$\begin{align} \sigma_{i_\tau a_\tau} = \langle \Phi(\kappa_0)|[[\hat{H}_0,( \hat{a}^\dagger_a \hat{a}_i - \hat{a}^\dagger_i \hat{a}_a )],\hat{t}_1]| \Phi(\kappa_0)\rangle \end{align}$$
The following code generates einsum expressions for $\sigma_{i_\tau a_\tau}$.
import pdaggerq
from pdaggerq.parser import contracted_strings_to_tensor_terms
pq = pdaggerq.pq_helper("fermi")
# we need to use the antihermitian cluster operator from UCC
pq.set_unitary_cc(True)
for ham in ['f', 'v']:
pq.add_double_commutator(1.0, [ham], ['a*(a)', 'a(i)'], ['t1'])
pq.add_double_commutator(-1.0, [ham], ['a*(i)', 'a(a)'], ['t1'])
pq.simplify()
for s1 in ['a', 'b']:
spin_labels = {
'i' : s1,
'a' : s1,
}
terms = pq.strings(spin_labels = spin_labels)
terms = contracted_strings_to_tensor_terms(terms)
for term in terms:
print("%s" % (term.einsum_string(update_val='sigma_' + s1 + s1, output_variables=('a', 'i'))))
pq.clear()
sigma_aa += -1.00 * einsum('ij,aj->ai', f_aa[oa, oa], t1_aa)
sigma_aa += 1.00 * einsum('ba,bi->ai', f_aa[va, va], t1_aa)
sigma_aa += -1.00 * einsum('ji,aj->ai', f_aa[oa, oa], t1_aa)
sigma_aa += 1.00 * einsum('ab,bi->ai', f_aa[va, va], t1_aa)
sigma_aa += -1.00 * einsum('jiab,bj->ai', g_aaaa[oa, oa, va, va], t1_aa)
sigma_aa += 1.00 * einsum('ijab,bj->ai', g_abab[oa, ob, va, vb], t1_bb)
sigma_aa += 1.00 * einsum('ibaj,bj->ai', g_aaaa[oa, va, va, oa], t1_aa)
sigma_aa += 1.00 * einsum('ibaj,bj->ai', g_abab[oa, vb, va, ob], t1_bb)
sigma_aa += 1.00 * einsum('jabi,bj->ai', g_aaaa[oa, va, va, oa], t1_aa)
sigma_aa += 1.00 * einsum('ajib,bj->ai', g_abab[va, ob, oa, vb], t1_bb)
sigma_aa += -1.00 * einsum('abji,bj->ai', g_aaaa[va, va, oa, oa], t1_aa)
sigma_aa += 1.00 * einsum('abij,bj->ai', g_abab[va, vb, oa, ob], t1_bb)
sigma_bb += -1.00 * einsum('ij,aj->ai', f_bb[ob, ob], t1_bb)
sigma_bb += 1.00 * einsum('ba,bi->ai', f_bb[vb, vb], t1_bb)
sigma_bb += -1.00 * einsum('ji,aj->ai', f_bb[ob, ob], t1_bb)
sigma_bb += 1.00 * einsum('ab,bi->ai', f_bb[vb, vb], t1_bb)
sigma_bb += 1.00 * einsum('jiba,bj->ai', g_abab[oa, ob, va, vb], t1_aa)
sigma_bb += -1.00 * einsum('jiab,bj->ai', g_bbbb[ob, ob, vb, vb], t1_bb)
sigma_bb += 1.00 * einsum('bija,bj->ai', g_abab[va, ob, oa, vb], t1_aa)
sigma_bb += 1.00 * einsum('ibaj,bj->ai', g_bbbb[ob, vb, vb, ob], t1_bb)
sigma_bb += 1.00 * einsum('jabi,bj->ai', g_abab[oa, vb, va, ob], t1_aa)
sigma_bb += 1.00 * einsum('jabi,bj->ai', g_bbbb[ob, vb, vb, ob], t1_bb)
sigma_bb += 1.00 * einsum('baji,bj->ai', g_abab[va, vb, oa, ob], t1_aa)
sigma_bb += -1.00 * einsum('abji,bj->ai', g_bbbb[vb, vb, ob, ob], t1_bb)
The code below defines a function implementing these expressions.
def Hx(t1_aa, t1_bb):
"""
evaluate contraction of Hessian matrix against a t1-shaped operator
:param t1_aa: the alpha-spin part of t1 (the response vector)
:param t1_bb: the beta-spin part of t1 (the response vector)
:return sigma_aa
:return sigma_bb
"""
sigma_aa = -1.00 * einsum('ij,aj->ai', f_aa[oa, oa], t1_aa)
sigma_aa += 1.00 * einsum('ba,bi->ai', f_aa[va, va], t1_aa)
sigma_aa += -1.00 * einsum('ji,aj->ai', f_aa[oa, oa], t1_aa)
sigma_aa += 1.00 * einsum('ab,bi->ai', f_aa[va, va], t1_aa)
sigma_aa += -1.00 * einsum('jiab,bj->ai', g_aaaa[oa, oa, va, va], t1_aa)
sigma_aa += 1.00 * einsum('ijab,bj->ai', g_abab[oa, ob, va, vb], t1_bb)
sigma_aa += 1.00 * einsum('ibaj,bj->ai', g_aaaa[oa, va, va, oa], t1_aa)
sigma_aa += 1.00 * einsum('ibaj,bj->ai', g_abab[oa, vb, va, ob], t1_bb)
sigma_aa += 1.00 * einsum('jabi,bj->ai', g_aaaa[oa, va, va, oa], t1_aa)
sigma_aa += 1.00 * einsum('ajib,bj->ai', g_abab[va, ob, oa, vb], t1_bb)
sigma_aa += -1.00 * einsum('abji,bj->ai', g_aaaa[va, va, oa, oa], t1_aa)
sigma_aa += 1.00 * einsum('abij,bj->ai', g_abab[va, vb, oa, ob], t1_bb)
sigma_bb = -1.00 * einsum('ij,aj->ai', f_bb[ob, ob], t1_bb)
sigma_bb += 1.00 * einsum('ba,bi->ai', f_bb[vb, vb], t1_bb)
sigma_bb += -1.00 * einsum('ji,aj->ai', f_bb[ob, ob], t1_bb)
sigma_bb += 1.00 * einsum('ab,bi->ai', f_bb[vb, vb], t1_bb)
sigma_bb += 1.00 * einsum('jiba,bj->ai', g_abab[oa, ob, va, vb], t1_aa)
sigma_bb += -1.00 * einsum('jiab,bj->ai', g_bbbb[ob, ob, vb, vb], t1_bb)
sigma_bb += 1.00 * einsum('bija,bj->ai', g_abab[va, ob, oa, vb], t1_aa)
sigma_bb += 1.00 * einsum('ibaj,bj->ai', g_bbbb[ob, vb, vb, ob], t1_bb)
sigma_bb += 1.00 * einsum('jabi,bj->ai', g_abab[oa, vb, va, ob], t1_aa)
sigma_bb += 1.00 * einsum('jabi,bj->ai', g_bbbb[ob, vb, vb, ob], t1_bb)
sigma_bb += 1.00 * einsum('baji,bj->ai', g_abab[va, vb, oa, ob], t1_aa)
sigma_bb += -1.00 * einsum('abji,bj->ai', g_bbbb[vb, vb, ob, ob], t1_bb)
return sigma_aa, sigma_bb
Now, we can solve the response equations without the need to explicitly store the entire hessian matrix. Note that we have rewritten this solver to use containers that are the same shape as that used for t1 in p$^\dagger$q.
# dipole integrals shaped as [v, o]
mu_a = []
mu_a.append(mua_x[va, oa])
mu_a.append(mua_y[va, oa])
mu_a.append(mua_z[va, oa])
mu_b = []
mu_b.append(mub_x[vb, ob])
mu_b.append(mub_y[vb, ob])
mu_b.append(mub_z[vb, ob])
# response vectors shaped as [v, o]
kappa_a = []
kappa_b = []
old_kappa_a = []
old_kappa_b = []
for xyz in range (0, 3):
kappa_a.append(np.zeros_like(mu_a[xyz]))
kappa_b.append(np.zeros_like(mu_b[xyz]))
old_kappa_a.append(np.zeros_like(mu_a[xyz]))
old_kappa_b.append(np.zeros_like(mu_b[xyz]))
# diagonal of hessian shaped as [v, o]
h_diag_a = h_diag[:noa*nva].reshape(noa, nva).T
h_diag_b = h_diag[noa*nva:].reshape(nob, nvb).T
print('')
print(' ==> First-Order Response Equations <==')
print('')
print(' iter ||dk||')
for i in range (0, 100):
conv = 0.0
for xyz in range (0, 3):
sigma_a, sigma_b = Hx(kappa_a[xyz], kappa_b[xyz])
new_kappa_a = 1.0 / h_diag_a * ( mu_a[xyz] - sigma_a + h_diag_a * kappa_a[xyz])
new_kappa_b = 1.0 / h_diag_b * ( mu_b[xyz] - sigma_b + h_diag_b * kappa_b[xyz])
kappa_a[xyz] = damp * old_kappa_a[xyz] + (1.0 - damp) * new_kappa_a
kappa_b[xyz] = damp * old_kappa_b[xyz] + (1.0 - damp) * new_kappa_b
conv += np.linalg.norm(kappa_a[xyz] - old_kappa_a[xyz])
conv += np.linalg.norm(kappa_b[xyz] - old_kappa_b[xyz])
old_kappa_a[xyz] = kappa_a[xyz].copy()
old_kappa_b[xyz] = kappa_b[xyz].copy()
print(' %4i %12.8f' % (i, conv))
if conv < 1e-8 :
break
if i == 99 :
print('')
print(' warning! response equations did not converge!')
print('')
else:
print('')
print(' response equations converged!')
print('')
cart = ['x', 'y', 'z']
for i in range (0, 3):
for j in range (i, 3):
alpha = 4 * np.einsum('ai,ai->', mu_a[i], kappa_a[j])
alpha += 4 * np.einsum('ai,ai->', mu_b[i], kappa_b[j])
print(' alpha(%s%s) = %8.4f' % (cart[i], cart[j], alpha))
print('')
==> First-Order Response Equations <==
iter ||dk||
0 2.01923099
1 0.83637955
2 0.39947896
3 0.19987308
4 0.10366042
5 0.05510283
6 0.02985827
7 0.01642840
8 0.00915160
9 0.00514943
10 0.00292122
11 0.00166823
12 0.00095787
13 0.00055244
14 0.00031979
15 0.00018568
16 0.00010808
17 0.00006305
18 0.00003685
19 0.00002157
20 0.00001264
21 0.00000742
22 0.00000436
23 0.00000256
24 0.00000151
25 0.00000089
26 0.00000052
27 0.00000031
28 0.00000018
29 0.00000011
30 0.00000006
31 0.00000004
32 0.00000002
33 0.00000001
34 0.00000001
response equations converged!
alpha(xx) = 3.0444
alpha(xy) = -0.0000
alpha(xz) = -0.0000
alpha(yy) = 6.6932
alpha(yz) = -0.0000
alpha(zz) = 4.9785