\documentclass[twocolumn, letterpaper, superscriptaddress, prl, showpacs]{revtex4}


\usepackage{ifthen}
\usepackage{graphicx}
\usepackage{amsmath}


\def\A{{\bf A}}
\def\B{{\bf B}}
\def\Be{{\bf B}^{\text{ext}}}
\def\Bi{{\bf B}^{\text{ind}}}
\def\bi{B^{\text{ind}}}
\def\be{B^{\text{ext}}}
\def\D{{\bf D}}
\def\F{{\bf F}}
\def\G{{\bf G}}
\def\H{{\bf H}}
\def\J{{\bf J}}
\def\j{{\bf j}}
\def\k{{\bf k}}
\def\M{{\bf M}}
\def\m{{\bf m}}
\def\n{{\bf n}}
\def\p{{\bf p}}
\def\r{{\bf r}}
\def\R{{\bf R}}
\def\bra#1{\langle#1\vert}
\def\ket#1{\vert#1\rangle}
\def\me#1#2#3{\langle#1\vert#2\vert#3\rangle}
\def\emac{\cal E}
\def\Emac{\boldsymbol{\cal E}}






\begin{document}




\title{A converse approach to the calculation of NMR shielding tensors}
\author{T. Thonhauser}
\affiliation{Department of Physics and Astronomy, Rutgers University,
Piscataway, New Jersey 08854, USA.}
\affiliation{Department of Materials Science and Engineering,
MIT, Cambridge, Massachusetts 02139, USA.}
\author{Arash A. Mostofi}
\altaffiliation{Current address: Department of Materials, Imperial College
London, SW7 2AZ, United Kingdom.}
\affiliation{Department of Materials Science and Engineering,
MIT, Cambridge, Massachusetts 02139, USA.}
\author{Nicola Marzari}
\affiliation{Department of Materials Science and Engineering,
MIT, Cambridge, Massachusetts 02139, USA.}
\author{R. Resta}
\affiliation{Dipartimento di Fisica Teorica Universit\`a di Trieste and
INFM-DEMOCRITOS, Trieste, Italy.}
\author{David Vanderbilt}
\affiliation{Department of Physics and Astronomy, Rutgers University,
Piscataway, New Jersey 08854, USA.}
\date{\today}

\begin{abstract}  
We introduce an alternative approach to the first-principles calculation of NMR
shielding tensors. These are obtained from the derivative of the orbital 
magnetization with respect to the application of a microscopic, localized
magnetic dipole. The approach is simple, general, and can be applied to either
isolated or periodic systems. Calculated results show very good agreement with 
established methods and experimental results.
\end{abstract}

\pacs{
71.15.-m, % Methods of electronic structure calculations
71.15.Mb, % Density functional theory, LDA, gradient and other corrections
75.20.-g, % Diamagnetism, paramagnetism, and superparamagnetism
76.60.Cq  % Chemical and Knight shifts
}


\maketitle


Nuclear magnetic resonance (NMR) measures the transition frequencies for the
reorientation of nuclear magnetic moments in an applied magnetic field. Since
the local magnetic field differs from the external one as a result of electronic
screening, NMR spectroscopy~\cite{NMR_encyclopedia} has been recognized since
1938~\cite{Rabi} to be a powerful experimental probe of local chemical
environments, including structural and functional information on molecules,
liquids, and increasingly, on solid-state systems.

First-principles calculations of NMR spectra were first developed in the quantum
chemistry community~\cite{Kutzelnigg_90} and applied to molecules and clusters,
but applications to extended crystalline systems were hindered by the difficulty
of including macroscopic magnetic fields, which require a non-periodic vector
potential and therefore destroy Bloch symmetry. In 1996, Mauri \emph{et al.}\
developed a linear-response approach for calculating NMR shieldings in periodic
crystals  based on the long-wavelength limit of a periodic modulation of the
applied magnetic field~\cite{Mauri_96}, while in 2001 Sebastiani and Parrinello
used the Wannier representation to derive an alternative linear-response
approach based on the application of an infinitesimal uniform magnetic
field~\cite{Sebastiani_01}.  More recently, attention has focused on the
development  of these approaches in the context of pseudopotential and
projector-augmented wave (PAW)
descriptions~\cite{Pickard_Mauri_01,Yates_07,Truflandier_07}, leading to a
growing use of these methods in combination with modern plane-wave
pseudopotential codes \cite{paratec, castep}. Despite these advances,  existing
methods for computing NMR shifts in crystalline systems remain complex, in that
they require a linear-response implementation with significant extra coding. In
principle, the long-wavelength approach~\cite{Mauri_96} could be implemented
using standard ground-state calculations, but only at the cost of introducing
prohibitively large supercells. Furthermore, constant magnetic fields can be
dealt with in supercell calculations---such as e.g.\ in
Ref.~\cite{Cai_07}---provided they are commensurate, with an integer number of
flux quanta per supercell. The drawback is that, even for large supercells, the
lower limit of commensurate fields is much larger than the field values actually
achievable in a laboratory.

In this Letter, we reformulate the problem of computing NMR shielding tensors in
such a way that the need for a linear-response framework is circumvented. 
For clarity, the previous formulations shall be referred to here as
\emph{direct} approaches, in that a magnetic field is applied and the local
field at the nucleus is computed. Our alternative, \emph{converse} approach
obtains the NMR shifts instead from the macroscopic magnetization induced by
magnetic point dipoles placed at the nuclear sites of interest. We demonstrate
the method by a first application to small molecules,  finding excellent
agreement with experiment and with calculations using the direct approach.  Our
new method is both simple and general, provides a straightforward alternative
avenue to the computation of NMR shifts in complex materials that requires very
little coding, and  can be extended to electronic-structure methods where a
linear-response formulation is impractical.

Let us start by considering a sample to which a constant external magnetic field
$\Be$ is applied. The field induces a current that, in turn, induces a magnetic
field $\Bi(\r)$ such that the total magnetic field is $\B(\r) = \Be  +
\Bi(\r)$.  In NMR experiments the applied fields are small compared to the
typical electronic scales;  the absolute chemical shielding tensor
$\tensor{\sigma}$ is then defined via the linear relationship
%
\begin{equation}\label{equ:sigma}
\Bi_s=-\tensor{\sigma}_s\cdot\;\Be\;,
\quad \sigma_{s,\alpha\beta} =
-\frac{\partial\bi_{s,\alpha}}{\partial \be_\beta}\;.
\end{equation}
%
The index $s$ indicates that the corresponding quantity is to be taken at
position $\r_s$, i.e., the site of nucleus $s$. NMR experiments usually report
the isotropic shielding $\sigma_s=\frac{1}{3}{\text{Tr}}[\,\tensor{\sigma}_s\,]$
via a chemical shift that is defined by convention as $-[\sigma_s -
\sigma_{\text{ref}}]$. Here $\sigma_{\text{ref}}$ is the  isotropic shielding of
a reference compound, e.g., tetramethylsilane for $^1$H and $^{13}$C NMR.

As mentioned above, direct approaches
\cite{Mauri_96,Sebastiani_01,Pickard_Mauri_01,Yates_07,Truflandier_07} calculate
the chemical shielding from the current response of the system to an external
magnetic field using perturbation theory.  We propose a fundamentally different
approach that circumvents the need to apply an external magnetic field, relying
instead on a thermodynamic relationship between mixed partial derivatives. Using
$B_{s,\alpha}=\be_\alpha + \bi_{s,\alpha}$, Eq.~(\ref{equ:sigma}) becomes
$\delta_{\alpha\beta} - \sigma_{s,\alpha\beta} = \partial
B_{s,\alpha}/\partial\be_\beta$. For the moment, we assume that $\Be$ can be
replaced by the total macroscopic $B$-field in the denominator of this equation,
thus neglecting the {\it macroscopic} induced field.  (This restriction,
appropriate for normal components in a slab geometry, will be relaxed shortly.)
The numerator may be written as $B_{s,\alpha} = -\partial E/\partial
m_{s,\alpha}$, where $E$ is interpreted either as the energy of a virtual
magnetic dipole $\m_s$ at one nuclear center $\r_s$ in the field $\B$ for a
finite system, or as the energy per cell of a periodic lattice of such dipoles;
we adopt the latter view.  Then, writing the macroscopic magnetization as
$M_\beta=-\Omega^{-1}\,\partial E/\partial B_\beta$ (where $\Omega$ is the cell
volume),
%
\begin{equation}\label{equ:converse}
\delta_{\alpha\beta} - \sigma_{s,\alpha\beta} =
   - \frac{\partial}{\partial B_\beta}
   \frac{\partial E}{\partial m_{s,\alpha}} =
   - \frac{\partial}{\partial m_{s,\alpha}}
   \frac{\partial E}{\partial B_\beta} =
   \Omega\frac{\partial M_\beta}{\partial m_{s,\alpha}}\;.
\end{equation}
%
Thus, $\tensor{\sigma}_s$ accounts for the shielding contribution to the
macroscopic magnetization induced by a magnetic point dipole $\m_s$ sitting at
nucleus $\r_s$ and all of its periodic replicas. In other words, instead of
applying a constant (or long-wavelength) field $\Be$ to an infinite periodic
system and calculating the induced field at all equivalent $s$ nuclei, we apply
an infinite array of magnetic dipoles to all equivalent sites $s$, and calculate
the change in magnetization. Since the perturbation is now periodic, it can
simply be computed using finite differences of ground-state calculations;
isolated system can be handled straightforwardly in the spirit of the supercell
approximation. This is our principal result, and it is viable for periodic
systems only because of the availability of the recently introduced modern
theory of magnetization \cite{Resta_05,Thonhauser_06,Ceresoli_06,Niu}. Note that
$\M=\m_s/\Omega + \M^{\text{ind}}$, where the first term is present merely
because we have included magnetic dipoles by hand. It follows that the shielding
is related to the true induced magnetization via $\sigma_{s,\alpha\beta}=
-\Omega\,\partial M^{\text{ind}}_\beta/\partial m_{s,\alpha}$.

It is useful to pause here and consider the analogy of the Born~\cite{Born} (or
dynamical) effective charge tensor $Z^*_{s,\alpha\beta}$, which may be regarded
as (i) the component of the force $\F_s$ in direction $\alpha$ on site $\r_s$ by
a unit macroscopic electric field $\Emac$ in direction $\beta$ (at zero nuclear
displacement), or alternatively as (ii) the $\beta$-component of the 
macroscopic electric polarization ${\bf P}$ linearly induced by a unit
displacement of nucleus $s$ and its periodic replicas in direction $\alpha$ in
vanishing macroscopic electric field.  Since the force on  nucleus $s$ is given
by $F_{s,\alpha}=-\partial  E/\partial r_{s,\alpha}$, (i) and (ii) are related
by
%
\begin{equation}\label{equ:analogy}
Z^*_{s,\alpha\beta}
= -\frac{\partial}{\partial\emac_\beta}
   \frac{\partial E}{\partial r_{s,\alpha}} 
= -\frac{\partial}{\partial r_{s,\alpha}}
   \frac{\partial E}{\partial\emac_\beta}
=  \Omega\frac{\partial P_\beta}{\partial r_{s,\alpha}}\;,
\end{equation}
%
in close analogy with Eq.~(\ref{equ:converse}). Note that, in order to comply
with the Born definition, one must choose the \emph{lattice-periodical} solution
of Poisson's equation, corresponding to vanishing macroscopic electric field.
Other choices are possible and lead to other kinds of effective charges, which
are all related to each other by means of the dielectric
constant~\cite{Ghosez98}. By comparing Eq.~(\ref{equ:analogy}) to
Eq.~(\ref{equ:converse}) we notice that the genuine analogue to $\tensor{Z^*_s}$
is  $1 - \tensor{\sigma}_s$ (and {\it not} $\tensor{\sigma}_s$), as indeed the
names ``effective'' vs. ``shielding'' imply. The analogy between the electric
and magnetic cases would be more direct if only the nucleus carried a
magnetic monopole charge; the NMR shielding would then be related to the
magnetic force on site $s$ induced by a macroscopic magnetic field,
just as $Z^*$ is the electric force induced at $\r_s$ by a macroscopic electric
field.

As in the electrical case~\cite{Ghosez98}, the choice of magnetic boundary
conditions corresponds to a choice of the {\it shape} of a macroscopic finite
sample. Following Ref.~\cite{Landau1}, the shape effects are embedded in the
depolarization coefficients ${n}_\alpha$ (with $\sum_\alpha {n}_\alpha=1$),
whose special cases are the sphere ($n_x$=$n_y$=$n_z$=1/3), the cylinder along
$z$ ($n_x$=$n_y$=1/2, $n_z$=0), and the slab normal to $z$ ($n_x$=$n_y$=0,
$n_z$=1). The main relationship for the macroscopic fields in Gaussian units may
be written as $B_\alpha = \be_\alpha + 4 \pi (1 - n_\alpha) M_\alpha$. It can be
seen that for the slab the normal component of $\B$ coincides with the one of
$\Be$. Hence our computed $\sigma_{s,zz}$ are suitable for direct comparison
with measurements of the normal component performed on a slab-shaped sample.
Assuming non-magnetic media with small, isotropic susceptibility $\chi$, it can
be shown that the shielding for a general shape is related to our  calculated
one by $\sigma_{s,\alpha\beta}^{\text{shape}}\simeq\sigma_{s,\alpha\beta} - 
\delta_{\alpha\beta}\,4\pi\chi(1-n_\beta)$. For the special case of a spherical
sample we have $\sigma_{s,\alpha\beta}^{\text{sphere}}\simeq\sigma_{s,\alpha
\beta} - (8\pi/3)\,\chi\,\delta_{\alpha\beta}$.

In order to calculate the shielding tensor of nucleus $s$ using
Eq.~(\ref{equ:converse}), it is necessary to calculate the induced orbital
magnetization due to the presence of an array of point magnetic dipoles $\m_s$
at all equivalent sites $\r_s$. The vector potential of a single dipole in
Gaussian units is given by~\cite{Jackson_dipole}
%
\begin{equation}\label{equ:A_sr}
\A^0(\r) = \frac{\m_s \times (\r - \r_s)}
{ |\r - \r_s|^3} \;.
\end{equation}
%
For an array of magnetic dipoles $\A(\r)=\sum_{\R}\A^{0}(\r-\R)$, where $\R$ is
a lattice vector. Since $\A$ is periodic, the average of its magnetic field
$\nabla\times\A$ over the unit cell vanishes; thus, the eigenstates of the
Hamiltonian remain Bloch-representable. In the Fourier representation
$\A(\r)=\sum_{\G\neq0}\tilde{\A}(\G)\,{\rm e}^{i\G\cdot\r}$ with
%
\begin{equation}
\tilde{\A}(\G) = -\frac{4 \pi i}{\Omega} \frac{\m_s \times
  \G}{G^2}\,{\rm e}^{-i \G \cdot \r_s},
\end{equation}
%
where the reciprocal lattice vector $\G=0$ may be excluded from the sum with no
loss of generality. Note that we have implicitly chosen the transverse gauge
$\nabla\cdot\A=0$, which is apparent from $\G\cdot(\m_s\times\G)=0$. The
periodic vector potential $\A(\r)$ can now be included in the Hamiltonian with
the usual substitution for the momentum operator $\p\rightarrow
\p-\frac{e}{c}\A$. As a result, the kinetic energy operator becomes
%
\begin{equation}\label{equ:new_terms}
\frac{p^2}{2m_{e}} \longrightarrow
 \frac{p^2}{2m_{e}} - \frac{e}{m_ec}\A\cdot\p + 
 \frac{e^2}{2m_ec^2}A^2\;, 
\end{equation}
%
where $m_{e}$ is the electronic mass and $c$ is the speed of light. Due to our
choice of gauge, $\p$ and $\A$ commute.  We can now calculate the shielding
according to Eq.~(\ref{equ:converse}) by solving for the ground state with the
additional terms of Eq.~(\ref{equ:new_terms}) included in the Hamiltonian, and
then calculating the resulting change in orbital magnetization \cite{A2}.

As an initial test, we have applied our converse approach to the calculation of
the NMR chemical shifts for small molecules in a supercell geometry, since these
can be directly compared with previous results obtained by Mauri \emph{et al.}\
using the direct method \cite{Mauri_96}.  Because we use a large supercell with
vacuum between molecules, the magnetic susceptibility of the system effectively
vanishes, and our computed $\tensor{\sigma}_s$ can be compared to the
experimental one without any shape correction.  For the same reason, the
recently-developed theory of orbital magnetization for extended systems
\cite{Resta_05, Thonhauser_06, Ceresoli_06, Niu} or its  single k-point
extension \cite{Ceresoli_07} is not needed, and we compute the induced
magnetization directly from the current distribution (paramagnetic plus
diamagnetic), which is
%
\begin{equation}
\label{equ:current}
\j(\r)= \sum_{n=1}^{\text{occ}}\Big[\frac{e\hbar}{m_e}\text{Im}
  \big[\psi_n^{*}(\r)\nabla\psi_n(\r)\big]
 -\frac{e^2}{m_ec}\A(\r)|\psi_n(\r)|^2\Big]\;.
\end{equation}
%
We implemented this approach into {\sc pwscf}, part of the {\sc
Quantum-ESPRESSO} distribution \cite{pwscf}. Here, the simplicity of the
converse method becomes obvious---the implementation  only required additional
coding of the order of 100 lines.  We used a local-density exchange-correlation
functional as in Ref.~\cite{Mauri_96}, and the structural parameters have been
taken from footnote [22] therein, in order to allow for a close comparison
between the direct and converse methods. We chose a dipole perturbation $|\m_s|$
of 1$\mu_\mathrm{B}$, although the results are numerically insensitive to 
dipoles in the range from $10^{-6}$ to $10^3$$\mu_\mathrm{B}$.

First, we used Coulombic potentials (i.e., a plane-wave all-electron approach)
with a high plane-wave energy cutoff of 500~Ry to calculate the $^1$H NMR
chemical shielding of several small molecules. Our  results are given in
Table~\ref{tab:results} together with experimental and theoretical values. It is
immediately obvious that the direct and converse methods give almost identical
results, validating our approach.


\begin{table}
\caption{\label{tab:results} Hydrogen NMR chemical shielding $\sigma$ in ppm for
several different molecules. Values for the direct method are taken from
Ref.~\cite{Mauri_96}. The pseudopotential results for the direct
method---according to footnote [21] of Ref.~\cite{Mauri_96}---have been
corrected by 0.3~ppm in order to allow a better comparison between the direct
and converse methods.}
\begin{tabular*}{\columnwidth}{@{}l@{\extracolsep{\fill}}ccccr@{}}\hline\hline
          &  exp.         &  \multicolumn{2}{c}{all-electron} &\multicolumn{2}{c}{pseudopotential}\\
          &               &    direct    &  converse          &  direct &  converse \\ \hline
H$_2$     &   26.26$\,^a$ &    25.9      &  25.9              &  25.6   &  25.6     \\
HF        &   28.51$\,^a$ &    28.4      &  28.5              &  28.1   &  28.1     \\
CH$_4$    &   30.61$\,^a$ &    30.7      &  30.7              &  30.4   &  30.4     \\
C$_2$H$_2$&   29.26$\,^b$ &    28.6      &  29.2              &  28.3   &  28.2     \\
C$_2$H$_4$&   25.43$\,^b$ &    24.5      &  25.4              &  24.2   &  24.0     \\
C$_2$H$_6$&   29.86$\,^b$ &    29.7      &  29.7              &  29.4   &  29.4     \\ \hline\hline
\end{tabular*}
\raggedright
$^a\,$Reference \cite{Raynes}.\\
$^b\,$Reference \cite{Schneider}.
\end{table}

Second, we performed these calculations in a pseudopotential approximation,
since plane-wave codes use pseudopotentials  to allow the use of a lower
kinetic-energy cutoff and, hence, greater computational efficiency.  Two
complications need to be taken into account in this case. First, a Bl\"ochl-like
PAW reconstruction needs to be performed, as shown by Pickard and Mauri
\cite{Pickard_Mauri_01}, for all elements beyond the first row. Since, at the
moment, we are only concerned with hydrogen shielding, we have neglected this issue. A
second and more subtle point concerns pseudopotentials that contain nonlocal
projectors (as it is usually the case when the Kleinman-Bylander separable form
is used \cite{kleinman-bylander}). This becomes apparent if the effect of
applying a constant vector potential $\A_0$ is considered. Such a gauge change
should not produce any observable effects, and indeed does not because of
cancellations between the substitution $\p\rightarrow \p-\frac{e}{c}\A_0$ for
the momentum operator and $\psi\rightarrow\psi e^{(ie/\hbar c)\A_0\cdot\r}$ for
all wavefunctions. However, if the pseudopotential is nonlocal, a similar
cancellation is ensured only if the nonlocal projectors $\ket{\beta_{lm}}$ are
phase-twisted according to \cite{Pickard_Mauri_01,ICL2001}
%
\begin{equation}\label{equ:phasetwist}
\sum_{lm}v_l\ket{\beta_{lm}}\bra{\beta_{lm}}\rightarrow
\sum_{lm}v_l\,e^{\frac{ie}{\hbar c}\A_0\cdot\r}\ket{\beta_{lm}}\bra{\beta_{lm}}
e^{-\frac{ie}{\hbar c}\A_0\cdot\r}\;,
\end{equation}
%
where the indices $l$ and $m$ correspond to orbital and azimuthal angular
momenta and $v_{l}$ is the Kleinman-Bylander coefficient. When the vector
potential is not constant, as in Eq.~(\ref{equ:A_sr}), the correct choice of
phase twist is generally not obvious. Nevertheless, it is a reasonable first
approximation to adopt Eq.~(\ref{equ:phasetwist}) with $\A_0$ replaced by
$\A(\r_i)$ for all nonlocal projectors located on the atom at $\r_i$. That is,
we assume that $\A(\r)$ is constant throughout the small core region in which
the projector differs from zero.  Because this approximation is not perfect,
however, the current density $\j(\r)$ still has an erroneous contribution
$\j_\text{err}(\r)$ inside the core regions. We have found that we can largely
correct for this remaining error by identifying $\j_\text{err}(\r)$ through its
unphysical property $\nabla\cdot\j_\text{err}(\r)\ne 0$. More precisely, we
calculate the first moment of the current divergence
%
\begin{equation}
\Delta\J_i = \int\limits_{\text{core }i} \!\! d^3r\;
\r\;\big(\nabla\cdot\j(\r)\big)\;,
\end{equation}
%
and add to the total current of Eq.~(\ref{equ:current}) a ``current
$\delta$-function'' of strength $-\Delta\J_i$ located at $\r_i$ before computing
the magnetization.

Taking all the above into account, we have again calculated the hydrogen NMR
chemical shieldings of the same molecules as before. Our pseudopotentials were
of the Troullier-Martins type \cite{TM}; for H we used a local pseudopotential,
and for all other atoms we included a nonlocal $s$-projector. We find
convergence of the shielding for a cutoff of 70~Ry, and our results for the
isotropic shielding are reported in Table~\ref{tab:results}. Again, it can be
seen that the direct and converse approaches give almost identical results.

We turn now to a discussion of the relative advantages and disadvantages of the
converse method compared with the direct method.  A potential disadvantage of
the converse method is that in order to calculate the shielding tensor for $N$
atoms, we need to perform $3N$ calculations. Using the direct method, on the
other hand, one obtains the shieldings of all $N$ atoms at once by only
performing three calculations. However, in many cases one is not interested in
the shielding of all atoms in the system, but rather only in a few, e.g., near a
binding site of a large molecule. In this case, the methods may be competitive.
But the main advantage of the converse method is really the simplicity of its
implementation, in that it works via finite differences of ground-state
calculations and does not require a linear-response implementation.  This is
likely to be a significant advantage for future applications in conjunction with
more complex forms of exchange-correlation functionals such as DFT+U, exact
exchange, or hybrid functionals.

In conclusion, we have derived an alternative first-principles method for
calculating NMR chemical shielding tensors.  The new approach is considerably
simpler, since difficulties concerning the choice of the gauge origin are
avoided and the need for a linear-response implementation is circumvented. We
have demonstrated the correctness and viability of our approach by calculating
the shieldings of simple molecular systems and finding excellent agreement with
previous theoretical and experimental results. Applications to more complex
systems, including crystalline systems for which the orbital magnetization will
need to be computed \cite{Resta_05,Thonhauser_06,Ceresoli_06,Niu}, are currently
in progress.

We are grateful for extensive discussions with F.~Mauri and D.~Hamann. This work
was supported by NSF grant DMR-0549198, ONR grant N00014-07-1-1095, and the
DOE/SciDAC project on Quantum Simulation of Materials and Nanostructures. 




\begin{thebibliography}{99}

\bibitem{NMR_encyclopedia} {\it Encyclopedia of NMR}, edited by D.M. Grant
   and R.K. Harris (Wiley, London, 1996).
  
\bibitem{Rabi} I. I. Rabi, J. R. Zacharias, S. Millman, and P. Kusch,
   Phys. Rev. {\bf 53}, 318 (1938).

\bibitem{Kutzelnigg_90} W. Kutzelnigg, U. Fleischer, and M. Schindler,
   \emph{NMR Basic Principles and Progress} (Springer, Berlin, 1990).

\bibitem{Mauri_96} F. Mauri, B. G. Pfrommer, and S. G. Louie,
   Phys. Rev. Lett. {\bf 77}, 5300 (1996).
   
\bibitem{Sebastiani_01} D. Sebastiani and M. Parrinello,
   J. Phys. Chem. A {\bf 105}, 1951 (2001).

\bibitem{Pickard_Mauri_01} C. J. Pickard and F. Mauri, Phys. Rev. B
   {\bf 63}, 245101 (2001); Phys. Rev. Lett. {\bf 91}, 196401 (2003).
   
\bibitem{Yates_07} J. R. Yates, C. J. Pickard, and F. Mauri, Phys. Rev. B
   {\bf 76}, 024401 (2007).

\bibitem{Truflandier_07} L. Truflandier, M. Paris, and F. Boucher,
   Phys. Rev. B {\bf 76}, 035102 (2007).
   
\bibitem{paratec} B. G. Pfrommer, J. Demmel, and H. Simon,
   J. Comp. Phys. {\bf 150}, 287 (1999); B. G. Pfrommer, M. Cot\'e, 
   S. G. Louie, and M. L. Cohen, J. Comp. Phys. {\bf 131}, 233 (1997).
   
\bibitem{castep} S. J. Clark, M. D. Segall, C. J. Pickard, P. J. Hasnip, M. J.
   Probert, K. Refson, M. C. Payne, Zeitschrift f\"ur Kristallographie
   {\bf 220(5-6)}, 567 (2005).
   
\bibitem{Cai_07} W. Cai and G. A. Galli, Phys. Rev. Lett. {\bf 92}, 186402
   (2004); E. Lee, W. Cai, and G. A. Galli, J. Chem. Phys., in print
   (2007).
   
\bibitem{Resta_05} R. Resta, D. Ceresoli, T. Thonhauser, and
   David Vanderbilt, ChemPhysChem {\bf 6}, 1815 (2005). 

\bibitem{Thonhauser_06} T. Thonhauser, D. Ceresoli, D.
   Vanderbilt, and R. Resta, Phys. Rev. Lett. {\bf 95}, 137205 (2005).

\bibitem{Ceresoli_06} D. Ceresoli, T. Thonhauser, D.
   Vanderbilt, and R. Resta, Phys. Rev. B {\bf 74}, 24408 (2006).
   
\bibitem{Niu} Di Xiao, J. Shi, and Q. Niu, Phys. Rev. Lett. {\bf 95},\
   137204 (2005); J. Shi, G. Vignale, Di Xiao, and Qian Niu,
   arXiv:0704.3824v1, unpublished.

\bibitem{Born} R. Resta, Rev. Mod. Phys. {\bf 66}, 899 (1994).

\bibitem{Ghosez98} Ph. Ghosez, J.-P. Michenaud, and X. Gonze, Phys. Rev. B
   {\bf 58}, 6224 (1998).

\bibitem{Landau1} L. D. Landau and E. M. Lifshitz, {\it Electrodynamics of
   Continuous Media} (Pergamon Press, Oxford, 1960).

\bibitem{Jackson_dipole} J. D. Jackson, {\it Classical Electrodynamics,
   2$^{\rm nd}$ ed.} (Wiley, New York, 1975).

\bibitem{A2} Insofar as our goal is to extrapolate the finite-difference
   results to reproduce the linearly-induced magnetization, the $A^2$ term
   truly has no influence, and may be included or not at one's convenience.

\bibitem{Ceresoli_07} D. Ceresoli and R. Resta, Phys. Rev. B {\bf 76},
   012405 (2007).

\bibitem{pwscf} S. Baroni \emph{et al.},
   % A. Dal Corso, S. de Gironcoli, P. Giannozzi, C.
   % Cavazzoni, G. Ballabio, S. Scandolo, G. Chiarotti, P. Focher, A.
   % Pasquarello, K. Laasonen, A. Trave, R. Car, N. Marzari, A. Kokalj,
   http://www.quantum-espresso.org.
   
\bibitem{Raynes} W. T. Raynes, \emph{Nuclear Magnetic Resonance} (The Chemical
   Society, London, 1977), Vol.\ 7.

\bibitem{Schneider} Values of $\sigma$ relative to CH$_4$ are reported in
   W. G. Schneider, H. J. Bernstein, and J. A. Pople,
   J. Chem. Phys. {\bf 28}, 601 (1958). Reference~\cite{Raynes} can then
   be used to convert to the absolute scale.

\bibitem{kleinman-bylander}
L. Kleinman and D. M. Bylander, Phys. Rev. Lett. {\bf 48}, 1425 (1982).

\bibitem{ICL2001}
  S. Ismail-Beigi, E. K. Chang, and S. G. Louie, Phys. Rev. Lett. {\bf
  87}, 087402 (2001).

\bibitem{TM} N. Troullier and J. L. Martins, Phys. Rev. B {\bf 43},
   1993 (1991).

\end{thebibliography}




\end{document}
