Molecular Hessians - eChem book

This section deals with the calculation of second derivatives of the electronic energy with respect to nuclear coordinates. While the numerical calculation of the Hessian based on an analytical gradient follows the same procedure as for the numerical gradient, the analytical calculation is in general much more involved than for first derivatives. After some special aspects of the numerical calculation of second derivatives, the calculation of the analytic Hessian at the Hartree–Fock (HF) or self-consistent field (SCF) level of theory is described as an example, followed by the required terms for methods based on density functional theory (DFT). Furthermore, the calculation of other second derivatives of the energy is described briefly, taking the example of electric field perturbations.

Numerical Hessians¶

If the analytical gradient is not available (and both first and second derivatives are needed), it is more efficient to resort to an equation for the second derivative that depends on the energy $E$ only, rather than employing one of the equations here with the energy itself replaced by its numerical gradient. Assuming the energy depends on (at least) two distinct nuclear coordinates $x$ and $y$ , the second derivative of $E$ can be calculated as

\begin{align} \frac{\mathrm{d}^2 E}{\mathrm{d} x \mathrm{d} y} &\approx \frac{E(x + x_0, y + y_0) - E(x + x_0, y) - E(x, y + y_0)}{2 x_0 y_0}\\ & + \frac{+ 2 E(x, y) - E(x - x_0, y) - E(x, y - y_0) + E(x - x_0, y - y_0)}{2 x_0 y_0} \end{align}

(1)

where $x_0$ and $y_0$ are some small changes in the respective variable. For the case that the derivative is taken with respect to the same nuclear coordinate, the formula simplifies to

\begin{align} \frac{\mathrm{d}^2 E}{\mathrm{d} x^2} \approx \frac{E(x + x_0, y) - 2 E(x, y) + E(x - x_0, y)}{x_0^2} \end{align}

(2)

Analytical Hessians¶

The calculation of analytical second derivatives of the energy is much more involved than for first derivatives, even for variational methods. We restrict ourselves to SCF methods here, i.e., Hartree–Fock (HF) and density functional theory (DFT). The difficulty comes from the fact that the derivative of the MO coefficients $\mathbf{C}$ with respect to the perturbation $x$ cannot be avoided anymore, leading to the so-called coupled-perturbed SCF equations, as shown in the following.

HF¶

The analytic gradient of the HF energy $E_{\mathrm{HF}}$ with respect to a nuclear coordinate $x$ , which was previously given in MO basis, can equivalently be written in AO basis as Pople et al. (1979)

\frac{\mathrm{d} E_{\text{HF}}}{\mathrm{d} x} = \sum_{\mu \nu} D_{\mu \nu} h_{\mu \nu}^{x} + \frac12 \sum_{\mu \nu \kappa \lambda} D_{\mu \nu} D_{\kappa \lambda} \langle \mu \kappa || \nu \lambda \rangle^x + \sum_{\mu \nu} \omega_{\mu \nu} S_{\mu \nu}^{x} + \frac{\mathrm{d} V_{nn}}{\mathrm{d} x}

(3)

where the HF density matrix $\mathbf{D}$ and energy-weighted density matrix $\boldsymbol{\omega}$ in a real AO basis are defined as

\begin{align} D_{\mu \nu} &= \sum_{i} C_{\mu i} C_{\nu i} \\ \omega_{\mu \nu} &= - \sum_{i} \varepsilon_i C_{\mu i} C_{\nu i} \end{align}

(4)

Straightforward differentiation of the above equation with respect to another nuclear coordinate $y$ gives second derivatives of the HF energy

\begin{align} \frac{\mathrm{d}^2 E_{\text{HF}}}{\mathrm{d} y \mathrm{d} x} &= \sum_{\mu \nu} D_{\mu \nu} h_{\mu \nu}^{x y} + \frac12 \sum_{\mu \nu \lambda \sigma} D_{\mu \nu} D_{\lambda \sigma} \langle \mu \lambda || \nu \sigma \rangle^{x y} + \sum_{\mu \nu} \omega_{\mu \nu} S_{\mu \nu}^{x y} + \frac{\mathrm{d}^2 V_{nn}}{\mathrm{d} y \mathrm{d} x} \\ &+ \sum_{\mu \nu} \frac{\mathrm{d} D_{\mu \nu}}{\mathrm{d} y} h_{\mu \nu}^{x} + \sum_{\mu \nu \lambda \sigma} \frac{\mathrm{d} D_{\mu \nu}}{\mathrm{d} y} D_{\lambda \sigma} \langle \mu \lambda || \nu \sigma \rangle^{x} + \sum_{\mu \nu} \frac{\mathrm{d} \omega_{\mu \nu}}{\mathrm{d} y} S_{\mu \nu}^{x} \end{align}

(5)

Now, the calculation of the perturbed density, which boils down to the calculation of the perturbed MO coefficients $C_{\mu p}^{x}$ , can no longer be avoided Pople et al. (1979). This leads to the so-called coupled-perturbed Hartree–Fock (CPHF) equations, as will be described in the following.

CPHF equations¶

To derive the coupled-perturbed equation, the starting point is HF

\sum_{\nu} F_{\mu \nu} (\mathbf{C}) C_{\nu i} = \varepsilon_i \sum_{\nu} S_{\mu \nu} C_{\nu i}

(6)

which can be differentiated and rearranged, such that one obtains Neese (2009)

\sum_{\nu} \bigg[ \frac{\mathrm{d} F_{\mu \nu} (\mathbf{C})}{\mathrm{d} x} C_{\nu i} + \Big( F_{\mu \nu} (\mathbf{C}) - \varepsilon_i S_{\mu \nu} \Big) \frac{\mathrm{d} C_{\nu i}}{\mathrm{d} x} \bigg] = \frac{\mathrm{d} \varepsilon_i}{\mathrm{d} x} \sum_{\nu} S_{\mu \nu} C_{\nu i} + \varepsilon_i \sum_{\nu} \frac{\mathrm{d} S_{\mu \nu}}{\mathrm{d} x} C_{\nu i}

(7)

Next, an ansatz for the perturbed MO coefficients $C_{\mu p}^{x} = \frac{\mathrm{d} C_{\mu p}}{\mathrm{d} x}$ needs to be made. Since the unperturbed MOs span the same space as the AOs, but form an orthonormal set, the perturbed MO coefficients are expanded in the basis of the unperturbed ones Pople et al. (1979)

C_{\mu p}^{x} = \sum_{q} C_{\mu q} U_{q p}^{x}

(8)

where the matrix $\mathbf{U}^x$ contains the unknown expansion coefficients. By taking the derivative of the orthonormality condition $\mathbf{C}^\text{T} \mathbf{SC} = \mathbf{1}$ , or in subscript notation Neese (2009)

\sum_{\mu \nu} C_{\mu p} S_{\mu \nu} C_{\nu q} = \delta_{pq}

(9)

it can be shown that the occupied-occupied block of $\mathbf{U}^x$ is related to the partial derivative of the overlap matrix

U_{ij}^{x} = - \frac12 S_{ij}^{(x)}

(10)

The virtual-virtual block $U_{ab}^{x}$ is not required, and hence only the occupied-virtual block $U_{ai}^x$ needs to be determined. The total derivative of the Fock matrix $F_{\mu \nu}^{x} = \frac{\mathrm{d} F_{\mu \nu} (\mathbf{C})}{\mathrm{d} x}$ is given by

F_{\mu \nu}^{x} = h_{\mu \nu}^{x} + \sum_{\kappa \lambda} P_{\kappa \lambda} \langle \mu \kappa || \nu \lambda \rangle^{x} + \sum_{qi} \sum_{\kappa \lambda} (U_{qi}^{x} C_{\kappa q} C_{\lambda i} + U_{qi}^{x} C_{\kappa i} C_{\lambda q} ) \langle \mu \kappa || \nu \lambda \rangle

(11)

where it can be seen that the derivative of the Fock matrix also depends on the unknowns $\mathbf{U}^{x}$ . The coupled-perturbed Hartree--Fock (CPHF) equations are obtained by multiplying the derivative of the HF equations from the left by $C_{\mu a}$ and summing over all $\mu$ Neese (2009). The CPHF equations, whose solution yields the $U_{ai}^{x}$ , can be written as Deglmann et al. (2002)

\begin{align*} (\varepsilon_i - \varepsilon_a) U_{ai}^{x} - 2 G_{ai}[U_{bj}^{x}] = R_{ai}^{x} \end{align*}

(12)

where the right-hand side (RHS), which depends on the perturbation $x$ , is given by

R_{ai}^{x} = F_{ai}^{(x)} - \varepsilon_{i} S_{ai}^{(x)} - G_{ai}[S_{jk}^{(x)}]

with the partial derivative of the Fock matrix being given by

F_{pq}^{(x)} = h_{pq}^{(x)} + \sum_{j} \langle pj||qj \rangle^{(x)}

(13)

and the matrices $\mathbf{G}[M_{rs}^{x}]$ are defined for arbitrary $M_{rs}^{x}$ as Deglmann et al. (2002)

G_{pq}[M_{rs}^{x}] = \sum_{\mu \nu} C_{\mu p} C_{\nu q} \bigg[ \sum_{\kappa \lambda} \langle \mu \kappa || \nu \lambda \rangle \Big( \sum_{rs} C_{\kappa r} M_{rs}^{x} C_{\lambda s} \Big) \bigg]

(14)

It should be noted that in contrast to the first derivative of the MP2 or CIS energy, where only a single orbital-response equation had to be solved, a response equation needs to be solved for every perturbation $x$ , i.e., the $3N$ Cartesian nuclear coordinates of a molecule containing $N$ atoms, to obtain the second derivative of the HF energy.

Perturbed density¶

Having solved for the CPHF coefficients $\mathbf{U}^{x}$ , the perturbed density matrix $\mathbf{D}^x$ can be calculated by using the chosen ansatz as

\begin{align} \frac{\mathrm{d} D_{\mu \nu}}{\mathrm{d} x} &= \sum_{i} \frac{\mathrm{d}}{\mathrm{d} x} (C_{\mu i} C_{\nu i}) = \sum_{i} ( C_{\mu i}^x C_{\nu i} + C_{\mu i} C_{\nu i}^x ) \\ &= - \frac12 \sum_{ij} ( C_{\mu i} C_{\nu j} S_{ij}^{(x)} + C_{\mu j} C_{\nu i} S_{ij}^{(x)} ) + \sum_{ia} ( C_{\mu a} C_{\nu i} + C_{\mu i} C_{\nu a} ) U_{ai}^{x} \\ \end{align}

(15)

As an alternative to this equation, the second derivatives of the HF energy can then be written in terms of the RHS $R_{ai}^{x}$ , the CPHF coefficients $U_{ai}^{x}$ , and partial derivatives of the Fock and overlap matrices as Deglmann et al. (2002)

\begin{align} \frac{\mathrm{d}^2 E_{\text{HF}}}{\mathrm{d} x \mathrm{d} y} &= \sum_{\mu \nu} D_{\mu \nu} h_{\mu \nu}^{x y} + \frac12 \sum_{\mu \nu \lambda \sigma} D_{\mu \nu} D_{\lambda \sigma} \langle \mu \lambda || \nu \sigma \rangle^{x y} + \sum_{\mu \nu} \omega_{\mu \nu} S_{\mu \nu}^{x y} + \frac{\mathrm{d}^2 V_{nn}}{\mathrm{d} x \mathrm{d} y} \\ &+ 2 \sum_{ia} R_{ai}^{x} U_{ai}^{y} - \sum_{ij} \Big( F_{ij}^{(x)} S_{ij}^{(y)} + F_{ij}^{(y)} S_{ij}^{(x)} - 2 \varepsilon_i S_{ij}^{(x)} S_{ij}^{(y)} - 2 G_{ij}[S_{kl}^{(x)}] S_{ij}^{(y)} \Big) \, , \end{align}

(16)

where the explicit construction of the perturbed $\mathbf{D}$ and $\boldsymbol{\omega}$ matrices is avoided.

DFT¶

Analogous to the case of the ground-state gradient, only some additional exchange-correlation (xc) contributions need to be considered for the analytic DFT Hessian Deglmann et al. (2002) compared to HF. First of all, the Fock matrix needs to be replaced by the Kohn--Sham matrix $\mathbf{F}$ , such that the partial derivative of the latter for a hybrid functional is given by

F_{pq}^{(x)} = h_{pq}^{(x)} + \sum_{j} \big( \langle pj|qj \rangle^{(x)} - c_{\text{x}} \langle pj|jq \rangle^{(x)} \big) + v_{pq}^{\text{xc}(x)}

(17)

where the derivative of the xc potential $v_{pq}^{\text{xc}}$ is given here, and the matrices $\mathbf{G}[M_{rs}^{x}]$ also need to include a contraction with the xc kernel $f^{\text{xc}}_{\mu\nu\kappa\lambda}$ in the same way as the Coulomb integral.

Finally, the second partial derivative of the xc energy contribution $E_\text{xc}$ with respect to two nuclear coordinates $x$ and $y$ is needed, which is given as

\begin{align} \frac{\partial^2 E_\text{xc}}{\partial y \partial x} &= \frac{\partial}{\partial y} \int \frac{\partial e_{\text{xc}}}{\partial n(\mathbf{r})} n(\mathbf{r})^{(x)} \mathrm{d} \mathbf{r} = \frac{\partial}{\partial \zeta} \int \frac{\partial e_{\text{xc}}}{\partial n(\mathbf{r})} \sum_{\mu\nu} D_{\mu\nu} \frac{\partial (\chi_{\mu} \chi_{\nu})}{\partial x} \mathrm{d} \mathbf{r} \\ &= \frac{\partial}{\partial y} \int \frac{\partial e_{\text{xc}}}{\partial n(\mathbf{r})} \sum_{\mu\nu} D_{\mu\nu} \frac{\partial^2 (\chi_{\mu} \chi_{\nu})}{\partial y \partial x} \mathrm{d} \mathbf{r} + \int \frac{\partial^2 e_{\text{xc}}}{\partial n(\mathbf{r}) \partial n'(\mathbf{r})} n(\mathbf{r})^{(x)} n(\mathbf{r})^{(y)} \mathrm{d} \mathbf{r} \end{align}

(18)

where we made use of the partial derivative of the ground-state density $n(\mathbf{r})^{(x)}$ as given here.

Second-order properties¶

Analogous to the case of first-order properties, some other time-independent molecular properties can be calculated as second derivatives of the energy $E$ . As an example, consider the static dipole polarizability $\boldsymbol{\alpha}$ , which corresponds to the derivative of the dipole moment $\boldsymbol{\mu}$ with respect to an external electric field $\boldsymbol{\mathcal{F}}$

\boldsymbol{\alpha} = \frac{\mathrm{d} \boldsymbol{\mu}}{\mathrm{d} \boldsymbol{\mathcal{F}}} \approx \frac{\boldsymbol{\mu}(\boldsymbol{\mathcal{F}}) - \boldsymbol{\mu}(-\boldsymbol{\mathcal{F}})}{2 \boldsymbol{\mathcal{F}}}

(19)

which can be obtained numerically by calculating the (analytical) dipole moment in presence of an external field in positive and negative direction. Since the dipole moment is already the first derivative of the energy $E$ with respect to the field, the polarizability is actually a second derivative,

\boldsymbol{\alpha} = \frac{\mathrm{d}^2 E}{\mathrm{d} \boldsymbol{\mathcal{F}}^2}

(20)

and as such can be referred to as a second-order property. The analytic calculation at the HF or DFT level procedes analogously to the nuclear Hessian, meaning that the CPHF equations need to be solved for the unknowns $U_{ia}^{\mathcal{F}}$ . However, since the basis functions do not depend on the electric field, the RHS is in this case simply given by a component $\{ x,y,z \}$ of the dipole integrals, $R_{ia}^{\mathcal{F}} = \mu_{ia}$ . As a matter of fact, this equation is then equivalent to the linear response equation with $\omega = 0$ . Having solved for the unknown CPHF coefficients, the perturbed density is obtained from

P_{\lambda \kappa}^{\mathcal{F}} = \frac{\mathrm{d} D_{\lambda \kappa}}{\mathrm{d} \mathcal{F}} = \sum_{ia} \big( C_{\lambda i} C_{\kappa a} + C_{\lambda a} C_{\kappa i} \big) U_{ia}^{\mathcal{F}}

(21)

and the polarizability components can be calculated as

\alpha = \sum_{\kappa \lambda} D_{\lambda \kappa}^{\mathcal{F}} \, \mu_{\kappa \lambda}

(22)

for each component of the perturbed density and dipole integrals. Another example of a second-order property is the derivative of the dipole moment with respect to nuclear coordinates, which is needed for IR intensities. This can be seen as a mixed second energy derivative (once with respect to nuclear coordinates, once with respect to an electric field), and is discussed in more detail here.

References¶

Pople, J. A., Krishnan, R., Schlegel, H. B., & Binkley, J. S. (1979). Derivative studies in Hartree–Fock and Møller–Plesset theories. Int. J. Quantum Chem., 16, 225–241. 10.1002/qua.560160825
Neese, F. (2009). Prediction of molecular properties and molecular spectroscopy with density functional theory: From fundamental theory to exchange-coupling. Coord. Chem. Rev., 253, 526–563. 10.1016/j.ccr.2008.05.014
Deglmann, P., Furche, F., & Ahlrichs, R. (2002). An efficient implementation of second analytical derivatives for density functional methods. Chem. Phys. Lett., 362, 511–518. 10.1016/S0009-2614(02)01084-9