# Ab initio Green-Kubo simulations of heat transport in solids: Method and
implementation

Florian Knoop  The NOMAD Laboratory at the FHI of the Max-Planck-Gesellschaft
and IRIS-Adlershof of the Humboldt-Universität zu Berlin, 14195 Berlin,
Germany  Theoretical Physics Division, Department of Physics, Chemistry and
Biology (IFM), Linköping University, SE-581 83 Linköping, Sweden  Matthias
Scheffler  The NOMAD Laboratory at the FHI of the Max-Planck-Gesellschaft and
IRIS-Adlershof of the Humboldt-Universität zu Berlin, 14195 Berlin, Germany
Christian Carbogno  The NOMAD Laboratory at the FHI of the Max-Planck-
Gesellschaft and IRIS-Adlershof of the Humboldt-Universität zu Berlin, 14195
Berlin, Germany

(September 2022)

###### Abstract

_Ab initio_ Green-Kubo (aiGK) simulations of heat transport in solids allow
for assessing lattice thermal conductivity in anharmonic or complex materials
from first principles. In this work, we present a detailed account of their
practical application and evaluation with an emphasis on noise reduction and
finite-size corrections in semiconductors and insulators. To account for such
corrections, we propose strategies in which all necessary numerical parameters
are chosen based on the dynamical properties displayed during molecular
dynamics simulations in order to minimize manual intervention. This paves the
way for applying the aiGK method in semi-automated and high-throughput
frameworks. The proposed strategies are presented and demonstrated for
computing the lattice thermal conductivity at room temperature in the mildly
anharmonic periclase MgO, and for the strongly anharmonic marshite CuI.

##  I Introduction

Heat transport is an important phenomenon in many branches of physics and
adjacent fields, be it materials science investigating technologically
relevant compounds [1, 2], or astrophysics and earth sciences, where
thermodynamic properties of planets are studied [3, 4]. In dielectric solids,
thermal transport is mostly determined by the conduction of heat energy in the
form of thermal nuclear motion (lattice thermal conductivity), and electronic
heat transport, photonic heat radiation, as well as convective contributions
due to mass transport can be neglected [5]. _Ab initio_ simulations of the
lattice thermal conductivity are typically performed in the framework of
phonon theory: The potential-energy surface (PES) is approximated by force
constants which can be obtained either as potential-energy derivatives or in a
renormalized, temperature-dependent fashion [6, 7, 8]. The equation of motion
is solved for the harmonic, second-order terms, which results in decoupled
phonon modes. Higher-order terms, up to fourth order [9], are included via
perturbation theory to compute phonon lifetimes. The Boltzmann transport
equation (BTE) [10, 11] then yields particle-like conduction contributions to
the thermal conductivity. Additional contributions from wave-like conduction
are accounted for in the Wigner transport formulation [12]. These
contributions are particularly important in complex crystals when the
individual phonon modes overlap.

However, thermal insulators with $\kappa<10$ W/mK of importance for
applications such as thermal barrier coatings in heat engines [2, 13] or
thermoelectric materials for waste-heat recovery [1], are often strongly
anharmonic [14], and the phonon picture underlying the Boltzmann or Wigner
transport equations likely breaks down, even when renormalized, temperature-
dependent force constants are used [15]. The Ioffe-Regel criterion [16] poses
a formal limit for the validity of the phonon picture and therefore
perturbative formalisms [12]. This intuitive criterion states that phonons are
only well-defined quasiparticles if their lifetimes exceed their oscillation
periods. Non-perturbative approaches based on Green-Kubo (GK) theory [17, 18,
19] do not suffer from this shortcoming since the heat flux is evaluated in
molecular dynamics simulations and all anharmonic effects are taken into
account. Accordingly, the GK method covers the transport mechanisms described
by the Boltzmann and Wigner transport equations, but also the regime beyond
the Ioffe-Regel criterion in which the phonon picture becomes invalid [12,
20]. Its extension to first principles frameworks, i. e., the _ab initio_
Green-Kubo (aiGK) technique, was introduced recently [21, 22]. By evaluating
the _ab initio_ heat flux along _ab initio_ molecular dynamics (aiMD)
trajectories, GK theory can be used to access the thermal conductivity in a
non-perturbative way on the basis of a fully _ab initio_ description of the
PES. Simplifying model assumptions about the PES such as the (quasi-)harmonic
approximation are therefore not needed. This makes aiGK a suitable tool for
the parameter-free study of materials of arbitrary anharmonic strength
whenever the nuclear dynamics can be described by aiMD simulations, i. e., at
temperatures where nuclear quantum effects can be neglected [23], and at which
the system is not close to a structural phase transition [24].

While aiGK offers an encompassing framework for first-principles heat
transport simulations, its practical implementation brings a set of challenges
that need to be addressed: The noise due to thermodynamic fluctuations in
small ensemble sizes of less than ten aiMD trajectories with tens of
picoseconds simulation time each, and finite-size effects when using
supercells of hundreds atoms only. These hurdles are of particular importance
in crystalline systems where finite-size effects can be significant because
readily accessible supercell sizes are not sufficient to describe all relevant
length scales [25, 26], as opposed to amorphous systems and liquids where
vibrations are more localized [27].

Several ideas to reduce noise [28, 29] and account for finite-size effects
[21] have been suggested in the literature. Their application to practical
calculations, however, requires educated choices for several numerical
parameters. Furthermore, no unified framework addressing both problems is
available to date, hindering broader application of aiGK methods for crystals.
In this work, we present and discuss such a framework and its numerical
implementation in FHI-vibes [30] for two test systems: Periclase magnesium
oxide (MgO), and marshite copper iodide (CuI). Both are simple binary, cubic
compounds, however, CuI is much more anharmonic [14], and BTE simulations
overestimate its thermal conductivity significantly [31].

We demonstrate the implementation in detail for the case of periclase MgO
which is well-known in the literature of first-principles heat transport
techniques [32, 33, 34]. We discuss the impact of noise-reduction and finite-
size-extrapolation techniques, and propose strategies to apply such
corrections without human intervention by choosing the required numerical
parameters using the available aiMD data: First, we present a new real-time
scheme to remove noise from the heat flux autocorrelation function (HFACF),
which enables to choose cutoff times in a numerically robust way based on a
“first dip” criterion, i. e., the time when the HFACF drops to zero for the
first time [35]. The two-step procedure comprises discarding non-contributing
terms from the flux based physical arguments [36], and subsequent noise
filtering that preserves the integrated thermal conductivity. Next, we discuss
a size extrapolation scheme for periodic solids, adapted from the one first
introduced in Ref. [21], which allows to correct for finite size effects of
simulation cells used in aiMD simulations. Finally, we discuss convergence in
the simulation times. This approach is then applied to the strongly anharmonic
CuI. Good agreement with the literature is obtained in both cases.

Both materials are studied at the level of the generalized gradient
approximation (GGA) using the PBEsol functional and light-default basis sets
in FHI-aims [37, 38]. Supercells are $3\times 3\times 3$ extensions of the
conventional, cubic unit cells, with 216 atoms each. The molecular dynamics
(MD) simulations are performed via FHI-vibes [30]. The aiGK method as
described here is implemented in FHI-vibes as well. Force constants for the
size extrapolation via harmonic mapping are obtained by regression from the MD
runs via the temperature dependent effective potentials code (TDEP) [39, 7].
The MD runs are thermalized using the pre-thermalization technique outlined in
Ref. [40] using finite-differences force constants obtained via phonopy [41].
Afterwards, a Langevin thermostat at the target temperature (300 K) is used to
perform NVT sampling. After an initial sampling period of 2.5 ps, the cell
parameters are adjusted such that thermal pressure is minimized to below $5$
kbar in order to account for thermal expansion [42]. Starting conditions for
the NVE simulations are chosen from an NVT run for the relaxed supercell at
least 2 ps apart. The time step for the MD simulation was chosen as 5 fs,
which corresponds to a tenth of the shortest period duration of the harmonic
spectrum of MgO ($\omega_{\rm max}\approx 20$ THz). The heat flux is sampled
less frequently since heat transport is dominated by the slow vibrations, and
all results are reported for a heat-flux sampling period of 20 fs. We have
checked that further decreasing the heat-flux sampling frequency does not
change the results significantly.

The work is organized as follows: In Sec. II, we review Green-Kubo theory in
order to highlight the steps necessary for the numerical implementation.
Section III presents our approach for noise reduction based on physical
arguments and real-time signal analysis which allows to truncate the GK time
integral in a numerically robust way. Section IV presents the updated version
of the size-extrapolation scheme first introduced in Ref. [21]. To complete
the method description, we discuss results for MgO in Sec. V, and compare to
available experimental and computational literature. After completing the
discussion of the method for MgO, we apply the scheme to CuI in Sec. VI, and
conclude with some remarks on simulation time convergence in Sec. VII.

##  II Green-Kubo theory

Let us start with a short summary of classical thermal transport in the
framework of GK theory [18, 19, 43]: The thermal conductivity tensor at
temperature $T$ is given as the canonical ensemble average defined by the
phase-space integral

$\displaystyle\kappa^{\alpha\beta}(T)=\frac{1}{\mathcal{Z}}\int{\rm
d}\Gamma~{}\kappa^{\alpha\beta}[\Gamma]~{}{\rm e}^{-\frac{1}{k_{\rm
B}T}\mathcal{H}[\Gamma]}~{},$ (1)

where $\Gamma\equiv({\bf R}_{1},\ldots{\bf R}_{N};{\bf P}_{1},\ldots{\bf
P}_{N})$ are phase-space configurations for $N$ atoms with positions ${\bf
R}_{I}$ and momenta ${\bf P}_{I}$. $\mathcal{H}[\Gamma]$ is the Hamiltonian of
the system with corresponding partition function $\mathcal{Z}$, $k_{\rm B}$ is
the Boltzmann constant, and $\alpha,\beta$ denote the Cartesian components of
the tensor. For each phase-space configuration $\Gamma$, the thermal
conductivity is computed as

$\displaystyle\kappa^{\alpha\beta}[\Gamma]$ $\displaystyle=\frac{V}{k_{\rm
B}T^{2}}\lim_{t_{\rm c}\to\infty}\int_{0}^{t_{\rm c}}{\rm
d}t~{}C_{JJ}^{\alpha\beta}[\Gamma](t)~{},$ (2)

with the HFACF,

$\displaystyle\begin{split}C^{\alpha\beta}_{JJ}[\Gamma](t)&=\lim_{t_{0}\to\infty}\frac{1}{t_{0}-t}\int_{0}^{t_{\rm
0}-t}{\rm d}s~{}J^{\alpha}[\Gamma({t+s})]J^{\beta}[\Gamma(s)]\\\
&\equiv\langle J^{\alpha}(t)J^{\beta}(0)\rangle~{},\end{split}$ (3)

where the phase-space points $\Gamma(t)$ in the trajectory are obtained from
the time evolution generated by the many-body Hamiltonian of the system,
$\mathcal{H}[\Gamma]$, by propagating the initial configuration
$\Gamma(0)\equiv\Gamma$ for a time $t$. $J^{\alpha}[\Gamma(t)]\equiv
J^{\alpha}(t)$ is the the heat flux component evaluated for the configuration
$\Gamma(t)$, and $\langle\cdot\rangle$ is the shorthand notation for the time
average in Eq. (3).

In order to evaluate these equations in finite simulations, the integrals need
to be discretized and truncated to finite domains. First, Eq. (1) is
approximated by taking a finite set of $M$ starting configurations
$\Gamma^{i}(0)\equiv\Gamma^{i}$, so that

$\displaystyle\kappa^{\alpha\beta}(T)\approx\frac{1}{M}\sum_{i=1}^{M}\kappa^{\alpha\beta}[\Gamma^{i}]~{},$
(4)

where the starting conditions $\Gamma^{i}$ are chosen from NVT MD simulations
for the thermodynamic conditions of interest. For each starting condition
$\Gamma^{i}$, NVE MD simulations are performed to generate the time evolution
of the system, $\Gamma^{i}(t)$, and evaluate the heat flux ${\bf J}(t)$ along
this trajectory. The simulation is performed for a total simulation time
$t_{0}$, thereby truncating the time integral in Eq. (3). This time needs to
be large enough to cover the time scales of the physical processes relevant
for heat transport. From the resulting autocorrelation function of finite
length, the thermal conductivity components are computed via Eq. (2). For each
component, a _cutoff time_ $t_{\rm c}<t_{0}$ is chosen to avoid integrating
parts of the HFACF after it has effectively decayed, since its tail can be
heavily affected by statistical fluctuations stemming from finite size and
time effects [44, 45, 46, 47], or tiny but systematic drifts that accumulate
in the long time limit, e. g., when the average flux $\langle{\bf
J}\rangle_{t}$ does not vanish exactly over the simulation time or appears
slightly skewed.

After computing the thermal conductivity tensor for each trajectory, the final
value is given by Eq. (4), i. e., by the _mean_ of the individual
trajectories. The statistical error due to the finite ensemble average is
estimated by the _standard error_ , i. e., the standard deviation of the mean,

$\displaystyle\Delta\kappa^{\alpha\beta}(T)=\frac{1}{\sqrt{N}}\sqrt{\frac{1}{N}\sum_{i}\left(\kappa^{\alpha\beta}(T)-\kappa^{\alpha\beta}[\Gamma^{i}]\right)^{2}}~{}.$
(5)

From the Cartesian components of the thermal conductivity
$\kappa^{\alpha\beta}(T)$, the scalar thermal conductivity $\kappa(T)$ is
obtained via

$\displaystyle\kappa(T)=\frac{1}{3}\sum_{\alpha}\kappa^{\alpha\alpha}(T)~{}.$
(6)

In principle, these equations can be evaluated as is, and convergence in size
and time can be checked by simply increasing the respective scales. While this
is computationally possible when using analytical force fields [48], this is
certainly not desirable in the _ab initio_ case, where the cost per time step
is considerably higher: Here, the accessible size and time scales are
typically orders of magnitude lower, and additional steps to increase the
amount of information that can be extracted from the comparatively short
simulations are pivotal.

##  III Cutoff time and noise reduction

For a robust identification of the cutoff time $t_{\rm c}$, we first reduce
noise from the HFACF as much as possible. This is achieved in two steps:
First, we re-define the _ab initio_ heat flux used in this work such that
terms not contributing to the thermal conductivity are discarded [36]. Second,
we filter remaining contributions from the HFACF that do not contribute to the
integrated thermal conductivity. This allows to determine the cutoff time even
in the presence of noise based on a “first dip” criterion, i. e., by choosing
the time when the HFACF drops below zero for the first time.

The raw _ab initio_ virial heat flux used in this work was introduced in Ref.
[21] and is given for a phase-space point $\Gamma(t)=\set{{\bf R}(t),{\bf
P}(t)}$ by

$\displaystyle{\bf J}_{\rm raw}[\Gamma(t)]\equiv{\bf J}_{\rm
raw}(t)=\frac{1}{V}\sum_{I}\sigma_{I}(t)\dot{\bf R}_{I}(t)~{},$ (7)

where $\sigma_{I}(t)\equiv\sigma_{I}[{\bf R}(t)]$ is the contribution of atom
$I$ to the virial stress tensor for the configuration ${\bf R}=({\bf
R}_{1},\ldots{\bf R}_{N})$ at the given time $t$ as derived and discussed in
Ref. [21, 49], and $\dot{\bf R}_{I}(t)={\bf P}_{I}(t)/M_{I}$ is the velocity
of atom $I$ with mass $M_{I}$. In this definition of the heat flux, convective
contributions that become important in liquids and gases are entirely
neglected [5, 50].

###  III.1 Discarding non-contributing terms

We split the raw flux in two parts,

$\displaystyle{\bf J}_{\rm
raw}(t)=\underset{1)}{\underbrace{\frac{1}{V}\sum_{I}\delta\sigma_{I}(t)\dot{\bf
R}_{I}(t)}}+\underset{2)}{\underbrace{\frac{1}{V}\sum_{I}\braket{\sigma_{I}}_{t}\dot{\bf
R}_{I}(t)}}~{},$ (8)

where $\braket{\sigma_{I}}_{t}$ is the time-averaged atomic virial, and
$\delta\sigma_{I}(t)$ is the time-dependent part. In the absence of diffusion,
the second term is the total time derivative of a bounded vector field,
$\sum_{I}\braket{\sigma_{I}}\dot{\bf R}_{I}(t)=\frac{{\rm d}}{{\rm
d}t}\sum_{I}\braket{\sigma_{I}}{\bf U}_{I}(t)$, where ${\bf U}_{I}(t)={\bf
R}_{I}(t)-{\bf R}_{I}$ is the displacement of atom $I$ from its reference
position in the lattice, ${\bf R}_{I}$. Contributions to the heat flux that
can be written in this form do not contribute to the integrated transport
coefficient, as can be elegantly shown using the “gauge theorem” discussed in
detail in Ref. [36, 51, 28]. We therefore discard the second term from the
flux, and proceed using the following gauge-fixed heat flux expression:

$\displaystyle{\bf J}(t)=\frac{1}{V}\sum_{I}\delta\sigma_{I}(t)\dot{\bf
R}_{I}(t)~{}.$ (9)

In compounds with two or more elements, the individual average virials
$\langle\sigma_{I}\rangle$ can be significant, and discarding the non-
contributing part from the raw heat flux reduces the noise in the HFACF
considerably, as shown for the case of MgO in Fig. 1 (red curves compared to
gray curves). We note that this amount of noise reduction is difficult to
achieve by means of mere filtering: The blue curves in Fig. 1 are obtained by
filtering the raw HFACF obtained with the flux defined in Eq. (8) with the
filter discussed below in Sec. III.2. It is apparent that leveraging the gauge
theorem by using the flux defined in Eq. (9) instead reduces the variance in
the HFACF much more significantly, and is furthermore physically rigorous.
Finally, we enforce a vanishing expectation of the flux to remove bias from
the resulting quantities due to the finite time of the simulation by removing
the finite-time average, ${\bf J}(t)\to\delta{\bf J}(t)={\bf
J}(t)-\braket{{\bf J}}_{t}$.

We note in passing that the above argument holds for any heat flux that can be
written in a virial-based form similar to Eq. (7), e. g., common heat flux
expressions for empirical or machine-learned force fields [52]. For such force
fields, a noise-reduction approach similar to the one presented here was very
recently developed independently and applied successfully by Pereverzev and
Sewell [53].

Figure 1: Heat flux autocorrelation function (HFACF) $C^{yy}_{JJ}(t)$ in MgO
for the $yy$-component as defined in Eq. (3), and its cumulative integral, i.
e., the thermal conductivity $\kappa^{yy}(t)$ defined in Eq. (10) as function
of integration time $t$. Shaded gray: $C^{yy}_{JJ}(t)$ and $\kappa^{yy}(t)$
obtained by using the raw flux as defined in Eq. (7). Purple: Using only the
filter discussed in Sec. III.2. Red: Using the gauge-fixed flux defined in Eq.
(8). Black dashed curves: After applying noise filtering to the gauge-fixed
quantities as explained in the main text. The cutoff time $t_{\rm c}$ is
chosen based on the “first dip” of the gauge-fixed and noise-filtered HFCAF.

###  III.2 Noise filtering

After obtaining the gauge-fixed heat flux by discarding the non-contributing
term, there is still a considerable level of noise in the HFACF which hinders
a robust identification of the time at which it is effectively decayed, i. e.,
the cutoff time $t_{\rm c}$. Available techniques to identify cutoff times,
such as the first avalanche method introduced in Ref. [35] typically require
system-dependent parameters, such as a tolerable signal-over-noise ratio or
window sizes for moving average computation. To overcome this issue, we
suggest an approach that does rely only on one single parameter which is
chosen based on the vibrational spectrum of the material: Motivated by the
fact that the _integrated_ HFACF, i. e., the cumulative thermal conductivity

$\displaystyle\kappa(t)=\frac{V}{k_{\rm B}T^{2}}\int_{0}^{t}{\rm
d}t^{\prime}~{}C_{JJ}(t^{\prime})~{},$ (10)

is already a much smoother function than the HFACF itself, we apply a moving
window average to $\kappa(t)$ instead of $C_{JJ}(t)$. The remaining parameter
is the window size for the filter. It is chosen based on the vibrational
spectrum of the material by taking the period length corresponding to the
slowest significant frequency, $t_{\rm window}=1/\omega_{\min}$, which is
chosen to be the first peak in the vibrational density of states (VDOS). To
ensure that $\kappa(t)$ vanishes identically at $t=0$, the filter is applied
to the cumulative thermal conductivity extended antisymmetrically to negative
times via $\kappa(-t)=-\kappa(t)$, a property which follows from the time
symmetry of $C_{JJ}(t)$ [54]. Thereby all noise and non-contributing parts of
higher frequency are effectively filtered from $\kappa(t)$, while all relevant
time integrals are preserved by construction. As required, the cumulative
kappa before (red curve) and after filtering (black dashed curve) lie right on
top of each other in the lower panel of Fig. 1. As in the previous section, we
stress that using the gauge-fixed heat flux is crucial to access the
integrated $\kappa(t)$, as shown in the lower panel of Fig. 1 in comparison to
the case where only the filter was used to smoothen $\kappa(t)$ obtained from
the raw flux (blue curve).

The filtering is carried over to the HFACF, $C_{JJ}(t)$, by numerically
differentiating the filtered cumulative thermal conductivity with respect to
time and applying the same filter on the numerical gradient of $\kappa(t)$.
The resulting HFACF $C_{JJ}(t)$ is shown as a black dashed curve in the upper
panel of Fig. 1, and the further reduced level of variance in the HFACF is
apparent. As seen in Fig. 1, $\kappa(t)$ reaches a plateau at approximately 8
ps that lasts for several picoseconds. After that, numerical noise dominates
and the accumulated numerical errors lead to a drop in $\kappa(t)$. The “first
dip” criterion is used to detect this plateau numerically in a reliable
fashion. For this purpose, a cutoff time $t_{\rm c}$ is chosen that
corresponds to the time when the signal-over-noise ratio vanishes, i. e., when
$C_{JJ}(t)$ drops to zero [35]. Note that also different numerical approaches
[55, 56, 29, 57] or a visual inspection of $\kappa(t)$ can be used to identify
this plateau. The main advantage of the “first-dip” criterion used here is
that no numerical parameters have to be chosen, which facilitates automatic
evaluation and the systematic comparison of aiGK simulations with different
trajectory lengths $t_{0}$. In turn, this results in smoothly converging
thermal conductivities with respect to $t_{0}$ as shown below in see Fig. 6
and 7. This reflects that the final, converged values for $\kappa(t)$ are
virtually independent from the details used to choose the cutoff time $t_{\rm
c}$. With the cutoff time $t_{\rm c}$, the resulting thermal conductivity for
a given component of the thermal conductivity tensor is given by the value
$\kappa=\kappa(t_{\rm c})$ as indicated by the horizontal line in Fig. 1. The
presented scheme will be used for all reported values of thermal conductivity
in the following.

We note that this filter corresponds to a low-pass filter in Fourier space
[29]. However, since we found the real-time noise-reduction scheme as
presented above sufficient to obtain robust results, we did not investigate
additional processing steps that involve discrete (inverse) Fourier transforms
with further numerical parameters and potential aliasing problems when only a
few thousand data points are available.

##  IV Size extrapolation for crystals

After we have seen how the GK formula is used to compute thermal
conductivities from the _ab initio_ heat flux evaluated along aiMD
trajectories, we discuss an update to the size-correction scheme for crystals
first introduced in Ref. [21]. The aim of this extrapolation is to correct for
size effects occurring in aiMD simulations because phonon modes of longer
wavelength than the supercell dimensions are not included. This is of
particular importance in crystals with periodic long-range order where these
modes can contribute significantly to heat transport. The correction works by
computing the harmonic contribution to the thermal conductivity $\kappa_{\rm
ha}$ within the supercell via

$\displaystyle\kappa_{\rm ha}^{\alpha\beta}=Vk_{\rm B}\sum_{b{\bf
q}}v^{\alpha}_{b{\bf q}}v^{\beta}_{b{\bf q}}{\tau}^{\phantom{\dagger}}_{b{\bf
q}}~{},$ (11)

where $V$ is the system volume, $k_{\rm B}$ is the classical heat capacity per
phonon mode, $v_{b{\bf q}}$ is the group velocity of a phonon mode with band
index $b$ and wave vector $\bf q$, and $\tau_{b{\bf q}}$ is the lifetime of
the mode extracted from the aiMD trajectory. This contribution is first
computed for the wave vectors $\bf q$ commensurate with the supercell,
$\kappa_{\rm ha-supercell}$, and then _extrapolated_ to bulk limit by
_interpolating_ the lifetimes to denser $\bf q$-meshes in the Brillouin zone,
$\kappa_{\rm ha-bulk}$. The resulting size-corrected thermal conductivity is
obtained as

$\displaystyle\kappa_{\rm corrected}=\kappa+\underset{\delta\kappa_{\rm ha-
correction}}{\underbrace{\kappa_{\rm ha-bulk}-\kappa_{\rm ha}}}~{}.$

The necessary steps to compute $\kappa_{\rm ha}$, $\kappa_{\rm ha-bulk}$, and
therefore $\delta\kappa_{\rm ha-correction}$ are presented below.

In passing, we like to mention the main differences to the original approach
[21]. Space-group symmetries are now systematically exploited and all
reciprocal-space quantities are generated from the irreducible part of the
Brillouin zone. This guarantees unambiguous branch matching during the
interpolation. Further minor changes are that lifetimes are extracted from
fitting an exponential decay, cf. Sec. IV.2, and scaled lifetimes are
interpolated linearly instead of using Fourier interpolation, cf. Sec. IV.3.

###  IV.1 Harmonic mapping

In order to map the real-space dynamics to the phonon picture which allows for
interpolating in reciprocal space, we first define an auxiliary harmonic model
determined by the real-space dynamical matrix

$\displaystyle
D_{{IJ}}^{{\alpha\beta}}=\frac{1}{\sqrt{M_{I}M_{J}}}\Phi_{{IJ}}^{{\alpha\beta}}~{},$
(12)

where $\Phi_{{IJ}}^{{\alpha\beta}}$ are the $(\alpha,\beta)$ components of
harmonic force constants in the supercell between atom pairs $(I,J)$. As noted
earlier, we obtain the force constants via the TDEP method to account for
finite-temperature renormalization of phonon frequencies, eigenvectors, and
group velocities [39, 7].

Using the crystal periodicity, the Fourier-transformed dynamical matrix reads

$\displaystyle D_{{\bf q},ij}^{\alpha\beta}=\sum_{L}{\rm e}^{{\rm i}{\bf
q}\cdot({\bf R}_{i}-{\bf R}_{j}-{\bf R}_{L})}D_{i0,jL}^{\alpha\beta}~{},$ (13)

where $\set{{\bf R}_{i},{\bf R}_{j}}$ denote reference positions in the unit
cell, ${\bf R}_{L}$ is a Bravais lattice vector, and $\bf q$ is a commensurate
wave vector fulfilling ${\bf q}\cdot{\bf R}_{L}=2\pi n$ with an integer $n$.
The dynamical matrix yields real eigenvalues $\omega^{2}_{b{\bf q}}$ and
complex eigenvectors ${\bf e}_{{b{\bf q}},i}$ via the eigenvalue equation

$\displaystyle\sum_{j\beta}D_{{\bf q},ij}^{\alpha\beta}~{}e_{{b{\bf
q}},j}^{\beta}=\omega^{2}_{b{\bf q}}e_{{b{\bf q}},i}^{\alpha}~{},$ (14)

where the band index $b$ was introduced to discern branches of solutions. To
directly translate between real space coordinates $I=(i,L)$ and reciprocal
space coordinates $b,{\bf q}$, we define the generalized eigenvector

$\displaystyle{\bf e}_{{b{\bf q}},I}\equiv\frac{1}{\sqrt{N_{\bf q}}}{\rm
e}^{-{\rm i}{\bf q}\cdot{\bf R}_{I}}\,{\bf e}_{{b{\bf q}},i}$ (15)

with ${\bf R}_{I}={\bf R}_{i}+{\bf R}_{L}$, which diagonalizes the real-space
dynamical matrix $D_{IJ}$ defined in Eq. (12), where $N_{\bf q}$ is the number
of lattice points or commensurate wave vectors in the supercell, respectively.

Using the generalized eigenvector defined in Eq. (15), we define normal
coordinates $u_{b{\bf q}}$ and $p_{b{\bf q}}$ as

$\displaystyle\begin{split}u_{b{\bf q}}(t)&=\sum_{I}\sqrt{M_{I}}\,{\bf
e}_{{b{\bf q}},I}\cdot{\bf U}_{I}(t)~{},\\\ p_{b{\bf
q}}(t)&=\sum_{I}\frac{1}{\sqrt{M_{I}}}\,{\bf e}_{{b{\bf q}},I}\cdot{\bf
P}_{I}(t)~{},\end{split}$ (16)

where ${\bf U}_{I}(t)={\bf R}_{I}(t)-{\bf R}_{I}$ is the instantaneous
displacement of atom $I$ from its reference position ${\bf R}_{I}$, and ${\bf
P}_{I}(t)$ is its momentum as before. From here, the time-dependent complex
mode amplitude $a_{b{\bf q}}(t)$ follows [58],

$\displaystyle a_{b{\bf q}}(t)$
$\displaystyle=\frac{1}{\sqrt{2}}\left(u_{b{\bf q}}(t)+\frac{{\rm
i}}{\omega_{b{\bf q}}}p_{b{\bf q}}(t)\right)~{},$ (17)

from which the time-dependent mode-resolved energy is obtained,

$\displaystyle E_{b{\bf q}}(t)=\omega_{b{\bf q}}^{2}\,{a}^{\dagger}_{b{\bf
q}}(t){a}^{\phantom{\dagger}}_{b{\bf q}}(t)~{},$ (18)

where we note in passing that the harmonic energy expression familiar from
quantum mechanics is recovered by substituting $a_{b{\bf
q}}\to\sqrt{\hbar/\omega_{b{\bf q}}}\,a_{b{\bf q}}$.

Using the mode-resolved energy, the harmonic heat flux can be defined as [59]

$\displaystyle{\bf J}_{\rm ha}(t)=\frac{1}{V}\sum_{b{\bf q}}E_{b{\bf
q}}(t){\bf v}_{b{\bf q}}~{}.$ (19)

For a classical harmonic system, we can use that $\langle E^{2}_{b{\bf
q}}\rangle=(k_{\rm B}T)^{2}$ and that cross correlations between different
modes $(b,{\bf q})\neq(b^{\prime},{\bf q}^{\prime})$ vanish, so that the
harmonic thermal conductivity follows by application of Eq. (2)-(3),

$\displaystyle\kappa^{\alpha\beta}_{\rm ha}=Vk_{\rm B}\sum_{b{\bf
q}}v^{\alpha}_{b{\bf q}}v^{\beta}_{b{\bf q}}\int_{0}^{\infty}G_{b{\bf
q}}(t)~{},$ (20)

with the mode-energy autocorrelation function

$\displaystyle G_{b{\bf q}}(t)=\frac{\braket{E_{b{\bf q}}(t)E_{b{\bf
q}}(0)}}{\braket{E_{b{\bf q}}^{2}}}~{}.$ (21)

The lifetime $\tau_{b{\bf q}}$ in Eq. (11) is therefore defined by the
integral

$\displaystyle\tau_{b{\bf q}}\equiv\int_{0}^{\infty}G_{b{\bf q}}(t)~{}.$ (22)

###  IV.2 Lifetime extraction

For a purely harmonic system, the mode-energy autocorrelation function
$G_{b{\bf q}}(t)$ does not decay, and the lifetime given by Eq. (22) diverges.
In the presence of phonon-phonon interactions due to anharmonicity however,
$G_{b{\bf q}}(t)$ decays and the resulting lifetime is finite.

In order to evaluate the lifetime $\tau_{b{\bf q}}$ via Eq. (22), it is again
necessary to integrate an autocorrelation function. As already mentioned in
the introduction and Sec. III, this typically requires to choose an
appropriate cutoff time to obtain numerically stable results. In the case of
the lifetimes $\tau_{b{\bf q}}$, we exploit the advantage that we can
integrate distinct phonon modes individually, in contrast to Sec.III in which
the flux for the whole system was processed at once. In perturbation theory,
the leading contribution to the decay of each phonon mode can be approximated
via [60]

$\displaystyle G_{b{\bf q}}(t)\approx{\rm e}^{-2{\rm Im}\Sigma_{b{\bf
q}}t}\equiv{\rm e}^{-t/\tau_{b{\bf q}}}~{},$ (23)

where ${\rm Im}\Sigma_{b{\bf q}}$ is the imaginary part of the phonon self
energy, and $\tau_{b{\bf q}}=1/2{\rm Im}\Sigma_{b{\bf q}}$ is the
corresponding lifetime. Under this approximation, the integral can be
performed analytically based on the early decay of $G_{b{\bf q}}(t)$, hence
allowing to capture also those long-lived, long-wavelength modes that are not
guaranteed to be accessible via brute-force integration on the simulation time
scales accessible in aiMD simulations. We compute $G_{b{\bf q}}(t)$ for each
mode in the supercell, and obtain the corresponding lifetime by fitting Eq.
(23) for times where $G_{b{\bf q}}(t)>0.1$, in order to avoid fitting noise
when $G_{b{\bf q}}(t)$ is effectively decayed. A comparison of numeric
correlation functions via Eq. (21), and the respective analytic correlation
functions given by Eq. (23), is shown in Fig. 2. We find that the analytic
behavior, i. e., exponential decay of $G_{b{\bf q}}(t)$, is indeed observed
for many modes, in particular those with long lifetime. Some modes deviate
more strongly from the exponential decay, for example the mode highlighted in
Fig. 2 where $G_{b{\bf q}}(t)$ shortly increases after about 2 ps, before
dropping off again. The integrated correlation function, i. e., the lifetime,
is however only mildly affected from wiggles like this.

Figure 2: Fit of mode lifetimes for MgO at 300 K. Simulation performed with a
timestep of 5 fs for a simulation time of 60 ps. Top: Normalized mode-energy
autocorrelation function $G_{b{\bf q}}(t)$ as obtained from the simulation by
Eq. (21). Bottom: Analytic expression given by Eq. (23) after fitting mode
lifetimes $\tau_{b{\bf q}}$. The correlation functions for the six modes at
${\bf q}=(-1/6,-1/6,0)$ are highlighted in red for comparison for times
$t<3\tau_{b{\bf q}}$. The y-axis is logarithmic such that exponential
functions appear as straight lines.

###  IV.3 Lifetime interpolation

For a given simulation $\set{\Gamma^{i}(t)}$, the lifetimes $\tau_{b{\bf q}}$
are evaluated for all commensurate $\bf q$-points, and projected to the
symmetry-inequivalent points in the Brillouin zone determined by the space
group operations of the system to improve the statistics: The irreducible
q-points in the Brillouin zone are obtained by iteratively reducing the given
grid with the available symmetry operations for the system obtained by the
spglib package [61]. To avoid band-index matching problems between different
$\bf q$-points, the eigenvectors for the full grid of commensurate $\bf
q$-points are created by solving the eigenvalue problem in Eq. (14) on the
irreducible grid, and transforming the eigenvectors to the reducible points
according to the transformation rules given in Ref. [62].

In the next step, the lifetimes $\tau_{b{\bf q}}$ are interpolated to denser
$\bf q$-point meshes. For this purpose, the fully anharmonic lifetimes
$\tau_{b{\bf q}}$ at the commensurate $\bf q$-points are used to define one
function $\lambda_{b}({\bf q})$ for each branch $b$ such that

$\displaystyle{{\tau}}^{\phantom{\dagger}}_{b}({\bf
q})={\lambda}^{\phantom{\dagger}}_{b}({\bf q})\omega_{b}^{-2}({\bf q})~{}.$
(24)

The frequency scaling ensures that $\lambda_{b}({\bf q})$ is only weakly ${\bf
q}$-dependent, which facilitates a linear interpolation of the lifetimes to
arbitrary values $\tilde{\bf q}$ in the Brillouin zone via Eq. (24). For the
acoustic modes at ${\bf q}=\Gamma=0$, where $\omega({\bf q\to 0})\to 0$, the
value for $\lambda_{b}(\Gamma)$ is obtained by averaging over values at the
surrounding $\bf q$-points. The scaling of lifetimes with
$\omega_{b}^{-2}({\bf q})$ used here is rooted in basic phonon theory as
discussed in detail by Pomeranchuk and Herring, but it is not universal [63,
64]. Other scaling laws can be obtained by using different limiting
assumptions that can, for example, depend on the crystal structure [64].
However, the quadratic scaling used here is generally the strongest possible
variation consistent with non-diverging thermal conductivities in the limit of
dense Brillouin zone sampling irrespective of further limiting assumptions
[63]. Therefore, it leads to a firm upper bound for the effect of size
extrapolation. Also due to the fact that this interpolation scheme is fully
mode- and q-resolved via the function $\lambda_{b}({\bf q})$ and incorporates
the fully anharmonic lifetimes at commensurate ${\bf q}$-points, no systematic
errors associated to the scaling and interpolation procedure have been
observed in extended validation calculations covering several different
crystal structures [65].

For sampling the interpolated points $\set{\tilde{\bf q}}$, we use even-
numbered Monkhorst-Pack grids as implemented in phonopy [66, 41], with a
maximum grid size of $20\times 20\times 20$. The symmetrized lifetimes
$\tau_{b{\bf q}}$ obtained from fitting Eq. (23), as well as the interpolated
lifetimes denoted by $\tilde{\tau}_{b}(\tilde{\bf q})$ are displayed as
scattering rates ($\propto\tau^{-1}$) in Fig. 3.

Figure 3: Scattering rates vs. frequency in MgO at 300 K after fitting Eq.
(23) and symmetrizing using space-group operations as explained in the main
text (black dots), and after interpolation to $20\times 20\times 20$ grid (red
dots). Gray dashed: Ioffe-Regel limit $\tau^{-1}=\omega$ with 50 % margins.
Figure 4: Scattering rates vs. frequency in MgO at 300 K. Comparisons of the
lifetimes obtained via interpolation, cf. Fig. 3, to those obtained using
third-order perturbation theory with the TDEP code [39, 7].

It is apparent that long-lived contributions stemming from modes with
$\omega<5$ THz are introduced through the interpolation (red dots),
effectively capturing modes with lifetimes $>100$ ps which is beyond the
simulation time of 60 ps. On the other hand, lifetimes for modes $\omega>5$
THz are already well-captured by the commensurate modes (black dots). This is
further substantiated in Fig. 4, which compares the lifetimes obtained by the
described interpolation procedure with those obtained using third-order
perturbation theory with TDEP force constants fitted to our MD simulations
[7]. The difference observed between the perturbative and the interpolated
scattering rates in the 7–12 THz range is attributed to higher-order
anharmonic scattering, in line with the findings discussed for MgO at higher
temperatures in Ref. [67] using a similar harmonic mapping procedure in larger
supercells without interpolation. We also note that all scattering rates are
well below the Ioffe-Regel limit $\tau^{-1}=\omega$ [68]. This limit can be
taken as a qualitative rule for estimating the validity of the phonon
quasiparticle picture: Since the scattering rate is defined as the width of
the phonon spectral function via Eq. (23), $\tau_{b{\bf q}}^{-1}=2{\rm
Im}\Sigma_{b{\bf q}}$, rates below this limit signify sharply peaked, well-
defined quasiparticles [12, 20].

###  IV.4 Thermal conductivity extrapolation

For the new, denser grid, an interpolated value,

$\displaystyle\kappa_{\rm ha-int}^{\alpha\beta}(N_{\tilde{\bf q}})=Vk_{\rm
B}\frac{N_{\bf q}}{N_{\tilde{\bf q}}}\sum_{b,\tilde{\bf
q}}v^{\alpha}_{b}(\tilde{\bf q})v^{\beta}_{b}(\tilde{\bf
q}){\tilde{\tau}}^{\phantom{\dagger}}_{b}(\tilde{\bf q})~{},$ (25)

can be obtained, where $N_{\tilde{\bf q}}$ is the number of points in the new
grid, and the factor $N_{\bf q}/N_{\tilde{\bf q}}$ accounts for the increased
number points. The bulk limit of Eq. (25) is obtained by computing
interpolated values for an increasing density of $\bf q$-points. The
convergence of Eq. (25) is approximately linear in $N_{\tilde{\bf
q}}^{-1/3}\equiv 1/n_{q}$, where $n_{q}$ is number of $\bf q$-points per
Cartesian direction. The slope of this curve can therefore be used to
extrapolate the value of $\kappa_{\rm ha}$ to bulk limit, as shown in Fig. 5.

Figure 5: Size extrapolation correction to bulk limit computed from Eq. (25)
assuming linear convergence in $1/n_{q}$, where $n_{q}$ is the number of
q-points per Cartesian direction. The offset between $\kappa_{\rm ha}$ and
$\kappa_{\rm ha-int}$ arises because even grids are used for the
extrapolation, whereas a $3\times 3\times 3$ supercell and respective grid of
commensurate ${\bf q}$-points are used in the simulation.

With the extrapolated value $\kappa_{\rm ha-bulk}$, a correction can be
obtained via

$\displaystyle\delta\kappa_{\rm ha-correction}=\kappa_{\rm ha-
bulk}-\kappa_{\rm ha}~{},$ (26)

from which the final result for the thermal conductivity is obtained via

$\displaystyle\kappa^{\alpha\beta}_{\rm
corrected}=\kappa^{\alpha\beta}+\delta\kappa^{\alpha\beta}_{\rm ha-
correction}~{},$ (27)

where $\kappa^{\alpha\beta}$ is the value from the aiGK simulation. The
interpolation scheme effectively subtracts harmonic contributions to the
thermal conductivity from vibrations commensurate with the supercell, and
extrapolates them to the bulk limit, thereby including long-range
contributions otherwise not present in the simulation cell. The size-corrected
contributions are subsequently added back to the total thermal conductivity.

We note that several approximations are involved in the scheme outlined above,
such as the assumption of exponential decay of the mode-energy autocorrelation
function in Eq. (23), or neglecting mode cross-correlations with $(b,{\bf
q})\neq(b^{\prime},{\bf q}^{\prime})$ in Eq. (11). However, the dominant
contribution to the size correction in Eq. (26) can be expected to come from
low-frequency, long-lived phonons missing in the simulation cell, as shown in
Fig. 3, for which the approximations listed above are well justified. The
shorter-range contributions of modes that interact more strongly are fully
captured on the _ab initio_ level, for which the lifetimes are only weakly
frequency-dependent, as seen in Fig. 3 when focusing on the regime where
$\omega>5\,{\rm THz}$.

We also note in passing that the force constants used for size extrapolation
could also be obtained via finite differences as, e. g., implemented in
phonopy [69, 41, 70]. However, especially for quite harmonic materials such as
MgO, we didn’t find the extrapolation scheme to be sensitive to subtle
differences in the force constants used to describe the auxiliary (effective)
harmonic model. Furthermore, effective harmonic models extend the
applicability of phonon theory to dynamically unstable high-temperature
phases, for example in SrTiO3 [71], PdH [8], or ZrO2 [72, 21]. A detailed
account of the feasibility of the size-extrapolation scheme presented here for
these systems is, however, beyond the scope of the current work.

##  V Results for magnesium oxide

After we have seen how the cutoff time $t_{\rm c}$ in Eq. (2) can be obtained,
and finite-size errors can be corrected, we discuss the convergence of
presented scheme as a function of the simulation time $t_{0}$ in Eq. (3). We
do this for the case of MgO for three independent trajectories of 60 ps length
each. We truncate every trajectory in 10 % steps down to a length of 6 ps, and
apply the workflow presented in the previous sections to each of the truncated
trajectories.

Figure 6: Thermal conductivity $\kappa$ as function of the simulation time
$t_{0}$ as defined in Eq. 29. Values are given as the ensemble average over
three independent trajectories. The error bars are computed according to Eq.
(5) as the standard error of the ensemble average. The blue curve is a
logistic curve defined in Eq. (28) fitted to the $\kappa$ values, the dashed
blue curve is the infinite time limit of the fitted function. Gray dots
represent the thermal conductivity as given by the simulation without the
size-correction scheme. Please note that the values for $\kappa$ shown here
cannot be directly compared to the value displayed in Fig. 1 or 5, because the
latter only show single components of single runs, which can vary
substantially from the total average.

Figure 6 shows that the thermal conductivity converges to a plateau after
about 40 ps, where the value of $\kappa$ stays constant within the error bars.
The overall shape of the curve can be described as follows: Simulations
shorter than 20 ps sample the early decay of the HFACF which contribute about
30 W/mK to the total thermal conductivity After a simulation time of 25 ps,
the late decay of the HFACF is sampled, contributing more than double the
amount to the total thermal conductivity of $68.8\pm 6.1$ W/mK after the total
simulation time. In the plot, this two-step behavior is approximated by a
logistic function

$\displaystyle f(t)=\frac{L}{1+\exp\left(-\frac{(t-t_{\rm
inflection})}{\tau}\right)}+f_{0}~{},$ (28)

which captures the second super-linear increase in $\kappa$ at $t_{\rm
inflection}\simeq 29\,{\rm ps}$, and models the correct asymptotic behavior
for long times. The asymptotic value of $\kappa(t_{0}\to\infty)=68.5$ W/mK
agrees very well with the value after 60 ps of $\kappa(60\,{\rm ps})=68.8\pm
6.1$ W/mK. As the largest lifetime in the simulation corresponds to
$\tau\approx 15.5$ ps, as highlighted in Fig. 2, we do not expect a
significant increase of $\kappa(t_{0})$ after this simulation time. We
therefore conclude that the simulation time of $t_{0}=60$ ps can be considered
converged, and that Eq. (28) can be used to model the late increase of
$\kappa(t_{0})$, in line with the division into short and long processes
commonly discussed in the literature [73, 74, 75]. Furthermore, we note that
the size extrapolation increases the value from $\kappa=45.5$ W/mK to
$\kappa_{\rm corrected}=68.8$ W/mK, i. e., the value increases significantly
by about 50 %.

We like to point out that we did not discard an initial time from the NVE
simulations to allow for further equilibrating after the thermostat is
switched off, as is common practice in the literature on GK simulations [75].
We did not find this procedure to be necessary: The truncation of simulation
times displayed in Fig. 6 was performed such that the _early_ time steps in
the simulation were discarded, i. e., a simulation time of $t_{0}=54$ ps
corresponds to discarding the first 6 ps from the trajectory. Since discarding
6 ps or more did not change the result significantly, its effect can be
assumed to be minor for aiGK simulations.

###  V.1 Comparison to literature

We conclude the discussion for MgO by comparing to available experimental and
theoretical references. These references are listed in Tab. 1.

Reference | Thermal conductivity
---|---
| at 300 K (W/mK)
Experiment [76, 77, 78, 79, 80, 81, 82] | 50-65
de Koker 2010 (LDA) [83] | $\approx 75^{\,\dagger}$
Stackhouse et al. 2010 (LDA) [33] | $58\pm 6^{\,\dagger}$
Tse et al. 2018 (PBE) [84] | $70.3\pm 8.9$
Dekura and Tsuchiya 2017 (LDA) [85] | $\approx 54$
Plata et al. 2017 (PBE) [86] | $54.06$
Xia et al. 2020 (PBE) [15] | $50.1-58.7$
This work | $68.8\pm 6.1$
Table 1: Reference values for the thermal conductivity of periclase MgO at
ambient conditions. Values marked with $\dagger$ are extrapolated values using
data from higher temperatures using Eq. (17) in Ref. [83] and Eq. (5) in Ref.
[33], respectively. See also discussion in Ref. [87]. For the _ab initio_
studies, the level of theory is indicated in parentheses: local-density
approximation (LDA) [88], or the GGA parametrized by Perdew, Burke, and
Ernzerhof (PBE) [89].

The given references show a significant spread. On the experimental side, the
main source of uncertainty arises from different sample quality, different
measurement techniques, and the fact that thermal conductivity is usually not
the direct observable, but inferred from thermal diffusivity measurements and
other material parameters which may imply additional sources of error [82].
The agreement between our aiGK simulation and experiment is satisfactory.
While our thermal conductivity is larger than that of the listed experiments,
this is to be expected since isotope effects are not included in our
simulations. The listed measurements on the other hand are not performed for
isotopically pure MgO. Neglecting these isotope scattering effects can lead to
increases of thermal conductivity. Using TDEP with force constants fitted to
our PBEsol trajectories [39, 7], we obtain an increase of 35% in the single-
mode relaxation time approximation and of 46% when solving the full linearized
BTE, in line with the 46 % increase reported in Ref. [34] using the LDA
functional. Corrected by this factor, our result would be $47.1\pm 4.2$ W/mK,
only slightly below and within error of the most recent experimental results
of 50.1 W/mK obtained by Hofmeister [82].

On the theoretical side, we compare to three other approaches based on aiMD
simulations by de Koker [32, 83], Stackhouse, Stixrude, and Karki [33], and
Tse et al.[84]. The quantitative agreement with de Koker and Stackhouse et al.
is acceptable, given that different xc functionals and otherwise computational
settings were used. We note that the higher values found by de Koker are a
little surprising, given that smaller supercells were used, and no size
extrapolation. However, the value of $\kappa\approx 75$ W/mK listed in Tab. 1
is an extrapolation from higher temperatures, where finite-size effects are
likely less important. The agreement with the study by Tse et al. [84] based
on the Einstein relation introduced in Ref. [90] is very good.

The other theoretical works are based on perturbative BTE approaches [10]. The
listed references use three-phonon scattering to compute phonon lifetimes, the
lowest value reported by Xia et al. [15] is obtained by additionally including
fourth-order scattering which further reduces lifetimes in MgO [9]. The BTE
approaches based on third-order scattering listed here account for isotope
scattering and are therefore consistently lower than our aiGK value, since
isotope scattering is more pronounced than higher-order phonon-phonon
scattering in MgO.

Given the comparatively large uncertainty inherent to thermal conductivity
measurement and simulation, the agreement between aiGK and the available
literature can therefore be considered satisfactory. The discussion for MgO
further shows that aiGK can be used for mostly harmonic materials with
considerable phonon lifetimes _when_ a suitable extrapolation scheme is
employed.

##  VI Results for copper iodide

Next, we apply the scheme as presented above to marshite CuI, a strongly
anharmonic material which becomes a superionic conductor above 643 K [91, 92].

The final ensemble-averaged thermal conductivity as function of the simulation
time is displayed in Fig. 7.

Figure 7: Thermal conductivity $\kappa$ for CuI as function of the simulation
time $t_{0}$. Curves and symbols are equivalent to those defined in Eq. 6.

Fitting the logistic function in Eq. (28) to $\kappa(t)$ as before, we can
pinpoint the second superlinear increase in $\kappa(t)$ to $t_{\rm
inflection}\simeq 18\,{\rm ps}$, although this increase is visually less
pronounced than in MgO. The early increase is dominated by the considerably
stronger scattering and therefore shorter lifetimes of modes in the optical
range, as shown in Fig. 8.

Figure 8: Scattering rates vs. frequency in CuI at 300 K after fitting Eq.
(23) and symmetrizing using space-group operations as explained in the main
text (black dots), and after interpolation to $20\times 20\times 20$ grid (red
dots). Gray dashed: Ioffe-Regel limit $\tau^{-1}=\omega$ with 50 % margins.

As the lifetimes in the lower frequency range are comparable to those of MgO,
we conclude that the simulation time of 60 ps is sufficient for CuI as well.
The size extrapolation scheme increases the thermal conductivity from
$\kappa=1.12$ W/mK to $\kappa_{\rm corrected}=1.38$ W/mK, i. e., by about 23
%, which is in line with the assumption that finite-size effects become less
important in strongly anharmonic materials.

As CuI is less abundant than MgO, it is less frequently covered in the
literature. The available reference is summarized in Tab. 2.

Reference | Thermal conductivity
---|---
| at 300 K (W/mK)
CRC Handbook [93] (experiment, bulk) | 1.68
Yang et al. [94] (experiment, thin film) | 0.55
Togo et al. [31] (theory) | 6.55–7.22
This work | $1.38\pm 0.14$
Table 2: Experimental values and one theoretical reference for the thermal
conductivity of marshite CuI at ambient conditions. The value from Yang et al.
marked by $\dagger$ is from a thin film experiment, and therefore can be
regarded as a lower bound of the bulk thermal conductivity [94].

We slightly underestimate the CRC Handbook reference of $1.68$ W/mK [93], but
clearly above the thin-film reference of about $0.55$ W/mK as reported by Yang
and coworkers [94], which we take as a firm lower limit to the intrinsic
thermal conductivity of CuI due to boundary scattering. It is noteworthy that
an earlier computational investigation based on Boltzmann transport theory by
Togo and coworkers reaches a much higher value of thermal conductivity in CuI
of about 7 W/mK [31]. We propose the following explanation based on the
findings for other strongly anharmonic zincblende compounds presented by Xia
and coworkers in Ref. [15]: The authors showed that higher-order phonon
scattering can limit the thermal conductivity in zincblende compounds
considerably, and including only third-order scattering can overestimate
$\kappa$ by up to 450 % in the case of HgTe, a compound which is less
anharmonic than CuI according to the quantification scheme presented in Ref.
[14]. Since the aiGK method is non-perturbative, anharmonic scattering up to
arbitrary order is naturally included, which explains the variance with BTE
results using third-order scattering only. This is further supported by
comparing the scattering rates displayed in Fig. 8 to the Ioffe-Regel
introduced earlier [68]: Nearly all modes in the optical part of the spectrum
$>2.5$ THz approach or exceed the range $\tau^{-1}=\omega$. Several modes in
the range 1-2 THz approach this limit. As discussed in Sec. IV.3 for MgO, a
scattering rate in this range signifies a strong broadening of the phonon
spectral function beyond the perturbative regime which is defined for
$\tau^{-1}\ll\omega$ [12, 20].

##  VII Remark on simulation times

It is clear that the minimal necessary simulation time $t_{0}$ is material-
dependent and needs to be checked in each study. It may therefore come as a
surprise that 60 ps turned out to sufficient both for MgO and CuI, although
their dynamical and anharmonic properties are quite different, and one might
_a priori_ expect much longer simulation times to be necessary for the more
harmonic MgO. However, when inspecting the vibrational properties of both
materials, one can infer that the _effective_ simulation time for MgO is
indeed much longer than for CuI: We define the dimensionless effective
simulation length via

$\displaystyle\tilde{t}_{0}=t_{0}\cdot\bar{\omega}_{\rm low}~{},$ (29)

where $t_{0}$ is the simulation time, and $\bar{\omega}_{\rm low}$ is a
characteristic frequency for the slow degrees of freedom of the system,
motivated by the fact that heat transport is usually dominated by these slow
processes. We choose $\bar{\omega}_{\rm low}$ as the mean frequency of the
lowest 20 % of the vibrational spectrum as shown in Fig. 9, but emphasize that
the argument is not sensitive to this somewhat arbitrary choice.

Figure 9: Significant portion of the vibrational density of states (VDOS)
[95], and representative low frequency $\bar{\omega}_{\rm low}$ for MgO (upper
panel) and CuI (lower panel). Note the different limits on the x-axis.

Using the effective simulation length as defined in Eq. (29), we see that the
simulation for MgO ($\tilde{t}_{0}=442.2$) is effectively 6.5 times longer
than that for CuI ($\tilde{t}_{0}=67.8$), due to the stiffer bonding and
therefore faster vibrations present in MgO. Noting that lifetimes tend to
decrease with increasing frequency as argued earlier [64], this explains why
the longest lifetimes in the simulation cell are of the same order of
magnitude in MgO and CuI as shown in Fig. 3 and 8, despite the much stronger
anharmonic character of CuI.

For novel materials, we therefore suggest to estimate the necessary simulation
times based on the vibrational spectrum and lifetime estimates according to
Eq. (23). Special care must be taken for slow but harmonic materials, for
which the necessary simulation times might easily be much longer than those
reported here, consistent with GK studies based on empirical force fields for
harmonic materials like silicon [96].

##  VIII Conclusion and outlook

We have presented and applied an implementation of aiGK simulations based on
the first-principles heat flux introduced in Ref. [21]. Systematically
removing noise from the HFACF allows to estimate cutoff times in a numerically
robust way. We also presented a detailed account of our updated size-
extrapolation scheme based on mapping the aiMD trajectories to an (effective)
harmonic model. The scheme was applied to MgO and CuI, two simple binary
systems of quite different harmonic character: MgO is an example for a
stiffly-bonded, harmonic material with quite high thermal conductivity for a
rock salt compound, whereas CuI is a strongly anharmonic compound that
dynamically destabilizes and becomes superionic conducting at higher
temperatures [14, 91]. Good agreement with the available literature is found
for CuI, and for MgO when correcting for isotope effects.

The presented scheme and its implementation in FHI-vibes makes performing,
post-processing, and analyzing aiGK simulations much more straightforward than
previously possible. The reduced human intervention enables to study heat
transport in materials across materials space, in particular for strongly
anharmonic or complex materials, as discussed in Ref. [65]. Furthermore,
approaches such as those based on BTE can be systematically benchmarked
against non-perturbative results obtained from aiGK in the future.

We note in passing that the presented approach can equally be applied to GK
studies based on empirical or machine-learned force fields [97, 98, 99, 100,
101, 102, 103]. Also in these cases, the presented strategies that rely on
physically motivated approaches and parameters can be helpful to obtain more
stable and systematic data from potentially noisy simulations.

## Data and code availability

The presented method is implemented and available in the open source package
FHI-vibes [30]. The package builds on the atomic simulation environment (ASE)
[104]. The data and scripts used to create the plots are made available via
figshare [105]. The raw DFT calculations including input and output files are
accessible via NOMAD. [106].

## Acknowledgments

This project was supported by the NOMAD Center of Excellence (European Union’s
Horizon 2020 research and innovation program, grant agreement No. 951786), the
ERC Advanced Grant TEC1p (European Research Council, grant agreement No.
740233), and the North-German Supercomputing Alliance (HLRN). F. K.
acknowledges support from the Swedish Research Council (VR) program
2020-04630, and the Swedish e-Science Research Centre (SeRC). F. K. would like
to thank Stefano Baroni and Federico Grasselli for inspiration and fruitful
discussions related to GK theory, and Marcel Langer for useful feedback on the
method and manuscript.

## References

  * Snyder and Toberer [2008] G. J. Snyder and E. S. Toberer, Nature Materials 7, 105 (2008).
  * Perepezko [2009] J. H. Perepezko, Science 326, 1068 (2009).
  * Nimmo  _et al._ [2007] F. Nimmo, J. R. Spencer, R. T. Pappalardo, and M. E. Mullen, Nature 447, 289 (2007).
  * Grasselli  _et al._ [2020] F. Grasselli, L. Stixrude, and S. Baroni, Nature Communications 11, 3605 (2020), 2003.12557 .
  * Vogelsang  _et al._ [1987] R. Vogelsang, C. Hoheisel, and G. Ciccotti, The Journal of Chemical Physics 86, 6371 (1987).
  * Esfarjani and Stokes [2008] K. Esfarjani and H. T. Stokes, Physical Review B 77, 144112 (2008).
  * Hellman  _et al._ [2013] O. Hellman, P. Steneteg, I. A. Abrikosov, and S. I. Simak, Physical Review B 87, 104111 (2013), 1303.1145 .
  * Paulatto  _et al._ [2015] L. Paulatto, I. Errea, M. Calandra, and F. Mauri, Physical Review B 91, 054304 (2015), 1411.5628 .
  * Feng  _et al._ [2017] T. Feng, L. Lindsay, and X. Ruan, Physical Review B 96, 161201(R) (2017).
  * Broido  _et al._ [2007] D. A. Broido, M. Malorny, G. Birner, N. Mingo, and D. A. Stewart, Applied Physics Letters 91, 231922 (2007).
  * Simoncelli  _et al._ [2019] M. Simoncelli, N. Marzari, and F. Mauri, Nature Physics 15, 809 (2019), 1901.01964 .
  * Simoncelli  _et al._ [2022] M. Simoncelli, N. Marzari, and F. Mauri, Physical Review X 12, 041011 (2022), 2112.06897 .
  * Clarke  _et al._ [2012] D. R. Clarke, M. Oechsner, and N. P. Padture, MRS Bulletin 37, 891 (2012).
  * Knoop  _et al._ [2020a] F. Knoop, T. A. R. Purcell, M. Scheffler, and C. Carbogno, Physical Review Materials 4, 083809 (2020a), 2006.14672 .
  * Xia  _et al._ [2020] Y. Xia, V. I. Hegde, K. Pal, X. Hua, D. Gaines, S. Patel, J. He, M. Aykol, and C. Wolverton, Physical Review X 10, 041029 (2020).
  * Ioffe and Regel [1960] A. Ioffe and A. Regel, Prog. Semicond 4, 237 (1960).
  * Green [1954] M. S. Green, The Journal of Chemical Physics 22, 398 (1954).
  * Kubo [1957] R. Kubo, Journal of the Physical Society of Japan 12, 570 (1957).
  * Kubo  _et al._ [1957] R. Kubo, M. Yokota, and S. Nakajima, Journal of the Physical Society of Japan 12, 1203 (1957).
  * Caldarelli  _et al._ [2022] G. Caldarelli, M. Simoncelli, N. Marzari, F. Mauri, and L. Benfatto, Physical Review B 106, 024312 (2022), 2202.02246 .
  * Carbogno  _et al._ [2017] C. Carbogno, R. Ramprasad, and M. Scheffler, Physical Review Letters 118, 175901 (2017), 1608.06917 .
  * Marcolongo  _et al._ [2016] A. Marcolongo, P. Umari, and S. Baroni, Nature Physics 12, 80 (2016).
  * Sutherland  _et al._ [2021] B. J. Sutherland, W. H. D. Moore, and D. E. Manolopoulos, The Journal of Chemical Physics 154, 174104 (2021), 2104.09566 .
  * Chun  _et al._ [2021] H.-M. Chun, Q. Gao, and J. M. Horowitz, Physical Review Research 3, 043172 (2021), 2103.09288 .
  * Schelling  _et al._ [2002] P. K. Schelling, S. R. Phillpot, and P. Keblinski, Physical Review B 65, 144306 (2002).
  * Sellan  _et al._ [2010] D. P. Sellan, E. S. Landry, J. E. Turney, A. J. H. McGaughey, and C. H. Amon, Physical Review B 81, 214305 (2010).
  * Allen and Feldman [1993] P. B. Allen and J. L. Feldman, Physical Review B 48, 12581 (1993).
  * Marcolongo  _et al._ [2020] A. Marcolongo, L. Ercole, and S. Baroni, Journal of Chemical Theory and Computation 16, 3352 (2020).
  * Ercole  _et al._ [2017] L. Ercole, A. Marcolongo, and S. Baroni, Scientific Reports 7, 15835 (2017), 1706.01381 .
  * Knoop  _et al._ [2020b] F. Knoop, T. Purcell, M. Scheffler, and C. Carbogno, Journal of Open Source Software 5, 2671 (2020b).
  * Togo  _et al._ [2015] A. Togo, L. Chaput, and I. Tanaka, Physical Review B 91, 094306 (2015), 1501.00691 .
  * deKoker [2009] N. deKoker, Physical Review Letters 103, 125902 (2009).
  * Stackhouse  _et al._ [2010] S. Stackhouse, L. Stixrude, and B. B. Karki, Physical Review Letters 104, 208501 (2010).
  * Tang and Dong [2010] X. Tang and J. Dong, Proceedings of the National Academy of Sciences 107, 4539 (2010).
  * Chen  _et al._ [2010a] J. Chen, G. Zhang, and B. Li, Physics Letters A 374, 2392 (2010a), 1005.4478 .
  * Ercole  _et al._ [2016] L. Ercole, A. Marcolongo, P. Umari, and S. Baroni, Journal of Low Temperature Physics 185, 79 (2016), 1601.00158 .
  * Perdew  _et al._ [2008] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Physical Review Letters 100, 136406 (2008), 0711.0156 .
  * Blum  _et al._ [2009] V. Blum, R. Gehrke, F. Hanke, P. Havu, V. Havu, X. Ren, K. Reuter, and M. Scheffler, Computer Physics Communications 180, 2175 (2009).
  * Hellman  _et al._ [2011] O. Hellman, I. A. Abrikosov, and S. I. Simak, Physical Review B 84, 180301(R) (2011), 1103.5590 .
  * West and Estreicher [2006] D. West and S. K. Estreicher, Physical Review Letters 96, 115504 (2006).
  * Togo and Tanaka [2015] A. Togo and I. Tanaka, Scripta Materialia 108, 1 (2015).
  * Knoop [2022] F. Knoop, “Heat transport in strongly anharmonic solids from first principles,”  (2022).
  * Baroni  _et al._ [2018] S. Baroni, R. Bertossa, L. Ercole, F. Grasselli, and A. Marcolongo, arXiv , 1 (2018), 1802.08006 .
  * Zwanzig and Ailawadi [1969] R. Zwanzig and N. K. Ailawadi, Physical Review 182, 280 (1969).
  * Jones and Mandadapu [2012a] R. E. Jones and K. K. Mandadapu, The Journal of Chemical Physics 136, 154102 (2012a).
  * Oliveira and Greaney [2017] L. S. Oliveira and P. A. Greaney, Physical Review E 95, 023308 (2017).
  * Wang  _et al._ [2017] Z. Wang, S. Safarkhani, G. Lin, and X. Ruan, International Journal of Heat and Mass Transfer 112, 267 (2017).
  * He  _et al._ [2012] Y. He, I. Savić, D. Donadio, and G. Galli, Physical Chemistry Chemical Physics 14, 16209 (2012).
  * Knuth  _et al._ [2015] F. Knuth, C. Carbogno, V. Atalla, V. Blum, and M. Scheffler, Computer Physics Communications 190, 33 (2015).
  * Ercole [2018] L. Ercole, “Ab initio Simulation of Heat Transport in Silica Glass,” (2018).
  * Isaeva  _et al._ [2019] L. Isaeva, G. Barbalinardo, D. Donadio, and S. Baroni, Nature Communications 10, 3853 (2019), 1904.02255 .
  * Fan  _et al._ [2015] Z. Fan, L. F. C. Pereira, H.-Q. Wang, J.-C. Zheng, D. Donadio, and A. Harju, Physical Review B 92, 094301 (2015), 1503.06565 .
  * Pereverzev and Sewell [2022] A. Pereverzev and T. Sewell, International Journal of Heat and Mass Transfer 188, 122647 (2022).
  * Coretti  _et al._ [2018] A. Coretti, S. Bonella, L. Rondoni, and G. Ciccotti, Molecular Physics 116, 1 (2018).
  * Chen  _et al._ [2010b] J. Chen, G. Zhang, and B. Li, Physics Letters A 374, 2392 (2010b).
  * Oliveira and Greaney [2016] L. S. Oliveira and P. A. Greaney, Physical review. E 95, 023308 (2016).
  * Jones and Mandadapu [2012b] R. E. Jones and K. K. Mandadapu, The Journal of Chemical Physics 136, 154102 (2012b).
  * Born and Huang [1954] M. Born and K. Huang,  _Dynamical theory of crystal lattices_ (Clarendon Press, 1954).
  * Hardy [1963] R. J. Hardy, Physical Review 132, 168 (1963).
  * Negele and Orland [1988] J. W. Negele and H. Orland,  _Quantum many-particle systems_ (Addison-Wesley, 1988).
  * Togo and Tanaka [2018] A. Togo and I. Tanaka, arXiv (2018), 1808.01590 .
  * Maradudin and Vosko [1968] A. A. Maradudin and S. H. Vosko, Reviews of Modern Physics 40, 1 (1968).
  * Pomeranchuk [1941] I. Pomeranchuk, Physical Review 60, 820 (1941).
  * Herring [1954] C. Herring, Physical Review 95, 954 (1954).
  * Knoop  _et al._ [2022] F. Knoop, T. A. R. Purcell, M. Scheffler, and C. Carbogno, arXiv (2022), 2209.12720 .
  * Monkhorst and Pack [1976] H. J. Monkhorst and J. D. Pack, Physical Review B 13, 5188 (1976).
  * Puligheddu  _et al._ [2019] M. Puligheddu, Y. Xia, M. Chan, and G. Galli, Physical Review Materials 3, 085401 (2019), 1902.08260 .
  * Sheng  _et al._ [1994] P. Sheng, M. Zhou, and Z.-Q. Zhang, Physical Review Letters 72, 234 (1994).
  * Parlinski  _et al._ [1997] K. Parlinski, Z. Q. Li, and Y. Kawazoe, Physical Review Letters 78, 4063 (1997).
  * Carreras  _et al._ [2017] A. Carreras, A. Togo, and I. Tanaka, Computer Physics Communications 221, 221 (2017), 1708.03435 .
  * Tadano and Tsuneyuki [2015] T. Tadano and S. Tsuneyuki, Physical Review B 92, 054301 (2015), 1506.01781 .
  * Carbogno  _et al._ [2014] C. Carbogno, C. G. Levi, C. G. VandeWalle, and M. Scheffler, Physical Review B 90, 144109 (2014).
  * Ladd  _et al._ [1986] A. J. C. Ladd, B. Moran, and W. G. Hoover, Physical Review B 34, 5058 (1986).
  * Kaburaki  _et al._ [1998] H. Kaburaki, J. Li, and S. Yip, MRS Proceedings 538, 503 (1998).
  * McGaughey and Kaviany [2004] A. McGaughey and M. Kaviany, International Journal of Heat and Mass Transfer 47, 1783 (2004).
  * Slack [1962] G. A. Slack, Physical Review 126, 427 (1962).
  * Touloukian  _et al._ [1 01] Y. S. Touloukian, R. W. Powell, C. Y. Ho, and P. G. Klemens, (1971-01).
  * MacPherson and Schloessin [1983] W. MacPherson and H. Schloessin, High Temperatures. High Pressures (Print) 15, 495 (1983).
  * Andersson and Bäckström [1986] S. Andersson and G. Bäckström, Review of Scientific Instruments 57, 1633 (1986).
  * Katsura [1997] T. Katsura, Physics of the Earth and Planetary Interiors 101, 73 (1997).
  * Dalton  _et al._ [2013] D. A. Dalton, W.-P. Hsieh, G. T. Hohensee, D. G. Cahill, and A. F. Goncharov, Scientific Reports 3, 2400 (2013).
  * Hofmeister [2014] A. M. Hofmeister, Physics and Chemistry of Minerals 41, 361 (2014).
  * Koker [2010] N. d. Koker, Earth and Planetary Science Letters 292, 392 (2010).
  * Tse  _et al._ [2018] J. S. Tse, N. J. English, K. Yin, and T. Iitaka, The Journal of Physical Chemistry C 122, 10682 (2018).
  * Dekura and Tsuchiya [2017] H. Dekura and T. Tsuchiya, Physical Review B 95, 184303 (2017).
  * Plata  _et al._ [2017] J. J. Plata, P. Nath, D. Usanmaz, J. Carrete, C. Toher, M. d. Jong, M. Asta, M. Fornari, M. B. Nardelli, and S. Curtarolo, npj Computational Materials 3, 45 (2017).
  * Haigis  _et al._ [2012] V. Haigis, M. Salanne, and S. Jahn, Earth and Planetary Science Letters 355, 102 (2012), 1210.7456 .
  * Ceperley and Alder [1980] D. M. Ceperley and B. J. Alder, Physical Review Letters 45, 566 (1980).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Physical Review Letters 77, 3865 (1996).
  * Kinaci  _et al._ [2012] A. Kinaci, J. B. Haskins, and T. Çağın, The Journal of Chemical Physics 137, 014106 (2012), 1201.6684 .
  * Boyce  _et al._ [1980] J. Boyce, T. Hayes, J. Mikkelsen, and W. Stutius, Solid State Communications 33, 183 (1980).
  * Boyce  _et al._ [1981] J. B. Boyce, T. M. Hayes, and J. C. Mikkelsen, Physical Review B 23, 2876 (1981).
  * Perry [2016] D. L. Perry,  _Handbook of Inorganic Compounds_, zeroth ed. (CRC Press, 2016).
  * Yang  _et al._ [7 12] C. Yang, D. Souchay, M. Kneiß, M. Bogner, H. M. Wei, M. Lorenz, O. Oeckler, G. Benstetter, Y. Q. Fu, and M. Grundmann, Nature Communications 8, 16076 (2017-12).
  * Dove [1993] M. T. Dove,  _Introduction to lattice dynamics_ (Cambridge University Press, 1993).
  * Volz and Chen [2000] S. G. Volz and G. Chen, Physical Review B 61, 2651 (2000).
  * Sosso  _et al._ [2012] G. C. Sosso, D. Donadio, S. Caravati, J. Behler, and M. Bernasconi, Physical Review B 86, 104301 (2012).
  * Korotaev  _et al._ [2019] P. Korotaev, I. Novoselov, A. Yanilkin, and A. Shapeev, Physical Review B 100, 144308 (2019).
  * Qian  _et al._ [2019] X. Qian, S. Peng, X. Li, Y. Wei, and R. Yang, Materials Today Physics 10, 100140 (2019).
  * Li  _et al._ [2020] R. Li, E. Lee, and T. Luo, Materials Today Physics 12, 100181 (2020).
  * Fan  _et al._ [2021] Z. Fan, Z. Zeng, C. Zhang, Y. Wang, K. Song, H. Dong, Y. Chen, and T. Ala-Nissila, Physical Review B 104, 104309 (2021), 2107.08119 .
  * Verdi  _et al._ [2021] C. Verdi, F. Karsai, P. Liu, R. Jinnouchi, and G. Kresse, npj Computational Materials 7, 156 (2021).
  * Langer  _et al._ [2022] M. Langer, F. Knoop, C. Carbogno, M. Scheffler, and M. Rupp, “Thermal conductivities with message-passing neural networks,” (2023), in preparation.
  * Larsen  _et al._ [2017] A. H. Larsen, J. J. Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dułak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C. Jennings, P. B. Jensen, J. Kermode, J. R. Kitchin, E. L. Kolsbjerg, J. Kubal, K. Kaasbjerg, S. Lysgaard, J. B. Maronsson, T. Maxson, T. Olsen, L. Pastewka, A. Peterson, C. Rostgaard, J. Schiøtz, O. Schütt, M. Strange, K. S. Thygesen, T. Vegge, L. Vilhelmsen, M. Walter, Z. Zeng, and K. W. Jacobsen, Journal of Physics: Condensed Matter 29, 273002 (2017).
  * [105] “https://dx.doi.org/10.6084/m9.figshare.20146847,” .
  * [106] “https://dx.doi.org/10.17172/NOMAD/2022.06.24-1,” .