Molecular Hessians

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{split} \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{split}\]

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

\[ \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} \]

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 [PKSB79]

\[%:label: eq:HF_gradient_in_ao \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} \]

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

\[\begin{split} 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{split}\]

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

\[\begin{split}%:label: eq:HF_Hessian_Pople \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{split}\]

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 [PKSB79]. 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} \]

which can be differentiated and rearranged, such that one obtains [Nee09]

\[%:label: eq:derivative_hf_equations \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} \]

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 [PKSB79]

\[%:label: eq:perturbed_mo_coefficients C_{\mu p}^{x} = \sum_{q} C_{\mu q} U_{q p}^{x} \]

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 [Nee09]

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

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

\[\begin{equation*} U_{ij}^{x} = - \frac12 S_{ij}^{(x)} \end{equation*}\]

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

\[%:label: eq:derivative_fock_matrix 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 \]

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\) [Nee09]. The CPHF equations, whose solution yields the \(U_{ai}^{x}\), can be written as [DFA02]

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

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

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

with the partial derivative of the Fock matrix being given by

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

and the matrices \(\mathbf{G}[M_{rs}^{x}]\) are defined for arbitrary \(M_{rs}^{x}\) as [DFA02]

\[\begin{equation*} 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] \, . \end{equation*}\]

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{split}%:label: eq:perturbed_density \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{split}\]

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 [DFA02]

\[\begin{split}%:label: eq:HF_Hessian_Furche \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{split}\]

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 [DFA02] 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)} \]

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{split} \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{split}\]

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

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

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

and the polarizability components can be calculated as

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

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.