The GHF magnetic inducibility
At the Hartree-Fock level, i.e. for the (G)HF wave function model \[\begin{equation} \require{physics} \ket{ \text{HF}(\boldsymbol{\kappa}) } = \exp(-\hat{\kappa}) \ket{\text{HF}} \thinspace , \end{equation}\] we can calculate the total current density as (at the converged, and current orbitals denoted by \(\boldsymbol{\kappa}_0\)): \[\begin{multline} \vb{j} ( \vb{r}; \boldsymbol{\kappa}_0, \vb{B}_{\text{ext}}, \vb{G} ) \\ = - \frac{i}{2} \sum_I^N \Bigg( \phi_I^\dagger(\vb{r}) \thinspace [ \grad{ \phi_I(\vb{r}) } ] - [ \grad{ \phi_I^\dagger(\vb{r}) } ] \thinspace \phi_I(\vb{r}) + i \thinspace \phi_I^\dagger(\vb{r}) \phi_I(\vb{r}) \thinspace \vb{B}_{\text{ext}} \cross (\vb{r} - \vb{G}) \Bigg) \end{multline}\] if we do not employ any field-dependent spinors. In the RHF case, this formula leads to the result we find in (Stevens, Pitzee, and Lipscomb 1963).
The magnetic inducibility has implicit and explicit contributions. Let us first tackle the calculation of the explicit contribution. The explicit contribution to the (first-order) magnetic inducibility can be calculated as follows: \[\begin{equation} \eval{ \pdv{ j_{i, \thinspace PQ} ( \vb{r}; \vb{B}_{\text{ext}}, \vb{G} ) }{ B_{\text{ext}, m} } }_{\vb{B}_{\text{ext}, \thinspace 0}} = \frac{1}{2} \phi^\dagger_P(\vb{r}) \phi_Q(\vb{r}) \thinspace \epsilon_{imk} [\vb{r} - \vb{G}] \thinspace , \end{equation}\] where \(\epsilon_{imk}\) is the Levi-Civita symbol and the Einstein summation convention is implied. From this, the explicit contribution to the (first-order) magnetic inducibility can be calculated as: \[\begin{align} &\mathcal{J}_{im, \thinspace \text{explicit}} ( \vb{r}; \vb{p}^\star(\vb{B}_{\text{ext}, \thinspace 0}), \vb{B}_{\text{ext}, \thinspace 0}, \vb{G} ) \\ &\hspace{24pt} = \frac{1}{2} \epsilon_{imk} [\vb{r} - \vb{G}]_k \sum_{PQ}^M \phi^\dagger_P(\vb{r}) \phi_Q(\vb{r}) D_{PQ}(\vb{p}^\star(\vb{B}_{\text{ext}, \thinspace 0})) \\ &\hspace{24pt} = \frac{1}{2} \rho(\vb{r}; \vb{p}^\star(\vb{B}_{\text{ext}, \thinspace 0})) \thinspace \epsilon_{imk} [\vb{r} - \vb{G}]_k \thinspace , \end{align}\] from which follows the explicit contribution to the first-order induced response current density, \(\vb{j}_{\text{ind}, \thinspace \text{explicit}}\): \[\begin{align} & j^{(1)}_{ \text{ind}, \thinspace i, \thinspace \text{explicit} } ( \vb{r}; \vb{p}^\star(\vb{B}_{\text{ext}, \thinspace 0}), \vb{B}_{\text{ext}, \thinspace 0}, \vb{B}_{\text{ext}}, \vb{G} ) \\ & \hspace{48pt} = \sum_m \mathcal{J}_{im, \thinspace \text{explicit}} ( \vb{r}; \vb{p}^\star(\vb{B}_{\text{ext}, \thinspace 0}), \vb{B}_{\text{ext}, \thinspace 0}, \vb{G} ) \thinspace B_{\text{ext}, \thinspace m} \\ & \hspace{48pt} = \frac{1}{2} \thinspace \rho(\vb{r}; \vb{p}^\star(\vb{B}_{\text{ext}, \thinspace 0})) \thinspace [ \vb{B}_{\text{ext}} \cross (\vb{r} - \vb{G}) ]_i \thinspace , \end{align}\] in which, at the Hartree-Fock level of theory, the density can be written as: \[\begin{equation} \rho(\vb{r}; \vb{p}^\star(\vb{B}_{\text{ext}, \thinspace 0})) = \sum_{I}^N \phi^\dagger_I(\vb{r}) \phi_I(\vb{r}) \thinspace . \end{equation}\] In the RHF case, we find a similar expression in (Stevens, Pitzee, and Lipscomb 1963) and (Keith and Bader 1993). Keith and Bader (Keith and Bader 1993) have called this the diamagnetic contribution but, as far as I am aware, this has no connection with the diamagnetic part of the scalar kinetic operator. In fact, following the derivation of the Pauli current density, it is actually the paramagnetic contribution to the scalar kinetic operator which gives rise to this term.
The implicit part of the first-order induced response current density, is given through a contraction of the first-order (response) derivative of the 1-DM and the field-free current density matrix elements. \[\begin{equation} \mathcal{J}_{im, \thinspace \text{implicit}} ( \vb{r}; \vb{p}^\star(\vb{B}_{\text{ext}, \thinspace 0}), \vb{B}_{\text{ext}, \thinspace 0}, \vb{G} ) = \sum_{PQ}^M j_{i, \thinspace PQ} ( \vb{r}; \vb{B}_{\text{ext}, \thinspace 0} ) \thinspace \qty( \eval{ \dv{ D_{PQ}(\vb{p}^\star(\vb{B}_{\text{ext}}, \vb{G})) }{ B_{\text{ext}, \thinspace m} } }_{\vb{B}_{\text{ext}, \thinspace 0}} ) \thinspace . \end{equation}\] Specifically, for (G)HF theory, we find for the first-order (response) derivative of the 1-DM: \[\begin{equation} \eval{ \dv{ D_{PQ}(\vb{p}^\star(\vb{B}_{\text{ext}}, \vb{G})) }{ B_{\text{ext}, \thinspace m} } }_{\vb{B}_{\text{ext}, \thinspace 0}} = \sum_{AI} \qty( \eval{ \pdv{ D_{PQ}(\boldsymbol{\kappa}) }{ \kappa_{AI} } }_{\boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0})} ) \vb{x}_{AI, m} \thinspace , \end{equation}\] where we have used a short-hand symbol for the first-order parameter response: \[\begin{equation} \qty[ \vb{x}(\vb{B}_{\text{ext}, \thinspace 0}, \vb{G}) ]_{AI,m} = \eval{ \pdv{ \kappa_{AI}^\star( \vb{B}_{\text{ext}}, \vb{G} ) }{ B_{\text{ext}, \thinspace m} } }_{ \vb{B}_{\text{ext}, \thinspace 0} } \thinspace . \end{equation}\] Since the (G)HF 1-DM \(D_{PQ}(\boldsymbol{\kappa})\) is given by \[\begin{equation} D_{PQ}(\boldsymbol{\kappa}) = \ev{ \exp(\hat{\kappa}) \hat{E}_{PQ} \exp(-\hat{\kappa}) }{\text{core}} \thinspace , \end{equation}\] its parameter (orbital) derivative can be calculated using the BCH expansion and is subsequently given by: \[\begin{equation} \eval{ \pdv{ D_{PQ}(\boldsymbol{\kappa}) }{ \kappa_{AI} } }_{\boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0})} = - \delta_{P I} \delta_{QA} \thinspace . \end{equation}\] Continuing our calculation of the implicit first-order magnetic inducibility by plugging in this intermediary result, we find: \[\begin{equation} \mathcal{J}_{im, \thinspace \text{implicit}} ( \vb{r}; \boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0}), \vb{B}_{\text{ext}, \thinspace 0}, \vb{G} ) = - \frac{i}{2} \sum_{AI} \Bigg( \phi_I^\dagger(\vb{r}) \thinspace [ \grad_i{ \phi_A(\vb{r}) } ] - [ \grad_i{ \phi_I^\dagger(\vb{r}) } ] \thinspace \phi_A(\vb{r}) \Bigg) (- \vb{x}_{AI,m}) \thinspace . \end{equation}\]
If we now contract this implicit part of the (first-order) magnetic inducibility with the external field, we obtain the implicit part of the first-order response current density: \[\begin{equation} j^{(1)}_{i, \thinspace \text{implicit}} ( \vb{r}; \boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0}), \vb{B}_{\text{ext}, \thinspace 0}, \vb{G} ) = - \frac{i}{2} \sum_I \Bigg( \phi_I^\dagger(\vb{r}) \thinspace [ \grad_i{ \phi^{(1)}_I(\vb{r}) } ] - [ \grad_i{ \phi_I^\dagger(\vb{r}) } ] \thinspace \phi^{(1)}_I(\vb{r}) \Bigg) \thinspace , \end{equation}\] if we define some sort of first-order perturbed spinors: Then, by making the connection \[\begin{equation} \phi_{I}^{(1)}(\vb{r}) = - \sum_A \qty( \sum_m \vb{x}_{AI,m} B_{\text{ext}, m} ) \phi_A(\vb{r}) \thinspace . \end{equation}\] The listed expression for the (implicit part of the) first-order response current density is analoguous to what is listed in (Stevens, Pitzee, and Lipscomb 1963), but (Keith and Bader 1993) lists a plus sign in the middle of the parentheses, where we have a minus sign. Furthermore, we can make the connection with this derivation and their papers, by observing that we can retrieve their expansion coefficients as: \[\begin{equation} C^{(1)}_{AI} = - \sum_m \vb{x}_{AI,m} B_{\text{ext}, m} \thinspace . \end{equation}\]
Finally, the (G)HF (first-order) magnetic inducibility can be written as: \[\begin{multline} \mathcal{J}_{im} ( \vb{r}; \vb{p}^\star(\vb{B}_{\text{ext}, \thinspace 0}), \vb{B}_{\text{ext}, \thinspace 0}, \vb{G} ) \\ = - \frac{i}{2} \sum_{AI} \Bigg( \phi_I^\dagger(\vb{r}) \thinspace [ \grad_i{ \phi_A(\vb{r}) } ] - [ \grad_i{ \phi_I^\dagger(\vb{r}) } ] \thinspace \phi_A(\vb{r}) \Bigg) (- \vb{x}_{AI,m}) + \frac{1}{2} \rho(\vb{r}; \vb{p}^\star(\vb{B}_{\text{ext}, \thinspace 0})) \thinspace \epsilon_{imk} [\vb{r} - \vb{G}]_k \end{multline}\] and the associated (first-order) induced response current density can be written as: \[\begin{multline} \vb{j}^{(1)} ( \vb{r}; \boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0}), \vb{B}_{\text{ext}, \thinspace 0}, \vb{G} ) \\ = - \frac{i}{2} \sum_I \Bigg( \phi_I^\dagger(\vb{r}) \thinspace [ \grad{ \phi^{(1)}_I(\vb{r}) } ] - [ \grad{ \phi_I^\dagger(\vb{r}) } ] \thinspace \phi^{(1)}_I(\vb{r}) \Bigg) + \frac{1}{2} \thinspace \rho(\vb{r}; \vb{p}^\star(\vb{B}_{\text{ext}, \thinspace 0})) \thinspace [ \vb{B}_{\text{ext}} \cross (\vb{r} - \vb{G}) ] \thinspace . \end{multline}\]
In order to find the first-order linear wave function response, we’ll have to solve the associated linear response equations. The response force constant for (G)HF is just the field-free (G)HF orbital Hessian: \[\begin{equation} \qty[ \vb{k}_{\boldsymbol{\kappa}} ]_{AI, BJ} = \eval{ \pdv{ E(\vb{B}_{\text{ext}, \thinspace 0}, \boldsymbol{\kappa}) }{\kappa_{AI}}{\kappa_{BJ}} }_{ \boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0}) } \end{equation}\] and the response force can be calculated as: \[\begin{equation} \qty[ \vb{F}_{\boldsymbol{\kappa}} ]_{AI, m} = \eval{ \pdv{ E(\vb{B}_{\text{ext}}, \vb{G}, \boldsymbol{\kappa}) }{\kappa_{AI}}{ B_{\text{ext}, \thinspace m} } }_{ \vb{B}_{\text{ext}, \thinspace 0}, \thinspace \boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0}) } \thinspace . \end{equation}\]