Linear response derivatives for RHF
In order to find second-order response properties using the RHF wave function model, we have to solve the linear response equations. In RHF theory, the optimizable parameters are the non-redundant orbital rotation parameters \(\{ \kappa_{ai} \}\), which can be parametrized by the two sets of real parameters \(\{ {}^{\text{R}} \kappa_{ai} \}\) and \(\{ {}^{\text{I}} \kappa_{ai} \}\): \[\begin{equation} \require{physics} \kappa_{ai} = {}^{\text{R}} \kappa_{ai} + i {}^{\text{I}} \kappa_{ai} \thinspace . \end{equation}\]
The linear response equations (in this case also called the CPHF, coupled perturbed Hartree-Fock, equations) then admit a blocked-out form: \[\begin{equation} \begin{pmatrix} {}^{\text{R}} \vb{F}_{\boldsymbol{\kappa}} \\ {}^{\text{I}} \vb{F}_{\boldsymbol{\kappa}} \end{pmatrix} + \begin{pmatrix} {}^{\text{RR}} \vb{k}_{\boldsymbol{\kappa}} & {}^{\text{RI}} \vb{k}_{\boldsymbol{\kappa}} \\ {}^{\text{IR}} \vb{k}_{\boldsymbol{\kappa}} & {}^{\text{II}} \vb{k}_{\boldsymbol{\kappa}} \\ \end{pmatrix} \begin{pmatrix} {}^{\text{R}} \vb{x} \\ {}^{\text{I}} \vb{x} \end{pmatrix} = \begin{pmatrix} \vb{0} \\ \vb{0} \end{pmatrix} \thinspace . \end{equation}\] In full, they can be written as: \[\begin{equation} \eval{ \pdv{ E(\boldsymbol{\eta}, \boldsymbol{\kappa}) }{ {}^{\text{R}} \kappa_{ai} }{\eta_m} }_{ \boldsymbol{\eta}_0, \boldsymbol{\kappa}^\star(\boldsymbol{\eta}_0) } + \sum_{bj} \qty( \eval{ \pdv{ E(\boldsymbol{\eta}_0, \boldsymbol{\kappa}) }{ {}^{\text{R}} \kappa_{ai} }{ {}^{\text{R}} \kappa_{bj} } }_{\boldsymbol{\kappa}^\star(\boldsymbol{\eta}_0)} ) \thinspace {}^{\text{R}} x_{bj} + \sum_{bj} \qty( \eval{ \pdv{ E(\boldsymbol{\eta}_0, \boldsymbol{\kappa}) }{ {}^{\text{R}} \kappa_{ai} }{ {}^{\text{I}} \kappa_{bj} } }_{\boldsymbol{\kappa}^\star(\boldsymbol{\eta}_0)} ) \thinspace {}^{\text{I}} x_{bj} = 0 \thinspace , \end{equation}\] and \[\begin{equation} \eval{ \pdv{ E(\boldsymbol{\eta}, \boldsymbol{\kappa}) }{ {}^{\text{I}} \kappa_{ai} }{\eta_m} }_{ \boldsymbol{\eta}_0, \boldsymbol{\kappa}^\star(\boldsymbol{\eta}_0) } + \sum_{bj} \qty( \eval{ \pdv{ E(\boldsymbol{\eta}_0, \boldsymbol{\kappa}) }{ {}^{\text{I}} \kappa_{ai} }{ {}^{\text{R}} \kappa_{bj} } }_{\boldsymbol{\kappa}^\star(\boldsymbol{\eta}_0)} ) \thinspace {}^{\text{R}} x_{bj} + \sum_{bj} \qty( \eval{ \pdv{ E(\boldsymbol{\eta}_0, \boldsymbol{\kappa}) }{ {}^{\text{I}} \kappa_{ai} }{ {}^{\text{I}} \kappa_{bj} } }_{\boldsymbol{\kappa}^\star(\boldsymbol{\eta}_0)} ) \thinspace {}^{\text{I}} x_{bj} = 0 \thinspace . \end{equation}\]
The elements of the response force are given by: \[\begin{align} {}^{\text{R}} F_{ai, m} (\boldsymbol{\kappa}^\star(\boldsymbol{\eta}_0)) = \eval{ \pdv{ E(\boldsymbol{\eta}, \boldsymbol{\kappa}) }{ {}^{\text{R}} \kappa_{ai} }{\eta_m} }_{ \boldsymbol{\eta}_0, \boldsymbol{\kappa}^\star(\boldsymbol{\eta}_0) } &= \ev{ \comm{ \hat{E}^-_{ai} }{\qty( \eval{ \pdv{ \hat{\mathcal{H}}(\boldsymbol{\eta}) }{\eta_m} }_{\boldsymbol{\eta}_0} )} }{\text{core}} \\ &= - 2 \qty[ {}^{i} F_{ia} \qty( \eval{ \pdv{ \hat{\mathcal{H}}(\boldsymbol{\eta}) }{\eta_m} }_{\boldsymbol{\eta}_0} ) + {}^{i} F_{ai} \qty( \eval{ \pdv{ \hat{\mathcal{H}}(\boldsymbol{\eta}) }{\eta_m} }_{\boldsymbol{\eta}_0} ) ] \thinspace , \end{align}\] and \[\begin{align} {}^{\text{I}} F_{ai, m} (\boldsymbol{\kappa}^\star(\boldsymbol{\eta}_0)) = \eval{ \pdv{ E(\boldsymbol{\eta}, \boldsymbol{\kappa}) }{ {}^{\text{I}} \kappa_{ai} }{\eta_m} }_{ \boldsymbol{\eta}_0, \boldsymbol{\kappa}^\star(\boldsymbol{\eta}_0) } &= i \ev{ \comm{ \hat{E}^+_{ai} }{\qty( \eval{ \pdv{ \hat{\mathcal{H}}(\boldsymbol{\eta}) }{\eta_m} }_{\boldsymbol{\eta}_0} )} }{\text{core}} \\ &= 2i \qty[ {}^{i}F_{ai} \qty( \eval{ \pdv{ \hat{\mathcal{H}}(\boldsymbol{\eta}) }{\eta_m} }_{\boldsymbol{\eta}_0} ) - {}^{i} F_{ia} \qty( \eval{ \pdv{ \hat{\mathcal{H}}(\boldsymbol{\eta}) }{\eta_m} }_{\boldsymbol{\eta}_0} ) ] \thinspace , \end{align}\] where we have done some simplifications due to the RHF Fockian and have assumed that the first-order perturbational partial derivative of the Hamiltonian is again Hermitian. The elements of the response force constant \(\vb{k}_{\boldsymbol{\kappa}}\) are exactly those of the RHF orbital Hessian.
Following (Ruud et al. 1993), if we may use real spin-orbitals for the perturbation-free optimization, and use an imaginary wave function response for the first-order perturbation, we can apply some simplifications. In the previous equations, this means that we set the real first-order response to zero: \[\begin{equation} {}^{\text{R}} x_{ai, m} = \qty( \eval{ \pdv{ {}^{\text{R}} \kappa_{ai}^\star(\boldsymbol{\eta}) }{\eta_m} }_{\boldsymbol{\eta}_0} ) = 0 \thinspace , \end{equation}\] which leads to \[\begin{equation} \begin{pmatrix} {}^{\text{R}} \vb{F}_{\boldsymbol{\kappa}} \\ {}^{\text{I}} \vb{F}_{\boldsymbol{\kappa}} \end{pmatrix} + \begin{pmatrix} {}^{\text{RI}} \vb{k}_{\boldsymbol{\kappa}} {}^{\text{I}} \vb{x} \\ {}^{\text{II}} \vb{k}_{\boldsymbol{\kappa}} {}^{\text{I}} \vb{x} \\ \end{pmatrix} = \begin{pmatrix} \vb{0} \\ \vb{0} \end{pmatrix} \thinspace , \end{equation}\] which can in turn be rewritten as: \[\begin{equation} \qty( {}^{\text{R}} \vb{F}_{\boldsymbol{\kappa}} -i {}^{\text{I}} \vb{F}_{\boldsymbol{\kappa}} ) + \qty( {}^{\text{RI}} \vb{k}_{\boldsymbol{\kappa}} -i {}^{\text{II}} \vb{k}_{\boldsymbol{\kappa}} ) {}^{\text{I}} \vb{x} = \vb{0} \thinspace . \end{equation}\] Upon closer inspection, we indeed end up with the equations listed in (Ruud et al. 1993), with \[\begin{align} \qty( {}^{\text{R}} \vb{F}_{\boldsymbol{\kappa}} -i {}^{\text{I}} \vb{F}_{\boldsymbol{\kappa}} )_{ai, m} &= 2 \ev{ \comm{ \hat{E}_{ai} }{\qty( \eval{ \pdv{ \hat{\mathcal{H}}(\boldsymbol{\eta}) }{\eta_m} }_{\boldsymbol{\eta}_0} )} }{\text{core}} \\ &= -4 {}^{i} F_{ia} \qty( \eval{ \pdv{ \hat{\mathcal{H}}(\boldsymbol{\eta}) }{\eta_m} }_{\boldsymbol{\eta}_0} ) \end{align}\] and \[\begin{align} \qty( {}^{\text{RI}} \vb{k}_{\boldsymbol{\kappa}} -i {}^{\text{II}} \vb{k}_{\boldsymbol{\kappa}} )_{aibj} &= 2i \ev*{ \comm*{ \hat{E}_{ai} }{ \comm*{ \hat{E}^+_{bj} }{ \hat{\mathcal{H}}(\boldsymbol{\eta}_0) } } }{ \text{core} } \\ &= -4i \qty[ \delta_{ij} \delta_{ab} (\epsilon_a - \epsilon_i) + g_{ibja} - g_{ijba} ] \thinspace . \end{align}\]
Working out these equations a little more, we find6: \[\begin{equation} {}^{i} F_{ai} \qty( \eval{ \pdv{ \hat{\mathcal{H}}(\boldsymbol{\eta}) }{\eta_m} }_{\boldsymbol{\eta}_0} ) - i \qty[ (\epsilon_a - \epsilon_i) {}^{\text{I}} x_{ai,m} + \sum_{bj} (g_{ibja} - g_{ijba}) {}^{\text{I}} x_{bj,m} ] = 0 \thinspace , \end{equation}\] where we have used the fact that in the basis of the canonical HF orbitals, the inactive Fock matrix is diagonal and we have used real permutational symmetries since we have emloyed real spin-orbitals in the perturbation-free optimization. If we now define the purely imaginary coefficients \[\begin{equation} C^{(1)}_{ai} (\boldsymbol{\eta}) = \sum_m \qty( -i {}^{\text{I}} x_{ai, m} ) \eta_m \thinspace , \end{equation}\] we may rewrite the CPHF equations in the following form: \[\begin{equation} \sum_m {}^{i} F_{ai} \qty( \eval{ \pdv{ \hat{\mathcal{H}}(\boldsymbol{\eta}) }{\eta_m} }_{\boldsymbol{\eta}_0} ) \eta_m + (\epsilon_a - \epsilon_i) C^{(1)}_{ai} (\boldsymbol{\eta}) + \sum_{bj} (g_{ibja} - g_{ijba}) C^{(1)}_{bj} (\boldsymbol{\eta}) = 0 \thinspace . \end{equation}\] In order to recognize the complete, exact form listed by (Stevens, Pitzee, and Lipscomb 1963) and (Keith and Bader 1993), we finally have to calculate the first term in the previous equation. First of all, we recognize the first-order perturbed Hamiltonian: \[\begin{equation} \mathcal{H}^{c, (1)}(\boldsymbol{\eta}) = \sum_m \qty( \eval{ \pdv{ \hat{\mathcal{H}}(\boldsymbol{\eta}) }{\eta_m} }_{\boldsymbol{\eta}_0} ) \eta_m \thinspace , \end{equation}\] and we then assume that the first-order perturbed Hamiltonian is a one-electron operator (i.e. the perturbed physical interaction only takes place in the one-electron part of the Hamiltonian). The inactive Fockian element then becomes: \[\begin{equation} \sum_m {}^{i} F_{ai} \qty( \eval{ \pdv{ \hat{\mathcal{H}}(\boldsymbol{\eta}) }{\eta_m} }_{\boldsymbol{\eta}_0} ) \eta_m = \matrixel{\phi_a}{ \mathcal{H}^{c, (1)}(\boldsymbol{\eta}) }{\phi_i} \thinspace . \end{equation}\] Plugging this in, we thus finally obtain the CPHF equations (Stevens, Pitzee, and Lipscomb 1963) (Keith and Bader 1993): \[\begin{equation} \matrixel{\phi_a}{ \mathcal{H}^{c, (1)}(\boldsymbol{\eta}) }{\phi_i} + (\epsilon_a - \epsilon_i) C^{(1)}_{ai} (\boldsymbol{\eta}) + \sum_{bj} (g_{ibja} - g_{ijba}) C^{(1)}_{bj} (\boldsymbol{\eta}) = 0 \thinspace . \end{equation}\]
Electric response derivatives
In the case of real RHF theory, we may discard any contributions arising from \({}^{\text{I}} \vb{x}\) and we may ignore the bottom half of the CPHF equations. Specifically, the (real) response force for an electric field \(\vb{F}\) is: \[\begin{equation} \qty( \vb{F}_{\boldsymbol{\kappa}} )_{ai, m} = \eval{ \pdv{ E(\vb{F}, \boldsymbol{\kappa}) }{\kappa_{ai}}{F_m} }_{\vb{F}_0, \boldsymbol{\kappa}^\star} = 4 \mu_{m, ai} \thinspace , \end{equation}\] where we have used the theory on perturbation-dependent Fock spaces and the Schrödinger Hamiltonian for a uniform electric field in order to find an expression for the electronic Hamiltonian in a uniform external electric field. We should note that because these gradients and Hessians use expressions that involve commutators of the total Hamiltonian, the scalar term can be dropped and we can work with the electronic Hamiltonian instead. The response force constant that we can use in this case is the (real) RHF orbital Hessian: \[\begin{equation} \vb{k}_{\boldsymbol{\kappa}} = {}^{\text{R} \text{R}} \vb{E}^{(2)}_{aibj} = 4 \qty[ \delta_{ij} {}^{i} F_{ab} - \delta_{ab} {}^{i} F_{ij} + 4 g_{aibj} - g_{abij} - g_{ajbi} ] \thinspace . \end{equation}\]
External magnetic fields without London orbitals
Let us discuss the case where the external perturbation is a uniform external magnetic field \(\vb{B}_{\text{ext}}\). If we are not employing London orbitals, we may treat the spin-orbital basis metric as perturbation-independent and can therefore proceed with the simplest form of the perturbation-dependent Hamiltonian. The linear response equations become \[\begin{equation} {}^{i} F_{ai} \qty( \eval{ \pdv{ \hat{\mathcal{H}}(\vb{B}_{\text{ext}}, \vb{G}) }{B_{\text{ext}, \thinspace m}} }_{\vb{B}_{\text{ext}, \thinspace 0}} ) - i \qty[ (\epsilon_a - \epsilon_i) {}^{\text{I}} x_{ai,m} (\vb{G}) + \sum_{bj} (g_{ibja} - g_{ijba}) {}^{\text{I}} x_{bj,m} (\vb{G}) ] = 0 \thinspace , \end{equation}\] which should be fulfilled for any choice of gauge origin \(\vb{G}\). From the Schrödinger Hamiltonian, we can find the one-electron perturbational terms due the presence of an external magnetic field. The first-order derivative operator that should be quantized can be calculated as (Keith and Bader 1993): \[\begin{align} \eval{ \pdv{ h^c(\vb{B}_{\text{ext}}, \vb{G}) }{B_{\text{ext}, m}} }_{\vb{B}_{\text{ext}, 0}} &= -\frac{i}{2} \qty[ (\vb{r} - \vb{G}) \cross \grad ]_m \\ &= \frac{1}{2} L^{c}_m(\vb{G}) \end{align}\] which is calculated from the perturbational derivative of the paramagnetic term and where \(L^{c}_m(\vb{G})\) is the orbital angular momentum operator about the reference point \(\vb{G}\). Its contraction with the perturbation \(\vb{B}_{\text{ext}}\) then yields the first-order perturbed Hamiltonian \[\begin{equation} \mathcal{H}^{(1)}(\vb{B}_{\text{ext}}) = \frac{1}{2} \vb{B}_{\text{ext}} \vdot \vb{L}^{c}(\vb{G}) \thinspace , \end{equation}\] which is for example listed in (Keith and Bader 1993).
In the Pauli Hamiltonian there are only one-electron terms affected by the uniform external magnetic field and since it only appears inside a commutator, we may use the following form of the first-order partial perturbational derivative of the Hamiltonian: \[\begin{equation} \eval{ \pdv{ \hat{\mathcal{H}}(\vb{B}_{\text{ext}}, \vb{G}) }{B_{\text{ext}, m}} }_{\vb{B}_{\text{ext}, 0}} = \sum_{PQ}^M \eval{ \pdv{ h_{PQ}(\vb{B}_{\text{ext}}, \vb{G}) }{B_{\text{ext}, m}} }_{\vb{B}_{\text{ext}, 0}} \hat{E}_{PQ} \thinspace , \end{equation}\] where the two-electron terms are not affected by the perturbation and the nuclear term can be ignored due to the aforementioned commutator-related discussion. The physical interactions that are included, therefore enter the calculations through the (derivative of the) one-electron integrals only. The derivative operator that should be quantized can be calculated as: \[\begin{equation} \eval{ \pdv{ h^c(\vb{B}_{\text{ext}}, \vb{G}) }{B_{\text{ext}, m}} }_{\vb{B}_{\text{ext}, 0}} = -\frac{i}{2} \qty[ (\vb{r} - \vb{G}) \cross \grad ]_m \vb{I}_2 + \frac{1}{2} \sigma_m \thinspace , \end{equation}\] where only the first part is considered in (Keith and Bader 1993) and (Ruud et al. 1993).
Gauge transformations (without London orbitals)
We will now parametrize the gauge origin \(\vb{G}\) as \[\begin{equation} \vb{G}(\vb{d}) = \vb{G}_0 + \vb{d} \thinspace , \end{equation}\] with \(\vb{G}_0\) a fixed point in space, i.e. \[\begin{equation} \vb{G}_0 = \vb{G}(\vb{d}_0) \end{equation}\] and \(\vb{d}\) a perturbation on the initial gauge origin \(\vb{G}_0\). By taking the partial derivative \[\begin{equation} \eval{ \pdv{d_n} \qty( \vphantom{\frac{1}{2}} \cdot ) }_{\vb{d}_0} \end{equation}\] of the linear response equations, we find: \[\begin{equation} {}^{i} F_{ai} \qty( \eval{ \pdv{d_n} \qty( \eval{ \pdv{ \hat{\mathcal{H}}(\vb{B}_{\text{ext}}, \vb{d}) }{B_{\text{ext}, \thinspace m}} }_{\vb{B}_{\text{ext}, \thinspace 0}} ) }_{\vb{d}_0} ) - i \qty[ (\epsilon_a - \epsilon_i) {}^{\text{I}} \varepsilon_{ai,mn} + \sum_{bj} (g_{ibja} - g_{ijba}) {}^{\text{I}} \varepsilon_{bj,mn} ] = 0 \thinspace , \end{equation}\] where we have defined a mixed second-order response property: \[\begin{equation} {}^{\text{I}} \varepsilon_{ai,mn} = \eval{ \pdv{ {}^{\text{I}} x_{ai, m} (\vb{d}) }{d_n} }_{\vb{d}_0} \thinspace . \end{equation}\] Since the operator that should be quantized in the inactive Fockian: \[\begin{equation} \eval{ \pdv{d_n} \qty( \eval{ \pdv{ h^c(\vb{B}_{\text{ext}}, \vb{d}) }{B_{\text{ext}, m}} }_{\vb{B}_{\text{ext}, 0}} ) }_{\vb{d}_0} = \frac{i}{2} \epsilon_{mnf} \partial_f \thinspace , \end{equation}\] contains the Levi-Civita tensor \(\epsilon_{mnf}\), if \(m=n\) (i.e. considering the same component for the magnetic field and gauge origin translation), the response force is zero such that the corresponding diagonal mixed-order responses vanish: \[\begin{equation} {}^{\text{I}} \varepsilon_{ai,mm} = 0 \thinspace . \end{equation}\] The contraction of the mixed-order perturbational derivative with \(\vb{d}\) and \(\vb{B}\) yields the mixed second-order perturbed Hamiltonian \[\begin{equation} \mathcal{H}^{(1, 1)}(\vb{B}_{\text{ext}}, \vb{d}) = \frac{i}{2} \vb{B}_{\text{ext}} \vdot \vb{d} \cross \grad \thinspace , \end{equation}\] which is listed in (Keith and Bader 1993) equation (10). The connection with equation (12) in (Keith and Bader 1993) is found by the following realization of their changes in first-order coefficients: \[\begin{equation} \delta C^{(1)}_{ai} (\vb{B}_{\text{ext}}, \vb{d}) = \sum_{m \neq n} (-i {}^{\text{I}} \varepsilon_{ai, mn}) \thinspace B_m d_n \thinspace . \end{equation}\]
References
or perhaps its complex conjugate↩︎