# Atomic forces by quantum Monte Carlo: Application to phonon dispersion
calculations

Kousuke Nakano1,2<EMAIL_ADDRESS>Tommaso Morresi3 Michele Casula3 Ryo
Maezono2 Sandro Sorella1<EMAIL_ADDRESS>1 International School for Advanced
Studies (SISSA), Via Bonomea 265, 34136, Trieste, Italy  2 Japan Advanced
Institute of Science and Technology (JAIST), Asahidai 1-1, Nomi, Ishikawa
923-1292, Japan  3 Institut de Minéralogie, de Physique des Matériaux et de
Cosmochimie (IMPMC), Sorbonne Université, CNRS UMR 7590, IRD UMR 206, MNHN, 4
Place Jussieu, 75252 Paris, France

###### Abstract

We report a successful application of the ab initio quantum Monte Carlo (QMC)
framework to a phonon dispersion calculation. A full phonon dispersion of
diamond is successfully calculated at the variational Monte Carlo (VMC) level,
based on the frozen-phonon technique. The VMC-phonon dispersion is in good
agreement with the experimental results, giving renormalized harmonic optical
frequencies very close to the experimental values, and improving upon previous
density functional theory estimates. Key to success for the QMC approach is
the statistical error reduction in the atomic force evaluation. We show that
this can be achieved by using well conditioned atomic basis sets and by
explicitly removing the basis-set redundancy, which reduces the statistical
error of forces by up to two orders of magnitude by combining it with the so-
called space-warp transformation algorithm. This leads to affordable and
accurate QMC-phonons calculations, which are up to $10^{4}$ times more
efficient than a bare force treatment, and paves the way to new applications,
particularly in correlated materials, where phonons have been poorly
reproduced so far.

The accurate description of phonons in a solid is one of the central research
topics in the field of condensed matter physics and materials science for
discussing phase stability (i.e., Gibbs-free energy), electron-phonon
interaction, and structural phase transitions of materials Martin (2004);
Sholl and Steckel (2011). Ab initio phonon calculations based on the Density
Functional Theory (DFT) Baroni  _et al._ (2001); Togo and Tanaka (2015) have
been successful for many compounds, but they often fail in strongly-correlated
materials. For example, DFT calculations severely underestimate the highest
frequency of the optical phonons of graphene at the $K$ point Mounet and
Marzari (2005); Lazzeri  _et al._ (2008), because the electron-electron
correlation is not taken into account with sufficient accuracy. Another
example is the elemental cerium, whose phonon dispersions measured by neutron
scattering strongly mismatch with the calculated DFT-PBE ones Krisch  _et al._
(2011). Interestingly, such failure was also seen in a high-$T_{c}$ cuprate
superconductor Reznik  _et al._ (2008) Some effort has been made to include
correlation in phonon calculations in the DFT+U framework Floris  _et al._
(2011) and also within the DMFT framework Leonov  _et al._ (2012, 2014). In
both cases, this requires modeling correlations by an empirical parameter,
though physically motivated (i.e., the Hubbard $U$). Indeed, a genuine ab
initio framework applicable to strongly correlated materials without any
empirical parameters remains, so far, a very important theoretical challenge.

The ab initio quantum Monte Carlo (QMC) framework, including variational
quantum Monte Carlo (VMC) and the diffusion quantum Monte Carlo (DMC) schemes,
is among the state-of-art numerical methods to obtain highly accurate many-
body wave functions Foulkes  _et al._ (2001), and cope with the electron
correlation more rigorously than DFT. It has been successfully applied to
challenging materials that DFT cannot tackle, such as cuprates Marchi  _et
al._ (2011), iron-arsenides Casula and Sorella (2013), and graphene Marchi
_et al._ (2011); Sorella  _et al._ (2018). So far, unfortunately, almost all
QMC applications are mainly based on energy and its first derivative (i.e.,
atomic force) calculations. Indeed, it is at present an open problem how to
evaluate, with an affordable computational effort, second and higher-order
derivatives, which are essential for computing various physical properties.

There are three routes to compute the second derivatives (i.e.,
$\frac{\partial^{2}E}{\partial{R_{\alpha}}\partial{R_{\beta}}}$), which are
needed for evaluating harmonic phonon properties, by ab initio calculations,
i.e., $\rm(\hskip 1.79993pti\hskip 1.79993pt)$ potential energy surface (PES)
fitting, $\rm(\hskip 0.80002ptii\hskip 0.80002pt)$ finite difference
expression based on atomic force evaluations , and $\rm(iii)$ direct
evaluation of second derivatives. All of the above attempts have been
successful for isolated molecular systems Zen  _et al._ (2012); Luo  _et al._
(2014); Liu  _et al._ (2019); Nakano  _et al._ (2019). On the other hand, for
solids, only strategy $\rm(\hskip 1.79993pti\hskip 1.79993pt)$ has been
successful so far within the QMC framework. For example, Maezono et al.
calculated Raman frequencies of diamond (phonons at $q=\Gamma$ point) Maezono
_et al._ (2007). However, QMC-phonon calculations of solids have been limited
to a single high-symmetry $q$-point. To the best of our knowledge, the full
($q$-resolved) phonon dispersion has been unaccessible so far at both VMC and
DMC levels.

In this paper, we report a successful phonon dispersion calculation of diamond
at the VMC level by adopting strategy $\rm(\hskip 0.80002ptii\hskip
0.80002pt)$ , the so-called frozen phonon technique Togo and Tanaka (2015).
The key to success is to reduce the statistical error of atomic forces. We
found that removing the nearly linear dependency of the basis set used for the
trial wave function parametrization Azadi  _et al._ (2010) dramatically lowers
the statistical error of forces. Its decrease reaches the order of $\times$
${\sim 10^{-2}}$, which corresponds to a speed-up of ${\sim 10^{4}}$ times in
a VMC computation. This drastic reduction enables us to construct a dynamical
matrix within an affordable computational time, and to eventually apply VMC-
phonon calculations to new interesting materials from first principles.

Table 1:  The highest harmonic phonon frequencies (THz) of diamond at high-symmetry $q$ points. All phonons are calculated using the experimental lattice parameter. In this table, (P) denotes the interpolated frequencies obtained by Phonopy, (F) denotes the phonon frequency at $\Gamma$ obtained by fitting forces of distorted structures along the Raman mode with a linear function at the VMC level, and (E) denotes the phonon frequency at $\Gamma$ obtained by fitting energies of undistorted and distorted structures along the Raman mode with a quadratic function at the VMC and LRDMC levels. The last column indicates the harmonic frequencies estimated by the raw experimental values: 39.938 THzLiu  _et al._ (2000), 35.299 THzSchwoerer-Böhning  _et al._ (1998), and 37.962 THzKulda  _et al._ (2002) at the $\Gamma$, $X$, and $L$ points, respectively. $q$ | DFT | Previous work | This work444These values include the one-body finite size corrections, i.e., $-$ 0.16 THz, $-$ 0.18 THz, and $-$ 0.23 THz for $q$ = $\Gamma$, $X$, and $L$, respectively. | Experiment
---|---|---|---|---
LDA-PZ | GGA-PBE | VMC333These values are taken from Ref. Maezono  _et al._ , 2007. | DMC333These values are taken from Ref. Maezono  _et al._ , 2007. | VMC(P) | VMC(F) | VMC(E) | LRDMC(E) | Harmonic (Estimated)
$\Gamma$ | 38.5511138.40 THz in the previous DFT study employing the LDA functional. See Ref. Maezono  _et al._ , 2007. | 38.8222238.73 THz in the previous DFT study employing the GGA-PBE functional. See Ref. Maezono  _et al._ , 2007. | 41.64(9) | 41.22(12) | 40.65(38) | 40.49(4) | 40.68(29) | 41.52(22) | 40.460555The anharmonic renormalization, $-$ 17.4 cm-1 = $-$ 0.522 THz was employed. See Ref. Vanderbilt  _et al._ , 1984., 40.349666The anharmonic renormalizations are $-$ 0.411 THz, $-$ 0.177 THz, and $-$ 0.280 THz for $\Gamma$, $X$, and $L$, respectively, which were estimated by molecular dynamics simulations performed in this work.
$X$ | 35.64 | 35.87 | - | - | 36.48(40) | - | - | - | 35.476666The anharmonic renormalizations are $-$ 0.411 THz, $-$ 0.177 THz, and $-$ 0.280 THz for $\Gamma$, $X$, and $L$, respectively, which were estimated by molecular dynamics simulations performed in this work.
$L$ | 37.31 | 37.47 | - | - | 38.01(31) | - | - | - | 38.242666The anharmonic renormalizations are $-$ 0.411 THz, $-$ 0.177 THz, and $-$ 0.280 THz for $\Gamma$, $X$, and $L$, respectively, which were estimated by molecular dynamics simulations performed in this work.
|  |  |  |  |  |  |  |  |
Figure 1: The phonon dispersion and Raman frequencies of diamond calculated
using a 2 $\times$ 2 $\times$ 2 conventional supercell at the VMC level, where
the anharmonic and one-body finite size corrections are included. The
experimental lattice parameter (i.e., 6.741 bohr Parrish (1960)) was employed.
The error bars were estimated by the jackknife method. Observed experimental
results are also plotted for comparison. Experimental data are taken from
Refs. Warren  _et al._ , 1965; *1967WAR; Liu  _et al._ , 2000; Kulda  _et al._
, 2002; Schwoerer-Böhning  _et al._ , 1998. The data in Ref. Schwoerer-Böhning
_et al._ , 1998 were digitalized using WebPlotDigitizer Rohatgi (2020).

All Variational Monte Carlo (VMC) and lattice regularized diffusion Monte
Carlo (LRDMC) Casula  _et al._ (2005); Nakano  _et al._ (2020a) calculations
in this study were performed by the TurboRVB Casula and Sorella (2003); Nakano
_et al._ (2020b) SISSA quantum Monte Carlo package. We employed the Jastrow
Slater-determinant(JSD) ansatz, i.e, $\Psi=\Psi_{\text{SD}}\times\exp J,$
where $\Psi_{\text{SD}}$ and $J$ are the Slater determinant and Jastrow terms,
respectively. The Slater determinant part is expressed in terms of molecular
orbitals $\phi_{k}(\mathbf{r})=\sum_{i=1}^{L}{c_{i,k}{\psi_{i}}(\mathbf{r})}$
expanded in a periodized Gaussian basis set $\\{\psi_{i}\\}$. The valence
triple-zeta (VTZ) basis set accompanying an energy-consistent pseudo potential
developed by Burkatzki et al. Burkatzki  _et al._ (2007) was employed for the
primitive Gaussian atomic orbitals (Table S-Ifoo ). The coefficients of atomic
orbitals (i.e., $c_{i,k}$) were obtained by a DFT calculation with the LDA-PZ
exchange-correlation functional Perdew and Zunger (1981) and were left
unchanged during the VMC optimization. The Jastrow factor was composed of
inhomogeneous one-, two- and three-body contributions
($J={J_{1}}^{inh}+{J_{2}}+{J_{3}}$). foo  The variational parameters in the
Jastrow terms were optimized by the stochastic reconfiguration Sorella  _et
al._ (2007) and/or the modified linear method Umrigar  _et al._ (2007); Nakano
_et al._ (2020b) implemented in TurboRVB. LRDMC calculations were performed by
the single-grid scheme Casula  _et al._ (2005) with a lattice space, $a$ = 0.2
bohr.

In this paper, we focus on diamond (Space group: $Fd\bar{3}m$) as a proof of
concept for the first example of a VMC-based phonon dispersion calculation. 2
$\times$ 2 $\times$ 2 conventional supercells (256 electrons / 64 carbon atoms
in a simulation cell) were used for most calculations, while 3 $\times$ 3
$\times$ 3 conventional supercells (864 electrons / 216 carbon atoms in a
simulation cell) were also used for several calculations to investigate the
finite size errors. $L$-twist (i.e., $k$ = $\pi$, $\pi$, $\pi$) was employed
for alleviating the so-called one-body finite-size effects Maezono  _et al._
(2007); Hennig  _et al._ (2010); Sorella  _et al._ (2011). Dynamical matrices
and the corresponding phonon dispersions were calculated based on the frozen-
phonon technique implemented in the Phonopy module Togo and Tanaka (2015),
where a 0.15 bohr displacement of carbon atoms was large enough to work with
reasonable signal/noise ratio in QMC forces. This displacement underestimates
harmonic frequencies only by $\sim$ 0.1 THz, as shown in Fig. S-3 foo . Error
bars in a phonon dispersion were estimated by the jackknife method Becca and
Sorella (2017). Phonon calculations based on DFT were performed using the
Quantum Espresso package Giannozzi  _et al._ (2009) with LDA-PZ Perdew and
Zunger (1981) and GGA-PBE Perdew  _et al._ (1996) exchange-correlation
functionals at the experimental lattice parameter. foo  The phonon dispersion
of diamond has already been studied using DFT calculations by many groups so
far, at the theoretical Pavone  _et al._ (1993); Kresse  _et al._ (1995) and
experimental Maezono  _et al._ (2007) lattice parameters, that makes this
system a very good testbed for any new methodological implementation of
phonons calculations. As shown later, our DFT calculations are consistent with
the previous study.

Table 2:  Vinet EOS parameters. Zero point energy (ZPE) and temperature effects (TE) are not included in these theoretical data. Parameter | DFT | QMC | Experiment
---|---|---|---
LDA-PZ | GGA-PBE | VMC111These values are taken from Ref. Maezono  _et al._ , 2007. Here, ZPE and TE are subtracted. | DMC111These values are taken from Ref. Maezono  _et al._ , 2007. Here, ZPE and TE are subtracted. | VMC | LRDMC | w/o ZPE, w/o TE | Observed
Lattice (Bohr) | 6.683 | 6.748 | 6.691(4) | 6.734(4) | 6.693(1) | 6.702(1) | 6.7193(5)222ZPE and TE are corrected, $-$ 0.37 Bohr3, $+$ 11 GPa, and $-$ 0.03 for $V_{0}$, $B_{0}$, and $B^{\prime}_{0}$, respectively. See Ref. Maezono  _et al._ , 2007. | 6.7410(5)333These values are taken from Ref. Occelli  _et al._ , 2003.
$V_{0}$ (bohr3) | 37.32 | 38.41 | 37.45(6) | 38.17(6) | 37.47(2) | 37.63(2) | 37.920(9)222ZPE and TE are corrected, $-$ 0.37 Bohr3, $+$ 11 GPa, and $-$ 0.03 for $V_{0}$, $B_{0}$, and $B^{\prime}_{0}$, respectively. See Ref. Maezono  _et al._ , 2007. | 38.290(9)333These values are taken from Ref. Occelli  _et al._ , 2003.
$B_{0}$ (GPa) | 465 | 433 | 483(4) | 448(3) | 476(6) | 463(5) | 457(1)222ZPE and TE are corrected, $-$ 0.37 Bohr3, $+$ 11 GPa, and $-$ 0.03 for $V_{0}$, $B_{0}$, and $B^{\prime}_{0}$, respectively. See Ref. Maezono  _et al._ , 2007., 453(5)222ZPE and TE are corrected, $-$ 0.37 Bohr3, $+$ 11 GPa, and $-$ 0.03 for $V_{0}$, $B_{0}$, and $B^{\prime}_{0}$, respectively. See Ref. Maezono  _et al._ , 2007. | 446(1)333These values are taken from Ref. Occelli  _et al._ , 2003., 442(5)444These values are taken from Refs. McSkimin and Andreatch Jr, 1972; Occelli  _et al._ , 2003.
$B^{\prime}_{0}$ | 3.65 | 3.70 | 3.8(1) | 3.7(1) | 4.0(6) | 4.9(6) | 3.0(1)222ZPE and TE are corrected, $-$ 0.37 Bohr3, $+$ 11 GPa, and $-$ 0.03 for $V_{0}$, $B_{0}$, and $B^{\prime}_{0}$, respectively. See Ref. Maezono  _et al._ , 2007., 4.0(7)222ZPE and TE are corrected, $-$ 0.37 Bohr3, $+$ 11 GPa, and $-$ 0.03 for $V_{0}$, $B_{0}$, and $B^{\prime}_{0}$, respectively. See Ref. Maezono  _et al._ , 2007. | 3.0(1)333These values are taken from Ref. Occelli  _et al._ , 2003., 4.0(7)444These values are taken from Refs. McSkimin and Andreatch Jr, 1972; Occelli  _et al._ , 2003.
|  |  |  |  |  |  |  |
Figure 2: Equation of states and the curves fitted by the Vinet exponential
function Vinet  _et al._ (1987). VMC and DMC results (Maezono et al.) are
taken from Ref. Maezono  _et al._ , 2007.  Figure 3: (a) VMC energies and
their statistical error bars and (b) VMC forces and their statistical error
bars vs. inverse of the condition number of the overlap matrix. The numerical
data is shown in Table S-II foo . The forces were calculated with and without
the space warp coordinate transformations (SWCT) Sorella and Capriotti (2010).
The Hellmann-Feynman and Pulay contributions to the forces are shown in Fig.
S-7 foo . The same VTZ atomic basis set and the pseudo potential Burkatzki
_et al._ (2007) as in the phonon and EOS calculations were adopted, while 1
$\times$ 1 $\times$ 1 conventional cell (8 carbon atoms in a simulation cell)
with $k$ = $\Gamma$ twist was employed. Only one carbon atom was displaced in
$x$ direction by 0.15 bohr (Table S-III) foo .

Fig. 1 shows the phonon dispersion obtained by our VMC calculations using the
conventional 2 $\times$ 2 $\times$ 2 diamond supercell 111 2 $\times$ 2
$\times$ 2 conventional supercell containing 64 atoms in the simulation cell
is large enough for obtaining a phonon dispersion almost consistent with
experimental ones, as shown in Fig. S-2 foo . This was confirmed by comparing
the phonon dispersion obtained by the finite displacement method and that
obtained by the linear-response method with very dense $k$ and $q$ grids at
the DFT level. . Observed experimental frequencies Warren  _et al._ (1967);
Liu  _et al._ (2000); Kulda  _et al._ (2002) are also plotted for comparison.
A phonon dispersion obtained by the finite displacement method does not
include anharmonic effects, which can decrease harmonic frequencies by up to
$\sim$ 5-10% for the lightest elements Harris (1995). Therefore, for
comparison with the experimental results, anharmonic corrections were added to
the VMC phonon dispersion in this study. The phonon dispersion before the
correction is shown in Fig. S-1 foo  The anharmonic renormalizations were
estimated in this study using path integral molecular dynamics
simulationsMouhat  _et al._ (2017); Morresi  _et al._ at the DFT level with
the PBE exchange-correlation function (Fig. S-5 foo ), giving $-$0.411 THz,
$-$0.177 THz, and $-$0.280 THz for $\Gamma$, $X$, and $L$, respectively (see
the SI for details). The value at $\Gamma$ is consistent with a reported
estimate of $-$17.4 cm-1 = $-$0.522 THz Vanderbilt  _et al._ (1984). Notice
that other possible sources of error were also considered. The phonon
dispersion in Fig. 1 also includes the one-body finite size corrections that
were estimated by DFT calculations (see Fig. S-2 foo ) The two-body finite-
size error is negligible because the average density and the volume do not
change in phonon calculations. Table 1 shows a detailed comparison of the
highest harmonic phonon frequencies at three $q$ points, i.e., $\Gamma$ =
(0,0,0), $X$ = (2$\pi$,0,0), and $L$ = ($\pi$,$\pi$,$\pi$). Raman frequencies
at $\Gamma$ obtained by a direct fit of VMC energies, VMC forces, and LRDMC
energies of the structures displaced along the corresponding eigenmode are
also plotted in Fig. 1. 222 In the Raman mode, two nearest-neighbor carbon
atoms are displaced in opposite directions by a distance $u$ from their high-
symmetry positions Maezono  _et al._ (2007). The distortion in the [1 0 0]
direction was employed in this study Vanderbilt  _et al._ (1984). We
calculated VMC forces for two distorted structures ($u=0.03$ Bohr and $u=0.05$
Bohr) and obtained the Raman frequency by fitting the forces with a linear
function. We also calculated VMC and LRDMC energies for the undistorted
structure and two distorted structures ($u=0.03$ Bohr and $u=0.05$ Bohr), then
obtained Raman frequencies by fitting the energies with a quadratic function.
In Table 1 and hereafter, (P) denotes the interpolated frequencies obtained by
Phonopy, (F) denotes the phonon frequency at $\Gamma$ obtained by force
fitting, and (E) denotes the phonon frequency at $\Gamma$ obtained by energy
fitting. The corresponding Raman frequencies are consistent within the error
bars, indicating that the Slater determinant obtained by DFT is almost optimal
also in the presence of the Jastrow factor, thus explaining why this
consistency is satisfied quite accurately Moroni  _et al._ (2014). In other
words, if the wavefunction is at its minimum, the consistency is a
consequence. If it is not at its minimum, the consistency may also be
satisfied by chance or good behavior of the used basis set.

All VMC(P), VMC(F), and VMC(D) calculations give the harmonic phonon frequency
at $\Gamma$ very close to the experimental value, considering the
renormalization of the anharmonic effect, i.e., the discrepancy is just $\sim$
0.3 THz. On the other hand, both DFT calculations with LDA-PZ and GGA-PBE
exchange-correlation functionals underestimate the highest frequency at
$\Gamma$ by $\sim$ 1.8 THz and $\sim$ 1.5 THz, respectively. Table 1 shows
that our VMC phonon dispersion calculation also gives accurate harmonic
frequencies at other $q$ points, i.e., at $X$ and $L$. Compared with the
previous VMC study Maezono  _et al._ (2007), our VMC frequency at $\Gamma$
point is closer to the experimental value, as reported in Table 1. 333Since
smaller displacements were employed in this study (i.e., 0.03 and 0.05 Bohr)
than the previous one, the error bars in phonon frequencies at $\Gamma$ point
are a bit larger. The improvement at the VMC level certainly derives from the
more flexible Jastrow factor employed in this study, while the one used in the
previous VMC study was much simpler Drummond  _et al._ (2004); Maezono  _et
al._ (2007).

It is intriguing that our LRDMC calculation gives a 1.2(2) THz higher Raman
frequency than the experimental one, as shown in Table 1. The previous DMC
study also overestimated the Raman frequency by 0.9(1) THz Maezono  _et al._
(2007). To discuss the origin of the discrepancy, we investigated the effect
of the lattice-space error in our LRDMC calculations. An extrapolation ($a\to
0$) with four lattice spaces (i.e., 0.20, 0.30, 0.40, and 0.50 bohr) yields
41.89(44) THz, suggesting that the lattice-space error is not the origin of
the discrepancy. We also suspected that the experimental lattice parameter
employed in the phonon calculation could be significantly different from the
theoretical one, but this is also not the origin, as shown later. Therefore,
the discrepancy should arise from the fixed-node error and this should be
alleviated by a nodal surface optimization Nakano  _et al._ (2020b), which
however is prohibitive in the 2 $\times$ 2 $\times$ 2 conventional supercells
due to the large number of variational parameters of the distorted structures.
A possible future work to study the nodal surface effect is the use of various
exchange-correlation functionals suitable for solids such as HSE06 Krukau  _et
al._ (2006) and SCAN Sun  _et al._ (2015)

Since the equilibrium lattice parameter also affects phonon frequencies, we
investigated the equation of states (EOSs) of diamond. Figure 2 shows plots of
internal energies vs. volumes fitted by the Vinet EOS foo  Previous VMC and
DMC results Maezono  _et al._ (2007), and the experimentally observed
equilibrium lattice parameter are also plotted in Fig. 2 and summarized in
Table 2. In Fig. 2, the zero point energy (ZPE) contributionMaezono  _et al._
(2007); Hao  _et al._ (2012) is subtracted, to make the comparison possible
with internal energies computed at $T=0$ and on static lattice. Table 2 shows
that our VMC calculation reproduces the previous VMC study, while our LRDMC
calculation gives a slightly smaller lattice parameter (6.702(1) bohr) than
the previous DMC study (6.734(4) bohr). This discrepancy likely derives from
the different nodal surfaces used in the two studies, namely the one
originating from the DFT-PBE orbitals in Ref. Maezono  _et al._ , 2007 and the
one coming from the LDA orbitals in our work. Notice that both one-body and
two-body finite size errors are negligible for the EOS calculations as shown
in Fig. S-6 foo  Table 2 indicates that the equilibrium lattice parameter of
diamond is $\sim$ 0.03 bohr ($\sim$ 0.02 bohr) smaller than the experimental
one at the VMC (LRDMC) level. Our DFT calculations (Fig. S-4 foo ) suggest
that the Raman frequency is implicitly proportional to the lattice parameter,
with a gradient of $-$0.18 THz/0.01 bohr. Therefore, if the equilibrium
lattice parameter (with ZPE) were employed instead of the experimental one,
VMC (LRDMC) calculations would give $\sim$ 0.5 THz ($\sim$ 0.4 THz) higher
Raman frequencies, while still staying close to the experimental values.

Reducing the statistical errors of atomic forces is key to a successful phonon
calculation. We found that alleviating the linear dependency of a localized
atomic basis set drastically decreases the statistical error. In general, a
basis set optimized for molecular systems is not suitable for solid state
calculations, (i.e., strongly linear dependent) due to the presence of
orbitals having small exponents (c.f., typically $<$ 0.1) Peintinger  _et al._
(2013). The quality of the basis set is systematically improved by a general
and efficient scheme implemented in TurboRVB Azadi  _et al._ (2010). Indeed,
the linear dependency of a localized atomic basis set ($\psi_{\mu}(\vec{r})$)
is characterized by the condition number, $\kappa({\bf S})$, where
$S_{\mu,\nu}=\braket{\psi_{\mu}}{\psi_{\nu}}$ is the overlap matrix
444$\kappa({\bf S})$ = the maximum eigenvalue / the minimum eigenvalue of the
overlap matrix ($S_{\mu,\nu}=\braket{\psi_{\mu}}{\psi_{\nu}}$). In TurboRVB, a
redundant basis set is converted into a well-conditioned one, by disregarding
small eigenvalues and the corresponding eigenvectors of the overlap matrix
${\bf S}$ Azadi  _et al._ (2010); Nakano  _et al._ (2020b). We note that a
well-conditioned basis set can also be constructed by simply removing orbitals
having small exponents, while the method employed in this work is more general
and systematic. Figure 3 shows the plots of VMC energies, VMC forces, and
their statistical errors vs. the inverse of the condition number ($\kappa({\bf
S})^{-1}$). Figure 3 (a) indicates that the statistical error of the energy is
independent of the condition number, at variance with the statistical error of
the force, which instead strongly depends on it (Fig. 3 (b)). The error bar in
the forces amounts to $\sim 3.4\times 10^{-2}$ (Ha/bohr) when the atomic basis
set is strongly linear dependent (i.e., $\kappa({\bf S})^{-1}=10^{-15}$), a
condition that also introduces some bias in the forces because, as we have
verified, they are no longer consistent with the finite difference energy
derivatives for $\kappa({\bf S})^{-1}<10^{-12}$ (see Fig. 3b) On the other
hand, the statistical error becomes much smaller $\sim 1.7\times 10^{-3}$
(Ha/bohr) by removing the linear dependency (i.e., $\kappa({\bf
S})^{-1}=10^{-7}$). The space warp coordinate transformation (SWCT) Sorella
and Capriotti (2010) is able to further reduce the statistical error. Indeed,
Fig. 3 (b) shows that the error bar of the force becomes $\sim 4.1\times
10^{-4}$ (Ha/bohr) by using the SWCT algorithm combined with a well-
conditioned basis set, corresponding to $\sim 10^{4}$ times more efficient
computation than a bare force treatment.

We analyze now in detail the reason of this behavior. TurboRVB evaluates
atomic forces in a differential form (i.e, by the algorithmic differentiation)
Attaccalite and Sorella (2008); Sorella and Capriotti (2010):

${F_{\alpha}}={\braket{\frac{{dE}}{{d{R_{\alpha}}}}}}_{\left|{\Psi}^{\bf{R}}\right|^{2}}\simeq{\braket{\frac{E\left({\bf{R}}+\Delta{\bf{R}_{\alpha}}\right)-E\left({\bf{R}}\right)}{{\Delta{R_{\alpha}}}}}}_{\left|{\Psi}^{\bf{R}}\right|^{2}},$

where ${\bf{R}}=(R_{1},\ldots,R_{\alpha},\ldots)$ and
$\Delta{\bf{R}_{\alpha}}=(0,\ldots,\Delta R_{\alpha},\ldots)$. This equation
suggests that the statistical error on forces depends on how much the
wavefunction changes after an atom is displaced. In other words, to minimize
the stochastic error, the overlap
$\braket{\Psi^{{\bf{R}}+\Delta{\bf{R_{\alpha}}}}}{\Psi^{\bf{R}}}$ should be
close to unity. To investigate the effect of the linear dependency on the
overlap, we calculated
$\braket{\Psi^{{\bf{R}}+\Delta{\bf{R}}_{\alpha}}}{\Psi^{\bf{R}}}$ with linear-
dependent and linear-independent basis sets, using correlated sampling
techniques Nakano  _et al._ (2020b), where only one carbon atom in the 1
$\times$ 1 $\times$ 1 conventional simulation cell was displaced in the $x$
direction by $\Delta R_{\alpha}=0.005$ bohr. We obtained 0.9999 and 0.9726 for
$\kappa({\bf S})^{-1}=10^{-7}$ and $\kappa({\bf S})^{-1}=10^{-15}$,
respectively. This clearly indicates that the linear dependency of the basis
set deteriorates the overlap
$\braket{\Psi^{{\bf{R}}+\Delta{\bf{R}}_{\alpha}}}{\Psi^{\bf{R}}}$, thus
increasing the statistical error of forces.

The deterioration is explained as follows: Here, a simple Slater wavefunction
without Jastrow factor is considered for the sake of clarity. In this case,
the overlap $\braket{\Psi^{{{\bf{R}}}}}{{\Psi^{{{\bf{R}}}}}}$ reads
$\braket{\Psi^{{{\bf{R}}}}}{{\Psi^{{{\bf{R}}}}}}={\rm
det}\langle\phi_{i}^{\bf{R}}|\phi_{j}^{\bf{R}}\rangle,$ where
$\phi_{i}^{\bf{R}}$ is the $i$-th molecular orbital depending on nuclear
positions R, defining the above $N\times N$ determinant matrix. The molecular
orbital is expanded over localized atomic orbitals, i.e.,
$\phi_{i}^{\bf{R}}=\sum_{a,l}{c_{i,\\{a,l\\}}}\psi_{\\{a,l\\}}^{R_{a}}$, where
$\psi_{\\{a,l\\}}^{R_{a}}$ are (periodized) atomic orbitals explicitly
dependent on a nuclear position $R_{a}$ $\\{a=1,\ldots,\alpha,\ldots\\}$,
while $c_{i,\\{a,l\\}}$ are nuclear position independent. We can readily
derive $\braket{\Psi^{{{\bf{R}}}}}{{\Psi^{{{\bf{R}}}}}}=1$ when the molecular
orbitals are orthonormalized (i.e.,
$\braket{\phi_{i}^{{{\bf{R}}}}}{\phi_{j}^{{{\bf{R}}}}}=\delta_{i,j}$). What
about
$\braket{\Psi^{{{\bf{R}}+\Delta{\bf{R_{\alpha}}}}}}{{\Psi^{{{\bf{R}}}}}}$? We
would like to show how the perturbation
$\frac{d}{d_{R_{\alpha}}}\phi_{i}^{\bf{R}}\equiv\sum_{a=\alpha,l}{c_{i,\\{a,l\\}}}\frac{d}{d_{R_{\alpha}}}\psi_{\\{a,l\\}}^{R_{a}}$
affects the overlap. In DFT calculations, it turns out that
$\left|c_{i,\\{a,l\\}}\right|\gg 1$ when the basis is redundant (e.g.,
$\kappa({\bf S})^{-1}=10^{-15}$), while $\left|c_{i,\\{a,l\\}}\right|\ll 1$
when the basis set is well-conditioned (e.g., $\kappa({\bf S})^{-1}=10^{-7}$).
When $\left|c_{i,\\{a,l\\}}\right|\ll 1$, the perturbation effect is rather
small, and the orthonormalization condition almost holds, while
$\left|c_{i,\\{a,l\\}}\right|\gg 1$ makes the perturbation effect significant,
and by consequence the orthonormalization condition is certainly deteriorated.
This is why the linear dependency of an atomic basis set deteriorates the
overlap and, thus, induces a large error bar in forces. Thus, a complete all-
electron and pseudopotential basis set database suitable for QMC solid state
calculations will be quite useful for the application of QMC-forces to the
calculations of phonons in realistic materials.

In summary, we report a VMC determination of the momentum-resolved phonon
dispersion of diamond. Our approach combines the ab initio quantum Monte Carlo
framework with the so-called frozen phonon technique. It gives results in very
good agreement with experiments and provides renormalized harmonic optical
frequencies consistent with the experimental findings. We estimated the purely
harmonic contribution to the phonon spectrum, by evaluating $q$-dependent
anharmonic corrections by means of a path integral molecular dynamics driven
by PBE forcesMorresi  _et al._ . After including these corrections, the VMC
phonon spectrum agrees very well with the experimental phonon dispersion. We
found that alleviating the atomic basis-set redundancy of the trial
wavefunction is key to reduce the statistical error of atomic forces and,
thus, to make the VMC phonons calculations feasible over the full Brillouin
zone. This achievement paves the way to new relevant applications, for
instance, in correlated materials and Van der Waals crystals (i.e., molecular
crystals Subedi and Boeri (2011); Casula  _et al._ (2012)), where sometimes
phonons are poorly reproduced within the DFT framework.

###### Acknowledgements.

K.N. is grateful for computational resources from the facilities of Research
Center for Advanced Computing Infrastructure at Japan Advanced Institute of
Science and Technology (JAIST). K.N. and S.S. are grateful for computational
resources from PRACE project No. 2019204934. T.M. and M.C. acknowledge that
this work was supported by French state funds managed by the ANR within the
Investissements d’Avenir programme under reference ANR-11-IDEX-0004-02, and
more specifically within the framework of the Cluster of Excellence MATISSE
led by Sorbonne University. M.C. is grateful to the French Grand équipement
national de calcul intensif (GENCI) for the computational time provided
through the Project No. 0906493. K.N., M.C., and S.S. are grateful for
computational resources of the supercomputer Fugaku provided by RIKEN through
the HPCI System Research Project (Project ID: hp200164). K.N. acknowledges a
support from the JSPS Overseas Research Fellowships and that from Grant-in-Aid
for Scientific Research on Innovative Areas (No. 16H06439). R.M. is grateful
for financial supports from MEXT-KAKENHI (19H04692 and 16KK0097), from Toyota
Motor Corporation, from the Air Force Office of Scientific Research (AFOSR-
AOARD/FA2386-17-1-4049;FA2386-19-1-4015), and from JSPS Bilateral Joint
Projects (with India DST). S.S. also acknowledges a financial support from
PRIN 2017BZPKSZ. The authors appreciate helpful comments by A. Zen.

##  Appendix A  Details of QMC calculations

All Variational Monte Carlo (VMC) and lattice regularized diffusion Monte
Carlo (LRDMC) Casula  _et al._ (2005); Nakano  _et al._ (2020a) calculations
in this study were performed by a SISSA quantum Monte Carlo package, called
TurboRVB Nakano  _et al._ (2020b). The package employs a many-body WF ansatz
$\Psi$ which can be written as the product of two terms, i.e.,
$\Psi=\Phi_{\text{AS}}\times\exp J\,,$ where the term $\exp J$ and
$\Phi_{\text{AS}}$ are conventionally called Jastrow and antisymmetric parts,
respectively. The antisymmetric part is denoted as the Antisymmetrized Geminal
Power (AGP) that reads:
${\Psi_{{\text{AGP}}}}\left({{{\mathbf{r}}_{1}},\ldots,{{\mathbf{r}}_{N}}}\right)={\hat{A}}\left[{\Phi\left({{\mathbf{r}}_{1}^{\uparrow},{\mathbf{r}}_{1}^{\downarrow}}\right)\Phi\left({{\mathbf{r}}_{2}^{\uparrow},{\mathbf{r}}_{2}^{\downarrow}}\right)\cdots\Phi\left({{\mathbf{r}}_{N/2}^{\uparrow},{\mathbf{r}}_{N/2}^{\downarrow}}\right)}\right],$
where ${\hat{A}}$ is the antisymmetrization operator, and
$\Phi\left({{\mathbf{r}}^{\uparrow},{\mathbf{r}}^{\downarrow}}\right)$ is
called the paring function Casula and Sorella (2003). The spatial part of the
geminal function is expanded over the Gaussian-type atomic orbitals:
${\Phi_{{\text{AGP}}}}\left({{{\mathbf{r}}_{i}},{{\mathbf{r}}_{j}}}\right)=\sum\limits_{l,m,a,b}{{f_{\left\\{{a,l}\right\\},\left\\{{b,m}\right\\}}}{\psi_{a,l}}\left({{{\mathbf{r}}_{i}}}\right){\psi_{b,m}}\left({{{\mathbf{r}}_{j}}}\right)}$
where ${\psi_{a,l}}$ and ${\psi_{b,m}}$ are primitive Gaussian atomic
orbitals, their indices $l$ and $m$ indicate different orbitals centered on
atoms $a$ and $b$, and $i$ and $j$ are coordinates of spin up and down
electrons, respectively, and
${{f_{\left\\{{a,l}\right\\},\left\\{{b,m}\right\\}}}}$ are the variational
parameters. In this study, a triple-zeta basis set (11$s$11$p$2$d$1$f$)
accompanying an energy-consistent pseudo potential Burkatzki  _et al._ (2007)
was employed for the atomic orbitals of the antisymmetric part (Table S-I).
The pairing function can be also written as
${\Phi_{{\text{AGPn}}}}\left({{{\mathbf{r}}_{i}},{{\mathbf{r}}_{j}}}\right)=\sum_{k=1}^{M}\lambda_{k}\phi_{k}(\mathbf{r}_{i})\phi_{k}(\mathbf{r}_{j})$
with $\lambda_{k}>0$, where $\phi_{k}(\mathbf{r})$ is a molecular orbital,
i.e., $\phi_{k}(\mathbf{r})=\sum_{i=1}^{L}{c_{i,k}{\psi_{i}}(\mathbf{r})}$.
When the paring function is expanded over $M$ molecular orbitals where $M$ is
equal to the half of the total number of electrons ($N/2$), the AGP coincides
with the Slater-Determinant ansatz Becca and Sorella (2017); Marchi  _et al._
(2009). In this study, we restricted ourselves to a Jastrow-Slater determinant
(JSD) by setting $M=\frac{1}{2}\cdot N$, wherein the coefficients of atomic
orbitals, i.e., ${c_{i,k}}$, were obtained by a Density Functional theory
(DFT) calculation, and were fixed during a VMC optimization.

Figure S-1:  The phonon dispersions of diamond calculated using 2 $\times$ 2 $\times$ 2 conventional supercell at the VMC level. “Raw data” denotes the phonon dispersion obtained by the finite-displacement method, while “Corrected” is the phonon dispersion where the anharmonic and one-body finite-size corrections are included. The error bars are not shown here for the sake of clarity. Observed experimental results are also plotted for comparison. Experimental data are taken from Refs. Warren  _et al._ , 1965; *1967WAR; Liu  _et al._ , 2000; Kulda  _et al._ , 2002; Schwoerer-Böhning  _et al._ , 1998. The data in Ref. Schwoerer-Böhning  _et al._ , 1998 were digitalized using WebPlotDigitizer Rohatgi (2020).  Table S-I:  The valence triple-zeta (VTZ) uncontracted basis set (11$s$11$p$2$d$1$f$) used for the determinant part in the present work. Orbital type | Exponent | Orbital type | Exponent
---|---|---|---
$s$ | 13.073594 | $p$ | 1.871016
$s$ | 6.541187 | $p$ | 0.935757
$s$ | 3.272791 | $p$ | 0.468003
$s$ | 1.637494 | $p$ | 0.376742
$s$ | 0.921552 | $p$ | 0.234064
$s$ | 0.819297 | $p$ | 0.126772
$s$ | 0.409924 | $p$ | 0.117063
$s$ | 0.205100 | $p$ | 0.058547
$s$ | 0.132800 | $p$ | 0.029281
$s$ | 0.102619 | $d$ | 1.141611
$s$ | 0.051344 | $d$ | 0.329486
$p$ | 7.480076 | $f$ | 0.773485
$p$ | 3.741035 |  |
|  |  |

The Jastrow term is composed of one-body, two-body and three/four-body factors
($J={J_{1}}{J_{2}}{J_{3/4}}$). The one-body and two-body factors are
essentially used to fulfill the electron-ion and electron-electron cusp
conditions, respectively, and the three/four-body factor is employed to
consider further electron-electron correlations (e.g., electron-nucleus-
electron). Since we employed an energy-consistent pseudo potential developed
by Burkatzki et al. Burkatzki  _et al._ (2007), we used only the inhomogeneous
part of ${J_{1}^{inh}}$, $J_{2}$ and $J_{3/4}$ in this study (i.e., no
electron-ion cusp corrections are needed). The inhomogeneous one-body Jastrow
factor ${J_{1}^{inh}}$ reads
${J_{1}^{inh}}\left({{{\mathbf{r}}_{1}},\ldots,{{\mathbf{r}}_{N}}}\right)=\sum_{i=1}^{N}\sum_{a=1}^{N_{\text{atom}}}\left({\sum\limits_{l}{M_{a,l}\chi_{a,l}\left({{{\mathbf{r}}_{i}}}\right)}}\right),$
where ${{{\mathbf{r}}_{i}}}$ are the electron positions,
${{{\mathbf{R}}_{a}}}$ are the atomic positions with corresponding atomic
number $Z_{a}$, $l$ runs over atomic orbitals $\chi_{a,l}$ (e.g., GTO)
centered on the atom $a$, ${N_{\text{atom}}}$ is the total number of atoms in
a system, and $\\{M_{a,l}\\}$ are variational parameters. The two-body Jastrow
factor is defined as:
${J_{2}}\left({{{\mathbf{r}}_{1}},\ldots{{\mathbf{r}}_{N}}}\right)=\exp\left({\sum\limits_{i<j}{v\left({\left|{{{\mathbf{r}}_{i}}-{{\mathbf{r}}_{j}}}\right|}\right)}}\right),$
where $v\left(r\right)$ is $v\left(r\right)=\frac{1}{2}r\cdot{\left({1-F\cdot
r}\right)^{-1}}$ and $F$ is a variational parameter. The three-body Jastrow
factor is:
${J_{3/4}}\left({{{\mathbf{r}}_{1}},\ldots{{\mathbf{r}}_{N}}}\right)=\exp\left({\sum\limits_{i<j}{{\Phi_{{\text{Jas}}}}\left({{{\mathbf{r}}_{i}},{{\mathbf{r}}_{j}}}\right)}}\right),$
and
${\Phi_{{\text{Jas}}}}\left({{{\mathbf{r}}_{i}},{{\mathbf{r}}_{j}}}\right)=\sum\limits_{l,m,a,b}{g_{a,l,m,b}\chi_{a,l}^{{\text{Jas}}}\left({{{\mathbf{r}}_{i}}}\right)\chi_{b,m}^{{\text{Jas}}}\left({{{\mathbf{r}}_{j}}}\right)},$
where the indices $l$ and $m$ again indicate different orbitals centered on
corresponding atoms $a$ and $b$. In this study, the coefficients of the
three/four-body Jastrow factor were set to zero for $a\neq b$ because it
significantly decreases the number of variational parameters while rarely
affects variational energies . A basis set (3$s$2$p$) was employed for the
atomic orbitals of the Jastrow part. The variational parameters in the Jastrow
factor were optimized by the so-called stochastic reconfiguration Sorella  _et
al._ (2007) and/or the modified linear method Umrigar  _et al._ (2007); Nakano
_et al._ (2020b) implemented in TurboRVB. Using the optimized wavefunction,
energies and forces are calculated at the VMC and the LRDMC levels. All LRDMC
calculations were performed by the original single-grid scheme Casula  _et
al._ (2005) with the discretization grid size $a=0.20$ Bohr. For the Raman
frequency calculation, several lattice spaces, i.e., $a$ = 0.20, 0.30, 0.35,
0.40, and 0.50 Bohr were also used to extrapolate the energies (i.e., $a\to
0$) with a quadratic function (i.e., $E(a)=E_{0}+k_{1}\cdot a^{2}+k_{2}\cdot
a^{4}$).

Figure S-1 shows the phonon dispersions and Raman frequencies of diamond
calculated using 2 $\times$ 2 $\times$ 2 conventional supercell at the VMC
level. “Raw data” denotes the phonon dispersion obtained by the finite-
displacement method, while “Corrected” is the phonon dispersion where the one-
body finite-size (see sec. B) and anharmonic (see sec. C) corrections are
included.

Figure S-2: Phonon dispersions of diamond obtained by different methods with
different computational conditions at the DFT level, wherein the experimental
lattice parameter was used. (a) A phonon dispersion calculated by the finite-
displacement method (FDM), or the so-called frozen-phonon method, implemented
in Phonopy Togo and Tanaka (2015) with 2 $\times$ 2 $\times$ 2 conventional
supercell at $L$ point. The displacement of 0.03 Bohr was employed for
calculating derivatives of atomic forces with respect to a nucleus position.
(b) The same as (b) except for $k$-grid. The shifted 6 $\times$ 6 $\times$ 6
Monkhorst-Pack grid Monkhorst and Pack (1976) was employed. (c) A phonon
dispersion calculated by the linear response method implemented Quantum
Espresso package Giannozzi  _et al._ (2009) with 1 $\times$ 1 $\times$ 1
primitive unit cell. The shifted (unshifted) 8 $\times$ 8 $\times$ 8
Monkhorst-Pack grid Monkhorst and Pack (1976) was employed for $k$-integration
($q$-interpolation). The DFT calculations were performed with GGA-PBE
functionals with Ultrasoft (US) pseudo potential taken from PS-Library Dal
Corso (2014).  Figure S-3: Phonon frequencies calculated by DFT at $q$ =
$\Gamma$ and $L$ points v.s. employed displacements in the finite displacement
method, where the GGA-PBE functional was employed. The phonon dispersions were
calculated using 2 $\times$ 2 $\times$ 2 conventional supercell with $k$ twist
= $L$ ($\pi$,$\pi$,$\pi$). Figure S-4: Plots of Raman Frequencies vs. lattice
parameters of diamond. They were calculated using the linear response theory
implemented in Quantum Espresso at a single $q$ point, i.e., $q$ = (0,0,0)
with the LDA-PZ exchange-correlation functional. The primitive lattice was
employed with $k$-grid = 8 $\times$ 8 $\times$ 8\.

##  Appendix B  Validation of computational conditions using DFT.

There are two major implementations to compute a harmonic phonon dispersion,
i.e., the finite displacement method (FDM) and the linear response method. The
former is also called frozen phonon method. Their results should be consistent
as far as employed hyperparameters are correct. Figure S-2 shows phonon
dispersions of diamond obtained by the two different methods at the DFT level,
wherein the experimental lattice parameter was used. (a) A phonon dispersion
calculated by the frozen-phonon method implemented in Phonopy Togo and Tanaka
(2015) with 2 $\times$ 2 $\times$ 2 conventional supercell at $L$ point. The
displacement of 0.03 Bohr was employed for calculating derivatives of atomic
forces with respect to a nucleus position. (b) The same as (a) except for
$k$-grid. The shifted 6 $\times$ 6 $\times$ 6 Monkhorst-Pack grid Monkhorst
and Pack (1976) was employed. (c) A phonon dispersion calculated by the linear
response method implemented Quantum Espresso package Giannozzi  _et al._
(2009) with 1 $\times$ 1 $\times$ 1 primitive unit cell. The 8 $\times$ 8
$\times$ 8 Monkhorst-Pack grids Monkhorst and Pack (1976) were employed for
$k$ and $q$ integrations. The DFT calculations were performed with GGA-PBE
functionals with Ultrasoft (US) pseudo potential taken from PS-Library Dal
Corso (2014). Comparison of (a)-red with (c)-green proves that the one-body
finite size error is very small in a phonon calculation thanks to the error
cancellation when $L$ twist for 2 $\times$ 2 $\times$ 2 conventional supercell
is employed. In detail, the one-body finite size errors are 0.16 THz, 0.18
THz, and 0.23 THz for $q$ = $\Gamma$, $X$, and $L$, respectively. Comparison
of (b)-blue with (c)-green, which are almost overlapped, suggests that the
finite difference method does not bias the result. As mentioned in the main
text, 0.15 Bohr was employed for the displacement of a carbon atom in the VMC
frozen-phonon calculation to decrease the statistical error. Figure. S-3 shows
the amplitude of 0.15 Bohr employed underestimates frequencies by $\sim$ 0.12
THz. A single-phonon phonon calculation at $\Gamma$ alleviates this
underestimation as mentioned in the main text. Fig. S-4 shows a plot of Raman
Frequencies vs. lattice parameters of diamond. They were calculated using the
linear response theory implemented in Quantum Espresso at a single $q$ point,
i.e., $q$ = (0,0,0) with the LDA-PZ exchange-correlation functional. The
primitive lattice was employed with 8 $\times$ 8 $\times$ 8.

##  Appendix C  Anharmonic renormalizations from path integral molecular
dynamics

The anharmonic renormalizations of phonon frequencies reported in Table 1 with
the superscript ’$f$’ were computed through the displacement-displacement
zero-time Kubo-correlator built upon path integral molecular dynamics (PIMD)
trajectories. The evaluation of the displacement-displacement correlator over
classical molecular dynamics trajectories is a standard method to predict
force constant matrices and phonon spectra, but it misses nuclear quantum
effects. The details of the extension of such a method to extract phonon
dispersions from path integral molecular dynamics can be found in Ref. Morresi
_et al._ , .

In Fig S-5, we report the whole spectrum that we have obtained using this new
approach compared to the linear response DFT curve.

Figure S-5: (a) Phonon spectrum at 300 K obtained from PIMD simulation versus
linear response (LR) DFT spectrum; (b) Zoom over the optical modes.

In particular, the ion equations of motion were integrated using the PIOUD
algorithm Mouhat  _et al._ (2017), while at each time step the electronic
potential energy surface was evaluated at DFT level using the Quantum Espresso
package Giannozzi  _et al._ (2009) and GGA-PBE Perdew  _et al._ (1996)
exchange-correlation functionals. The simulation was carried out at 300 Kelvin
for 34 picoseconds employing 12 beads. We used the same 64-atoms supercell as
for the QMC calculations, with a 2 $\times$ 2 $\times$ 2 $k$-grid for
electronic integration and a cutoff for wavefunctions equal to 60 Ry.

Figure S-6: Finite-size extrapolation for different lattice parameters. The differences of energies are constant regardless of the simulation cell size, indicating that the finite-size error is negligible for an EOS calculation thanks to the error cancellation. Since diamond is an insulator, raw QMC energies obtained at $L$ point can be smoothly extrapolated by a linear function. See. Ref. Hennig  _et al._ , 2010 Figure S-7: VMC forces, its Hellmann-Feynman, and its Pulay contributions v.s. inverse of the condition number of the overlap matrix. The forces were calculated with the space warp coordinate transformations (SWCT) Sorella and Capriotti (2010). The same VTZ atomic basis set and the pseudo potential Burkatzki  _et al._ (2007) as in the phonon and EOS calculations were adopted, while 1 $\times$ 1 $\times$ 1 conventional cell (8 carbon atoms in a simulation cell) with $k$ = $\Gamma$ twist was employed. Only one carbon atom was displaced in $x$ direction by 0.15 Bohr.  Table S-II:  Error bars of forces corresponding to Fig. 3(b) in the main text. $\kappa({\bf S})^{-1}$ denotes the inverse of the condition number of the overlap matrix ${\bf S}$. The used geometry is showin in Table S-III $\kappa({\bf S})^{-1}$ | with SWCT | without SWCT
---|---|---
Force ($F_{x}$ of C1) (Ha/bohr) | Error bar (Ha/bohr) | Force ($F_{x}$ of C1) (Ha/bohr) | Error bar (Ha/bohr)
$10^{-5}$ | -6.659 $\times$ $10^{-2}$ | 4.1 $\times$ $10^{-4}$ | -6.71 $\times$ $10^{-2}$ | 1.7 $\times$ $10^{-3}$
$10^{-6}$ | -6.622 $\times$ $10^{-2}$ | 4.1 $\times$ $10^{-4}$ | -6.67 $\times$ $10^{-2}$ | 1.4 $\times$ $10^{-3}$
$10^{-7}$ | -6.523 $\times$ $10^{-2}$ | 4.4 $\times$ $10^{-4}$ | -6.68 $\times$ $10^{-2}$ | 1.3 $\times$ $10^{-3}$
$10^{-8}$ | -6.635 $\times$ $10^{-2}$ | 4.6 $\times$ $10^{-4}$ | -6.77 $\times$ $10^{-2}$ | 1.5 $\times$ $10^{-3}$
$10^{-10}$ | -6.734 $\times$ $10^{-2}$ | 8.0 $\times$ $10^{-4}$ | -6.59 $\times$ $10^{-2}$ | 1.6 $\times$ $10^{-3}$
$10^{-12}$ | -13.87 $\times$ $10^{-2}$ | 1.8 $\times$ $10^{-3}$ | -13.86 $\times$ $10^{-2}$ | 2.4 $\times$ $10^{-3}$
$10^{-15}$ | -4.7 $\times$ $10^{-2}$ | 3.3 $\times$ $10^{-2}$ | -4.9 $\times$ $10^{-2}$ | 3.4 $\times$ $10^{-2}$
|  |  |  |

##  Appendix D  Finite-size errors in QMC calculations

The effect of the finite-size errors on an EOS calculation has been
investigated. Figure S-6 shows that finite-size extrapolations for two
different lattice parameters ($a$ = 6.6107 and 6.7407 Bohr). The differences
of the energies are constant regardless of the simulation cell size,
indicating that the finite-size error is negligible for an EOS calculation
thanks to the error cancellation. Since diamond is an insulator, raw QMC
energies obtained at $L$ point can be smoothly extrapolated by a linear
function. See. Ref. Hennig  _et al._ , 2010.

##  Appendix E  Vinet exponential function

In this study, the obtained energies were fitted by the Vinet exponential
function Vinet  _et al._ (1987) that reads:

$\displaystyle E(V)$ $\displaystyle=$ $\displaystyle
E_{0}+\frac{2B_{0}V_{0}}{\left(B_{0}^{\prime}-1\right)^{2}}\left\\{2-\left[5+3\left(\frac{V}{V_{0}}\right)^{\frac{1}{3}}(B_{0}^{\prime}-1)-3B_{0}^{\prime}\right]\right.$
$\displaystyle\times$
$\displaystyle\left.\exp\left[-\frac{3}{2}\left(B_{0}^{\prime}-1\right)\left[\left(\frac{V}{V_{0}}\right)^{\frac{1}{3}}-1\right]\right]\right\\},$

where $E(V)$ is the total energy per atom, $V$ is volume per atom, and
$E_{0}$, $V_{0}$, $B_{0}$, and $B_{0}^{\prime}$ are parameters. The non-linear
fittings were performed using SciPy module Virtanen  _et al._ (2020)
implemented in Python.

##  Appendix F Atomic forces

TurboRVB evaluates atomic forces in a differential expression Sorella and
Capriotti (2010):

${{\mathbf{F}}_{a}}=-\left\langle{\frac{d}{{d{{\mathbf{R}}_{a}}}}{e_{L}}}\right\rangle+2\left({\left\langle{{e_{L}}}\right\rangle\left\langle{\frac{d}{{d{{\mathbf{R}}_{a}}}}\ln\left({{{\mathcal{J}}^{\frac{1}{2}}}\Psi}\right)}\right\rangle-\left\langle{{e_{L}}\frac{d}{{d{{\mathbf{R}}_{a}}}}\ln\left({{{\mathcal{J}}^{\frac{1}{2}}}\Psi}\right)}\right\rangle}\right)$

where ${\mathcal{J}}$ is the Jacobian of the space warp coordinate
transformation (SWCT) Sorella and Capriotti (2010), ${e_{L}}$ is the local
energy, and the brackets indicates a Monte Carlo average over the trial WF.
The first term is called “Hellmann-Feynman” force (HF) and the second and
third terms are called “Pulay” force (PF). All the terms above can be written
by the partial derivatives of the local energy and those of the logarithm of
the WF Sorella and Capriotti (2010). These differential expressions are
efficiently computed in TurboRVB, by using the adjoint algorithmic
differentiation technique Sorella and Capriotti (2010). TurboRVB also employs
the so-called reweighting technique developed by Attaccalite and Sorella
Attaccalite and Sorella (2008) to avoid divergence of the derivatives in the
vicinity of the nodal surfaces. Fig. S-7 shows VMC forces, its Hellmann-
Feynman, and its Pulay contributions v.s. inverse of the condition number of
the overlap matrix, corresponding to Fig. 3 in the main text. Fig. S-7
clarifies that the erratic behavior of the forces at $\kappa({\bf
S})^{-1}=10^{-12}$ comes from its Pulay contribution, indicating that the
linear-dependency of a basis-set also deteriorates absolute values of forces.

Table S-III:  The geometry of the displaced diamond used to obtain forces in Table S-II. The lattice parameter is 6.741 Bohr. $x$, $y$, and $z$ are represented by fractional coordinates. Element | Label | $x$ | $y$ | $z$
---|---|---|---|---
C | C1 | 0.02225 | 0.00000 | 0.00000
C | C2 | 0.00000 | 0.50000 | 0.50000
C | C3 | 0.50000 | 0.00000 | 0.50000
C | C4 | 0.50000 | 0.50000 | 0.00000
C | C5 | 0.25000 | 0.25000 | 0.75000
C | C6 | 0.75000 | 0.75000 | 0.75000
C | C7 | 0.25000 | 0.75000 | 0.25000
C | C8 | 0.75000 | 0.25000 | 0.25000
|  |  |  |

##

  * Martin (2004) R. M. Martin,  _Electronic structure : basic theory and practical methods_ (Cambridge University Press, 2004).
  * Sholl and Steckel (2011) D. Sholl and J. A. Steckel,  _Density functional theory: a practical introduction_ (John Wiley & Sons, 2011).
  * Baroni  _et al._ (2001) S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Rev. Mod. Phys. 73, 515 (2001).
  * Togo and Tanaka (2015) A. Togo and I. Tanaka, Scr. Mater. 108, 1 (2015).
  * Mounet and Marzari (2005) N. Mounet and N. Marzari, Phys. Rev. B 71, 205214 (2005).
  * Lazzeri  _et al._ (2008) M. Lazzeri, C. Attaccalite, L. Wirtz, and F. Mauri, Phys. Rev. B 78, 081406(R) (2008).
  * Krisch  _et al._ (2011) M. Krisch, D. L. Farber, R. Xu, D. Antonangeli, C. M. Aracne, A. Beraud, T.-C. Chiang, J. Zarestky, D. Y. Kim, E. I. Isaev, R. Ahuja, and B. Johansson, Proc. Natl. Acad. Sci. U.S.A. 108, 9342 (2011).
  * Reznik  _et al._ (2008) D. Reznik, G. Sangiovanni, O. Gunnarsson, and T. Devereaux, Nature 455, E6 (2008).
  * Floris  _et al._ (2011) A. Floris, S. de Gironcoli, E. K. U. Gross, and M. Cococcioni, Phys. Rev. B 84, 161102(R) (2011).
  * Leonov  _et al._ (2012) I. Leonov, A. I. Poteryaev, V. I. Anisimov, and D. Vollhardt, Phys. Rev. B 85, 020401(R) (2012).
  * Leonov  _et al._ (2014) I. Leonov, V. I. Anisimov, and D. Vollhardt, Phys. Rev. Lett. 112, 146401 (2014).
  * Foulkes  _et al._ (2001) W. Foulkes, L. Mitas, R. Needs, and G. Rajagopal, Rev. Mod. Phys. 73, 33 (2001).
  * Marchi  _et al._ (2011) M. Marchi, S. Azadi, and S. Sorella, Phys. Rev. Lett. 107, 086807 (2011).
  * Casula and Sorella (2013) M. Casula and S. Sorella, Phys. Rev. B 88, 155125 (2013).
  * Sorella  _et al._ (2018) S. Sorella, K. Seki, O. O. Brovko, T. Shirakawa, S. Miyakoshi, S. Yunoki, and E. Tosatti, Phys. Rev. Lett. 121, 066402 (2018).
  * Zen  _et al._ (2012) A. Zen, D. Zhelyazov, and L. Guidoni, J. Chem. Theory Comput. 8, 4204 (2012).
  * Luo  _et al._ (2014) Y. Luo, A. Zen, and S. Sorella, J. Chem. Phys. 141, 194112 (2014).
  * Liu  _et al._ (2019) Y. Y. F. Liu, B. Andrews, and G. J. Conduit, J. Chem. Phys. 150, 034104 (2019).
  * Nakano  _et al._ (2019) K. Nakano, R. Maezono, and S. Sorella, J. Chem. Theory Comput. 15, 4044 (2019).
  * Maezono  _et al._ (2007) R. Maezono, A. Ma, M. D. Towler, and R. J. Needs, Phys. Rev. Lett. 98, 025701 (2007).
  * Azadi  _et al._ (2010) S. Azadi, C. Cavazzoni, and S. Sorella, Phys. Rev. B 82, 125112 (2010).
  * Liu  _et al._ (2000) M. S. Liu, L. A. Bursill, S. Prawer, and R. Beserman, Phys. Rev. B 61, 3391 (2000).
  * Schwoerer-Böhning  _et al._ (1998) M. Schwoerer-Böhning, A. T. Macrander, and D. A. Arms, Phys. Rev. Lett. 80, 5572 (1998).
  * Kulda  _et al._ (2002) J. Kulda, H. Kainzmaier, D. Strauch, B. Dorner, M. Lorenzen, and M. Krisch, Phys. Rev. B 66, 241202(R) (2002).
  * Vanderbilt  _et al._ (1984) D. Vanderbilt, S. G. Louie, and M. L. Cohen, Phys. Rev. Lett. 53, 1477 (1984).
  * Parrish (1960) W. Parrish, Acta Cryst. 13, 838 (1960).
  * Warren  _et al._ (1965) J. L. Warren, R. G. Wenzel, and J. L. Yarnell, in  _Inelastic Scattering of Neutrons_ , Vol. 1 (International Atomic Energy Agency, Vienna, 1965) pp. 361–371.
  * Warren  _et al._ (1967) J. L. Warren, J. L. Yarnell, G. Dolling, and R. A. Cowley, Phys. Rev. 158, 805 (1967).
  * Rohatgi (2020) A. Rohatgi, “Webplotdigitizer: Version 4.3,” https://automeris.io/WebPlotDigitizer (2020).
  * Casula  _et al._ (2005) M. Casula, C. Filippi, and S. Sorella, Phys. Rev. Lett. 95, 100201 (2005).
  * Nakano  _et al._ (2020a) K. Nakano, R. Maezono, and S. Sorella, Phys. Rev. B 101, 155106 (2020a).
  * Casula and Sorella (2003) M. Casula and S. Sorella, J. Chem. Phys. 119, 6500 (2003).
  * Nakano  _et al._ (2020b) K. Nakano, C. Attaccalite, M. Barborini, L. Capriotti, M. Casula, E. Coccia, M. Dagrada, C. Genovese, Y. Luo, G. Mazzola, _et al._ , J. Chem. Phys. 152, 204121 (2020b).
  * Burkatzki  _et al._ (2007) M. Burkatzki, C. Filippi, and M. Dolg, J. Chem. Phys. 126, 234105 (2007).
  * (35) See Supplemental Material for details of the DFT and QMC calculations, and for the supplemental Figures S-1-S-7 and Tables S-I-S-III.
  * Perdew and Zunger (1981) J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981).
  * Sorella  _et al._ (2007) S. Sorella, M. Casula, and D. Rocca, J. Chem. Phys. 127, 014105 (2007).
  * Umrigar  _et al._ (2007) C. J. Umrigar, J. Toulouse, C. Filippi, S. Sorella, and R. G. Hennig, Phys. Rev. Lett. 98, 110201 (2007).
  * Hennig  _et al._ (2010) R. G. Hennig, A. Wadehra, K. P. Driver, W. D. Parker, C. J. Umrigar, and J. W. Wilkins, Phys. Rev. B 82, 014101 (2010).
  * Sorella  _et al._ (2011) S. Sorella, M. Casula, L. Spanu, and A. Dal Corso, Phys. Rev. B 83, 075119 (2011).
  * Becca and Sorella (2017) F. Becca and S. Sorella,  _Quantum Monte Carlo approaches for correlated systems_ (Cambridge University Press, 2017).
  * Giannozzi  _et al._ (2009) P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, _et al._ , J. Phys. Condens. Matter. 21, 395502 (2009).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Pavone  _et al._ (1993) P. Pavone, K. Karch, O. Schütt, W. Windl, D. Strauch, P. Giannozzi, and S. Baroni, Phys. Rev. B 48, 3156 (1993).
  * Kresse  _et al._ (1995) G. Kresse, J. Furthmüller, and J. Hafner, EPL 32, 729 (1995).
  * Occelli  _et al._ (2003) F. Occelli, P. Loubeyre, and R. LeToullec, Nat. Mater. 2, 151 (2003).
  * McSkimin and Andreatch Jr (1972) H. McSkimin and P. Andreatch Jr, J. Appl. Phys. 43, 2944 (1972).
  * Vinet  _et al._ (1987) P. Vinet, J. R. Smith, J. Ferrante, and J. H. Rose, Phys. Rev. B 35, 1945 (1987).
  * Sorella and Capriotti (2010) S. Sorella and L. Capriotti, J. Chem. Phys. 133, 234111 (2010).
  * Note (1) 2 $\times$ 2 $\times$ 2 conventional supercell containing 64 atoms in the simulation cell is large enough for obtaining a phonon dispersion almost consistent with experimental ones, as shown in Fig. S-2 foo . This was confirmed by comparing the phonon dispersion obtained by the finite displacement method and that obtained by the linear-response method with very dense $k$ and $q$ grids at the DFT level.
  * Harris (1995) N. J. Harris, J. Phys. Chem. 99, 14689 (1995).
  * Mouhat  _et al._ (2017) F. Mouhat, S. Sorella, R. Vuilleumier, A. M. Saitta, and M. Casula, J. Chem. Theory Comput. 13, 2400 (2017).
  * (53) T. Morresi, L. Paulatto, R. Vuilleumier, and M. Casula, arXiv:2103.04094 .
  * Note (2) In the Raman mode, two nearest-neighbor carbon atoms are displaced in opposite directions by a distance $u$ from their high-symmetry positions Maezono  _et al._ (2007). The distortion in the [1 0 0] direction was employed in this study Vanderbilt  _et al._ (1984). We calculated VMC forces for two distorted structures ($u=0.03$ Bohr and $u=0.05$ Bohr) and obtained the Raman frequency by fitting the forces with a linear function. We also calculated VMC and LRDMC energies for the undistorted structure and two distorted structures ($u=0.03$ Bohr and $u=0.05$ Bohr), then obtained Raman frequencies by fitting the energies with a quadratic function.
  * Moroni  _et al._ (2014) S. Moroni, S. Saccani, and C. Filippi, J. Chem. Theory Comput. 10, 4823 (2014).
  * Note (3) Since smaller displacements were employed in this study (i.e., 0.03 and 0.05 Bohr) than the previous one, the error bars in phonon frequencies at $\Gamma$ point are a bit larger.
  * Drummond  _et al._ (2004) N. D. Drummond, M. D. Towler, and R. J. Needs, Phys. Rev. B 70, 235119 (2004).
  * Krukau  _et al._ (2006) A. V. Krukau, O. A. Vydrov, A. F. Izmaylov, and G. E. Scuseria, J. Chem. Phys. 125, 224106 (2006).
  * Sun  _et al._ (2015) J. Sun, A. Ruzsinszky, and J. P. Perdew, Phys. Rev. Lett. 115, 036402 (2015).
  * Hao  _et al._ (2012) P. Hao, Y. Fang, J. Sun, G. I. Csonka, P. H. T. Philipsen, and J. P. Perdew, Phys. Rev. B 85, 014111 (2012).
  * Peintinger  _et al._ (2013) M. F. Peintinger, D. V. Oliveira, and T. Bredow, J. Comput. Chem. 34, 451 (2013).
  * Note (4) $\kappa({\bf S})$ = the maximum eigenvalue / the minimum eigenvalue of the overlap matrix ($S_{\mu,\nu}=\mathinner{\langle{\psi_{\mu}|\psi_{\nu}}\rangle}$).
  * Attaccalite and Sorella (2008) C. Attaccalite and S. Sorella, Phys. Rev. Lett. 100, 114501 (2008).
  * Subedi and Boeri (2011) A. Subedi and L. Boeri, Phys. Rev. B 84, 020508 (2011).
  * Casula  _et al._ (2012) M. Casula, M. Calandra, and F. Mauri, Phys. Rev. B 86, 075445 (2012).
  * Marchi  _et al._ (2009) M. Marchi, S. Azadi, M. Casula, and S. Sorella, J. Chem. Phys. 131, 154116 (2009).
  * Monkhorst and Pack (1976) H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976).
  * Dal Corso (2014) A. Dal Corso, Comput. Mater. Sci. 95, 337 (2014).
  * Virtanen  _et al._ (2020) P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, Nat. Methods 17, 261 (2020).


# Enabling and Assessing Trust when Cooperating with Robots in Disaster
Response (EASIER)*

Laurent Frering1, Matthias Eder1, Bettina Kubicek2 Dietrich Albert2, Denis
Kalkofen3, Thomas Gschwandtner4,
Heimo Krajnz5 and Gerald Steinbauer-Wagner1 *This work was partially supported
by the Austrian Research Promotion Agency (FFG) with the project EASIER
(FO999886357).1Laurent Frering, Matthias Eder and Gerald Steinbauer-Wagner are
with the Institute of Software Technology, Graz University of Technology,
Graz, Austria. {laurent.frering, matthias.eder,
<EMAIL_ADDRESS>Kubicek and Dietrich Albert are with the
Institute of Psychology, University of Graz, Graz, Austria. {bettina.kubicek,
<EMAIL_ADDRESS>Denis Kalkofen is with the Institute of Computer
Graphics and Vision, Graz University of Technology, Graz, Austria.
<EMAIL_ADDRESS>Thomas Gschwandtner is with Rosenbauer International
AG, Linz, Austria<EMAIL_ADDRESS>Heimo Krajnz is with
Professional Fire Brigade Graz, Graz, Austria<EMAIL_ADDRESS>
###### Abstract

This paper presents a conceptual overview of the EASIER project and its scope.
EASIER focuses on supporting emergency forces in disaster response scenarios
with a semi-autonomous mobile manipulator. Specifically, we examine the
operator’s trust in the system and his/her cognitive load generated by its
use. We plan to address different research topics, exploring how shared
autonomy, interaction design, and transparency relate to trust and cognitive
load. Another goal is to develop non-invasive methods to continuously measure
trust and cognitive load in the context of disaster response using a
multilevel approach. This project is conducted by multiple academic partners
specializing in artificial intelligence, interaction design, and psychology,
as well as an industrial partner for disaster response equipment and end-users
for framing the project and the experiments in real use-cases.

##  I INTRODUCTION

Assistance robots are useful tools to support emergency forces in disaster
response situations. The robots can be used to collect information, operate
fixtures or manipulate dangerous objects. However, fully autonomous systems
are currently not feasible, both technically and due to a lack of acceptance
by emergency forces [1]. There are multiple reasons for this: firstly,
responders must have an appropriate level of trust in their equipment, which
needs to be highly reliable in the field. Secondly, automated skills are
currently still far from human-like performance in terms of complexity and
reliability. Thus, mainly fully remotely-operated systems with limited skills
are used in the field [2]. The main goal of using such robots is to keep
humans out of dangerous zones while still allowing them to perform their
tasks.

To increase the use of assistance robots among first responders, an
interdisciplinary team of roboticists, psychologists, user-interface
specialists, firefighters, and a renewed vendor of firefighting gears work on
mitigating the above drawbacks through an advanced Shared Autonomy concept and
a robust non-invasive assessment of trust and cognitive load as part of the
research project ”Enabling and Assessing Trust when Cooperating with Robots in
Disaster Response” (EASIER). The addressed use case which sets the frame for
the project is mobile manipulation in disaster response. An industrial-grade
mobile manipulator is used to conduct disaster-related tasks such as
reconnaissance of hazard areas, clearing of entrances from debris, operation
of valves and handles, or manipulation of dangerous goods. The research is
conducted in close cooperation with active firefighters to identify relevant
user requirements and conduct user studies in realistic disaster environments,
answering a recent call for more field studies in HRI (Human-Robot
Interaction) research [3]. Figure 1 shows an example firefighting situation
using the robot platform.

Figure 1: The used mobile platform in a firefighting situation. An additional
manipulator is developed for the project, which is placed on top of the mobile
base.

The research pursues two closely related goals. The first goal is to capture
trust in the robotic system and the cognitive load of using the system, as
both are crucial factors for human-robot interaction [4]. To assess trust and
cognitive load without disturbing the operator during operation, non-invasive
methods need to be developed, preferably through a multilevel approach that
combines the operator’s subjective evaluation with existing information from
the robotic system, the user interface, and information about the user’s
interaction with them (e.g. reaction time, gaze). The need for using non-
invasive methods results from the fact that mitigating a disaster is a very
demanding task and any additional intervention compromises the realistic
interaction. The second goal is to develop an advanced interaction concept
based on the triangle of user interface, shared autonomy level, and
transparency. A general overview of the concept is shown in Figure 2.

Figure 2: Conceptual overview of the proposed Human-Robot Interaction concept
in the area of disaster response.

Using a valid non-invasive measurement technique we can investigate the
effects of variations in the user interface, the shared autonomy level, and
the transparency on trust and cognitive load. For the former, we will adapt
the input and output of the user interface according to the actual task,
situation, as well as estimated trust and cognitive load; these adaptations
may happen offline or online and comprise abstraction or enhancement of the
visualization or changing the control interface. The latter two both relate to
Shared Autonomy as we aim at an optimized and transparent sharing of the
responsibility and task. Based on the trust of the user in the robot’s skills,
the level of autonomy can be adapted from full teleoperation over assistant
functions (similar to driver assistance systems) to full delegation of tasks.
We plan to increase transparency to support the operator’s situational
awareness using the situation awareness-based agent transparency (SAT) model
[5]. Based on the SAT model, we structure the information to enhance
perception, comprehension, and projection of the robot’s state, capabilities,
and decision-making processes, which aims to increase the operator’s trust in
the robot without increasing cognitive load.

Moreover, we envision that the robot is also able to estimate its reliability
and certainty level for various situations to reduce the autonomy level in
case of unforeseen circumstances the robot is not able to handle. Such an
adaptation builds trust without increasing cognitive load only if transparent
explanations are given by the robot. Such an interaction schema is rather
novel in the area of disaster response.

In sum, the main research questions of the presented project are the
following: 1) How can trust and cognitive load be measured in a non-invasive
way and their importance be evaluated in the context of disaster response? and
2) How can interaction design, shared autonomy, and transparency work together
to improve trust and reduce cognitive load in disaster response scenarios?

The remainder of this paper is structured as follows. First, related work on
the topic of shared autonomy, recent projects on disaster response, trust and
cognitive load in robotic systems is presented. Then, the conceptual overview
of the system design for answering the research questions is presented.
Furthermore, integration details on the system concept are discussed and the
next steps for the implementation of the proposed concept are presented.

Figure 3: System design highlighting the different interconnected components.

##  II RELATED WORK

This section discusses related work and developments regarding shared
autonomy, as well as trust and cognitive load in robot operation for disaster
response.

###  II-A Related Projects

Multiple recent projects have tackled the problem of human-robot collaboration
in disaster response.

NIFTi [6] was a project with close interaction with fire brigades, developing
adaptive robots and interfaces taking into account the cognitive load of the
users. The project notably experimented on multiple interfaces for mobile
first responders [7]. Building on NIFTi, TRADR [8] focused on persistent
models for perception, action, and collaboration in long-term missions. SHERPA
[9] was a project with goals relating to multimodal HRI and decision making in
the context of a multi-robot team commanded by one human in alpine scenarios
where the human could be more or less available depending on the situation.
Within S-HELP [10] a psychological framework to enable effective cognitive
processing in the design of emergency management information systems (EMIS)
focusing on decision-making processes, perception and information processing
as well as personal and contextual variables was established. In the current
project TRUST ME [11], an instrument for measuring individual trust in robots
will be developed and validated.

While similar in scope and research areas, these projects did not approach the
issue of assistance robots for disaster response by focusing on how trust and
cognitive load can be assessed using non-invasive methods and how the HRI can
be improved by investigating the effects of different configuration of the
level of autonomy, the interface design, and transparency.

###  II-B Shared Autonomy

While shared autonomy gained increasing research interest in recent years,
different definitions exist which define the term only roughly (see e.g. [12,
13]). While some definitions focus on the technical aspects [14], others
highlight the importance of Human-Robot Interaction [12]. As a result, the
focus in this area relates to (combined) aspects such as the optimization of
user input [14], the recognition of human/robot intentions [15], the
optimization of task distribution [13], the monitoring of automated systems
[16], or the general optimization of interaction [17]. Additionally, research
to design interfaces to mitigate performance issues in robot teleoperation
[18] is carried out. While all these definitions and fields of work do justify
the term of shared autonomy, we see the term of shared autonomy in the
presented project as a general interdisciplinary mean to improve assistance
robots with a focus on aspects such as interface optimization, switching the
level of assistant functionality, and transparent human/robot cooperation.

###  II-C Trust and Cognitive Load in Robot Operation

The operator’s trust in a robot system and the cognitive load caused by using
the system are decisive for the acceptance of an assistance system by rescue
teams [19]. Both, the level of trust and cognitive load depend on human-
related, robot-related, and contextual factors [4]. For example, user studies
showed a positive association between the level of transparency and the
overall trust in the system [20]. Moreover, the robot’s reliability affects
the operator’s trust [21] and trust repair [22]. Operators’ trust and workload
are also influenced by the level of autonomy [23].It is also shown that
explanations regarding robots behaviour can increase trust in robots [24].

Using a non-invasive multilevel approach including subjective data,
psychophysiological measurements, performance variables, reconstructional
video, and observation techniques, we will assess trust and cognitive load
without interfering with the operator during the current mission [25].

##  III System Overview

One important aspect of the project is the design and construction of the
robot. This is done in close collaboration with an industrial partner
experienced in robotic devices for assistance in firefighting, who will
provide the mobile base and the robotic arm and integrate the required
sensors. The robot will have two main skills, the navigation of the mobile
platform and the control of the robotic arm to manipulate objects in the
environment. For navigation, the robot will be equipped with two LiDARs and
with 2 camera arrays composed of a stereo camera and a depth camera to
increase the situational awareness of the robot and to improve the robot
operation. One camera array will be attached to the front of the robot and one
attached to the arm, to provide visual feedback for the autonomy module as
well as to the operator. The cameras are to be used for object detection for
save navigation and manipulation, both either when using teleoperation or in
the assistance mode with shared autonomy. Furthermore, a 360-degree stereo
camera array will be mounted on the back of the robot to improve the
operator’s situational awareness during teleoperation.

As the necessary backbone to start answering the aforementioned research
questions, the system design revolves around 4 components: 1) the user
interface (UI) providing a continuous flux of information to the user, 2) the
robot mentioned above allowing for dynamically switching between levels of
autonomy for multiple skills 3) the central intelligence computing reliability
and certainty levels for different skills, triggering autonomy level switches
and UI adaptation, and 4) the sensors creating a world representation. An
overview is given in Figure 3.

The main workflow is that given a certain state of the world, its observation
leads to a world representation that allows for skills to be performed at
specific levels of autonomy with estimated reliability and certainty levels.
This information is reflected on the user interface, and the user inputs
determine the specific skill and autonomy level to be used. The UI then
provides a specific interface tailored for each skill and autonomy level.

During skill execution, the central intelligence also updates a confidence
metric for the given task. This metric is defined task by task, for example by
estimating a collision probability during navigation [26]. This can be used
both to better inform the user and to trigger or advise for autonomy level
switches, for example when a failure is detected or the task is too complex
for the robot. In parallel to this, explanations for these decisions are also
computed for transparency purposes.

As a first approach, we focus on two different autonomy levels, namely
teleoperation and semi-autonomous operation. Teleoperation will be supported
by assistance functions (e.g. inverse kinematics, obstacle avoidance) and
semi-autonomous operation will make use of different tools to support robot
skills (e.g. waypoint navigation, grasping target selection).

The available skills and their specific interaction functionalities are
reflected on the user interface alongside reliability and certainty
information and a world representation devised from the sensors. The interface
is designed to be used either using traditional displays or virtual reality,
the design process is done alongside end-users to guarantee the best possible
applicability of the system.

Explanations for the decisions made by the central intelligence (e.g.
triggering an autonomy level switch, estimating reliability and certainty for
a task) will also be computed and displayed. One of the goals for future
experiments is to have a better understanding of the amount and content of the
explanations to guarantee enough transparency to build trust without impairing
the usability of the system due to information overload.

The presented concept aims to develop psychological methods and models to
determine the operator’s trust and cognitive load, focusing on a non-invasive
multilevel approach to not impair the control and communication with the robot
in case of a disaster by unavoidable disturbances.

This project is done in cooperation with end-users, with a continuous feedback
loop to improve the methodology and the individual components.

##  IV Current work and Future Plans

As a first step to improve the HRI, a fast method to provide detailed
explanations of failures in robot motion planning was proposed in [27]. By
identifying failures origin from planning constraints, possible failure
explanations can be provided in logarithmic time, also for large constraint
sets. A study to investigate the impact on HRI was also conducted on a
research platform and shows significant improvements in trust and task
performance when explanations of the errors are provided.

In the current state of the project, the next goal is to finish the
integration of the mobile manipulator with its sensors and implement basic
navigation and manipulation skills, as well as providing a first iteration of
the interaction design. Moreover, the concept of the central intelligence and
its interaction with the user, the world model, and the skills needs to be
refined.

The software development will be fully ROS-compatible and includes a simulated
environment to facilitate experiment design, and in a second stage help in
investigating assistance functions and skill-related UI.

Another experiment is planned to measure trust and cognitive load in a
navigation task with varying difficulty (e.g. faulty sensors, different
obstacle setups) and autonomy level switch between teleoperation and semi-
autonomous navigation. We aim at integrating an estimation of motion planning
reliability and certainty for informing the user and triggering autonomy level
switches. The previously mentioned robot sensors will be used, as well as an
overall scene camera and cameras pointing towards the user for further
analysis.

For estimating the operator’s level of trust and cognitive load using non-
invasive methods, data collection studies will be conducted in these
experiments with experienced rescue operators as test participants. Research
question 1) will be addressed using data collected during the experiments and
with posthoc reconstruction. To answer research question 2) user studies
following a multi-level approach to investigate the influences of the defined
levers are to be conducted.

##  V Conclusion

This paper presented the current scope and plans for the EASIER project. The
main objective is to integrate shared autonomy, interaction design, and
transparency to provide a trustworthy system for disaster response scenarios
while limiting cognitive load when using said system. Multiple developments
are expected, from robot design to a better understanding of how trust and
cognitive load can be measured and improved on in the context of disaster
response. To this end, we specifically plan to explore topics relating to
shared autonomy, interaction design, and transparency. The ultimate goal is to
validate these developments in real-world scenarios designed directly with end
users.

## References

  * [1] Sethu Vijayakumar, “Shared autonomy for robots in dynamic environments: Advances in learning control and representations,” Boston, 2019.
  * [2] R. Murphy, _Disaster Robotics_. MIT Press, 2014.
  * [3] A. Weiss, A.-K. Wortmeier, and B. Kubicek, “Cobots in industry 4.0: A roadmap for future practice studies on human–robot collaboration,” _IEEE Transactions on Human-Machine Systems_ , vol. 51, no. 4, pp. 335–345, 2021.
  * [4] M. Ahmad, J. Bernotat, K. Lohan, and F. Eyssel, “Trust and cognitive load during human-robot interaction,” _AAAI Fall Symposium on Artificial Intelligence for Human-Robot Interaction_ , 2019.
  * [5] J. Y. Chen, K. Procci, M. Boyce, J. Wright, A. Garcia, and M. Barnes, “Situation awareness-based agent transparency,” Army research lab aberdeen proving ground md human research and engineering …, Tech. Rep., 2014.
  * [6] G.-J. M. Kruijff, I. Kruijff-Korbayová, S. Keshavdas, B. Larochelle, M. Janíček, F. Colas, M. Liu, F. Pomerleau, R. Siegwart, M. A. Neerincx, _et al._ , “Designing, developing, and deploying systems to support human–robot teams in disaster response,” _Advanced Robotics_ , vol. 28, no. 23, pp. 1547–1570, 2014.
  * [7] B. Larochelle, G.-J. M. Kruijff, N. J. J. M. Smets, and J. van Diggelen, “Experiences with usar mobile interfaces: The need for persistent geo-localized information,” in _2013 IEEE/RSJ International Conference on Intelligent Robots and Systems_. IEEE, 2013, pp. 5333–5338.
  * [8] I. Kruijff-Korbayová, F. Colas, M. Gianni, F. Pirri, J. de Greeff, K. Hindriks, M. Neerincx, P. Ögren, T. Svoboda, and R. Worst, “Tradr project: Long-term human-robot teaming for robot assisted disaster response,” _KI-Künstliche Intelligenz_ , vol. 29, no. 2, pp. 193–201, 2015.
  * [9] L. Marconi, C. Melchiorri, M. Beetz, D. Pangercic, R. Siegwart, S. Leutenegger, R. Carloni, S. Stramigioli, H. Bruyninckx, P. Doherty, _et al._ , “The sherpa project: Smart collaboration between humans and ground-aerial robots for improving rescuing activities in alpine environments,” in _2012 IEEE International Symposium on Safety, Security, and Rescue Robotics (SSRR)_. IEEE, 2012, pp. 1–4.
  * [10] C. M. Steiner, A. Nussbaumer, K. Neville, and D. Albert, “A psychological framework to enable effective cognitive processing in the design of emergency management information systems,” _Electronic Journal of Information Systems Evaluation_ , vol. 20, no. 1, pp. 39–54, 2017.
  * [11] B. Kubicek and R. Esterbauer, “Trust me: Vertrauen in humanoide Roboter im Arbeitskontext,” _Research project funded by the University of Graz_ , 1019\.
  * [12] C. Brooks and D. Szafir, “Balanced information gathering and goal-oriented actions in shared autonomy,” in _2019 14th ACM/IEEE International Conference on Human-Robot Interaction (HRI)_ , 2019, pp. 85–94.
  * [13] S. Hart, J. Kramer, S. Gee, and R. R. Burridge, “The pharaoh procedure execution architecture for autonomous robots or collaborative human-robot teams,” in _2016 25th IEEE International Symposium on Robot and Human Interactive Communication (RO-MAN)_ , 2016, pp. 888–895.
  * [14] S. Nikolaidis, “Mathematical models of adaptation in human-robot collaboration,” Ph.D. dissertation, Carnegie Mellon University, Pittsburgh, PA, 2017.
  * [15] S. Javdani, H. Admoni, S. Pellegrinelli, S. S. Srinivasa, and J. A. Bagnell, “Shared autonomy via hindsight optimization for teleoperation and teaming,” _International Journal of Robotics Research_ , vol. 37, no. 7, pp. 717–742, 2018.
  * [16] P. Marion, M. Fallon, R. Deits, A. Valenzuela, C. P. D’Arpino, G. Izatt, L. Manuelli, M. Antone, H. Dai, T. Koolen, J. Carter, S. Kuindersma, and R. Tedrake, “Director: A user interface designed for robot operation with shared autonomy,” _Springer Tracts in Advanced Robotics_ , vol. 121, pp. 237–270, 2018.
  * [17] S. J. Anderson, J. M. Walker, and K. Iagnemma, “Experimental performance analysis of a homotopy-based shared autonomy framework,” _IEEE Transactions on Human-Machine Systems_ , vol. 44, no. 2, pp. 190–199, 2014.
  * [18] J. Y. Chen, E. C. Haas, and M. J. Barnes, “Human performance issues and user interface design for teleoperated robots,” _IEEE Transactions on Systems, Man, and Cybernetics, Part C (Applications and Reviews)_ , vol. 37, no. 6, pp. 1231–1245, 2007.
  * [19] V. Groom and C. Nass, “Psychological benchmarks of human–robot interaction,” _Interaction Studies. Social Behaviour and Communication in Biological and Artificial Systems_ , vol. 8, no. 3, pp. 483–500, 2007.
  * [20] M. W. Boyce, J. Y. Chen, A. R. Selkowitz, and S. G. Lakhmani, “Effects of agent transparency on operator trust,” in _Proceedings of the Tenth Annual ACMIEEE International Conference on Human-Robot Interaction Extended Abstracts_ , J. A. Adams, Ed. New York, NY: ACM, 2015, pp. 179–180.
  * [21] M. Desai, M. Medvedev, M. Vázquez, S. McSheehy, S. Gadea-Omelchenko, C. Bruggeman, A. Steinfeld, and H. Yanco, “Effects of changing reliability on trust of robot systems,” in _Proceedings of the seventh annual ACMIEEE international conference on Human-Robot Interaction_ , H. Yanco, Ed. New York, NY: ACM, 2012, p. 73.
  * [22] A. L. Baker, E. K. Phillips, D. Ullman, and J. R. Keebler, “Toward an understanding of trust repair in human-robot interaction,” _ACM Transactions on Interactive Intelligent Systems_ , vol. 8, no. 4, pp. 1–30, 2018\.
  * [23] L. Onnasch, C. D. Wickens, H. Li, and D. Manzey, “Human performance consequences of stages and levels of automation: An integrated meta-analysis,” _Human factors_ , vol. 56, no. 3, pp. 476–488, 2014.
  * [24] N. Wang, D. V. Pynadath, E. Rovira, M. J. Barnes, and S. G. Hill, “Is it my looks? or something i said? the impact of explanations, embodiment, and expectations on trust and performance in human-robot teams,” in _International Conference on Persuasive Technology_. Springer, 2018, pp. 56–69.
  * [25] T. R. Colin, N. J. Smets, T. Mioch, and M. A. Neerincx, “Real time modeling of the cognitive load of an urban search and rescue robot operator,” in _2014 RO-MAN: the 23rd IEEE International Symposium on Robot and Human Interactive Communication_. Piscataway, NJ: IEEE, 2014, pp. 874–879.
  * [26] B. Axelrod, L. P. Kaelbling, and T. Lozano-Pérez, “Provably safe robot navigation with obstacle uncertainty,” _The International Journal of Robotics Research_ , vol. 37, no. 13-14, pp. 1760–1774, 2018.
  * [27] M. Eder and G. Steinbauer-Wagner, “A fast method for explanations of failures in optimization-based robot motion planning,” in _Advances in Service and Industrial Robotics_ , A. Müller and M. Brandstötter, Eds. Cham: Springer International Publishing, 2022, pp. 114–121.


${h_{1}\omega_{0}\eta}{\hat{\delta}}^{2}/2$ occurs at the minimum of the
“potential energy”, namely at $\phi=0$ where $U(\phi)=0$. Then

$H_{0}=\frac{h_{1}\omega_{0}\eta}{2}{\hat{\delta}}^{2}=\frac{e\omega_{0}}{2\pi
E_{b}\beta^{2}}U(\hat{\phi}_{-})=\frac{e\omega_{0}}{2\pi
E_{b}\beta^{2}}U(\hat{\phi}_{+})$ (72)

so that $H_{0}$ expresses the amplitude of the synchrotron motion in terms of
$\hat{\delta}$, ${\hat{\phi}}_{-}$ or ${\hat{\phi}}_{+}$.

Eq. (VII.3) can be written as

$\displaystyle\frac{h_{1}\omega_{0}\eta}{2}\delta^{2}=H_{0}-\frac{e\omega_{0}}{2\pi
E_{0}\beta^{2}}U(\phi))$

Then the synchrotron oscillation period can be obtained by using this to
integrate the first equation of Eq. (VII.3):

$T_{z}=\oint
d\phi\left(2h_{1}\omega_{0}\eta\left[H_{0}-\frac{e\omega_{0}}{2\pi
E_{0}\beta^{2}}U(\phi)\right]\right)^{-1/2}$ (73)

This can be numerically evaluated for a specified $H_{0}$, so that we can then
obtain the synchrotron tunes, $\nu_{z}=\frac{2\pi}{T_{z}\omega_{0}}$, plotted
in Fig. 9a. For each $H_{0}$, a measure of the maximum longitudinal distance
of a particle from the center of the bunch is
$z_{H_{0}}\equiv\frac{({\hat{\phi}}_{+}-{\hat{\phi}}_{-})R}{2h_{1}}$ and this,
being more recognisable, is used instead of $H_{0}$ as the horizontal axis in
Fig. 9(a).

In the presence of radiation damping and quantum excitation, the longitudinal
phase-space distribution settles to an equilibrium given by the static
solution of the Fokker-Planck equation [76], namely,
$\tilde{\psi}(\delta,\phi)\propto\exp(-\rm{constant}\times H(\phi,\delta))$.
This distribution is approximately Gaussian in $\delta$ with an rms relative
energy spread $\sigma_{\delta}$, while the distribution of $\phi$ is

$\tilde{\rho}(\phi)=A_{\rho}\exp\left(-\frac{eU(\phi)}{2\pi
E_{0}h_{1}\eta\sigma_{\delta}^{2}}\right)$ (74)

where $A_{\rho}$ is a normalization factor chosen so that
$\int_{\phi_{\textrm{ufp}}}^{\phi_{\textrm{sep}}}\tilde{\rho}(\phi)d\phi=1$.
The first-order and second-order beam moments can be numerically evaluated by

$\displaystyle\bar{\phi}$ $\displaystyle=$
$\displaystyle\int_{\phi_{\textrm{ufp}}}^{\phi_{\textrm{sep}}}\phi\tilde{\rho}(\phi)d\phi$
$\displaystyle\Sigma_{\phi}$ $\displaystyle=$
$\displaystyle\int_{\phi_{\textrm{ufp}}}^{\phi_{\textrm{sep}}}(\phi-\bar{\phi})^{2}\tilde{\rho}(\phi)d\phi$
(75)

The longitudinal distributions, $\tilde{\rho}(\phi)$, are shown in Fig. 9(b),
and are obtained by numerically evaluating the r.h.s. of Eq. (74). The
horizontal axis in Fig. 9(b) is the longitudinal coordinate $z=-\beta c\Delta
t=\frac{\phi R}{h_{1}}$. The rms bunch length is then

$\sigma_{z}=\frac{\sqrt{\Sigma_{\phi}}R}{h_{1}}$ (76)

Note that the second-order beam moments at equilibrium are usually calculated
using the beam-envelope matrix as in [77]. That method is still applicable for
calculation of the rms energy spread, but it fails to deliver the correct rms
bunch length in the presence of a double-RF system, because $U(\phi)$ now
contains components higher than second order in $\phi$, which are not properly
taken into account by that method. In this case, one can use Eq. (75) and Eq.
(76) to numerically evaluate $\sigma_{z}$, or use Monte-Carlo simulation to
obtain $\sigma_{z}$.

In any case, Fig. 15 shows the probability density functions (PDF) of the beam
distributions in $\delta$ and $z$, for Case C at the ideal lengthening
condition, where the theoretical distributions match well with the histograms
from the final beam distribution of a Monte-Carlo simulation with localised
cavities and full synchro-betatron motion tracking with 4000 particles for 10
damping times (about 25000 turns). The beam energy is set to 45.72 GeV
($a\gamma_{0}=103.76$) and $V_{1}=103.372$ MV.

(a) Distribution of relative energy deviation $\delta$.

(b) Distribution of relative longitudinal positron $z$.

Figure 15:  The longitudinal beam distributions of Case C at the ideal
lengthening condition. The red curves are the theoretical probability density
function (PDF), the histograms are obtained from the final beam distribution
of a Monte-Carlo simulation, where 4000 particles are tracked for 10 damping
times (about 25000 turns). The beam energy is set to 45.72 GeV
($a\gamma_{0}=103.76$) and $V_{1}=103.372$ MV.