Perturbation-dependent Fock spaces

Since we would like to work in a second-quantized framework, we would naturally like to employ a Fock space that can be used to describe the system for every value of that perturbation. Unfortunately, when working with a finite basisset, we cannot assume that a basisset that is able to characterize a system at one value of the perturbation, is adequate for another. This means that, in general, we should work with a perturbation-dependent Fock space. However, following (T. U. Helgaker and Almlöf 1984, Helgaker2012), we can with little effort disguise the perturbation-dependent Fock space as a small complication in the Hamiltonian using spinor connections.

Let us start by summarizing the situation at zero perturbation, denoted by \(\require{physics} \boldsymbol{\eta}_0\). We expand an orthonormal set of \(M = K_\alpha + K_\beta\) spinors in its underlying scalar bases \(\qty{ \chi^\alpha, \chi^\beta }\), as is done in section : \[\begin{equation} \phi_P(\vb{r}; \boldsymbol{\eta}_0) % = % \begin{pmatrix} \displaystyle \sum_\mu^{K_\alpha} % \chi^\alpha_\mu(\vb{r}; \boldsymbol{\eta}_0) % C^{\alpha}_{\mu P}(\boldsymbol{\eta}_0) \\ \displaystyle \sum_\mu^{K_\beta} % \chi^\beta_\mu(\vb{r}; \boldsymbol{\eta}_0) % C^{\beta}_{\mu P}(\boldsymbol{\eta}_0) \end{pmatrix} \thinspace , \end{equation}\] for which orthonormality is expressed as an identity overlap matrix/metric: \[\begin{align} S_{PQ}(\boldsymbol{\eta}_0) % &= \int \dd{\vb{r}} \phi^*_P(\vb{r}; \boldsymbol{\eta}_0) % \phi_Q(\vb{r}; \boldsymbol{\eta}_0) \\ &= \vb{C}^\dagger(\boldsymbol{\eta}_0) % \begin{pmatrix} \vb{S}^{\alpha \alpha} % & \vb{0} \\ \vb{0} % & \vb{S}^{\beta \beta} \\ \end{pmatrix} \vb{C}(\boldsymbol{\eta}_0) \\ &= \delta_{PQ} % \thinspace . \end{align}\]

The perturbation affects the metric

In this case, we assume that the perturbation affects the metric of the underlying scalar bases, which is for example the case when a geometric distortion is applied, or when a magnetic field is applied when using London orbitals. This means that the old spinors \(\qty{ \phi_P(\vb{r}; \boldsymbol{\eta}_0) }\) are no longer orthonormal when the system is perturbed: \[\begin{align} S_{PQ}(\boldsymbol{\eta}) % &= % \int \dd{\vb{r}} % \phi_P^*(\vb{r}, \boldsymbol{\eta}) % \phi_Q(\vb{r}, \boldsymbol{\eta}) \\ &\neq \delta_{PQ} % \thinspace , \end{align}\] or equivalently stated: \[\begin{align} \vb{S}(\boldsymbol{\eta}) % &= \vb{C}^\dagger(\boldsymbol{\eta}) % \begin{pmatrix} \vb{S}^{\alpha \alpha}(\boldsymbol{\eta}) % & \vb{0} \\ \vb{0} % & \vb{S}^{\beta \beta}(\boldsymbol{\eta}) \\ \end{pmatrix} \vb{C}(\boldsymbol{\eta}) \\ &\neq \vb{I} % \thinspace . \end{align}\]

Since we consider the metric \(S^{\sigma \sigma}_{\mu \nu}\) of the underlying scalar bases to be perturbation-dependent, the expansion coefficients of the spinors \(\vb{C}\) must depend on the perturbation in order to maintain the flexibility to construct orthonormal spinors at every value of the perturbation \(\boldsymbol{\eta}\). Using a so-called connection matrix \(\vb{T}(\boldsymbol{\eta})\) (T. Helgaker et al. 2012) that satisfies the relation \(\eqref{eq:transformation_orthonormalization_relation}\): \[\begin{equation} \vb{T}^\dagger(\boldsymbol{\eta}) % \vb{S}(\boldsymbol{\eta}) % \vb{T}(\boldsymbol{\eta}) = \vb{I} % \thinspace , \end{equation}\] we can proceed to orthonormalizing the spinor basis at every value of the perturbation \(\boldsymbol{\eta}\). Different conditions on \(\vb{T}(\boldsymbol{\eta})\) lead to different connection schemes. (T. Helgaker et al. 2012) Since we require the spinors to be orthonormal at no perturbation \(\boldsymbol{\eta}_0\), the spinor connection matrix at zero perturbation must be unitary: \[\begin{equation} \vb{T}(\boldsymbol{\eta}_0)^\dagger % \vb{T}(\boldsymbol{\eta}_0) % = \vb{I} \end{equation}\] and we choose to make it satisfy \[\begin{equation} \vb{T}(\boldsymbol{\eta}_0) = \vb{I} % \thinspace . \end{equation}\] Using the spinor connection \(\vb{T}(\boldsymbol{\eta})\) as a transformation matrix always leads to an orthonormal spinor basis, which we will call the orthonormalized spinors: \[\begin{equation} \tilde{\phi}_P(\vb{r}, \boldsymbol{\eta}) % = % \sum_Q^M % \phi_Q(\vb{r}; \boldsymbol{\eta}) % T_{QP}(\boldsymbol{\eta}) % \thinspace , \end{equation}\] which are the spinors that we will use to construct a Fock space at every value of the perturbation. Since the orthonormalized spinors are by construction orthonormal, we have usual the anticommutation relations at each specific value of the perturbation \(\boldsymbol{\eta}\): \[\begin{align} & % \comm{ % \hat{\tilde{a}}_P(\boldsymbol{\eta}) % }{ % \hat{\tilde{a}}^\dagger_Q(\boldsymbol{\eta}) % }_+ % = \delta_{PQ} \\ & % \comm{ % \hat{\tilde{a}}^\dagger_P(\boldsymbol{\eta}) % }{ % \hat{\tilde{a}}^\dagger_Q(\boldsymbol{\eta}) % }_+ = 0 \thinspace . \label{eq:perturbation_dependent_overlap_matrix} \end{align}\] We should now emphasize that the values of these anticommutators (i.e. the right-hand terms) are in no way dependent of the perturbation if we do not mix up different perturbation values, so we can actually regard the elementary second quantization operators as perturbation-independent! Therefore, we can write the (perturbation-dependent) Hamiltonian as \[\begin{equation} \label{eq:perturbation_dependent_Hamiltonian} \hat{\mathcal{H}}(\boldsymbol{\eta}) % = % \sum_{PQ}^M % \tilde{h}_{PQ}(\boldsymbol{\eta}) % \hat{E}_{PQ} % + \frac{1}{2} \sum_{PQRS}^K % \tilde{g}_{PQRS}(\boldsymbol{\eta}) \thinspace % \hat{e}_{PQRS} % + h_{\text{nuc}}(\boldsymbol{\eta}) % \thinspace , \end{equation}\] with the transformed one-electron integrals \[\begin{equation} \tilde{h}_{PQ}(\boldsymbol{\eta}) % = \sum_{RS}^K % T^*_{RP}(\boldsymbol{\eta}) \thinspace % h_{RS}(\boldsymbol{\eta}) \thinspace % T_{SQ}(\boldsymbol{\eta}) % \thinspace . \end{equation}\] and the transformed two-electron integrals \[\begin{equation} \tilde{g}_{PQRS}(\boldsymbol{\eta}) % = \sum_{TUVW}^K % T^*_{TP}(\boldsymbol{\eta}) \thinspace % T^*_{VR}(\boldsymbol{\eta}) \thinspace % g_{TUVW}(\boldsymbol{\eta}) \thinspace % T_{UQ}(\boldsymbol{\eta}) \thinspace % T_{WS}(\boldsymbol{\eta}) % \thinspace . \end{equation}\] Keeping in mind that for molecular properties, we will need various first-order and second-order derivatives of these integrals, we introduce the following function \[\begin{equation} h'_{PQ}(\boldsymbol{\eta}, t) % = \sum_{RS}^M % \qty[ % \exp( % t \ln{\vb{T}(\boldsymbol{\eta})} % ) % ]^*_{RP} % \thinspace h_{RS}(\boldsymbol{\eta}) \thinspace % \qty[ % \exp( % t \ln{\vb{T}(\boldsymbol{\eta})} % ) % ]_{SQ} % \thinspace . \end{equation}\] At first sight, this may seem a strange construction, but evaluating this function at \(t=0\) yields the original spinors: \[\begin{equation} h'_{PQ}(\boldsymbol{\eta}, t=0) % = h_{PQ}(\boldsymbol{\eta}) \end{equation}\] and evaluating at \(t=1\) yields the orthonormalized spinors: \[\begin{equation} h'_{PQ}(\boldsymbol{\eta}, t=1) % = \tilde{h}_{PQ}(\boldsymbol{\eta}) % \thinspace . \end{equation}\] All in all, this means that we can expand \(h'_{PQ}(\boldsymbol{\eta}, t)\) around \(t=0\) in order to find its value at \(t=1\), leading to: \[\begin{equation} \label{eq:h'_expanded_0_1} h'_{PQ}(\boldsymbol{\eta}, 1) % = h'_{PQ}(\boldsymbol{\eta}, 0) % + % \eval{ % \dv{h'_{PQ}(\boldsymbol{\eta}, t)}{t} }_{t=0} % \cdot 1 % + \frac{1}{2!} % \eval{ % \dv[2]{h'_{PQ}(\boldsymbol{\eta}, t)}{t} }_{t=0} \cdot 1^2 % \thinspace , \end{equation}\] up to second order. I believe it is instructional to go through the full derivation, so we will start with calculating the first derivative of \(h'_{PQ}\) with respect to \(t\): \[\begin{equation} \begin{split} \dv{h'_{PQ}(\boldsymbol{\eta}, t)}{t} % = &\sum_{RS}^M % \qty[ % \exp( % t \ln{\vb{T}(\boldsymbol{\eta})} % ) % \ln{\vb{T}(\boldsymbol{\eta})} ]^*_{RP} % \thinspace h_{RS}(\boldsymbol{\eta}) \thinspace % \qty[ % \exp( % t \ln{\vb{T}(\boldsymbol{\eta})} % ) % ]_{SQ} \\ &+ \sum_{RS}^M % \qty[ % \exp( % t \ln{\vb{T}(\boldsymbol{\eta})} % ) % ]^*_{RP} % \thinspace h_{RS}(\boldsymbol{\eta}) \thinspace % \qty[ % \exp( % t \ln{\vb{T}(\boldsymbol{\eta})} % ) % \ln{\vb{T}(\boldsymbol{\eta})} ]_{SQ} % \thinspace . \end{split} \end{equation}\] Evaluating this derivative at \(t=0\) leads to \[\begin{align} \eval{ % \dv{h'_{PQ}(\boldsymbol{\eta}, t)}{t} % }_{t=0} % = \qty{ % \vb{h}(\boldsymbol{\eta}) % , % \ln{\vb{T}(\boldsymbol{\eta})} % }_{PQ} % \thinspace , \end{align}\] in which we have defined the following one-index transformation for the one-electron integrals \(\vb{f}\) with a transformation matrix \(\vb{T}\): \[\begin{equation} \qty{ \vb{f}, \vb{T} }_{PQ} % = \sum_R^M T^*_{RP} f_{RQ} % + \sum_R^M f_{PR} T_{RQ} % \thinspace . \end{equation}\] This one-index transformation is additive (but not linear) in its arguments: \[\begin{align} & \qty{ \vb{h} + \vb{f}, \vb{T} } % = \qty{ \vb{h}, \vb{T} } % + \qty{ \vb{f}, \vb{T} } % \label{eq:one_index_transformation_additive_left} \\ % & \qty{ \vb{f}, \vb{T} + \vb{R} } % = \qty{ \vb{f}, \vb{T} } % + \qty{ \vb{f}, \vb{R} } % \label{eq:one_index_transformation_additive_right} \thinspace , \end{align}\] and furthermore it is Hermitian if the one-electron integrals \(\vb{f}\) are Hermitian: \[\begin{equation} \qty{ \vb{f}, \vb{T} }^\dagger% = \qty{ \vb{f}, \vb{T} } % \thinspace , \end{equation}\] so the one-index transformation inherits the symmetry properties of the underlying operator.

We then find for the second derivative evaluated at \(t=0\): \[\begin{equation} \eval{ % \dv[2]{h'_{PQ}(\boldsymbol{\eta}, t)}{t} % }_{t=0} = % \qty{ % \qty{ % \vb{h}(\boldsymbol{\eta}) % , % \ln{\vb{T}(\boldsymbol{\eta})} % } % , % \ln{\vb{T}(\boldsymbol{\eta})} % }_{PQ} % \thinspace . \end{equation}\] Helgaker introduces a symmetrized double one-index transformation (T. Helgaker et al. 2012), but we have no need for that here.

Collecting the first and second derivatives, and recognizing the special cases for \(t=0\) and \(t=1\), equation \(\eqref{eq:h'_expanded_0_1}\) then reduces to: \[\begin{equation} \tilde{h}_{PQ}(\boldsymbol{\eta}) % = h_{PQ}(\boldsymbol{\eta}) % + \qty{ % \vb{h}(\boldsymbol{\eta}) % , % \ln{\vb{T}(\boldsymbol{\eta})} % }_{PQ} % + \frac{1}{2!} \qty{ % \qty{ % \vb{h}(\boldsymbol{\eta}) % , % \ln{\vb{T}(\boldsymbol{\eta})} % } % , % \ln{\vb{T}(\boldsymbol{\eta})} % }_{PQ} % \thinspace . \end{equation}\]

We now proceed by expanding the orthonormalized one-electron integrals \(\tilde{h}(\boldsymbol{\eta})\) around \(\boldsymbol{\eta} = \boldsymbol{\eta}_0\) up to second order: \[\begin{equation} \label{eq:expansion_h_tilde_perturbation} \tilde{h}_{PQ}(\boldsymbol{\eta}) % = \tilde{h}_{PQ}(\boldsymbol{\eta}_0) % + \sum_m % \eval{ % \pdv{ % \tilde{h}_{PQ}(\boldsymbol{\eta}) % }{\eta_m} % }_{\boldsymbol{\eta}_0} % \eta_m + \frac{1}{2!} \sum_{mn} % \eval{ % \pdv{ % \tilde{h}_{PQ}(\boldsymbol{\eta}) % }{\eta_m}{\eta_n} % }_{\boldsymbol{\eta}_0} % \eta_m \eta_n \thinspace , \end{equation}\] in which we will now tackle the unknown terms that appear, but in order to do so, we have to derive some more properties concerning the spinor connections. Expanding \(\ln{\vb{T}(\boldsymbol{\eta})}\) around \(\vb{T}(\boldsymbol{\eta}) = \vb{I}\): \[\begin{equation} \ln{\vb{T}(\boldsymbol{\eta})} % = \qty( \vb{T}(\boldsymbol{\eta}) - \vb{I} ) % - \frac{1}{2} % \qty( \vb{T}(\boldsymbol{\eta}) - \vb{I} )^2 % + \frac{1}{3} % \qty( \vb{T}(\boldsymbol{\eta}) - \vb{I} )^3 % - \cdots % \thinspace , \end{equation}\] we can express perturbational derivatives of \(\ln{\vb{T}(\boldsymbol{\eta})}\) in terms of those of related to \(\vb{T}(\boldsymbol{\eta})\). The zeroth-order term vanishes: \[\begin{equation} \ln{\vb{T}(\boldsymbol{\eta}_0)} % = \vb{0} % \thinspace . \end{equation}\] For the first-order perturbational derivative, we have: \[\begin{align} \eval{ % \pdv{ % [ \ln{\vb{T}(\boldsymbol{\eta})} ]_{PQ} % }{\eta_m} % }_{\boldsymbol{\eta}_0} % & = % \eval{ % \pdv{ % T_{PQ}(\boldsymbol{\eta}) % }{\eta_m} % }_{\boldsymbol{\eta}_0} \\ & \equiv % \qty[ % \vb{T}_m(\boldsymbol{\eta}_0) % ]_{PQ} % \end{align}\] and for the second-order perturbational derivative, we find: \[\begin{align} \eval{ % \pdv{ % [ \ln{\vb{T}(\boldsymbol{\eta})} ]_{PQ} % }{\eta_m}{\eta_n} % }_{\boldsymbol{\eta}_0} % & = % \eval{ % \pdv{ % T_{PQ}(\boldsymbol{\eta}) % }{\eta_m}{\eta_n} % }_{\boldsymbol{\eta}_0} % - \frac{1}{2} (1 + P_{mn}) % \sum_R^M % \eval{ % \pdv{ % T_{PR}(\boldsymbol{\eta}) % }{\eta_m} % }_{\boldsymbol{\eta}_0} % \eval{ % \pdv{ % T_{RQ}(\boldsymbol{\eta}) % }{\eta_n} % }_{\boldsymbol{\eta}_0} \\ & \equiv % \qty[ % \vb{T}_{mn}(\boldsymbol{\eta}_0) % ]_{PQ} % \thinspace , \end{align}\] where we have introduced short-hand notation for the first and second derivatives \(\vb{T}_m(\boldsymbol{\eta}_0)\) and \(\vb{T}_{mn}(\boldsymbol{\eta}_0)\).

The zeroth order term in equation \(\eqref{eq:expansion_h_tilde_perturbation}\) is then just given by the original integrals: \[\begin{equation} \tilde{h}_{PQ}(\boldsymbol{\eta}_0) % = h_{PQ}(\boldsymbol{\eta}_0) % \thinspace . \end{equation}\] The first-order term can be calculated as \[\begin{equation} \label{eq:first_order_derivative_one_electron_integrals} \eval{ % \pdv{ % \tilde{h}_{PQ}(\boldsymbol{\eta}) % }{\eta_m} % }_{\boldsymbol{\eta}_0} % = % \eval{ % \pdv{ % h_{PQ}(\boldsymbol{\eta}) % }{\eta_m} % }_{\boldsymbol{\eta}_0} % + % \qty{ % \vb{h}(\boldsymbol{\eta}_0) % , % \vb{T}_m(\boldsymbol{\eta}_0) % }_{PQ} \end{equation}\] by using the property \[\begin{equation} \pdv{ % \qty{ % \vb{f}(x) % , % \vb{T}(x) % } % }{x} % = \qty{ % \vb{f}(x) % , % \pdv{ % \vb{T}(x) % }{x} % } % + \qty{ % \pdv{ % \vb{f}(x) % }{x} % , % \vb{T}(x) % } % \thinspace , \label{eq:one_index_transformation_derivative_property} \end{equation}\] We can then work out the second-order term to become: \[\begin{equation} \begin{split} \eval{ % \pdv{ % \tilde{h}_{PQ}(\boldsymbol{\eta}) % }{\eta_m}{\eta_n} % }_{\boldsymbol{\eta}_0} % = % & \eval{ % \pdv{ % h_{PQ}(\boldsymbol{\eta}) % }{\eta_m}{\eta_n} % }_{\boldsymbol{\eta}_0} % + % \qty{ % \vb{h}(\boldsymbol{\eta}_0) % , % \vb{T}_{mn}(\boldsymbol{\eta}_0) % }_{PQ} \\ & + (1 + P_{mn}) % \qty{ % \eval{ % \qty{ % \vb{h}(\boldsymbol{\eta}) % , % \vb{T}(\boldsymbol{\eta}) % }_m % }_{\boldsymbol{\eta}_0} % , % \vb{T}_n(\boldsymbol{\eta}) % }_{PQ} % \thinspace , \end{split} \label{eq:second_order_derivative_one_electron_integrals} \end{equation}\] where we have introduced a kind of symmetrized derivative of a combined one-index transformation \[\begin{equation} \qty{ % \vb{h}(\boldsymbol{\eta}) % , % \vb{T}(\boldsymbol{\eta}) % }_m % = % \pdv{ % \vb{h}(\boldsymbol{\eta}) % }{\eta_m} % + \frac{1}{2!} % \qty{ % \vb{h}(\boldsymbol{\eta}) % , % \vb{T}_m(\boldsymbol{\eta}) % } % \thinspace . \end{equation}\]

We can proceed analogously for the two-electron integrals \(\tilde{g}_{PQRS}(\boldsymbol{\eta})\) by defining \[\begin{equation} \begin{split} g'_{PQRS}(\boldsymbol{\eta}) % = \sum_{TUVW}^K % & \qty[ % \exp( % t \ln{\vb{T}(\boldsymbol{\eta})} % ) % ]^*_{TP} \thinspace % \qty[ % \exp( % t \ln{\vb{T}(\boldsymbol{\eta})} % ) % ]^*_{VR} \thinspace \\ & g_{TUVW}(\boldsymbol{\eta}) \thinspace % \qty[ % \exp( % t \ln{\vb{T}(\boldsymbol{\eta})} % ) % ]_{UQ} \thinspace % \qty[ % \exp( % t \ln{\vb{T}(\boldsymbol{\eta})} % ) % ]_{WS} \end{split} \end{equation}\] and the combined one-index transformation for the two-electron integrals \(\vb{g}\) with a transformation matrix \(\vb{T}\): \[\begin{equation} \begin{split} \qty{ % \vb{g} % , % \vb{T} % }_{PQRS} % = &\sum_T^M % T^*_{TP} % g_{TQRS} % + \sum_T^M % T^*_{TR} % g_{PQTS} \\ &+ \sum_T^M % g_{PTRS} % T_{TQ} % + \sum_T^M % g_{PQRT} % T_{TS} % \thinspace . \end{split} \end{equation}\] This combined one-index transformation has the same properties (cfr. additivity \(\eqref{eq:one_index_transformation_additive_left}\) and \(\eqref{eq:one_index_transformation_additive_right}\), and a property related to derivatives \(\eqref{eq:one_index_transformation_derivative_property}\)) as the for the one-electron integrals, and also inherits the symmetries of the underlying two-electron integrals \(\vb{g}\): \[\begin{align} & \qty{ % \vb{g} % , % \vb{T} % }^*_{PQRS} % = % \qty{ % \vb{g} % , % \vb{T} % }^*_{QPSR} % \\ & \qty{ % \vb{g} % , % \vb{T} % }_{RSPQ} % = % \qty{ % \vb{g} % , % \vb{T} % }_{PQRS} % \thinspace . \end{align}\] We then obtain the same functional form for the first derivative of \(g'_{PQRS}(\boldsymbol{\eta})\) as we did for the first derivative of the corresponding one-electron integral function: \[\begin{equation} \eval{ % \dv{ % g'_{PQRS}(\boldsymbol{\eta}, t) % }{t} % }_{t=0} % = % \qty{ % \vb{g}(\boldsymbol{\eta}) % , % \ln{\vb{T}(\boldsymbol{\eta})} % }_{PQRS} % \thinspace , \end{equation}\] so we can find expressions for the first-order perturbational derivative of the two-electron integrals \[\begin{equation*} \eval{ % \pdv{ % \tilde{g}_{PQRS}(\boldsymbol{\eta}) % }{\eta_m} % }_{\boldsymbol{\eta}_0} % \end{equation*}\] by replacing \(\tilde{h}\) with \(\tilde{g}\) in equation \(\eqref{eq:first_order_derivative_one_electron_integrals}\) and second-order perturbational derivatives of the two-electron integrals \[\begin{equation} \eval{ % \pdv{ % \tilde{g}_{PQRS}(\boldsymbol{\eta}) % }{\eta_m}{\eta_n} % }_{\boldsymbol{\eta}_0} % \end{equation}\] by replacing \(\tilde{h}\) with \(\tilde{g}\) in equation \(\eqref{eq:second_order_derivative_one_electron_integrals}\) and changing the appropriate indices.

We are now at the point where we have calculated the first- and second-order perturbational partial derivative of the one- and two-electron integrals that compromise the Hamiltonian. We can then bring all of the previous derivations together, to finally find the first-order perturbational derivative of the Hamiltonian: \[\begin{equation} \begin{split} \eval{ % \pdv{ % \hat{\mathcal{H}}(\boldsymbol{\eta}) % }{\eta_m} % }_{\boldsymbol{\eta}_0} % = % & \sum_{PQ}^M % \eval{ % \pdv{% h_{PQ}(\boldsymbol{\eta}) % }{\eta_m} % }_{\boldsymbol{\eta}_0} % \hat{E}_{PQ} % + \frac{1}{2} \sum_{PQRS}^M % \eval{ % \pdv{ % g_{PQRS}(\boldsymbol{\eta}) % }{\eta_m} % }_{\boldsymbol{\eta}_0} % \hat{e}_{PQRS} % + % \eval{ % \pdv{ % h_{\text{nuc}}(\boldsymbol{\eta}) % }{\eta_m} % }_{\boldsymbol{\eta}_0} \\ &+ \sum_{PQ}^M % \qty[ % \vb{T}^\dagger_m(\boldsymbol{\eta}_0) % ]_{PQ} % \hat{\mathscr{F}}_{QP} \qty( % \hat{\mathcal{H}}_{\text{elec}} ( % \boldsymbol{\eta}_0 % ) % ) % + \sum_{PQ}^M \qty[ % \vb{T}_m(\boldsymbol{\eta}_0) % ]_{PQ} % \hat{\mathscr{F}}^\dagger_{QP} \qty( % \hat{\mathcal{H}}_{\text{elec}} ( % \boldsymbol{\eta}_0 % ) % ) % \thinspace , \end{split} \label{eq:first_order_perturbation_derivative_Hamiltonian} \end{equation}\] in which we have used the following properties: \[\begin{equation} \sum_{PQ}^M % \qty{ % \vb{f} % , % \vb{T} % }_{PQ} % \hat{E}_{PQ} % = % \sum_{PQ}^M % \qty[ \vb{T}^\dagger ]_{PQ} % \hat{\mathscr{F}}_{QP} ( \hat{f} ) % + % \sum_{PQ}^M % T_{PQ} % \hat{\mathscr{F}}^\dagger_{QP} ( \hat{f} ) % % \end{equation}\] and \[\begin{equation} \sum_{PQRS}^M % \qty{ % \vb{g} % , % \vb{T} % }_{PQRS} % \hat{e}_{PQRS} % = % \sum_{PQ}^M % \qty[ \vb{T}^\dagger ]_{PQ} % \hat{\mathscr{F}}_{QP} ( \hat{g} ) % + % \sum_{PQ}^M % T_{PQ} % \hat{\mathscr{F}}^\dagger_{QP} ( \hat{g} ) % % \end{equation}\] for Fockians (cfr. section \(\ref{sec:fockians}\)) of one-electron operators \(\hat{f}\) and Hermitian two-electron operators \(\hat{g}\) (including the prefactor \(1/2\)). For the second-order perturbational derivative of the Hamiltonian, we find: \[\begin{equation} \begin{split} \eval{ % \pdv{ % \hat{\mathcal{H}}(\boldsymbol{\eta}) % }{\eta_m}{\eta_n} % }_{\boldsymbol{\eta}_0} % \hspace{-3pt} = % & \sum_{PQ}^M % \eval{ % \pdv{% h_{PQ}(\boldsymbol{\eta}) % }{\eta_m}{\eta_n} % }_{\boldsymbol{\eta}_0} % \hat{E}_{PQ} % + \frac{1}{2} \sum_{PQRS}^M % \eval{ % \pdv{ % g_{PQRS}(\boldsymbol{\eta}) % }{\eta_m}{\eta_n} % }_{\boldsymbol{\eta}_0} % \hat{e}_{PQRS} % + % \eval{ % \pdv{ % h_{\text{nuc}}(\boldsymbol{\eta}) % }{\eta_m}{\eta_n} % }_{\boldsymbol{\eta}_0} \\ &+ \sum_{PQ}^M % \qty[ % \vb{T}^\dagger_{mn}(\boldsymbol{\eta}_0) % ]_{PQ} % \hat{\mathscr{F}}_{QP} \qty( % \hat{\mathcal{H}}_{\text{elec}} ( % \boldsymbol{\eta}_0 % ) % ) \\ &+ \sum_{PQ}^M \qty[ % \vb{T}_{mn}(\boldsymbol{\eta}_0) % ]_{PQ} % \hat{\mathscr{F}}^\dagger_{QP} \qty( % \hat{\mathcal{H}}_{\text{elec}} ( % \boldsymbol{\eta}_0 % ) % ) \\ &+ (1 + P_{mn}) \Bigg[ % \sum_{PQ}^M % \qty[ % \vb{T}^\dagger_n(\boldsymbol{\eta}_0) % ]_{PQ} % \hat{\mathscr{F}}_{QP} \qty( % \eval{ % \qty{ % \widehat{ % \hat{\mathcal{H}}_{\text{elec}} ( % \boldsymbol{\eta} % ) % , % \vb{T}(\boldsymbol{\eta}) % } % }_m % }_{\boldsymbol{\eta}_0} % ) \\ &\hspace{72pt} + \sum_{PQ}^M % \qty[ % \vb{T}_n(\boldsymbol{\eta}_0) % ]_{PQ} % \hat{\mathscr{F}}^\dagger_{QP} \qty( % \eval{ % \qty{ % \widehat{ % \hat{\mathcal{H}}_{\text{elec}} ( % \boldsymbol{\eta} % ) % , % \vb{T}(\boldsymbol{\eta}) % } % }_m % }_{\boldsymbol{\eta}_0} % ) % \Bigg] \thinspace , \end{split} \label{eq:second_order_perturbation_derivative_Hamiltonian} \end{equation}\] in which we have introduced the notation \[\begin{equation*} \qty{ % \widehat{ % \hat{O} % , \vb{T} % } % } % \end{equation*}\] to denote the resulting operator after \(\hat{O}\)’s one- and two-electron integrals have been one-index-transformed using the transformation matrix \(\vb{T}\).

Using the symmetric or Löwdin orthonormalization matrix as orbital connection, i.e. \[\begin{equation} \vb{T}(\boldsymbol{\eta}) % = \vb{S}^{-1/2}(\boldsymbol{\eta}) % \thinspace , \end{equation}\] we can use the formulas listed in this section, by replacing the first-order derivatives of the connection \(\vb{T}(\boldsymbol{\eta})\): \[\begin{equation} \eval{ % \pdv{ % \vb{T}(\boldsymbol{\eta}) % }{\eta_m} % }_{\boldsymbol{\eta}_0} % = \eval{ % \pdv{ % \vb{S}^{-1/2}(\boldsymbol{\eta}) % }{\eta_m} % }_{\boldsymbol{\eta}_0} % = - \frac{1}{2} % \eval{ % \pdv{ % \vb{S}(\boldsymbol{\eta}) % }{\eta_m} % }_{\boldsymbol{\eta}_0} \end{equation}\] and the second-order derivatives of the connection \(\vb{T}(\boldsymbol{\eta})\): \[\begin{align} \eval{ % \pdv{ % \vb{T}(\boldsymbol{\eta}) % }{\eta_m}{\eta_n} % }_{\boldsymbol{\eta}_0} % &= \eval{ % \pdv{ % \vb{S}^{-1/2}(\boldsymbol{\eta}) % }{\eta_m}{\eta_n} % }_{\boldsymbol{\eta}_0} \\ & = % - \frac{1}{2} % \eval{ % \pdv{ % \vb{S}(\boldsymbol{\eta}) % }{\eta_m}{\eta_n} % }_{\boldsymbol{\eta}_0} % + \frac{3}{8} (1 + P_{mn}) % \qty( % \eval{ % \pdv{ % \vb{S}(\boldsymbol{\eta}) % }{\eta_m} % }_{\boldsymbol{\eta}_0} % ) % \qty( % \eval{ % \pdv{ % \vb{S}(\boldsymbol{\eta}) % }{\eta_n} % }_{\boldsymbol{\eta}_0} % ) % \thinspace . \end{align}\] We have obtained these results by expanding (T. U. Helgaker and Almlöf 1984) the overlap matrix of the spinor basis \(\vb{S}(\boldsymbol{\eta})\) as: \[\begin{equation} \vb{S}(\boldsymbol{\eta}) % = \vb{I} % + \boldsymbol{\varepsilon}(\boldsymbol{\eta}) % \thinspace , \end{equation}\] where \[\begin{equation} \boldsymbol{\varepsilon}(\boldsymbol{\eta}_0) = \vb{0} \thinspace , \end{equation}\] and by using the Taylor expansion \[\begin{equation} \vb{S}^{-1/2}(\boldsymbol{\eta}) = % \vb{I} % - \frac{1}{2} \boldsymbol{\varepsilon}(\boldsymbol{\eta}) % + \frac{3}{8} \boldsymbol{\varepsilon}^2(\boldsymbol{\eta}) % + \order{\boldsymbol{\varepsilon}^3(\boldsymbol{\eta})} % \thinspace . \end{equation}\]

Since the coefficient matrix \(\vb{C}(\boldsymbol{\eta}\) is only implicitly dependent on the external perturbation \(\boldsymbol{\eta}\), the partial derivatives of the overlap matrix in spinor basis that are required can then by found as: \[\begin{equation} \begin{split} \eval{ % \pdv{ % S_{PQ}(\boldsymbol{\eta}) % }{\eta_m} % }_{\boldsymbol{\eta}_0} % = % & \sum_{\mu \nu}^{K_\alpha} % C^{\alpha *}_{\mu P}(\boldsymbol{\eta}_0) % C^{\alpha}_{\nu Q}(\boldsymbol{\eta}_0) % \eval{ % \pdv{ % S^{\alpha \alpha}_{\mu \nu}(\boldsymbol{\eta}) % }{\eta_m} % }_{\boldsymbol{\eta}_0} \\ & + \sum_{\mu \nu}^{K_\beta} % C^{\beta *}_{\mu P}(\boldsymbol{\eta}_0) % C^{\beta}_{\nu Q}(\boldsymbol{\eta}_0) % \eval{ % \pdv{ % S^{\beta \beta}_{\mu \nu}(\boldsymbol{\eta}) % }{\eta_m} % }_{\boldsymbol{\eta}_0} % \thinspace , \end{split} \end{equation}\] or, equivalently: \[\begin{equation} \eval{ % \pdv{ % \vb{S}(\boldsymbol{\eta}) % }{\eta_m} % }_{\boldsymbol{\eta}_0} % = % \vb{C}^\dagger(\boldsymbol{\eta}_0) % \begin{pmatrix} \eval{ % \dfrac{ % \partial \vb{S}^{\alpha \alpha}(\boldsymbol{\eta}) % }{ % \partial \eta_m % } % }_{\boldsymbol{\eta}_0} % & \vb{0} \\ \vb{0} % & \eval{ % \dfrac{ % \partial \vb{S}^{\beta \beta}(\boldsymbol{\eta}) % }{ % \partial \eta_m % } % }_{\boldsymbol{\eta}_0} % \end{pmatrix} \vb{C}(\boldsymbol{\eta}_0) \end{equation}\] for the first-order perturbational derivative. The second-order perturbational derivative can then be calculated using \[\begin{equation} \begin{split} \eval{ % \pdv{ % S_{PQ}(\boldsymbol{\eta}) % }{\eta_m}{\eta_n} % }_{\boldsymbol{\eta}_0} % = % & \sum_{\mu \nu}^{K_\alpha} % C^{\alpha *}_{\mu P}(\boldsymbol{\eta}_0) % C^{\alpha}_{\nu Q}(\boldsymbol{\eta}_0) % \eval{ % \pdv{ % S^{\alpha \alpha}_{\mu \nu}(\boldsymbol{\eta}) % }{\eta_m}{\eta_n} % }_{\boldsymbol{\eta}_0} \\ & + \sum_{\mu \nu}^{K_\beta} % C^{\beta *}_{\mu P}(\boldsymbol{\eta}_0) % C^{\beta}_{\nu Q}(\boldsymbol{\eta}_0) % \eval{ % \pdv{ % S^{\beta \beta}_{\mu \nu}(\boldsymbol{\eta}) % }{\eta_m}{\eta_n} % }_{\boldsymbol{\eta}_0} % \thinspace , \end{split} \end{equation}\] or equivalently: \[\begin{equation} \eval{ % \pdv{ % \vb{S}(\boldsymbol{\eta}) % }{\eta_m}{\eta_n} % }_{\boldsymbol{\eta}_0} % = % \vb{C}^\dagger(\boldsymbol{\eta}_0) % \begin{pmatrix} \eval{ % \dfrac{ % \partial^2 \vb{S}^{\alpha \alpha}(\boldsymbol{\eta}) % }{ % \partial \eta_m \partial \eta_n% } % }_{\boldsymbol{\eta}_0} % & \vb{0} \\ \vb{0} % & \eval{ % \dfrac{ % \partial^2 \vb{S}^{\beta \beta}(\boldsymbol{\eta}) % }{ % \partial \eta_m \partial \eta_n % } % }_{\boldsymbol{\eta}_0} % \end{pmatrix} \vb{C}(\boldsymbol{\eta}_0) % \thinspace . \end{equation}\]

The perturbation does not affect the scalar metric}

If, however, the metric of the underlying scalar bases does not change by turning on the perturbation, the overlap matrix remains the identity for every value of the perturbation: \[\begin{align} S_{PQ}(\boldsymbol{\eta}) % &= \int \dd{\vb{r}} \phi^*_P(\vb{r}; \boldsymbol{\eta}) % \phi_Q(\vb{r}; \boldsymbol{\eta}) \\ &= \vb{C}^\dagger(\boldsymbol{\eta}) % \begin{pmatrix} \vb{S}^{\alpha \alpha} % & \vb{0} \\ \vb{0} % & \vb{S}^{\beta \beta} \\ \end{pmatrix} \vb{C}(\boldsymbol{\eta}) \\ &= \delta_{PQ} % \thinspace . \end{align}\] which means that the initial set of spinors remains orthonormal throughout the perturbation. In this case, the Hamiltonian is written as \[\begin{equation} \hat{\mathcal{H}}(\boldsymbol{\eta}) % = % \sum_{PQ}^M % h_{PQ}(\boldsymbol{\eta}) % \hat{E}_{pq} % + \frac{1}{2} \sum_{PQRS}^K % g_{PQRS}(\boldsymbol{\eta}) % \thinspace \hat{e}_{PQRS} % + h_{\text{nuc}}(\boldsymbol{\eta}) % \thinspace , \end{equation}\] with the regular one- and two-electron integrals \(h_{PQ}\) and \(g_{PQRS}\), i.e. without the tildes that refer to transformations with an orbital connection. Therefore, the first-order perturbation partial derivative of the Hamiltonian become particularly simple: \[\begin{equation} \eval{ % \pdv{ % \hat{\mathcal{H}}(\boldsymbol{\eta}) % }{\eta_m} % }_{\boldsymbol{\eta}_0} % = \sum_{PQ}^M % \eval{ % \pdv{% h_{PQ}(\boldsymbol{\eta}) % }{\eta_m} % }_{\boldsymbol{\eta}_0} % \hat{E}_{PQ} % + \frac{1}{2} \sum_{PQRS}^M % \eval{ % \pdv{ % g_{PQRS}(\boldsymbol{\eta}) % }{\eta_m} % }_{\boldsymbol{\eta}_0} % \hat{e}_{PQRS} % + % \eval{ % \pdv{ % h_{\text{nuc}}(\boldsymbol{\eta}) % }{\eta_m} % }_{\boldsymbol{\eta}_0} % \end{equation}\] and so does the second-order partial derivative: \[\begin{equation} \eval{ % \pdv{ % \hat{\mathcal{H}}(\boldsymbol{\eta}) % }{\eta_m}{\eta_n} % }_{\boldsymbol{\eta}_0} % = \sum_{PQ}^M % \eval{ % \pdv{% h_{PQ}(\boldsymbol{\eta}) % }{\eta_m}{\eta_n} % }_{\boldsymbol{\eta}_0} % \hat{E}_{PQ} % + \frac{1}{2} \sum_{PQRS}^M % \eval{ % \pdv{ % g_{PQRS}(\boldsymbol{\eta}) % }{\eta_m}{\eta_n} % }_{\boldsymbol{\eta}_0} % \hat{e}_{PQRS} % + % \eval{ % \pdv{ % h_{\text{nuc}}(\boldsymbol{\eta}) % }{\eta_m}{\eta_n} % }_{\boldsymbol{\eta}_0} % \thinspace . \end{equation}\]

References

Helgaker, Trygve U., and Jan Almlöf. 1984. A second‐quantization approach to the analytical evaluation of response properties for perturbation‐dependent basis sets.” International Journal of Quantum Chemistry. https://doi.org/10.1002/qua.560260211.
Helgaker, Trygve, Sonia Coriani, Poul Jørgensen, Kasper Kristensen, Jeppe Olsen, and Kenneth Ruud. 2012. Recent Advances in Wave Function-Based Methods of Molecular-Property Calculations.” Chemical Reviews 112 (1): 543–631. https://doi.org/10.1021/cr2002239.