{}

$missing\text{\,}\mathrm{missing}$

# Towards quantum utility for NMR quantum simulation on a NISQ computer

Artemiy Burov  School of Life Sciences, University of Applied Sciences
Northwestern Switzerland (FHNW), Hofackerstrasse 30, CH-4132 Muttenz,
Switzerland  École Polytechnique Fédérale de Lausanne (EPFL), CH-1015
Lausanne, Switzerland  MatterDecoder, Lausanne, CH-1015, Switzerland  Oliver
Nagl  School of Life Sciences, University of Applied Sciences Northwestern
Switzerland (FHNW), Hofackerstrasse 30, CH-4132 Muttenz, Switzerland  ETH
Zurich, CH-8093 Zurich, Switzerland  Clément Javerzac-Galy
<EMAIL_ADDRESS>School of Life Sciences, University of Applied
Sciences Northwestern Switzerland (FHNW), Hofackerstrasse 30, CH-4132 Muttenz,
Switzerland  MatterDecoder, Lausanne, CH-1015, Switzerland

###### Abstract

While the recent demonstration of accurate computations of classically
intractable simulations on noisy quantum processors brings quantum advantage
closer, there is still the challenge of demonstrating it for practical
problems. Here we investigate the application of noisy intermediate-scale
quantum devices for simulating nuclear magnetic resonance (NMR) experiments in
the high-field regime. In this work, the NMR interactions are mapped to a
quantum device via a product formula with minimal resource overhead, an
approach that we discuss in detail. Using this approach, we show the results
of simulations of liquid-state proton NMR spectra on relevant molecules with
up to $11$ spins, and up to a total of $47$ atoms, and compare them with real
NMR experiments. Despite current limitations, we show that a similar approach
will eventually lead to a case of _quantum utility_ , a scenario where a
practically relevant computational problem can be solved by a quantum computer
but not by conventional means. We provide an experimental estimation of the
amount of quantum resources needed for solving larger instances of the problem
with the presented approach. The polynomial scaling we demonstrate on real
processors is a foundational step in bringing practical quantum computation
closer to reality.

##  I Introduction

The demonstration of a useful application of the noisy intermediate-scale
quantum (NISQ) technology is of substantial interest to the quantum technology
industry. It is still an open question whether the utility of quantum
computing can be achieved in the pre-fault-tolerant era. In this work, we show
a path towards this possibility. As recently outlined by Kim, Eddins et al.
Kim  _et al._ (2023), quantum advantage can be approached in two ways.
Firstly, by engineering artificial classically hard problems that can be
efficiently solved on a quantum device, as was demonstrated in Kim  _et al._
(2023). Secondly, by finding practical problems with the same requirement:
classically hard problems that can be efficiently solved on a quantum device.
While the first way has been demonstrated, albeit it was recently shown that
the exact problem statement considered in Kim  _et al._ (2023) is in fact not
classically sufficiently hard to construct a bulletproof argument Begušić  _et
al._ (2024), we focus on the second one.

Here we discuss the problem of Heisenberg Hamiltonian simulation, an important
building block for a set of problems revolving around spin dynamics of nuclear
magnetic resonance (NMR) (figure 1 a and b) previously discussed in the
context of zero-field NMR Sels  _et al._ (2020); Seetharam  _et al._ (2023);
O’Brien  _et al._ (2022). Compared to some other quantum chemistry tasks this
application offers a particularly favorable polynomial scaling, as illustrated
in figure 1 c. Simulation of fermionic systems described by second
quantization Hamiltonians is interesting in particular for understanding the
mechanisms behind high-energy superconductivity or in many electronic quantum
chemistry problems but requires a costly encoding (for example the Bravyi-
Kitaev transformation) to be mapped to a quantum circuit. The task of finding
low energy eigenvalues of molecules is another practically interesting problem
statement, but this is a _BQNP-hard_ problem. BQNP (Bounded-error Quantum Non-
deterministic Polynomial-time, closely related to well-known QMA) is the
quantum counterpart of the classical NP complexity class as described in
Kitaev  _et al._ (2002), and resources required for finding solutions to such
BQNP-hard problems scale superpolynomially with the size of the instance. This
type of quantum computation is to our knowledge beyond the capabilities of
quantum computers available today and in the next few years for practically
interesting problem instances.

Figure 1: (a) shows the main characters of this work: either the NMR
spectrometer or the quantum processor generates spectra of a given quantum
spin system. (b) shows the pathway of NMR Hamiltonian simulation (quantum or
classical). The signal is calculated for a set of time points by applying the
time propagator $\exp{\frac{-i\hat{H}t}{\hbar}}$ to the initial state
$\ket{\psi(0)}$. The propagator can be applied either analytically or by a
quantum calculation. (c) shows the scaling of various mathematical objects
involved in the discussion. The Hamiltonian scales exponentially - indicating
the need for an exponential amount of classical resources required for the
simulation. Required quantum resources scale only polynomially since the
product formula quantum circuit scales polynomially. (d) shows a part of the
mapping between a molecule (3-Amino-5-ethyl-1H-pyrazole) and its quantum
circuit for the FID signal calculation. The spin-field couplings are
represented by single qubit rotations around the $Z$ axis, and the spin-spin
couplings are represented by rotations around $XX$, $YY$, and $ZZ$ axes of the
corresponding two-qubit Hilbert space.

In the proposed nuclear magnetic spin dynamics simulations, the objects of
interest are spatially localized fermions, which can be described with similar
mathematics as qubits. This removes any encoding overhead.

Furthermore, NMR is in fact the grandfather of quantum computers, ultimately
going into early retirement as a quantum computing architecture facing
scalability and noise issues, but leaving us with a whole array of spin
manipulation techniques that are used in quantum computing today Vandersypen
and Chuang (2005). It is no surprise that a system that can be used as a
quantum computer is itself well suited to be simulated by a quantum computer.

Apart from this conceited interest from the quantum computing community, we
believe that the present set of problems is also of profound interest to the
NMR community, in particular for applications such as chemical structure
verification and elucidation via NMR spectroscopy. NMR spectroscopy is widely
used in many fields of science and industry of Japan (2017), such as organic
chemistry, medicine, agriculture, food chemistry, and environmental sciences.
This work could also ultimately impact computationally hard NMR problems in
catalyst and battery development Rei (2021).

As experience and previous work Sels  _et al._ (2020); Seetharam  _et al._
(2023); O’Brien  _et al._ (2022) show, the simulation of NMR spectra is a
crucial building block for these practical applications. It can benefit from
quantum information processing, since this problem is equivalent to simulating
the action of a Hamiltonian on a quantum system, which is believed to be a
non-trivial problem, though among the easiest to address with early quantum
devices Childs  _et al._ (2018). The idea of simulating microscopic nature, a
natively quantum problem, on a quantum machine, was one of the initial
motivations behind the development of quantum computing Manin (1980); Feynman
(1982).

Figure 2: Illustrates the concepts of the spin-spin coupling matrix and the
spin-spin coupling graph. The presented molecules are Butyraldehyde (a, b, and
c) and Adenosine triphosphate (d, e, and f). The high connectivity of spins in
Butyraldehyde can be seen by the number of connections in the corresponding
graph c) and the number of off-diagonal elements in the matrix b). The number
of off-diagonal elements is proportional to the number of _two-qubit gates_
needed for the simulation of the spin cluster, an essential resource in
quantum computing. Block-diagonal matrices such as seen in e) indicate a spin
system that is separable into independent subsystems, each with independent
unitary dynamics.

In this work, we show the simulation of high-field NMR free induction decay
(FID) experiments on a quantum processor (QPU). To our knowledge, this is the
first such simulation. We compare the results with the output of NMR
experiments performed on a $\quantity{400}{}$ instrument.

We first explain the principles of the NMR FID experiment and discuss cases
when this problem becomes computationally complex. Without loss of generality,
we only consider 1H hydrogen spins in this work, in which case one speaks of
_proton_ NMR. The way to map this experiment on a QPU is then discussed. This
includes the discussion of the quantum algorithm and the technicalities of
implementing this algorithm on quantum hardware. We then show the results of
the quantum simulation for some of the hardest cases that we can currently
simulate with the proposed approach, around $10$ spins, depending on the
structure of the molecule. The last results in Seetharam  _et al._ (2023)
showed simulations for a molecule containing a methyl group with $3$ hydrogen
atoms that were probed in a zero-field experiment. We also cover some topics
that may be discussed in more detail in future work. In particular, the
resources needed for harder cases of the NMR FID simulation problem are
covered.

##  II NMR Theoretical Background

The NMR Hamiltonian describes the dynamics of the magnetic system of nuclear
spins. Two types of interactions define these dynamics: the _spin-field_
interaction that gives rise to _chemical shifts_ in the NMR spectrum and the
_spin-spin_ interaction that causes splittings in the chemical shifts - or at
least such is the heuristic point of view that is commonly assumed in
interpreting spectra obtained in _high field_ _liquid-state_ NMR experiments.
In more general cases, such as _near-zero field_ NMR Ledbetter  _et al._
(2011), _solid-state_ NMR Rei (2021), and also in the case of _strongly
correlated_ high-field liquid-state NMR Appelt  _et al._ (2010), this
heuristic point of view and the corresponding heuristic methods of NMR spectra
analysis break down.

It is in general necessary to calculate the action of the Hamiltonian on the
Hilbert space of the spin system to get information about its dynamics. This
calculation is computationally hard - in general, it scales exponentially with
the size of the spin system, as does the Hamiltonian matrix, see figure 1 c.
It is mathematically equivalent to applying an exponential of a $2^{N}$ by
$2^{N}$ sparse matrix to a vector, where $N$ is the size of the spin system.
In many cases, numerical methods can be employed to perform this computation
efficiently Kuprov (2023). This is especially true for the cases where $N$ is
small or the spins in the molecule are weakly connected.

In the most difficult cases, such as molecules with highly interconnected
spins, the numerical methods fail Karabanov  _et al._ (2011). A good example
where exponential complexity is unavoidable is the NMR of single crystals,
famously calcium fluoride Parker and Lado (1974). For this reason, the
Hamiltonian simulation was proposed as a solution Sels  _et al._ (2020);
Seetharam  _et al._ (2023); O’Brien  _et al._ (2022).

### Spin connectivity.

To visualize the connectivity of the molecule, it helps to introduce the
concepts of the spin-spin coupling matrix and the spin-spin coupling graph. In
figure 2, two molecules are presented as examples, both with $8$ 1H atoms in
the considered spin systems.

The first example, Butyraldehyde, is a relatively small molecule with a highly
connected spin system. This can be seen in the number of off-diagonal elements
in the spin-spin coupling matrix, each of them corresponds to a spin-spin
coupling. It also can be seen in the highly connected spin-spin coupling
graph.

The second molecule, adenosine triphosphate (ATP), has more 1H atoms than the
$8$ in the shown spin system. Often, a spin will be decoupled from other
spins, for example, due to being in a bond with an oxygen atom. This results
in a 1H atom heavily shielded from its environment, which frequently exchanges
with solvent 1H atoms at higher temperatures Exarchou  _et al._ (2002). The
dynamics of the whole system can be considered as the set of separate unitary
evolutions of non-interacting subsystems of the whole system. For the NMR
spectra, it effectively means that they can be produced as a sum of the
spectra of all the non-interacting (or weakly interacting) subsystems of the
whole spin system - the spin clusters.

Two of the non-interacting atoms are included in the considered spin system.
This can be seen in the spin-spin coupling matrix, as it is block-
diagonalizable. The blocks in this representation correspond to the
interacting subsystems. It also can be seen in the disconnected spin-spin
coupling graph.

When the spin system together with its interactions is mapped to a quantum
circuit, the amount of the most expensive computational resource needed (in
terms of noise) - the two-qubit quantum gates - is directly proportional to
the number of the non-zero off-diagonal elements. In the worst-case scenario
of the fully connected spin system, the number of the two-qubit quantum gates
scales as the number of elements in the upper half of the spin-spin coupling
matrix, so it is quadratic in the size $N$ of the spin system. This is rare in
liquid-state NMR.

### NMR FID experiment.

In the NMR FID experiment, the signal is produced by the magnetization of
nuclear spins of a molecular system of interest precessing around a magnetic
field Levitt (2013). Usually, only one species of nuclei is excited, for
example, the 1H of the molecule in the case of proton NMR, as considered in
this work without loss of generality.

This FID signal allows us to gain information about the structure of the
studied sample. Mathematically the structure of the sample is contained in the
Hamiltonian $\hat{H}$ of the studied system. The Hamiltonian also defines the
dynamics of the studied system via the Schrödinger’s equation, provided for
the case of a pure state $\ket{\psi}$ by

$i\hbar\derivative{t}\ket{\psi(t)}=\hat{H}\ket{\psi(t)}.$ (1)

For the case of time-independent Hamiltonian, the formal solution for the time
evolution can be written simply as

$\ket{\psi(t)}=\exp{\frac{-i\hat{H}t}{\hbar}}\ket{\psi(0)}.$ (2)

The exponential of the Hamiltonian matrix in eq. (2) is hard to compute. It is
this computational barrier that hopefully can be traversed using quantum
computing.

Once the state $\ket{\psi(t)}$ is stored in a register of a quantum computer,
a projective quantum measurement of the total magnetization observable
$\hat{M}_{X}=\sum^{N}_{k=1}I^{X}_{k}$ can be performed to obtain the value of
magnetization at time $t$, where $I^{X}_{k}$ is the $X$ component of the spin
operator for the $k$-th spin. The measurement collapses the quantum state into
one of the eigenspaces of the observable, and the information about its
previous evolution is lost. For this reason, the evolution of the state has to
be performed for each recorded time point of the simulated FID. This seemingly
large overhead is however only constant, that is, it does not depend on the
size of the simulated spin system.

The NMR Hamiltonian is often described as a Heisenberg Hamiltonian Levitt
(2013). It is commonly written as

$\hat{H}=-\sum^{N}_{k=1}\gamma\overline{I}_{k}(1-\sigma_{k})\overline{B}_{0}+2\pi\sum_{k<l}\overline{I}_{k}J_{kl}\overline{I}_{l}$
(3)

for a liquid-state NMR experiment with spin $\frac{1}{2}$ nuclei, which is the
case considered in this work. Term-by term it can be understood as
$\gamma\overline{I}_{k}\overline{B}_{0}=\omega_{0}\overline{I}_{k}$ describes
the interaction of unshielded nucleus with the magnetic field, where $\gamma$
is the gyromagnetic ratio of the studied nuclei, hydrogens in our case.
$\overline{I}_{k}$ is the vector of spin operators. $\overline{B}_{0}$ is the
magnetic field vector. $-\gamma\overline{I}_{k}\sigma_{k}\overline{B}_{0}$
describes the modification to this interaction due to the shielding by the
electrons of the molecule, the _chemical shielding_ , where $\sigma_{k}$ is
the chemical shift. $2\pi\sum_{k<l}\overline{I}_{k}J_{kl}\overline{I}_{l}$
describes the spin-spin interaction, where $J_{kl}$ is the J-coupling. We
write $I_{k}$ here and in the following without the characteristic operator
hat to avoid unnecessary cluttering of text. Nevertheless, it is to be
understood as an operator, and $\overline{I}_{k}$ as a vector of operators. In
general, both $\sigma_{k}$ and $J$ can be a tensor, in which case the
interaction is called _anisotropic_. In this work, we consider only isotropic
interactions commonly encountered in liquid-state NMR.
$\overline{I}_{k}J_{kl}\overline{I}_{l}$ is a bilinear form, and in the
isotropic case it can be rewritten simply as
$J_{kl}\overline{I}_{k}\overline{I}_{l}$ since $J_{kl}$ is a scalar.

We will now rewrite this Hamiltonian in a format more suitable for the desired
simulations. Without loss of generality, the field $\overline{B}_{0}$ can be
considered to point along the $Z$ direction. Also, as in this work, we assume
interactions to be isotropic, both $\sigma_{k}$ and $J$ are scalars. With
these considerations, we can rewrite the Hamiltonian as

$\hat{H}=\sum^{N}_{k=1}\omega_{0k}I^{Z}_{k}+2\pi\sum_{k<l}J_{kl}\overline{I}_{k}\overline{I}_{l}$
(4)

with _shifted_ Larmor frequency of $k$-th nucleus given by $\omega_{0k}$.

Finally, for many applications, including the simulation presented in this
work, it helps to consider the Hamiltonian in the rotating frame, where
$\omega_{k}$ now is the shifted Larmor frequency in the rotating frame:

$\hat{H}^{r}=\sum^{N}_{k=1}\omega_{k}I^{Z}_{k}+2\pi\sum_{k<l}J_{kl}\overline{I}_{k}\overline{I}_{l}$
(5)

Indeed, since the generator of rotations around the $Z$ axis given by
$\sum_{i}I^{Z}_{i}$ commutes with all terms containing the $Z$ angular
momentum operators $I^{Z}_{k}$ of the Hamiltonian $\hat{H}$, and additionally
that it also commutes with all other terms that come up in this Hamiltonian,

$\begin{split}[\sum_{i}I^{Z}_{i},I^{X}_{k}I^{X}_{l}+I^{Y}_{k}I^{Y}_{l}]=[I^{Z}_{k}+I^{Z}_{l},I^{X}_{k}I^{X}_{l}+I^{Y}_{k}I^{Y}_{l}]=\\\
[I^{Z}_{k},I^{X}_{k}I^{X}_{l}]+[I^{Z}_{l},I^{X}_{k}I^{X}_{l}]+[I^{Z}_{k},I^{Y}_{k}I^{Y}_{l}]+[I^{Z}_{l},\hat{\sigma}^{Y}_{k}I^{Y}_{l}]=\\\
iI^{Y}_{k}I^{X}_{l}+iI^{X}_{k}I^{Y}_{l}-iI^{X}_{k}I^{Y}_{l}-iI^{Y}_{k}I^{X}_{l}=0,\end{split}$
(6) $\begin{split}[\sum_{i}I^{Z}_{i},I^{Z}_{k}I^{Z}_{l}]=0\end{split}$ (7)

we conclude that the Hamiltonian is invariant under rotations around the $Z$
axis. Unsurprisingly, since the unitary evolutions of all connected subsystems
of the spin system can be considered separately, the Hamiltonian of the full
spin systems also commutes with a transformation that involves different
rotations for each subsystem. This can also be seen from eq. (6) by grouping
on the right side of the commutator $[\sum_{i}I^{Z}_{i},\hat{H}^{r}]$ all
$2\pi\sum_{k<l}J_{kl}\overline{I}_{k}\overline{I}_{l}$ terms that belong to
the same connected subsystem. In that case on the left side of the commutator,
the whole sum over all spins is not necessary - it is enough to sum over the
spins of that subcluster to get the commutator to be equal to $0$. This sum
gives a generator of rotations of that subcluster, and following this
procedure for other subclusters we receive a generator of rotations that
commutes with the Hamiltonian of the full system for each subcluster.

In this work we consider the initial state of the spin system to be a pure
state. The spins are initially aligned with the $Z$ axis and then turned $90$
degrees around the $Y$ axis to point along the $X$ axis with a $\frac{\pi}{2}$
pulse.

The FID signal is obtained as an ensemble measurement of the total
magnetization operators along $X$ and $Y$ axes at different points in time.
These operators can be written as $\hat{M}_{X}=\sum^{N}_{k=1}I^{X}_{k}$ and
$\hat{M}_{Y}=\sum^{N}_{k=1}I^{Y}_{k}$, respectively.

Figure 3: Shows the cosine distance between the spectra obtained by (a) Lie-
Trotter and (b) Suzuki-Trotter product formulas and by exact Hamiltonian
exponentiation methods. The calculation is done on $4$-spin molecules from the
GISSMO Dashti  _et al._ (2018) dataset. The grey bar plots in (a) show the
scaling of the depth of two-qubit gates in the quantum circuit needed for the
FID simulation, this can be thought of as the computational cost of the
simulation. Figure 4: a) The left y-axis (red and blue data points, arrows
pointing to the corresponding axis for visibility) shows the total number and
depth of two-qubit echoed cross-resonance (ECR) gates needed to simulate a
single repetition of a product formula for a given number of spins. The
molecules are taken from the GISSMO database Dashti  _et al._ (2018). The
quantum circuits are transpired to the IBM ”brisbane” quantum device (Eagle
r3). The transpilation algorithm is probabilistic, but in the considered cases
the deviations in the number of gates for single molecules are low, in the
order of $10$. The deviations in the box plots come mostly from the fact that
the number of required gates depends on the connectivity of the molecule, and
there are multiple molecules with the same number of hydrogen spins. In the
extreme case where the molecule is not connected at all, the gate depth can be
as low as $2$. The right y-axis (grey data points, arrow pointing to the
corresponding axis for visibility) shows an estimation of the required runtime
on a classical computer, calculated by taking the matrix exponential of the
Hamiltonian, as a visual guide. We want to stress that for weakly correlated
spin clusters these computations can be sped up by employing numerical methods
to take the exponential of the Hamiltonian matrix, but in general cases, the
computational complexity of this task scales up exponentially (see figure 1).
b) Examples of molecules with 1H atoms highlighted in blue. These molecules
are chosen as representative examples with the average number of atoms for a
given number of spins, taken from the GISSMO database Dashti  _et al._ (2018).
Molecules considered even in this regime with low numbers of hydrogen atoms
can get quite large. Data for larger molecules is, at least in the GISSMO
database, scarce. Figure 5: Shows the $\quantity{400}{}$ FID signals and
spectra of two spin clusters obtained from IBM Eagle r3 ”brisbane” quantum
device (red). The spin clusters are (a, b, c) ATP ($8$ spins, low
connectivity) and (d, e, f) 4-methylvaleric acid ($11$ spins, high
connectivity). The spectra show correspondence to the reference spectra
obtained experimentally (grey). In some cases, as shown in the inset in (f),
the J-coupling structure can be reproduced in the simulation. Despite that,
the noise floor prevents obtaining the signals from all peaks, and the
artifacts from using a product formula are also visible. The FID signals
consist of $4096$ temporally equidistant points. The expectation value of the
$X$ magnetization operator is calculated by taking $4000$ shots of the quantum
circuit. The shown FID is an average of over $5$ runs on the quantum device.

##  III Implementation and results

In this work, we used one repetition of the Lie-Trotter product formula as the
algorithm for Hamiltonian simulation. To our knowledge, this method requires
the least amount of quantum operations - _gates_ \- and is easy to implement,
but it is also a very rough approximation. It is not common to use product
formulas in this way, usually the time evolution is discretized with some time
step $\delta t$. This time step is chosen by weighing the cost of introducing
smaller time steps versus the benefit in the accuracy that the smaller time
step provides. The reason for our choice of taking only _one_ time step is
that the main goal of this work is to obtain NMR signals from some of the
largest spin systems that can be simulated on quantum processors today. The
model error stemming from the choice of the product formula is secondary to
this work.

Another, quantitative justification for this approach is detailed in figure 3.
Here we compare two product formulas, the Lie-Trotter and Suzuki-Trotter
formulas. We analyzed a subset of the GISSMO database Dashti  _et al._ (2018)
consisting of $4$-spin molecules with non-trivial connectivity to get an idea
of the scaling of the model error due to the product formula approximation. In
our investigations, we have found no reason to believe that in this regime of
low number of repetitions increasing the number of repetitions leads to a
better model.

As a metric for the error, we used the cosine distance between the spectra
approximated by a product formula and the spectra calculated from first
principles by taking the exponential of the Hamiltonian. The cosine distance
treats spectra as high-dimensional vectors and is defined as
$1-\frac{\overline{x}\cdot\overline{y}}{\absolutevalue{\overline{x}\cdot\overline{y}}}$
for $\overline{x}$ and $\overline{y}$ \- two sets of points that define the
spectra. This metric was chosen because it allows us not to worry about
scaling the spectra to the same magnitude, which in this case is not a well-
defined operation since the spectra may have quite different shapes. The
cosine distance is bound by 1 for low similarity and 0 for high similarity of
spectra.

We found that the error increases initially up to $4$ repetitions for the Lie-
Trotter formula (fig. 3 a) and up to the fourth order for the Suzuki-Trotter
formula (fig. 3 b) for the considered subset of the GISSMO database. Given
that the number of required gates and the gate depth grows linearly with the
number of repetitions (fig. 3 a), and we already are operating at the limits
of the NISQ hardware available today, we concluded that in the current regime,
it makes sense to limit the number of Trotter repetitions to just one. Recent
results on optimizing Trotter parameters could lead to better results in the
future Mansuroglu  _et al._ (2023); Mc Keever and Lubasch (2023); Trenkwalder
_et al._ (2023).

### Product formula.

The order of gates in the product formula is chosen as follows:

$\begin{split}e^{-it(\sum^{N}_{k=1}\omega_{k}I^{Z}_{k}+2\pi\sum_{k<l}J_{kl}\overline{I}_{k}\overline{I}_{l})}=e^{-2it\pi\sum_{k<l}J_{kl}I^{Z}_{k}I^{Z}_{l}}*\\\
e^{-it\sum^{N}_{k=1}\omega_{k}I^{Z}_{k}}*\\\
e^{-2it\pi\sum_{k<l}J_{kl}I^{Y}_{k}I^{Y}_{l}}*\\\
e^{-2it\pi\sum_{k<l}J_{kl}I^{X}_{k}I^{X}_{l}}\end{split}$ (8)

The reasoning here is largely heuristic. If all Hamiltonian terms would
commute with each other, the product formula would be not an approximation but
an identity. This is not the case here, instead, the Hamiltonian terms that
_do_ commute with each other are grouped into $4$ groups, as seen in equation
8. Additionally, $\sum^{N}_{k=1}\omega_{k}I^{Z}_{k}$ always commutes with
$2\pi\sum_{k<l}J_{kl}I^{Z}_{k}I^{Z}_{l}$ and sometimes commutes with
$2\pi\sum_{k<l}J_{kl}I^{Y}_{k}I^{Y}_{l}+2\pi\sum_{k<l}J_{kl}I^{X}_{k}I^{X}_{l}$,
see equation 6.

### Gate counts - quadratic scaling.

The size of the quantum circuit required to implement this product formula is
displayed in figure 4. The circuits were transpiled to IBM’s $127$ qubit
architecture quantum chip and optimized with a standard Qiskit Qiskit
contributors (2023) pipeline, which includes some basic optimizations such as
combining sequences of one-qubit gates into one gate and rewriting rules based
on commutation relations of quantum operations.

The chemical shift and spin-spin coupling data was taken from the GISSMO
Dashti  _et al._ (2018) database. It consists of $1324$ entries of molecules,
their (1H) spin-spin coupling matrices, and experimental and simulated
spectral data. The spin clusters present in this database are mostly in the
range of up to $20$ spins.

The calculation in figure 4 is done based on the complete GISSMO database, the
box plots show the distribution of required gate counts for molecules
consisting of up to $20$ spins, and several more data points are available for
larger molecules. The dispersions in values for gate counts come from the fact
that the structure of the quantum circuit needed for FID simulation depends on
the structure of the simulated molecule, and there are many molecules with the
same number of spins in the GISSMO database, for instance, there are $93$
molecules with $10$ proton spins. The dispersion from different compilations
of the same quantum circuit is negligible. The higher the connectivity of the
spin cluster of the molecule, the higher the depth of the needed quantum
circuit. The presented numbers of gates are calculated for the IBM ”brisbane”
quantum device (Eagle r3 architecture).

The two-qubit gate requirements for the NMR FID simulation for molecules from
the GISSMO database were analyzed. Figure 4 a shows the total number and depth
of two-qubit echoed cross-resonance (ECR) gates needed to simulate a single
repetition of a product formula for a given number of spins. A quadratic fit
was used, matching the data quite well. We do expect the quadratic behavior
for the total number of two-qubit gates: it is proportional to the number of
off-diagonal elements in the spin-spin coupling matrix, which is quadratic in
the number of spins $N$.

The classically intractable regime starts as early as $N=20$. From the figure
4 we can see that in this regime the expected depth of two-qubit gates is
around $100$, with more connected spin clusters requiring larger gate depth,
potentially up to $300$ gates. The successful execution of quantum circuits
with this number of gates is feasible to expect in the next few years. The
larger instance of the two simulated in this work has the two-qubit gate depth
around $150$ with a small deviation in the order $10$. In our investigations,
we found that no significant signal can be achieved reliably for larger
instances. Given the estimated quadratic scaling of the quantum resources
needed, we believe that the simulations of Heisenberg Hamiltonian that are
inaccessible via classical computation have a good potential to be solvable on
quantum devices when quantum devices can execute circuits of depth in the
order of $1000$ two-qubit gates.

It is worth mentioning that while the regime relevant for quantum utility
starts at around $20$ (nuclear) spins, the actual molecules are often much
bigger than $20$ atoms. Examples of molecules with the average number of atoms
for several indicated spin counts are given in figure 4 b as a reference. This
shows that molecules chemically relevant for useful NMR applications could be
simulated already in the NISQ era.

### Classical exponential barrier.

As a guide to the eye, an estimation of classical computation times is given
in figure 4 a in grey. This estimation was made by a straightforward no-tricks
exponentiation of the Hamiltonian matrix on an Intel i7-1165G7 2.80GHz
computer. It is important to keep in mind that for the presented cases the
simulation can be done much faster with many spin dynamics simulation
packages, for example with SPINACH Hogben  _et al._ (2011); Kuprov (2023). But
the exponential scaling still holds for large and highly connected spin
clusters - for these cases numerical methods fail.

The gate counts in figure 4 were estimated by optimizing the circuits using
the same pipeline also used for optimizing and transpiling the circuits to the
quantum hardware for the presented QPU simulations.

### QPU experiment.

In this work, we have simulated the NMR spectra of two spin clusters on a QPU
as shown in figure 5. In the circuit model, the qubits are usually initialized
to point along the $Z$ axis. For the sake of this work, we assume this also to
be the case here, which means that as a first step, we need to initialize the
qubit register to correspond to the state of the spins after the
$\frac{\pi}{2}$ pulse. This step can be skipped by relabeling the coordinate
system in such a way as to exchange the $X$ and $Y$ axes. Then we can assume
that the spins point along the $X$ axis immediately after initialization. We
kept this step in to keep the correspondence to the actual NMR experiment
since it only changes the gate depth by one.

After that, the gates that simulate the interactions of spins with the field
and between each other are added according to the chosen product formula.
Finally, the transverse magnetization operator $\sum^{N}_{k=1}I^{X}_{k}$ is
measured along any direction in the transverse plane, we chose the $X$
direction.

The mapping can be seen schematically in figure 1 d for the example of
3-Amino-5-ethyl-1H-pyrazole.

The quantum circuits for the experiments presented in this work were run on
the IBM Eagle r3 ”brisbane” quantum chip. It consists of $127$ fixed-frequency
transmon qubits with heavy-hex connectivity and median $T_{1}=222\mu s$ and
$T_{2}=145\mu s$.

The FID signal (fig. 5 b and e) is obtained by taking $4096$ equidistant
points of time, calculating the evolution of the quantum state, and measuring
the expectation value of the total $X$ magnetization operator for each of
those points. The expectation value is calculated by taking $4000$ shots. The
FID signal is time-averaged over $5$ such runs on the quantum hardware.
Overall, repeating the experiment $N$ times, one expects an improvement in the
signal-to-noise ratio by a factor of $\sqrt{N}$.

The choice of the sampling interval between the time points is guided by the
range of frequencies that need to be captured for the NMR signal in the
rotating frame according to the Nyquist theorem. In the case of proton NMR
with the magnetic field equal to {400} this amounts to the frequency range
$[0,4000]$ Hz. This gives the sampling rate of $8000$ Hz and the sampling
interval of $0.125$ ms. In NMR, the magnetic field is commonly and
conveniently measured in terms of the Larmor frequency $\omega_{0}=\gamma
B_{0}$ of a given nucleus type.

The spectrum (fig. 5 c and f) is obtained by applying the Fourier
transformation to the FID signal.

We have simulated the NMR spectra of two spin clusters. The larger spin
cluster is the 4-methylvaleric acid C6H12O2 (fig. 5 a), an organic acid with
$11$ proton spins. The corresponding NMR FID simulation needs $11$ qubits and
the two-gate depth in the range of $140-160$, where the depth of gates is
probabilistic since the optimization algorithm for the quantum circuit is also
probabilistic.

The other spin cluster is ATP (fig. 5 a), containing $8$ spins. The two-qubit
gate depth needed for the simulation is in the range of $35-40$.
Interestingly, ATP is found in all known forms of life and is often referred
to as the ”molecular unit of currency” of intracellular energy transfer.
Moreover with the chemical formula C10H16N5O13P3, ATP is composed of 47 atoms,
showing that already within the NISQ era, one can simulate large and relevant
molecules from the standpoint of NMR.

The demonstrated simulations required around $2$ hours of quantum computing
run time according to the IBM Quantum Platform web interface.

### Verification with NMR experimental results.

The practical application considered here offers the advantage of being
efficiently verifiable since one can easily conduct an NMR experiment to
benchmark the results from the quantum computer. For this work, we have been
able to compare the spectra calculated on the QPU with real experimental
spectra, thus validating this advantage towards quantum utility. The
experimental spectra for ATP and 4-methylvaleric acid (fig. 5 c and f) were
obtained on a Brucker-Avance DPX FT-NMR $\quantity{400}{}$ instrument. The
solvents used are D2O for ATP and CDCl3 containing $0.03\%$ Tetramethylsilane
(resulting in a peak at $\quantity{7.26}{ppm}$) for 4-Methylvaleric acid.

The simulated spectra show a correspondence to the experimental ground truth
spectra when it comes to the general location of the peaks, artifacts often
cover the finer details such as peak splittings due to the product formula
approximation. In some cases, the splitting structure corresponds to the one
in the experimental, as can be seen in the inset of figure 5 f. Smaller peaks
cannot be faithfully represented since they fall below the noise level.

##  IV Conclusion

We have demonstrated that the nuclear spin dynamics problem can benefit from
the NISQ-based approaches. Our investigation into the application of NISQ
devices for simulating high-field regime NMR experiments marks a significant
step towards realizing quantum advantage for practical problems. The three
properties of a quantum utility experiment Aaronson (2023) were demonstrated:
its in-principle quantum advantage, the possibility to implement it with NISQ
resources, and crucially its efficient verifiability by benchmarking against
NMR experiments. By employing a product formula with minimal resource overhead
to map NMR interactions onto a QPU, we have successfully simulated liquid-
state proton NMR spectra for relevant molecules with up to $11$ spins and a
total of $47$ atoms. Through comparisons with real NMR experiments, despite
current limitations, our findings lay the groundwork for achieving quantum
utility, where computationally relevant problems can be solved by quantum
computers but not by conventional means. Additionally, we provided an
experimental estimation of the quantum resources required for solving larger
instances of the problem using our approach. The demonstrated polynomial
scaling on real processors represents a foundational step in advancing
practical quantum computation towards reality. We believe that in the next few
years new QPUs, combined with more advanced circuit optimization and error
mitigation methods will achieve useful classically intractable simulations of
NMR FID or other experiments.

##  V Data availability

The datasets generated and analyzed during this study are available at Zenodo.

##  VI Acknowledgment

We acknowledge the use of IBM Quantum services for this work. The views
expressed are those of the authors and do not reflect the official policy or
position of IBM or the IBM Quantum team. We acknowledge Christelle Jablonski
for NMR support, and Flavio Rump, Ilya Kuprov, and Julien Baglio for useful
discussions. We acknowledge support from the National Centre of Competence in
Research (NCCR) SPIN, funded by the Swiss National Science Foundation (grant
number 51NF40-180604), and from uptownBasel. We thank QuantumBasel for the
access to the IBM Quantum services and support. This research used the
computational cluster resource provided by FHNW HLS.

## References

  * Kim  _et al._ (2023) Y. Kim, A. Eddins, S. Anand, K. X. Wei, E. van den Berg, S. Rosenblatt, H. Nayfeh, Y. Wu, M. Zaletel, K. Temme, and A. Kandala, Nature 618, 500 (2023).
  * Begušić  _et al._ (2024) T. Begušić, J. Gray, and G. K.-L. Chan, Sci Adv 10, eadk4321 (2024).
  * Sels  _et al._ (2020) D. Sels, H. Dashti, S. Mora, O. Demler, and E. Demler, Nature Machine Intelligence 2, 396 (2020).
  * Seetharam  _et al._ (2023) K. Seetharam, D. Biswas, C. Noel, A. Risinger, D. Zhu, O. Katz, S. Chattopadhyay, M. Cetina, C. Monroe, E. Demler, and D. Sels, Science Advances 9, eadh2594 (2023), https://www.science.org/doi/pdf/10.1126/sciadv.adh2594 .
  * O’Brien  _et al._ (2022) T. E. O’Brien, L. B. Ioffe, Y. Su, D. Fushman, H. Neven, R. Babbush, and V. Smelyanskiy, PRX Quantum 3, 030345 (2022).
  * Kitaev  _et al._ (2002) A. Kitaev, A. Shen, and M. Vyalyi,  _Classical and Quantum Computation_, Graduate studies in mathematics (American Mathematical Society, 2002).
  * Vandersypen and Chuang (2005) L. M. K. Vandersypen and I. L. Chuang, Rev. Mod. Phys. 76, 1037 (2005).
  * of Japan (2017) T. N. M. R. S. of Japan,  _Experimental Approaches of NMR Spectroscopy: Methodology and Application to Life Science and Materials Science_ (Springer Nature Singapore, 2017).
  * Rei (2021) Nature Reviews Methods Primers 1, 4 (2021).
  * Childs  _et al._ (2018) A. M. Childs, D. Maslov, Y. Nam, N. J. Ross, and Y. Su, Proceedings of the National Academy of Sciences 115, 9456 (2018), https://www.pnas.org/doi/pdf/10.1073/pnas.1801723115 .
  * Manin (1980) Y. I. Manin,  _Computable and Uncomputable_ (Soviet Radio, 1980).
  * Feynman (1982) R. P. Feynman, International Journal of Theoretical Physics 21, 467 (1982).
  * Ledbetter  _et al._ (2011) M. P. Ledbetter, T. Theis, J. W. Blanchard, H. Ring, P. Ganssle, S. Appelt, B. Blümich, A. Pines, and D. Budker, Phys. Rev. Lett. 107, 107601 (2011).
  * Appelt  _et al._ (2010) S. Appelt, F. W. Häsing, U. Sieling, A. Gordji-Nejad, S. Glöggler, and B. Blümich, Phys. Rev. A 81, 023420 (2010).
  * Kuprov (2023) I. Kuprov,  _Spin: From Basic Symmetries to Quantum Optimal Control_ (Springer International Publishing, 2023).
  * Karabanov  _et al._ (2011) A. Karabanov, I. Kuprov, G. T. P. Charnock, A. van der Drift, L. J. Edwards, and W. Köckenberger, The Journal of Chemical Physics 135, 084106 (2011), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.3624564/15439618/084106_1_online.pdf .
  * Parker and Lado (1974) G. Parker and F. Lado, Phys. Rev. B 9 (1974), 10.1103/PhysRevB.9.22.
  * Exarchou  _et al._ (2002) V. Exarchou, A. Troganis, I. P. Gerothanassis, M. Tsimidou, and D. Boskou, Tetrahedron 58, 7423 (2002).
  * Levitt (2013) M. Levitt,  _Spin Dynamics: Basics of Nuclear Magnetic Resonance_ (Wiley, 2013).
  * Dashti  _et al._ (2018) H. Dashti, J. R. Wedell, W. M. Westler, M. Tonelli, D. Aceti, G. K. Amarasinghe, J. L. Markley, and H. R. Eghbalnia, Analytical Chemistry 90, 10646 (2018).
  * Mansuroglu  _et al._ (2023) R. Mansuroglu, F. Fischer, and M. J. Hartmann, Phys. Rev. Res. 5, 043035 (2023).
  * Mc Keever and Lubasch (2023) C. Mc Keever and M. Lubasch, Phys. Rev. Res. 5, 023146 (2023).
  * Trenkwalder  _et al._ (2023) L. M. Trenkwalder, E. Scerri, T. E. O’Brien, and V. Dunjko, “Compilation of product-formula hamiltonian simulation via reinforcement learning,” (2023), arXiv:2311.04285 [quant-ph] .
  * Qiskit contributors (2023) Qiskit contributors, “Qiskit: An open-source framework for quantum computing,”  (2023).
  * Hogben  _et al._ (2011) H. Hogben, M. Krzystyniak, G. Charnock, P. Hore, and I. Kuprov, Journal of Magnetic Resonance 208, 179 (2011).
  * Aaronson (2023) S. Aaronson,  _The Future of Quantum Supremacy Experiments_ (Q2B23 SV, 2023).


# Electric-field gradients at the nuclei from all-electron, four-component
relativistic density-functional theory using Gaussian-type orbitals

Marc Joosten<EMAIL_ADDRESS>Hylleraas Center for Quantum Molecular
Sciences, Department of Chemistry, University of Tromsø—The Arctic University
of Norway, 9037 Tromsø, Norway  Michal Repisky<EMAIL_ADDRESS>Hylleraas Center for Quantum Molecular Sciences, Department of Chemistry,
University of Tromsø—The Arctic University of Norway, 9037 Tromsø, Norway
Department of Physical and Theoretical Chemistry, Faculty of Natural Sciences,
Comenius University, Ilkovicova 6, SK-84215 Bratislava, Slovakia  Marius Kadek
<EMAIL_ADDRESS>Hylleraas Center for Quantum Molecular Sciences,
Department of Chemistry, University of Tromsø—The Arctic University of Norway,
9037 Tromsø, Norway  Pekka Pyykkö<EMAIL_ADDRESS>Department of
Chemistry, P.O. Box 55, FIN-00014, University of Helsinki, Finland  Kenneth
Ruud<EMAIL_ADDRESS>Hylleraas Center for Quantum Molecular Sciences,
Department of Chemistry, University of Tromsø—The Arctic University of Norway,
9037 Tromsø, Norway  Norwegian Defence Research Establishment, Instituttveien
20, 2027 Kjeller, Norway

(August 28, 2024)

###### Abstract

We present an all-electron, four-component relativistic implementation of
electric field gradients (EFGs) at the nuclei using Gaussian-type orbitals and
periodic boundary conditions. This allows us to include relativistic effects
variationally, which is important for compounds containing heavy elements and
for a property dependent the electronic structure close to the nuclei. The
all-electron approach ensures an accurate treatment of both core and valence
orbitals, as both are important in the evaluation of EFGs. Computational
efficiency is achieved through the use of a recent implementation of density
fitting in combination with quaternion algebra and restricted kinetic balance.
We use the relativistic approach to calculate the EFGs in different arsenic,
antimony and bismuth halides and oxyhalides, and explore the importance of
relativistic effects on EFGs in solids and compare these with results obtained
for molecular species. Our calculations contribute to establishing a reliable
estimate for the nuclear quadrupole moment of ${}^{209}\mathrm{Bi}$, for which
our best estimate is $-428(17)$ mb, in excellent agreement both with molecular
data and a recent reevaluation of the nuclear quadrupole moment obtained from
atomic data and ab initio calculations. Our results suggest that there is a
need to revisit the experimental data for the EFGs of several bismuth
oxyhalides.

##  I Introduction

Nuclear electric quadrupole moments (NQMs) are important parameters providing
fundamental information about nuclear structure, in particular whether the
nuclear structure is spherically symmetric or not. A prerequisite for the
existence of a non-vanishing NQM is a nuclear spin larger than $1/2$. The
nuclear quadrupole moments interact with the electric field gradients (EFGs)
at the nuclei, leading to a splitting of atomic or molecular energy levels
that can be observed for instance in nuclear quadrupole resonance spectroscopy
[1] or rotational spectroscopy [2, 3] and correspond to excitation energies in
the microwave region. Through measurements of the coupling of NQMs and EFGs,
one can also obtain information about the electron density close to the
nuclei, provided that accurate values of the NQMs are available. Much effort
has therefore been devoted to determine NQMs [4].

One way of determining NQMs is to combine accurate ab initio calculations of
EFGs with highly accurate transition energy data obtained from different
experiments, and derive the NQMs using the relation between the EFGs and
transition energies. A number of NQMs have been determined with very high
accuracy using this hybrid experimental/computational approach. For further
information on various methods of determining NQMs, their expected accuracy,
and the most up-to-date list of NQMs, we recommend the review by Pyykkö [4].

EFGs probe the electron density close to the nuclei, making relativistic
effects crucial, particularly for heavier nuclei [5]. The calculation of EFGs
is relatively straightforward, since the EFGs are first-order molecular
properties. They can be evaluated either as expectation values for variational
wave functions or by employing Z-matrix or Lagrangian techniques for highly
correlated, non-variational wave functions [6, 7], such as coupled-cluster or
configuration interaction methods. The incorporation of relativistic effects
is also fairly straightforward, even at the fully relativistic level of theory
using highly correlated wave functions [8, 9]. However, when relativistic
effects are approximated using a transformed Hamiltonian, significant picture-
change effects must be accounted for [10]. We also note that special care must
be taken if a finite-field approach is applied as an alternative way of
calculating EFGs. [11]

Compounds containing heavy elements often exist as solids. In the field of
computational materials modeling, a plane-wave description of the valence
electrons is often combined with the pseudopotential technique where the
oscillatory core electronic orbitals are effectively substituted with pseudo-
wave-functions that mimic the effect of the core electrons on the valence
electrons without the intricate oscillatory pattern. The use of
pseudopotentials allows relativistic effects on the core electrons to easily
be included, but the approximate description of the core electron density
makes them unsuited for evaluating properties such as the EFGs. The projector
augmented wave method allows for calculations of EFGs at the nuclei [12, 13],
but reaching accurate results can require that the EFGs are corrected using
molecular calculations at a more advanced theory level (_e.g._ with a hybrid
exchange-correlation functional) [14]. Alternative relativistic approaches
include the use of four-component numerical [15, 16, 17] or two-component
Slater-type atomic orbitals [5, 18] or techniques that combine plane waves
with auxiliary functions for explicit treatment of core electrons, while
potentially also including relativistic effects [19, 12, 4]. Nevertheless,
such hybrid approaches are inherently implementationally complex and
computationally intensive, necessitating meticulous convergence settings and a
profound understanding of the methods to ensure the accuracy of the calculated
results in complex systems.

We have adopted an alternative, efficient all-electron approach that
incorporates relativistic effects variationally at the four-component level of
theory and which utilizes Gaussian-type basis sets with periodic boundary
conditions, enabling systematic improvements in the results as the basis set
is expanded. Computational efficiency is ensured by the use of quaternion
algebra and restricted kinetic balance [20, 21], as well as the resolution-of-
the-identity approximation [22, 23] for the Coulomb integrals in combination
with multipole expansion of the far-field contributions [24]. Our previous
work has shown that our method reliably achieves accurate band gaps already
with double-zeta quality basis sets, and that smooth convergence of the
density is assured without issues even when diffuse functions are incorporated
into the basis set [21]. Indeed, we have established that retaining diffuse
functions is essential for rapid basis set convergence; conversely, the more
commonly employed “reduced” basis sets (where the diffuse functions are
removed) [25, 26, 27] yield inferior results and slower basis set convergence.

We note here that the use of Gaussian-type orbitals (GTOs) for computational
studies of periodic solids has an extensive history at the non-relativistic
level of theory, pioneered by the work of Pisani and Dovesi [28, 29], followed
by several more recent implementations [30, 31, 32, 33, 34, 35]. Some of these
GTO-based codes also reported calculations of EFGs for solid-state systems
[36, 37, 38, 39, 40], albeit without the full consideration of relativistic
effects.

Here, we extend our relativistic all-electron approach based on GTOs to the
calculation of EFGs in solids and use the first-principles calculations in
combination with experimental data to determine the NQM of As, Sb, and Bi from
the EFGs in their halides and oxyhalides. Furthermore, we investigate the
discrepancy between the atomic and molecular values of the NQM of 209Bi by
using solid-state Nuclear Quadrupole Resonance (NQR) data.

The remainder of the paper is organized as follows: In Section II, we provide
the basic theory for the calculation of EFGs in solids at the four-component
level of theory using Gaussian-type orbitals. In Section III, we give the
computational details for the calculations presented in this work. In Section
IV, we explore the importance of relativistic effects and choice of basis sets
on the EFGs in a series of pnicogen halides and oxyhalides, with a particular
focus on bismuth-containing systems. We obtain an estimate of the NQM of 209Bi
and compare it with atomic and molecular data. Finally, in Section V, we give
some concluding remarks and an outlook.

##  II Theory

NQR utilizes the interaction of the nuclear charge density with the
surrounding electric field. The quadrupole interaction is the leading
contribution, _i.e._ the interaction between the nuclear quadrupole moment and
the electric field gradient tensor. A detailed discussion can be found in Ref.
[41]. This leads to the NQR Hamiltonian in the eigensystem of the EFG, given
by

$H_{NQR}=\frac{Q}{2I(2I-1)}\left(V^{xx}I^{x}+V^{yy}I^{y}+V^{zz}I^{z}\right)=\frac{qQ}{4I(2I-1)}\left(3I^{z}-I^{2}+\frac{\eta}{2}\left(I^{2}_{+}+I^{2}_{-}\right)\right).$
(II.1)

We have here introduced the eigenvalues of the electric field gradient
components defined according to $|V^{zz}|>|V^{yy}|>|V^{xx}|$, the asymmetry
parameter $\eta=(V^{xx}-V^{yy})/V^{zz}$ and the nuclear spin $I$. $Q$ denotes
the nuclear quadrupole moment (NQM) and $q$ is the largest eigenvalue $V^{zz}$
and is often referred to as just the EFG. Diagonalization then leads to
different spin energy levels, with the magnitude depending on $q$ and $Q$,
while the general structure depends on $\eta$ and the total spin of the
nucleus $I$. In the following, only the traceless part of the EFG is
calculated; the traced part leads to a constant shift in the spectra [42] that
is not considered in this work.

The EFG tensor $V^{ij}$ is defined as the second derivative of the electric
potential at the position of the reference nucleus $K$ (using the notation
$\partial^{i}\equiv\partial/\partial r^{i}$):

$V_{K}^{ij}=\left.\partial^{i}\partial^{j}\int
d\mathbb{r^{\prime}}\rho(\mathbb{r^{\prime}})\frac{1}{\lvert\mathbb{r^{\prime}}-\mathbb{r}\rvert}\right|_{\mathbb{r}=\mathbb{r}_{K}}=\int
d\mathbb{r^{\prime}}\rho(\mathbb{r^{\prime}})v^{ij}\left(\mathbb{r^{\prime}}-\mathbb{r}_{K}\right),$
(II.2)

with

$v^{ij}(\mathbb{r})=\frac{3r^{i}r^{j}-\delta^{ij}r^{2}}{r^{5}},$ (II.3)

and the total charge density of the infinite periodic system

$\rho(\mathbb{r})=\sum_{\mathbb{m}}\rho_{\mathbb{m}}(\mathbb{r})=\sum_{\mathbb{m}}\rho_{\mathbb{0}}\left(\mathbb{r}-\mathbb{r}_{\mathbb{m}}\right),$
(II.4)

where $\mathbb{m}$ labels the unit cell centered at $\mathbb{r}_{\mathbb{m}}$.
In applying the derivative in Eq. (II.2), the Poisson term is neglected. This
has no effect on the traceless part of $V$ [42]. Equation (II.2) can be
evaluated as stated by keeping the nucleus $K$ fixed and summing over the
charge densities of every unit cell. Alternatively, the summation can be
shifted to the nucleus by using translation symmetry

$V_{K}^{ij}=\int
d\mathbb{r^{\prime}}\rho_{\mathbb{0}}(\mathbb{r^{\prime}})\sum_{\mathbb{m}}v^{ij}\left(\mathbb{r^{\prime}}-\mathbb{r}_{K}+\mathbb{r}_{\mathbb{m}}\right)=\sum_{\mathbb{m}}\int
d\mathbb{r^{\prime}}\rho_{\mathbb{0}}(\mathbb{r^{\prime}})v_{\mathbb{m},K}^{ij}\left(\mathbb{r^{\prime}}\right).$
(II.5)

The EFG consists of contributions from the electronic charge density and from
all other non-reference nuclei.
Expressing the electron density using the one-electron reduced density matrix
$D$ and Gaussian-type basis functions $\chi_{\mu\mathbb{n}}$ centered on each
unit cell $\mathbb{n}$, the electronic part of the EFG tensor is obtained as

$V_{e,K}^{ij}=-\sum_{\mathbb{m}}\sum_{\mu\nu\mathbb{n}}\int
d\mathbb{r}\text{Tr}\left[\chi^{\dagger}_{\mu\mathbb{0}}(\mathbb{r})v_{\mathbb{m},K}^{ij}\left(\mathbb{r}\right)\chi_{\nu\mathbb{n}}(\mathbb{r})D^{\nu\mathbb{n},\mu\mathbb{0}}\right]=-\sum_{\mathbb{m}u}\text{Tr}\left[I_{\mathbb{m},K,u}^{ij}D^{\bar{u}}\right],$
(II.6)

where $u\equiv(\mu\nu\mathbb{n})$, $\bar{u}\equiv(\nu\mathbb{n}\mu)$, and the
EFG integrals are labeled as $I$. Note here that the basis functions
$\chi_{\nu\mathbb{n}}$ are four-component Dirac bispinors, and the trace Tr
needs to be introduced to sum over the four components, _i.e._ four diagonal
matrix elements in the product of $I$ and $D$. More details about the
implementation of the calculation of the relativistic electronic structure of
solids in the ReSpect program is found in Ref. [21]. Similarly, the nuclear
part is given by

$V_{n,K}^{ij}=\sum_{\mathbb{m}}\sum_{A}Z_{A}\frac{3\left(\mathbb{r}_{A}-\mathbb{r}_{\mathbb{m},K}\right)^{i}\left(\mathbb{r}_{A}-\mathbb{r}_{\mathbb{m},K}\right)^{j}-\left|\mathbb{r}_{A}-\mathbb{r}_{\mathbb{m},K}\right|^{2}\delta^{ij}}{\left|\mathbb{r}_{A}-\mathbb{r}_{\mathbb{m},K}\right|^{5}}.$
(II.7)

In the summations over all atoms $A$ and unit cells $\mathbb{m}$, the
reference nucleus term $A=K$ must be excluded for $\mathbb{m}=\mathbb{0}$.

For sufficiently large values of $\mathbb{m}$, the charge density of the unit
cell $\mathbb{m}=\mathbb{0}$ has no significant overlap with the reference
nuclei in the cell $\mathbb{m}$ 111Due to the use of GTOs, the charge density
decays exponentially.. Therefore the lattice sum in Eq. (II.5) can be split
into near and far-field (NF and FF) parts. The NF is calculated directly using
Eqs. (II.6) and (II.7). The FF terms can be expanded using a multipole
expansion. This is based on

$\partial^{i}\partial^{j}\frac{1}{\left|\mathbb{r^{\prime}}-\mathbb{r}_{K}-\mathbb{r}_{\mathbb{m}}\right|}=\partial^{i}\partial^{j}\frac{1}{\left|\left(\mathbb{r^{\prime}}-\mathbb{r}_{C}\right)-\left(\mathbb{r}_{K}+\mathbb{r}_{\mathbb{m}}-\mathbb{r}_{C}\right)\right|}\\\
=\sum_{lm}\partial^{i}\partial^{j}R_{lm}\left(\mathbb{r^{\prime}}-\mathbb{r}_{C}\right)I^{*}_{lm}\left(\mathbb{r}_{K}+\mathbb{r}_{\mathbb{m}}-\mathbb{r}_{C}\right)=\sum_{lm}R_{lm}\left(\mathbb{r^{\prime}}-\mathbb{r}_{C}\right)\partial^{i}\partial^{j}I^{*}_{lm}\left(\mathbb{r}_{K}+\mathbb{r}_{\mathbb{m}}-\mathbb{r}_{C}\right).$
(II.8)

The center of the unit cell is chosen as the center of the multipole expansion
$\mathbb{r}_{C}$ and $R$ and $I$ refer to the scaled regular and irregular
solid harmonics as defined in Ref. [44]. The derivative in (II.8) acting on
$\mathbb{r^{\prime}}$ is equivalent to a derivative on
$\mathbb{r}_{\mathbb{m}}$ up to a minus sign. This can be used to shift the
derivatives between $R$ and $I$. Recursive expressions for the derivatives are
defined in Ref. [45], albeit using a different sign convention. The solid
harmonics defined above are complex functions. However, in practical
calculations, it is more convenient to use real harmonics [44]. Thus, the FF
contribution to the EFG becomes

$V^{ij}_{K,FF}=\sum_{lm}2\frac{1}{2^{\delta_{m0}}}Q_{lm}\left(\mathbb{r}_{C}\right)\partial^{i}\partial^{j}\sum_{\mathbb{m}\in
FF}I_{lm}\left(\mathbb{r}_{K}+\mathbb{r}_{\mathbb{m}}-\mathbb{r}_{C}\right)=\sum_{lm}Q_{lm}\left(\mathbb{r}_{C}\right)I^{ij}_{K,FF,lm}\left(\mathbb{r}_{C}\right),$
(II.9)

with the total multipole moments of the unit cell defined as

$Q_{lm}\left(\mathbb{r}_{C}\right)=\int
d\mathbb{r}\rho_{\mathbb{0}}(\mathbb{r})R_{lm}\left(\mathbb{r}-\mathbb{r}_{C}\right).$
(II.10)

Eq. (II.9) contains both electronic and nuclear contributions and is evaluated
in a similar manner as the NF terms in Eqs. (II.6) and (II.7). The additional
factor compensates for the switch to real harmonics and is absorbed into $I$.
The multipole moments are independent both of the lattice vector $\mathbb{m}$
and the specific nucleus $K$ and only need to be calculated once for a given
system. To evaluate the lattice sum, the renormalization scheme from Ref. [46]
is used, with additional details provided in Ref. [21]. This scheme only needs
to be employed once per system as it can be shifted to the position of a
different nucleus $K^{\prime}$ using the translation properties of $I$. Hence,
the computational cost of evaluating the FF terms is negligible, and the
entire computational effort comes from the direct calculation of the NF
contributions in Eq. (II.6).

The calculation of EFGs for molecules is significantly simpler. Eqs. (II.6)
and (II.7) can then be used directly by removing the lattice indices and the
summation over $\mathbb{m}$. The NF-FF split is in this case not needed. Our
implementation in ReSpect agrees for four-component molecular systems with
DIRAC[47] and for one-component solids with CRYSTAL[48].

##  III Computational details

We use the four-component Dirac–Coulomb Kohn–Sham level of theory together
with the generalized gradient approximation (GGA) for the nonrelativistic
exchange–correlation functional (PBE) [49]. The fully uncontracted
triple-$\zeta$ (TZ) Dyall basis sets are used unless otherwise noted [50, 51].
To accelerate the calculation of the Coulomb electrostatic interactions, we
employed the resolution of the identity approach [22, 23, 24]. A detailed
discussion of the implementation and accuracy of the RI-J approximation used
in this work will be provided elsewhere [24], but calculations on BiOCl lead
to errors of the order of 10-3 V/m2 in the EFG both at the non-relativistic
and relativistic levels of theory, a relative error that is less than 0.01%.

The momentum space was sampled with a $9\times 9\times 9$ mesh of $\mathbb{k}$
points for the bismuth oxyhalides and a $5\times 5\times 5$ mesh for the other
systems. The convergence with respect to the $\mathbb{k}$-point mesh was
assessed by comparing the EFG results obtained from calculations across a
progressively increasing sequence of meshes.

##  IV Results and discussion

In order to demonstrate the importance of relativistic effects for EFGs, we
will have a particular focus on 209Bi, as there until recently was a rather
significant discrepancy between the NQMs established based on comparing
experimental and computed data between molecular [52] and atomic studies [53].
This discrepancy was however recently resolved [54, 55], bringing the atomic
NQMs into very good agreement with the molecular NQMs. We also note that
nuclear calculations by Karayonchev et al. support the molecular NQM, with a
NQM of $-428$ mbarn [56]. It would be of interest to see whether our solid-
state implementation would confirm these predictions.

In order to better assess the accuracy of our calculated results, we have
performed benchmark calculations on molecular BiCl3 and solid BiOCl using
Dyall’s double (DZ), triple (TZ )and quadrupole-$\zeta$ (QZ) basis sets. [50,
51]. We have also included truncated versions of these basis sets for the
solid, in which we apply the commonly used approximation of removing all basis
functions with exponents $<0.1$ when calculating properties of periodic
systems using GTOs [25, 27, 26].

Table 1: Basis set convergence of Bi EFGs. Results for basis $X$ reported as EFG($X$)/EFG(QZ). For details, see text. System | tDz | DZ | tTz | TZ | tQZ | QZ
---|---|---|---|---|---|---
BiCl3 | - | 0.979 | - | 0.997 | - | 1
BiOCl | 0.879 | 0.938 | 0.984 | 0.986 | 1.0003 | 1

The value of the Bi EFG for the different basis sets are reported in Table 1
relative to the Bi EFG obtained using the QZ basis set. BiCl3 uses the
isolated molecular structure (vide infra), whereas for BiOCl we use the
experimental solid-state structure of Ref. [57].

For the molecular system, both the DZ and TZ results are very close to the QZ
value, being 97.9% and 99.7% of the QZ value, respectively. For the solid,
convergence is worse, with the DZ results clearly being of inferior quality,
with only 93.8% of the QZ value, whereas TZ is performing reasonably well with
98.6% of the QZ value. The Dyall basis sets are optimized for molecular
calculations, and a poorer basis set convergence in solids is not surprising.
The truncated DZ basis set performed poorly, whereas for the TZ basis, the
truncated result is only 0.2% off the full basis, and they are identical at
the QZ level.

We recently demonstrated that the removal of diffuse functions is not
necessary and inhibits basis-set convergence [21, 58] and prevents convergence
to the basis-set limit altogether for band gaps and structures. Even though
EFGs are a core-electron property, this is not the reason for the reduced
basis sets not affecting the quality of the calculated results. Instead, only
two diffuse basis functions are removed on Bi, whereas in our previous work,
as many as 6 basis functions were removed in the case of tungsten. This
highlights that truncation is highly dependent on the systems studied and the
basis sets used, and should in general be applied with care and ideally not at
all. In this work, all exponents of the original basis set are retained.

For BiOCl, we have also performed test calculations with the BLYP functional
[59, 60, 61]. The results are in good agreement with those obtained using the
PBE functional, being less than 3% larger than the PBE results both at the
nonrelativistic and relativistic levels of theory for the EFG at the Bi
nucleus using the TZ basis set.

Considering the possible uncertainty in the NQM of the 209Bi nucleus, we have
included a few arsenic and antimony halides and oxyhalides for which there
exist experimental data, in order to benchmark our results against reliable
experimental data. This also allows us to explore the relative importance of
relativistic corrections as we go down in the periodic table. All our results
are collected in Table 2. We report, in addition to the calculated EFG value
and the corresponding value for the asymmetry parameter $\eta$, also the
derived value for the NQM using the experimental resonance frequencies. The
reference values for the 75As, 121Sb, 123Sb and 209Bi NQMs are -311, -543,
-692 and -420 mbarn, respectively [62, 63, 64], and we also report the ratio
between the calculated and reference values for the NQM.

Table 2: Calculated EFGs [$10^{21}$ V/$m^{2}$], asymmetry parameters $\eta$ and resulting 209Bi quadrupole moments Q [mb]. Compound |  | Non-relativistic | Relativistic | Exp. | Ref.
---|---|---|---|---|---
| geom | EFG | $\eta$ | Q | Qcalc/Qref | EFG | $\eta$ | Q | Qcalc/Qref | $\eta_{\mathrm{exp}}$ |
AsCl3 | [65] | -19.11 | 0.0 | -341.7 | 1.10 | -18.47 | 0.0 | -353.5 | 1.14 | 0222Assumed to be zero. | [66]
AsBr3 | [67] | -14.47 | 0.076 | -349.7 | 1.12 | -13.65 | 0.082 | -370.7 | 1.19 | 0a | [66]
121SbCl3 | [68] | -26.87 | 0.151 | -590.6 | 1.09 | -25.60 | 0.170 | -619.9 | 1.14 | 0.184 | [69]
$\alpha$-121SbBr3 | [70] | -22.63 | 0.114 | -604.9 | 1.08 | -20.76 | 0.155 | -659.3 | 1.21 | 0.08 | [71]
$\beta$-121SbBr3333Exp data only available at 304K | [72] | -23.31 | 0.178 | -560.7 | 1.03 | -21.50 | 0.202 | -607.9 | 1.12 | 0.18 | [71]
123Sb4O5Cl2 | [73] | -32.81 | 0.261 | -805.0 | 1.16 | -33.37 | 0.302 | -791.5 | 1.14 | 0.307 | [57]
|  | -34.77 | 0.380 | -805.0 | 1.16 | -35.40 | 0.411 | -791.4 | 1.14 | 0.396 | [57]
BiCl3 | [74] | -31.95 | 0.362 | -421.3 | 1.00 | -27.36 | 0.565 | -492.1 | 1.17 | 0.583 | [75]
BiBr3 | [76] | -26.90 | 0.446 | -523.5 | 1.25 | -21.81 | 0.850 | -645.8 | 1.54 | 0.553 | [77]
BiBr3444Refitted experimental data, see text. | [76] |  |  | -387.1 | 0.92 |  |  | -477.6 | 1.14 | 0.840 |
BiOF | [57] | -29.90 | 0.000 | -410.1 | 0.98 | -26.70 | 0.000 | -459.3 | 1.09 | 0.0 | [57]
BiOCl | [57] | -15.39 | 0.000 | -424.8 | 1.01 | -10.62 | 0.001 | -615.8 | 1.47 | 0.0 | [57]
BiOBr | [57] | -12.56 | 0.001 | -410.5 | 0.98 | -8.22 | 0.002 | -627.1 | 1.49 | 0.0 | [57]
BiOI | [57] | -10.69 | 0.001 | -372.5 | 0.89 | -7.00 | 0.002 | -569.4 | 1.36 | 0.0 | [57]

For all As and Sb compounds, we observe good agreement with experimental data.
The differences are between 14 and 20% for the NQMs, which is in agreement
with the error for DFT calculations on EFGs in the literature [12]. The error
also includes temperature and other effects. For $\beta$-121SbBr3,
experimental data is only available at 300 K instead of at 77 K. If we assume
a similar temperature behavior as $\alpha$-121SbBr3, the ratio between
theoretical and experimental NQMs should increase to about 1.15. The agreement
for $\eta$ is generally good, although this quantity potentially benefits more
from error cancellation. The only exception is $\alpha$-121SbBr3, with a
measured value of 0.08 compared to the calculated value of 0.155.

Turning now to the bismuth compounds, the results for BiCl3 are in good
agreement with experiment. An overshoot of 17% is within expectations and
$\eta$ differs by only $0.02$. For BiBr3, neither the calculated EFG nor the
asymmetry parameter are in agreement with experiment. However, if we used a
different assignment of the experimental data, as discussed below, the
differences disappear and the agreement between theory and experiment is
comparable to that of BiCl3. Both BiCl3 and BiBr3, assuming the different
frequency assignment, support a value of -420 mb [64, 54] for the NQM of
209Bi.

The results for the bismuth oxyhalides are worse. The calculated EFGs decrease
for the heavier halides, following the same trend as observed for the other
systems. However, when combined with the experimental data, the resulting NQMs
for 209Bi are erratic. BiOF yields an expected error of 10%, but the other
three systems are much worse, without any discernible trend. Cl and Br have an
overshoot of almost 50%, while I is only off by 36%. These are beyond expected
variation both in relative accuracy as well as variation between the different
oxyhalides, and we are currently unable to resolve these discrepancies.
Whereas BiOF also supports a value of -420 mb for the NQM, it’s reliability is
questionable given the differences observed for the other three BiOXs.

As seen from Table 2, the agreement between the calculated and experimental
values for $\eta$ is good for all compounds, except for BiBr3, which shows a
large disagreement. This general agreement may be due to error cancellation
considering that $\eta$ involves a ratio between different EFG components.
Furthermore, across all other compounds the EFG decreases for heavier halides.
In contrast, the experimental data for BiBr3 suggest a small increase in the
EFG compared to BiCl3.

To investigate this conundrum, we have reoptimized the values for the nuclear
quadrupole coupling constant and $\eta$ from the experimentally measured
frequencies. For most systems, we obtain a better agreement with the measured
frequencies with a pure quadrupole spectra from this refitting process. While
the actual values of the EFGs do not change significantly, the refitting leads
to a reduction in the experimental uncertainties. For example, for BiCl3, the
agreement is below 0.1%. The only major change is for 121Sb in the
$\alpha$-structure, where the value changes from 343.95 MHz to 331.41 MHz,
which appears to be due to a mistake in one of the tables of Ref. [71].

The agreement is also excellent for BiBr3, however, only three out of the four
frequencies were observed. Therefore the frequency assignment is ambiguous,
and different choices are possible. We find a different assignment for
$\eta=0.840$. This value of $\eta$ agrees with our calculated value for
$\eta$, gives a good agreement for the NQM of 209Bi with the same overshoot as
seen for the arsenic and antimony compounds, and follows the trend of lowering
the EFG for heavier halides. However, this fit has a significantly higher
$\chi^{2}$ value of $4\cdot 10^{-1}$ compared to the experimental fit of
$6\cdot 10^{-4}$. For the other compounds, we obtain $\chi^{2}$ values ranging
from $2\cdot 10^{-2}$ to $4\cdot 10^{-6}$, making this fit an outlier. We are
at this point unable to provide an decisive conclusion on the origin of this
discrepancy.

Table 3: Results for EFGs [$10^{21}$ V/$m^{2}$] and asymmetry parameters $\eta$ in gas phase, for isolated molecular structures using solid-state geometries as well as for solids. For details, see text. Compound | EFG | $\eta$
---|---|---
| gas | ’molecule’ | solid | gas | ’molecule’ | solid
AsCl3 | -16.95 | -19.22 | -18.47 | 0 | 0.012 | 0.027
AsBr3 | -13.86 | -14.79 | -13.65 | 0 | 0.048 | 0.082
SbCl3 | -25.15 | -28.71 | -25.60 | 0 | 0.090 | 0.170
$\alpha$-SbBr3 | -21.36 | -24.99 | -20.76 | 0 | 0.10 | 0.155
$\beta$-SbBr3 | 555Only one configuration in gas phase. | -24.61 | -21.50 |  | 0.11 | 0.202
BiCl3 | -30.95 | -39.77 | -27.36 | 0 | 0.229 | 0.565
BiBr3 | -25.55 | -32.67 | -21.81 | 0 | 0.171 | 0.850

To estimate a value for the NQM of 209Bi based on the BiCl3 data, we apply a
scaling factor to the calculated value. The scaling factor is taken as the the
average ratios of calculated and experimental NQMs for AsCl3, AsBr3, 121SbCl3,
and the two sites in 123Sb4O5Cl2, as the temperature is the same for all these
studies, and good agreement between calculated and measured $\eta$ is
observed. We then obtain a NQM of -428(17) mbarn, with the number in
parenthesis representing the standard deviation, heavily influenced by the
result for AsBr3.

We have also implemented the calculation of EFGs in the molecular ReSpect
program. As we have a uniform level of theory for molecules and solids, and
many of the pnicogen halides are molecular solids, we have also explored the
effects of crystal packing as well as the geometry effects arising from
differences in bond lengths and bond angles going from the gas-phase structure
to the structure observed in the solid state. We have thus performed
calculations on individual molecules by taking the positions of one pnicogen
and its three nearest bonding halides from the respective solid-state
structures to study solid-state effects, as well as optimized the molecular
structure in the gas phase to explore the relaxation effects. All
computational details are the same, including XC functional and basis sets.
The calculated EFGs and asymmetry parameters are collected in Table 3.
There are significant effects both from changing the bond length going from
gas phase to the isolated molecule retaining the solid-state structure, as
well as due to neighboring atoms in the solid. Only the solid state
calculations give results for both EFGs and $\eta$ that are in agreement with
experiment, illustrating the need to use both a proper geometry as well as
include the solid-state effects. The increasing difference in EFGs from
isolated solid-state structure molecules to the solid phase for the heavier
compounds is not directly caused by relativity. This is rather due to the
closer crystal packing, leading to larger crystal-packing effects. Both the
pnicogen-pnicogen distances as well as the distance to the nearest non-bonding
halides are for these heavier congeners smaller. $\eta$ in the gas phase is
zero due to the trigonal pyramidal geometry, [78] whereas the value of $\eta$
for isolated solid-state structures and the solids show that it depends
strongly on the environment and is not only a results of the distorted
trigonal pyramid in solids.

##  V Concluding remarks

We have presented the theory for calculating EFGs in solids within the all-
electron, four-component framework using GTOs as implemented in the ReSpect
program. For As and Sb, these values are well established and we can use them
to determine the accuracy of EFGs calculated using our approach, and in
particular the inherent error of the DFT approximation applied in this work.
We have obtained good agreement for various systems, with errors in the range
of 14-20%, in line with other findings in the literature.

We have demonstrated that our solid-state approach for calculating electric
field gradients is reliable enough to determine the nuclear quadrupole moment
from observed nuclear quadrupole coupling constants. For Q(Bi), our current,
independent, result is -428(17) mb. It is comparable to the published ’World
average’ of -420(17) mb of Barzakh et al. [64] or the ’molecular’ values of
Teodoro and Haiduke (-420(8) mb) [52] or Dognon and Pyykkö (-422(3) mb) [54]
and others, as well as with recent nuclear calculations by Karayonchev et al.
(-428 mb) [56]. Note that the old atomic value of Bieron and Pyykkö (-516(15)
mb) [53] suffered from a too limited electron correlation space. Further
values and a new World average Q(Bi) should be proposed later.

However, the nuclear quadrupole moments derived from bismuth oxyhalides behave
erratically, even though the calculated EFGs appear reasonable. Although the
results for BiOF also supports the molecular value for the NQM of 209Bi, its
reliability is in question until the disagreements for the other oxyhalides
are resolved. We recommed that new experimental studies are performed for
bismuth oxyhalides.

##  VI Acknowledgment

This work was supported by the Research Council of Norway through its Centres
of Excellence scheme, project number 262695, its FRIPRO grant nr. 315822, and
the Mobility Grant scheme (project no. 301864). For the computational and data
storage resources, we acknowledge the support of Sigma2—the National
Infrastructure for High Performance Computing and Data Storage in Norway,
Grant Nos. NN4654K and NN14654K. In addition, MR acknowledges funding from the
European Union’s Horizon 2020 research and innovation program under the Marie
Skłodowska-Curie Grant Agreement No. 945478 (SASPRO2), and the Slovak Research
and Development Agency (Grant No. APVV-21-0497). PP is Professor Emeritus at
University of Helsinki.

## References

  * Smith [1971] J. A. S. Smith, Nuclear quadrupole resonance spectroscopy. General principles, J. Chem. Educ. 48, 39 (1971).
  * Puzzarini  _et al._ [2010] C. Puzzarini, J. F. Stanton, and J. Gauss, Quantum-chemical calculation of spectroscopic parameters for rotational spectroscopy, Int. Rev. Phys. Chem. 29, 273 (2010).
  * Puzzarini [2013] C. Puzzarini, Rotational spectroscopy meets theory, Phys. Chem. Chem. Phys. 15, 6595 (2013).
  * Pyykkö [2018] P. Pyykkö, Year-2017 nuclear quadrupole moments, Mol. Phys. 116, 1328 (2018).
  * van Lenthe and J. Baerends [2000] E. van Lenthe and E. J. Baerends, Density functional calculations of nuclear quadrupole coupling constants in the zero-order regular approximation for relativistic effects, J. Chem. Phys. 112, 8279 (2000).
  * Handy and Schaefer III [1984] N. C. Handy and H. F. Schaefer III, On the evaluation of analytic energy derivatives for correlated wave functions, J.Chem.Phys. 81, 5031 (1984).
  * Helgaker and Jørgensen [1979] T. Helgaker and P. Jørgensen, Configuration-interaction energy derivatives in a fully variational formulation, Theor. Chim. Acta 75, 111 (1979).
  * Visscher  _et al._ [1998] L. Visscher, T. Enevoldsen, T. Saue, and J. Oddershede, Molecular relativistic calculations of the electric field gradients at the nuclei in the hydrogen halides, J. Chem. Phys. 109, 9677 (1998).
  * De Jong  _et al._ [1998] W. De Jong, L. Visscher, and W. Nieuwpoort, On the bonding and the electric field gradient of the uranyl ion, J. Mol. Struct.: THEOCHEM 458, 41 (1998).
  * Autschbach  _et al._ [2012] J. Autschbach, D. Peng, and M. Reiher, Two-Component Relativistic Calculations of Electric-Field Gradients Using Exact Decoupling Methods: Spin-orbit and Picture-Change Effects, J.Chem.Theory Comput. 8, 4239 (2012).
  * Pernpointner  _et al._ [1998] M. Pernpointner, M. Seth, and P. Schwerdtfeger, A point-charge model for the nuclear quadrupole moment: Coupled-cluster, Dirac–Fock, Douglas–Kroll, and nonrelativistic Hartree–Fock calculations for the Cu and F electric field gradients in CuF, J.Chem.Phys. 108, 6722 (1998).
  * Petrilli  _et al._ [1998] H. M. Petrilli, P. E. Blöchl, P. Blaha, and K. Schwarz, Electric-field-gradient calculations using the projector augmented wave method, Phys. Rev. B 57, 14690 (1998).
  * Profeta  _et al._ [2003] M. Profeta, F. Mauri, and C. J. Pickard, Accurate first principles prediction of 17O NMR parameters in SiO2: assignment of the zeolite ferrierite spectrum, J. Am. Chem. Soc. 125, 541 (2003).
  * Hartman  _et al._ [2021] J. D. Hartman, A. Mathews, and J. K. Harper, Fast and accurate electric field gradient calculations in molecular solids with density functional theory, Front. Chem. 9, 751711 (2021).
  * Jeglič  _et al._ [2009] P. Jeglič, J.-W. Bos, A. Zorko, M. Brunelli, K. Koch, H. Rosner, S. Margadonna, and D. Arčon, Influence of the Nd 4 f states on the magnetic behavior and the electric field gradient of the oxypnictides superconductors NdFeAsO1-xFx, Phys. Rev. B 79, 094515 (2009).
  * Koch  _et al._ [2009] K. Koch, R. Kuzian, K. Koepernik, I. Kondakova, and H. Rosner, Analysis of the electric field gradient in the perovskites SrTiO3 and BaTiO3: Density functional and model calculations, Phys. Rev. B 80, 125113 (2009).
  * Koch [2009] K. Koch,  _Crystal structure, electron density and chemical bonding in inorganic compounds studied by the electric field gradient_ , Ph.D. thesis, Technische Universität Dresden (2009).
  * Autschbach  _et al._ [2010] J. Autschbach, S. Zheng, and R. W. Schurko, Analysis of electric field gradient tensors at quadrupolar nuclei in common structural motifs, Concepts Magn. Reson. A 36, 84 (2010).
  * Blaha  _et al._ [1985] P. Blaha, K. Schwarz, and P. Herzig, First-principles calculation of the electric field gradient of Li3N, Phys. Rev. Lett. 54, 1192 (1985).
  * Repisky  _et al._ [2020] M. Repisky, S. Komorovsky, M. Kadek, L. Konecny, U. Ekström, E. Malkin, M. Kaupp, K. Ruud, O. L. Malkina, and V. G. Malkin, Respect: Relativistic spectroscopy DFT program package, J. Chem. Phys. 152, 184101 (2020).
  * Kadek  _et al._ [2019] M. Kadek, M. Repisky, and K. Ruud, All-electron fully relativistic Kohn-Sham theory for solids based on the Dirac-Coulomb Hamiltonian and Gaussian-type functions, Physical Review B 99, 205103 (2019).
  * Eichkorn  _et al._ [1995] K. Eichkorn, O. Treutler, H. Öhm, M. Häser, and R. Ahlrichs, Auxiliary basis sets to approximate Coulomb potentials, Chem.Phys.Lett. 240, 283 (1995).
  * Eichkorn  _et al._ [1997] K. Eichkorn, F. Weigend, O. Treutler, and R. Ahlrichs, Auxiliary basis sets for main row atoms and transition metals and their use to approximate Coulomb potentials, Theor.Chem.Acc. 97, 119 (1997).
  * [24] M. Joosten, M. Kadek, K. Ruud, and M. Repisky, Manuscript in preparation.
  * Peralta  _et al._ [2005] J. E. Peralta, J. Uddin, and G. E. Scuseria, Scalar relativistic all-electron density functional calculations on periodic systems, J. Chem. Phys. 122, 084108 (2005).
  * Perry  _et al._ [2001] J. K. Perry, J. Tahir-Kheli, and W. A. Goddard III, Antiferromagnetic band structure of La2 CuO4: Becke-3–Lee-Yang-Parr calculations, Phys. Rev. B 63, 144510 (2001).
  * Usvyat  _et al._ [2011] D. Usvyat, B. Civalleri, L. Maschio, R. Dovesi, C. Pisani, and M. Schütz, Approaching the theoretical limit in periodic local MP2 calculations with atomic-orbital basis sets: the case of LiH, J. Chem. Phys. 134, 214105 (2011).
  * Pisani and Dovesi [1980] C. Pisani and R. Dovesi, Exact-exchange Hartree–Fock calculations for periodic systems. I. Illustration of the method, Int. J. Quantum Chem. 17, 501 (1980).
  * Dovesi  _et al._ [1983] R. Dovesi, C. Pisani, C. Roetti, and V. Saunders, Treatment of Coulomb interactions in Hartree-Fock calculations of periodic systems, Phys. Rev. B 28, 5781 (1983).
  * Lippert  _et al._ [1997] G. Lippert, J. Hutter, and M. Parrinello, A hybrid Gaussian and plane wave density functional scheme, Mol. Phys. 92, 477 (1997).
  * Boettger and Trickey [2000] J. Boettger and S. Trickey, Inclusion of relativistic effects in Gaussian-basis density functional calculations for extended systems, J. Mol. Struct. 501, 285 (2000).
  * Kudin and Scuseria [2000] K. N. Kudin and G. E. Scuseria, Linear-scaling density-functional theory with Gaussian orbitals and periodic boundary conditions: Efficient evaluation of energy and forces via the fast multipole method, Phys. Rev. B 61, 16440 (2000).
  * Łazarski  _et al._ [2015] R. Łazarski, A. M. Burow, and M. Sierka, Density functional theory for molecular and periodic systems using density fitting and continuous fast multipole methods, J. Chem. Theory Comput. 11, 3029 (2015).
  * Becker and Sierka [2019] M. Becker and M. Sierka, Density functional theory for molecular and periodic systems using density fitting and continuous fast multipole method: Stress tensor, J. Comput. Chem. 40, 2563 (2019).
  * Sun  _et al._ [2017] Q. Sun, T. C. Berkelbach, J. D. McClain, and G. K. Chan, Gaussian and plane-wave mixed density fitting for periodic systems, J. Chem. Phys. 147, 164119 (2017).
  * Palmer and Blair-Fish [1994] M. H. Palmer and J. A. Blair-Fish, Quadrupole Coupling Assignments in Inorganic Periodic Systems by ab initio Calculation of Electric Field Gradients, Z. Naturforsch. A 49, 137 (1994).
  * Stare  _et al._ [2020] J. Stare, A. Gradišek, and J. Seliger, Nuclear quadrupole resonance supported by periodic quantum calculations: A sensitive tool for precise structural characterization of short hydrogen bonds, Phys. Chem. Chem. Phys. 22, 27681 (2020).
  * Declerck  _et al._ [2006] R. Declerck, E. Pauwels, V. Van Speybroeck, and M. Waroquier, First-principles calculations of hyperfine parameters with the gaussian and augmented-plane-wave method: Application to radicals embedded in a crystalline environment, Phys. Rev. B 74, 245103 (2006).
  * Schmidt  _et al._ [2008] J. Schmidt, J. Hutter, H.-W. Spiess, and D. Sebastiani, Beyond isotropic tumbling models: nuclear spin relaxation in liquids from first principles, ChemPhysChem 9, 2313 (2008).
  * Sun  _et al._ [2020] Q. Sun, X. Zhang, S. Banerjee, P. Bao, M. Barbry, N. S. Blunt, N. A. Bogdanov, G. H. Booth, J. Chen, Z.-H. Cui, _et al._ , Recent developments in the PySCF program package, J. Chem. Phys. 153, 024119 (2020).
  * Seliger [2000] J. Seliger, Nuclear quadrupole resonance: Theory, in  _Encyclopedia of Spectroscopy and Spectrometry_ , edited by J. Lindon, G. Tranter, and J. Holmes (Academic Press, San Diego, 2000) pp. 1672–1680.
  * Thyssen  _et al._ [2001] J. Thyssen, P. Schwerdtfeger, M. Bender, W. Nazarewicz, and P. B. Semmes, Quadrupole and hexadecapole couplings for 127I in Li${}^{127}\mathrm{I}$, Phys. Rev. A 63, 022505 (2001).
  * Note [1] Due to the use of GTOs, the charge density decays exponentially.
  * Helgaker  _et al._ [2014] T. Helgaker, P. Jorgensen, and J. Olsen,  _Molecular Electronic-Structure Theory_ (Wiley, 2014).
  * Pérez-Jordá and Yang [1996] J. M. Pérez-Jordá and W. Yang, A concise redefinition of the solid spherical harmonics and its use in fast multipole methods, J. Chem. Phys. 104, 8003 (1996).
  * Kudin and Scuseria [2004] K. N. Kudin and G. E. Scuseria, Revisiting infinite lattice sums with the periodic fast multipole method, J. Chem. Phys. 121, 2886 (2004).
  * [47] DIRAC, a relativistic ab initio electronic structure program, Release DIRAC24 (2024), written by L. Visscher, H. J. Aa. Jensen, R. Bast, A. S. P. Gomes and T. Saue, with contributions from I. A. Aucar, V. Bakken, J. Brandejs, C. Chibueze, J. Creutzberg, K. G. Dyall, S. Dubillard, U. Ekström, E. Eliav, T. Enevoldsen, E. Faßhauer, T. Fleig, O. Fossgaard, L. Halbert, E. D. Hedegård, T. Helgaker, B. Helmich–Paris, J. Henriksson, M. van Horn, M. Iliaš, Ch. R. Jacob, S. Knecht, S. Komorovský, O. Kullie, J. K. Lærdahl, C. V. Larsen, Y. S. Lee, N. H. List, H. S. Nataraj, M. K. Nayak, P. Norman, A. Nyvang, G. Olejniczak, J. Olsen, J. M. H. Olsen, A. Papadopoulos, Y. C. Park, J. K. Pedersen, M. Pernpointner, J. V. Pototschnig, R. di Remigio, M. Repisky, K. Ruud, P. Sałek, B. Schimmelpfennig, B. Senjean, A. Shee, J. Sikkema, A. Sunaga, A. J. Thorvaldsen, J. Thyssen, J. van Stralen, M. L. Vidal, S. Villaume, O. Visser, T. Winther, S. Yamamoto and X. Yuan (available at https://doi.org/110.5281/zenodo.10680560, see also https://www.diracprogram.org).
  * Dovesi  _et al._ [2020] R. Dovesi, F. Pascale, B. Civalleri, K. Doll, N. M. Harrison, I. Bush, P. D’Arco, Y. Noël, M. Rérat, P. Carbonniere, _et al._ , The CRYSTAL code, 1976–2020 and beyond, a long story, J. Chem. Phys. 152, 204111 (2020).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Dyall [2002] K. G. Dyall, Relativistic and nonrelativistic finite nucleus optimized triple-zeta basis sets for the 4p, 5p and 6p elements, Theor. Chem. Acc. 108, 335 (2002).
  * Dyall [2006] K. G. Dyall, Relativistic quadruple-zeta and revised triple-zeta and double-zeta basis sets for the 4p, 5p, and 6p elements, Theor. Chem. Acc. 115, 441 (2006).
  * Teodoro and Haiduke [2013] T. Q. Teodoro and R. L. Haiduke, Nuclear electric quadrupole moment of bismuth from the molecular method, Phys.Rev.A 88, 052504 (2013).
  * Bieroń and Pyykkö [2001] J. Bieroń and P. Pyykkö, Nuclear Quadrupole Moments of Bismuth, Phys.Rev.Lett. 87, 133003 (2001).
  * Dognon and Pyykkö [2023] J.-P. Dognon and P. Pyykkö, Determining nuclear quadrupole moments of Bi and Sb from molecular data, Phys.Chem.Chem.Phys. 25, 2758 (2023).
  * Eliav  _et al._ [2024] E. Eliav, A. Borschevsky, A. Zaitsevskii, A. V. Oleynichenko, and U. Kaldor, Relativistic Fock-Space Coupled Cluster Method: Theory and Recent Applications, in  _Comprehensive Computational Chemistry_ , Vol. 3, edited by M. Yañes and R. J. Boyd (Elsevier, 2024) pp. 79–93.
  * Karayonchev  _et al._ [2019] V. Karayonchev, A. Blazhev, A. Esmaylzadeh, J. Jolie, M. Dannhoff, F. Diel, F. Dunkel, C. Fransen, L. M. Gerhard, R.-B. Gerst, L. Knafla, L. Kornwebel, C. Müller-Gatermann, J.-M. Régis, N. Warr, K. O. Zell, M. Stoyanova, and P. Van Isacker, Lifetimes in 211At and their implications for the nuclear structure above 208Pb, Phys. Rev. C 99, 024326 (2019).
  * Bastow and Whitfield [1977] T. Bastow and H. Whitfield, Nuclear quadrupole resonance of 209Bi in Bismuth oxyhalides, Journal of Magnetic Resonance (1969) 26, 461 (1977).
  * Kadek  _et al._ [2023] M. Kadek, B. Wang, M. Joosten, W.-C. Chiu, F. Mairesse, M. Repisky, K. Ruud, and A. Bansil, Band structures and $\mathbb{Z}_{2}$ invariants of two-dimensional transition metal dichalcogenide monolayers from fully relativistic Dirac-Kohn-Sham theory using Gaussian-type orbitals, Phys. Rev. Mater. 7, 064001 (2023).
  * Becke [1988] A. D. Becke, Density-functional exchange-energy approximation with correct asymptotic behavior, Phys. Rev. A 38, 3098 (1988).
  * Lee  _et al._ [1988] C. Lee, W. Yang, and R. G. Parr, Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density, Phys. Rev. B 37, 785 (1988).
  * Miehlich  _et al._ [1989] B. Miehlich, A. Savin, H. Stoll, and H. Preuss, Results obtained with the correlation energy density functionals of Becke and Lee, Yang and Parr, Chem.Phys.Lett. 157, 200 (1989).
  * Demovic  _et al._ [2010] L. Demovic, V. Kellö, and A. Sadlej, The quadrupole moment of the As nucleus from molecular microwave data and calculated relativistic electric field gradients, Chem.Phys.Lett. 498, 10 (2010).
  * Haiduke  _et al._ [2006] R. L. Haiduke, A. B. da Silva, and L. Visscher, The nuclear electric quadrupole moment of antimony from the molecular method, J.Chem.Phys. 125, 064301 (2006).
  * Barzakh  _et al._ [2021] A. Barzakh, A. N. Andreyev, C. Raison, J. G. Cubiss, P. Van Duppen, S. Péru, S. Hilaire, S. Goriely, B. Andel, S. Antalic, M. Al Monthery, J. C. Berengut, J. Bieroń, M. L. Bissell, A. Borschevsky, K. Chrysalidis, T. E. Cocolios, T. Day Goodacre, J.-P. Dognon, M. Elantkowska, E. Eliav, G. J. Farooq-Smith, D. V. Fedorov, V. N. Fedosseev, L. P. Gaffney, R. F. Garcia Ruiz, M. Godefroid, C. Granados, R. D. Harding, R. Heinke, M. Huyse, J. Karls, P. Larmonier, J. G. Li, K. M. Lynch, D. E. Maison, B. A. Marsh, P. Molkanov, P. Mosat, A. V. Oleynichenko, V. Panteleev, P. Pyykkö, M. L. Reitsma, K. Rezynkina, R. E. Rossel, S. Rothe, J. Ruczkowski, S. Schiffmann, C. Seiffert, M. D. Seliverstov, S. Sels, L. V. Skripnikov, M. Stryjczyk, D. Studer, M. Verlinde, S. Wilman, and A. V. Zaitsevskii, Large Shape Staggering in Neutron-Deficient Bi Isotopes, Phys. Rev. Lett. 127, 192501 (2021).
  * Galy  _et al._ [2002] J. Galy, R. Enjalbert, P. Lecante, and A. Burian, AsCl3: From the crystalline to the liquid state. XRD (176< T (K)< 250) and WAXS (295 k) studies, Inorg. Chem. 41, 693 (2002).
  * Barnes and Bray [1955a] R. Barnes and P. Bray, Nuclear Quadrupole Resonances of As75, J. Chem. Phys. 23, 407 (1955a).
  * Singh and Swaminathan [1967] A. Singh and S. Swaminathan, Refinement of the crystal structure of arsenic tribromide, Z. Kryst. 124, 375 (1967).
  * Lipka [1979] A. Lipka, An X-ray structure redetermination of antimony trichloride, Acta Crystallogr., Sect. B 35, 3020 (1979).
  * Barnes and Bray [1955b] R. Barnes and P. Bray, Nuclear Quadrupole Resonances of Sb121 and Sb123, J. Chem. Phys. 23, 1177 (1955b).
  * Cushen and Hulme [1964] D. W. Cushen and R. Hulme, 798. the crystal and molecular structure of antimony tribromide: $\alpha$-antimony tribromide, J. Chem. Soc. , 4162 (1964).
  * Ogawa [1958] S. Ogawa, Nuclear quadrupole resonances in arsenic and antimony trihalides and sulphur molecular addition compounds, J. Phys. Soc. Jpn. 13, 618 (1958).
  * Cushen and Hulme [1962] D. W. Cushen and R. Hulme, 427. the crystal and molecular structure of antimony tribromide: $\beta$-antimony tribromide, J. Chem. Soc. , 2218 (1962).
  * Paul  _et al._ [2022] S. Paul, B. Sen, N. Chakraborty, S. Das, S. Mondal, A. P. Chattopadhyay, and S. I. Ali, pH-regulated hydrothermal synthesis and characterization of Sb4O5X2 (X= Br/Cl) and its use for the dye degradation of methyl orange both with and without light illumination, RSC Adv 12, 8374 (2022).
  * Bartl [1982] H. Bartl, Meßzeitersparnis durch Profilauswertung bei der Registrierung von Neutroneneinkristallreflexen. Die Kristallstruktur von Wismuttrichlorid, BiCl3, Fresenius J Anal Chem 312, 17 (1982).
  * Robinson [1955] H. G. Robinson, Quadrupole Spectrum of BiCl3, Physical Review 100, 1731 (1955).
  * Benda [1980] H. v. Benda, Zur Polymorphie des Wismuttribromids, Z. Kryst. 151, 271 (1980).
  * Swiger  _et al._ [1965] E. D. Swiger, P. J. Green, G. L. McKown, and J. D. Graybeal, Nuclear Quadrupole Resonance of Bi209 in BiBr3, J. Phys. Chem. 69, 949 (1965).
  * Pyykkö and Elmi [2008] P. Pyykkö and F. Elmi, Deuteron quadrupole coupling in benzene: librational corrections using a temperature-dependent Einstein model, and summary. The symmetries of electric field gradients and conditions for $\eta$ = 1, Phys. Chem. Chem. Phys. 10, 3867 (2008).