# Fast and accurate quantum Monte Carlo for molecular crystals

Andrea Zen  Department of Physics and Astronomy, University College London,
Gower Street, London, WC1E 6BT, U.K.  Thomas Young Centre and London Centre
for Nanotechnology, 17–19 Gordon Street, London, WC1H 0AH, U.K.  Jan Gerit
Brandenburg  Department of Physics and Astronomy, University College London,
Gower Street, London, WC1E 6BT, U.K.  Thomas Young Centre and London Centre
for Nanotechnology, 17–19 Gordon Street, London, WC1H 0AH, U.K.  Jiří Klimeš
J. Heyrovský Institute of Physical Chemistry, Academy of Sciences of the Czech
Republic, Dolejškova 3, CZ-18223 Prague 8, Czech Republic and Department of
Chemical Physics and Optics, Faculty of Mathematics and Physics, Charles
University, Ke Karlovu 3, CZ-12116 Prague 2, Czech Republic  Alexandre
Tkatchenko  Physics and Materials Science Research Unit, University of
Luxembourg, L-1511 Luxembourg, Luxembourg  Dario Alfè  Thomas Young Centre and
London Centre for Nanotechnology, 17–19 Gordon Street, London, WC1H 0AH, U.K.
Department of Earth Sciences, University College London, Gower Street, London
WC1E 6BT, U.K.  Angelos Michaelides  angelos.michaelidesucl.ac.uk Department
of Physics and Astronomy, University College London, Gower Street, London,
WC1E 6BT, U.K.  Thomas Young Centre and London Centre for Nanotechnology,
17–19 Gordon Street, London, WC1H 0AH, U.K.

###### Abstract

Computer simulation plays a central role in modern day materials science. The
utility of a given computational approach depends largely on the balance it
provides between accuracy and computational cost. Molecular crystals are a
class of materials of great technological importance which are challenging for
even the most sophisticated _ab initio_ electronic structure theories to
accurately describe. This is partly because they are held together by a
balance of weak intermolecular forces but also because the primitive cells of
molecular crystals are often substantially larger than those of atomic solids.
Here, we demonstrate that diffusion quantum Monte Carlo (DMC) delivers sub-
chemical accuracy for a diverse set of molecular crystals at a surprisingly
moderate computational cost. As such, we anticipate that DMC can play an
important role in understanding and predicting the properties of a large
number of molecular crystals, including those built from relatively large
molecules which are far beyond reach of other high accuracy methods.

##  Significance

Computational approaches based on the fundamental laws of quantum mechanics
are now integral to almost all materials design initiatives in academia and
industry. If computational materials science is genuinely going to deliver on
its promises then an electronic structure method with consistently high
accuracy is urgently needed. We show that, thanks to recent algorithmic
advances and the strategy developed in our manuscript, quantum Monte-Carlo
(QMC) yields extremely accurate predictions for the lattice energies of
materials at a surprisingly modest computational cost. It is thus no longer a
technique that requires a world-leading computational facility in order to
obtain meaningful results. While we focus on molecular crystals the
significance of our findings extends to all classes of materials.

## Introduction

Computer simulations, in particular those based on the fundamental laws of
quantum mechanics, play a key role in modern materials science research. Such
electronic structure approaches serve as an indispensable complement to
experiment by helping to understand and predict the properties of materials as
well as guiding the development of new ones. Density functional theory (DFT)
is the leading electronic structure technique in materials science, thanks to
its favorable balance between accuracy and computational cost as well as
efficient implementation in modern codes. Indeed, DFT has been described as
one of the great success stories of modern science Burke (2012), with
widespread use in materials science and cognate disciplines Curtarolo  _et
al._ (2013); Marzari (2016); Sun  _et al._ (2016). However, DFT as generally
applied with standard exchange-correlation functionals, suffers from a number
of well-known deficiencies Cohen  _et al._ (2008); Peverati and Truhlar
(2014); Medvedev  _et al._ (2017); Hammes-Schiffer (2017). One notable
shortcoming is in the description of weak interactions such as London
dispersion forces. Although the last decade has seen impressive developments
with the incorporation of dispersion forces in the DFT framework Grimme  _et
al._ (2016); Klimeš and Michaelides (2012); Beran (2016); Hermann  _et al._
(2017), and new density functionals with improved accuracy have been developed
Sun  _et al._ (2015); Mardirossian and Head-Gordon (2014); Wang  _et al._
(2017), such approaches are not systematically improvable and their accuracy
for condensed phases is open to question. For example, DFT cannot be relied
upon to routinely deliver an accuracy below 4 kJ/mol for the lattice energy of
molecular crystals. This is the accuracy (so called chemical accuracy) that is
often needed to discriminate between different polymorphs of a given material,
and in some pharmaceutical molecules even more stringent accuracies (down to 1
kJ/mol) are needed Cruz-Cabeza  _et al._ (2015); A. M. Reilly et al. (2016).

Traditionally, highly accurate computations of interaction energies have been
based on quantum chemistry techniques, in particular, coupled cluster with
single, double and perturbative triple excitations (CCSD(T)). However, CCSD(T)
calculations on solids have been notoriously challenging. Great strides
forward have recently been made and CCSD(T) can now be used to study solids
through calculations in periodic boundary conditions Booth  _et al._ (2013),
or by employing approaches based on embedding and fragment decomposition Wen
_et al._ (2012); Bygrave  _et al._ (2012). The benzene crystal, for example,
was recently considered in a _tour de force_ study using fragment
decomposition Yang  _et al._ (2014). The recent introduction of local
approaches Werner and Schütz (2011) promises to extend the range of
applicability of CCSD(T) methods. However, the cost of CCSD(T) calculations
for large systems will remain high for the foreseeable future and their large-
scale application is additionally hindered by enormous memory requirements.
The random phase approximation (RPA) is emerging as a promising approach for
materials, in particular if singles corrections are introduced Schimka  _et
al._ (2010); Ren  _et al._ (2011); Klimeš (2016). Although it is less accurate
than CCSD(T) it is considerably more affordable and currently offers a very
good balance between accuracy and computational cost. Quantum Monte Carlo, in
particular within the fixed node diffusion Monte Carlo scheme Foulkes  _et
al._ (2001), is an established method for reference quality calculations of
molecular systems and condensed phases. Systematic studies in cases of non-
covalent bonding have shown that DMC has an accuracy comparable to CCSD(T)
Dubecký  _et al._ (2016). An advantage of DMC over traditional quantum-
chemical methods like CCSD(T) is that it is essentially unaffected by basis
set issues, thanks to the deployment of an efficient ground state projection
scheme and the use of B-splines. However, the enormous computational cost of
this approach means that DMC studies generally require a world-leading
computational facility and even then it remains highly challenging to obtain
demonstrably converged results. As a result, DMC studies are often restricted
to one-off benchmark calculations of specific systems.

Building upon a recent algorithmic development Zen  _et al._ (2016) which
significantly reduces the computational cost of DMC, we have developed a
scheme (described in the Methods section) for obtaining converged DMC energies
of molecular crystals in a fully periodic treatment, without the need for
fragment decompositions. The approach is effective for both small and large
molecules and yields lattice energies of periodic molecular crystals at a
computational cost comparable to RPA but with the accuracy of CCSD(T). For
example, a chemically accurate lattice energy of a benzene crystal can be
computed in as few as $10^{4}$ CPU hours. In addition, the scheme proposed
comes with relatively low manpower costs and can be seen as a step towards the
automated DMC treatment of molecular crystals. This opens up new horizons for
DMC as a tool for accurately and rapidly predicting the properties of
molecular crystals, including those built from pharmaceutical molecules. More
broadly the insights obtained from this study should also prove beneficial to
DMC studies of other crystalline materials and their surfaces.

## Results

The main quantity to consider in order to assess the stability of a crystal is
its lattice energy $E_{\text{latt}}$, which is the energy per molecule gained
upon assuming the crystal form with respect to the gas state. It can be
computed as

$E_{\textrm{latt}}={E_{\textrm{crys}}}-E_{\textrm{gas}}$ (1)

with $E_{\textrm{crys}}$ the energy per molecule in the crystal state and
$E_{\textrm{gas}}$ the energy of the isolated molecule. While the accurate
computation of $E_{\textrm{gas}}$ is straightforward, $E_{\textrm{crys}}$
involves an ideally infinite system. Since any simulated system will be
necessarily finite, the evaluation of $E_{\textrm{crys}}$ will be affected by
finite size errors (FSE). The most common approach to tackle crystals is to
exploit their periodicity and to simulate a cell within periodic boundary
conditions. In DFT the energy of the infinite system is recovered by
simulating the primitive cell and by sampling the first Brillouin zone. In DMC
the wave function is sampled by considering a number of electron
configurations, including configurations with local dipoles that in a small
simulation cell (such as the primitive cell) couple with the periodic images,
yielding substantial FSE. However, established procedures developed for DMC to
reduce or eliminate these errors are available, such as the model periodic
Coulomb (MPC) interaction Fraser  _et al._ (1996); the approach used in this
study. MPC treats the explicit two-body interactions within the minimum-image
convention in the simulation cell, while outside the cell the interaction is
with a pre-calculated charge density normally obtained from a DFT calculation.
In principle, MPC can miss some of the correlation energy. Thus, in DMC it is
good practice to simulate supercells of increasing size and to extrapolate the
energy of the infinite system. On the other hand, molecular crystals have
primitive cells that are already large (and related to the size of the
molecule) meaning that calculations in the larger cells needed for the
extrapolations can be computationally prohibitive. However, thanks to recent
developments Zen  _et al._ (2016) it is now possible to simulate large
supercells up to cases where FSE are negligible. This, in turn, allows us to
establish that accurate results from DMC can in practice be obtained in small
primitive cells if appropriate corrections for FSE are taken into account.

Typically the computation of $E_{\text{latt}}$ is performed at zero
temperature and considering only the electronic contribution, i.e. quantum
nuclear effects are neglected Beran  _et al._ (2016). The lattice energy is
not directly assessable experimentally, but it can be indirectly obtained from
experimental measures of the sublimation enthalpy $\Delta_{\textrm{sub}}H(T)$
at a given temperature $T$ by including a (theoretically evaluated) energy
contribution $\Delta_{\textrm{T\&QN}}(T)$ accounting for contributions from
thermal and quantum nuclear effects:

$\Delta_{\textrm{sub}}H(T)=-E_{\text{latt}}+\Delta_{\textrm{T\&QN}}(T).$ (2)

The evaluation of $\Delta_{\textrm{T\&QN}}(T)$ can be challenging, especially
for large molecules where anharmonic contributions are important Reilly and
Tkatchenko (2013). In the SI Appendix we provide details about the theoretical
evaluations of $\Delta_{\textrm{T\&QN}}(T)$ and the uncertainty associated
with experimental evaluations of $\Delta_{\textrm{sub}}H$. Since both
$\Delta_{\textrm{sub}}H(T)$ and $\Delta_{\textrm{T\&QN}}(T)$ are affected by
errors, accurate theoretical evaluations of $E_{\text{latt}}$ are of help for
comparison. To assess the accuracy of a method one needs to use a diverse test
set of systems.

Figure 1:  Molecular crystals considered in this work with DMC. Only the primitive cell is shown in each case. The systems treated are of considerable size and contain up to 144 molecules (ice VIII) or 1728 electrons (CO2).  Table 1:  Lattice energy [kJ/mol] for the molecular crystals under consideration in this work, computed with DMC compared to values from experimental measures of sublimation enthalpy. | DMC(lc)${}^{\text{[a]}}$ | DMC(sc)${}^{\text{[b]}}$ | Experiment${}^{\text{[c]}}$
---|---|---|---
Ice Ih | -59.3$\pm$0.5 | -59.2$\pm$0.2 | -58.8
Ice II | -59.1$\pm$0.6 | -59.0$\pm$0.3 | -58.8
Ice VIII | -57.3$\pm$0.6 | -57.4$\pm$0.1 | -57.4
Carbon dioxide | -28.2$\pm$1.3 | -28.5$\pm$0.4 | -28.4
Ammonia | -37.1$\pm$0.4 | -37.5$\pm$0.1 | -37.2
Benzene | -52.1$\pm$0.4 | -51.2$\pm$0.2 | -50.6
Naphtalene | -78.8$\pm$0.8 | -78.0$\pm$0.6 | -79.2
Anthracene | -105.5$\pm$1.7 | -103.9$\pm$1.0 | -105.8
${}^{\text{[a]}}$DMC using a large supercell. ${}^{\text{[b]}}$DMC using a
small supercell, additional values in SI Appendix. ${}^{\text{[c]}}$See SI
Appendix, Sec. S12 for details.

In this study we considered eight molecular crystals (see Fig. 1): carbon
dioxide (CO2), ammonia (NH3), benzene (C6H6), naphthalene (C10H8) and
anthracene (C14H10) crystals from the C21 test set of Otero-de-la-Roza and
Johnson Otero-de-la Roza and Johnson (2012), plus three polymorphs of ice: the
hexagonal ice Ih, ice II and a high-pressure phase ice VIII. This set of
molecular crystals comprises a diversity in intermolecular interactions
involving strong hydrogen bonds and London dispersion of saturated and
unsaturated molecules. A range of interactions such as this is a tough test
for any electronic structure method. It is a test that must be passed through
if a method is to be truly predictive and applicable to complex molecular
crystals, including those of industrial interest Price and Reutzel-Edens
(2016); Pulido  _et al._ (2017). DMC values for $E_{\text{latt}}$ of the eight
molecular crystals are summarized in Table 1. Two sets of DMC results are
reported: DMC(lc) and DMC(sc). The former is obtained using large supercells
(containing around a thousand valence electrons); the latter using smaller
cells and relying on the FSE correction via the MPC interaction (see Methods
and SI Appendix). We find that in all cases the DMC(sc) results are in
excellent agreement with the DMC(lc) results, as a confirmation of the quality
of the MPC approach for the correction of FSE. MPC was introduced two decades
ago Fraser  _et al._ (1996), but here using the algorithm of Ref. Zen  _et
al._ (2016) we have been able to explicitly demonstrate how well MPC performs
for complex systems such as molecular crystals. Table 1 also reports lattice
energies derived from experiments. The experimental values are obviously not
free from error and have an uncertainty coming from both the actual measure of
$\Delta_{\textrm{sub}}H$ and the computed term $\Delta_{\textrm{T\&QN}}$ on
the order of the chemical accuracy, $\sim 4$ kJ/mol. For naphthalene and
anthracene which have the largest values for $E_{\text{latt}}$ the
experimental uncertainties are likely to be larger. See the SI Appendix for a
detailed discussion on the experimental values. Upon comparing DMC to
experiment we find that both DMC(lc) and DMC(sc) always fall within the
accuracy of the experimental value. This is remarkable if we consider that
this accuracy is achieved over a large range of lattice energies,
$E_{\text{latt}}$ from 28 to more than 100 kJ/mol. DMC gets correct lattice
energies for hydrogen bonded, dispersion bonded and mixed bonded crystals. DMC
also predicts the correct relative energies of the ice polymorphs, yielding
slightly improved lattice energies over those reported in ref. Santra  _et
al._ (2011a).

## Discussion

A comparison of the results obtained from DMC to experiment and to other
reference-quality computational approaches is shown in Fig. 2. This includes
MP2 results for all systems Gillan  _et al._ (2013); Cutini  _et al._ (2016);
Wen and Beran (2011), and CCSD(T) for all molecules up to benzene Gillan  _et
al._ (2013); Yang  _et al._ (2014); Wen and Beran (2011). RPA and RPA with
singles corrections (RPA+GWSE) lattice energies for ice are calculated in this
work, and the other values are from Ref. Klimeš (2016). From this comparison
we notice that CCSD(T) and RPA+GWSE perform well, whereas RPA systematically
underbinds all systems and MP2 severely overbinds in systems with delocalized
electrons such as benzene, naphthalene and anthracene. Among the computational
approaches reported in Fig. 2, only CCSD(T) is acknowledged for an accuracy
comparable to DMC, and indeed they show excellent agreement. However, all
CCSD(T) (and MP2) values reported come from fragment decomposition approaches,
which involves the computation of many small contributions to the lattice
energy, all of which must be converged to high accuracy and typically the
correlation contribution from long-range fragments is computed at a lower
level of theory. This can be a painstaking process. Also the range of values
obtained from the widely studied benzene crystal (-50 to -56 kJ/mol) Yang  _et
al._ (2014); Wen and Beran (2011); Podeszwa  _et al._ (2008); Bludský  _et
al._ (2008); Ringer and Sherrill (2008), suggests that the decisions made in
carrying out the fragment decomposition can have a noticeable effect on the
final result. A big advantage of methods employing periodic boundary
conditions, such as DMC, is that the $E_{\text{latt}}$ is obtained from a
single calculation (provided that FSE are corrected for), which makes such
approaches more suitable for rapid screening. As an added bonus, methods such
as DMC also yield information on the electronic structure and electron density
on the full periodic system; information that can be used for the calculation
of experimental observables and to obtain deeper understanding of the system
under consideration.

Figure 2:  Accurate and fast DMC results for a range of molecular crystals.
(Top panel) Difference in the value of $E_{\text{latt}}$ between the
experimental value and several computational approaches often used as
reference methods. Here DMC(lc) and DMC(sc) indicates that large or small
supercells have been employed, respectively. RPA and RPA+GWSE values for ice
have been computed in this work, other values are from Ref. Klimeš (2016). MP2
and CCSD(T) values for ice are from Ref. Gillan  _et al._ (2013), benzene from
Ref. Yang  _et al._ (2014), MP2 values for naphthalene and anthracene from
Ref. Cutini  _et al._ (2016), other values are from Ref. Wen and Beran (2011).
(Bottom panel) Approximate computational cost for DMC(sc), DMC(lc), RPA and
RPA+GWSE (see SI Appendix for details). The DMC cost is intended for a
precision of 0.7 kJ/mol. Reported timings are intended only to provide an
indication; differences in the codes and computation facilities can yield very
different timings

Computational cost is of utmost importance when making comparisons of
computational methods. Whilst DMC(lc) and DMC(sc) produce almost equal values
for $E_{\text{latt}}$, each DMC(sc) is much cheaper than DMC(lc). For example,
as shown in the bottom panel of Fig. 2, DMC(sc) is typically 1-2 orders of
magnitude cheaper. Indeed, most of the DMC(sc) results take around $10^{4}$
core hours, and can be obtained in around a day on a few hundred processors.
This is much more affordable than CCSD(T), which is also only feasible for
relatively small molecules with the fragment decomposition approach or small
crystals in periodic boundary conditions. RPA+GWSE was so far providing a good
compromise between accuracy and computational cost. Fig. 2 shows that the cost
for DMC(sc), for a precision on $E_{\text{latt}}$ of around 1 kJ/mol, is in
general comparable to RPA.

The computational efficiency of the DMC simulations and the fact that we have
periodic boundary conditions makes it relatively straightforward to
investigate other properties beyond the lattice energy. For instance, we have
obtained the equation of state (EOS) for both ammonia and benzene; these are
crystals held together predominantly by hydrogen bonds and dispersion
interactions, respectively. The results of these simulations, along with fits
to the Murnaghan EOS, are reported in Fig. S2. From this we find that the
equilibrium volumes ($V_{0}$) predicted by DMC agree very well with
experiment, coming out $\sim$3% smaller than experiment for both crystals.
Slightly smaller DMC volumes are to be expected since our calculations do not
take into account anharmonic thermal expansion and quantum nuclear effects
present in experiment. The EOS calculations are also useful because they allow
us to test the sensitivity of our computed $E_{\text{latt}}$ to the volume
used in our calculations. The DMC $E_{\text{latt}}$ values listed in Table 1
have been obtained at experimentally measured densities. For the two crystals
reported in Fig. S2, the bias on $E_{\text{latt}}$ arising from the use of the
experimental volume appears very small, less than 0.2 kJ/mol. For the other
crystals reported in Table 1 we expect, on the basis of DFT tests Klimeš
(2016), a bias on $E_{\text{latt}}$ due to the volume on the order of 1 kJ/mol
or less. A second source of bias on the values of $E_{\text{latt}}$ reported
in Table 1 is due to the geometries used for the DMC calculations. Indeed,
DMC, CCSD(T), MP2 and RPA are typically too expensive for a geometry
optimization, which is often performed via DFT with a reliable functional. In
the SI Appendix we report the EOS of ammonia and benzene obtained using the
geometries from two different DFT functionals. The uncertainty on
$E_{\text{latt}}$ appears to be less than 1 kJ/mol both for ammonia and
benzene.

Figure 3:  Equation of state (EOS) for the ammonia and benzene crystals (with
zero-point motion not accounted for). In both cases, the DMC $E_{\text{latt}}$
values are calculated using the $2\times 2\times 2$ cell and FSE are corrected
for with MPC. The dashed line is the Murnaghan EOS fitting the DMC values,
which for ammonia yields a minimum $E_{0}$ of $-37.25\pm 0.05$ kJ/mol at a
volume $V_{0}$ of $32.6\pm 0.1$ Å3 per molecule and a bulk modulus $B_{0}$ of
$8.5\pm 0.5$ GPa, and for benzene $E_{0}$ is $-52.37\pm 0.06$ kJ/mol, $V_{0}$
is $115.7\pm 0.1$ Å3 per molecule and $B_{0}$ is $7.7\pm 0.2$ GPa. The value
of $B_{0}^{\prime}$ in the Murnaghan EOS is set to 4.

To conclude, we have demonstrated that DMC provides a route towards the fast
and accurate determination of the properties of molecular crystals. In
essence, the scheme makes use of the size-consistent DMC algorithm introduced
earlier Zen  _et al._ (2016) and an accurate approach for correcting for
finite size errors. We have applied this approach to a range of exemplar
systems held together with a range of intermolecular interactions (hydrogen
bonds to London dispersion). The calculations have confirmed previous results
on water-ice polymorphs but with minimal computational cost and with much more
control over the numerical accuracy of the results than before. Our results
also include EOS calculations for benzene – the “fruitfly” molecular crystal
in computational materials science – and anthracene, the largest molecule in
the C21 dataset. The consistently high accuracy demonstrated by DMC along with
its moderate computational cost suggests that DMC can play an increasingly
important role in studies of molecular crystals. In particular, DMC could
prove to be the method of choice in challenging polymorph prediction studies.
Similarly, molecules of direct pharmaceutical interest could now be tackled
with DMC; opening up their study with a high-level ab initio approach for the
first time. To provide full phase diagrams for molecular crystals, our
accurate lattice energies have to be combined with estimates of zero-point and
thermal effects. This is traditionally computed at the DFT level, where our
study will further provide an important benchmark to test and calibrate these
approximate methods. Looking further to the future we note that several steps
of the proposed methodology could be used in a full configuration interaction
QMC approach Booth  _et al._ (2013), which would yield essentially exact
solutions to the Schrödinger equation for molecular crystals. Finally, we note
that beyond molecular crystals the improved efficiencies and improved
understanding of finite size errors obtained here will also be of direct
relevance to DMC simulations on other classes of material, e.g. absorption in
metal organic frameworks and binding to surfaces.

## Materials and Methods

Geometries for the C21 crystals are taken from Klimeš (2016) (where the
geometries for molecules and crystals are optimized via DFT using the
optB88-vdW functional Klimeš  _et al._ (2010), and crystals are in the
experimental unit cell). For ice phases we took the geometries used in ref.
Santra  _et al._ (2011a). DMC simulations were carried out with the casino
code Needs  _et al._ (2010) in order to evaluate $E_{\textrm{crys}}$ and
$E_{\textrm{gas}}$. We used Dirac-Fock pseudopotentials Trail and Needs
(2005a, b) with the locality approximation Mitas  _et al._ (1991). The trial
wavefunctions were of the Slater-Jastrow type with single Slater determinants
and the single particle orbitals obtained from DFT-LDA plane-wave calculations
performed with pwscf pws  and re-expanded in terms of B-splines Alfè and
Gillan (2004). The Jastrow factor included electron-electron, electron-nucleus
and electron-electron-nucleus terms. Further details on the wave-function and
the optimization are provided in the SI Appendix, as well as some comparative
tests with the recently introduced correlated electron pseudopotentials Trail
and Needs (2013).

In the computation of $E_{\textrm{crys}}$, periodic boundary conditions are
employed. Simulations with DMC in periodic boundary conditions can be subject
to significant FSE, as previously discussed. In order to assess the converged
value of $E_{\text{latt}}$, for any molecular crystal several simulation cells
were considered as well as twist boundary conditions Lin  _et al._ (2001) for
the smallest cells. This has revealed that the $E_{\text{latt}}$ obtained from
the primitive cell can be overestimated by as much as 300% due to FSE.
However, correction schemes to reduce FSE are available in DMC, such as the
model periodic Coulomb (MPC) interaction Fraser  _et al._ (1996); Williamson
_et al._ (1997); Kent  _et al._ (1999), the correction proposed in ref. Chiesa
_et al._ (2006) and the one in ref. Kwee  _et al._ (2008). We have tested all
of them and observed that MPC provides the best results, as shown in the SI
Appendix. Here we report results obtained exclusively with MPC. A second and
smaller source of FSE in DMC stems from the use of single particle orbitals
obtained from a DFT calculation on a single point in the Brillouin zone
(typically the $\Gamma$-point). This error, called the independent particle
finite size error (IPFSE), can be easily estimated and corrected for by
performing few additional DFT calculations. Further computational details are
reported in the SI Appendix, including the atomic coordinates for each
molecular system studied.

The time step, $\tau$, is a key issue affecting the accuracy of DMC
calculations. In DMC a propagation according to the imaginary time Schrödinger
equation is performed in order to project out the exact ground state from a
trial wave function Foulkes  _et al._ (2001). A time step $\tau$ must be
chosen, keeping in mind that the efficiency of DMC is directly proportional to
$\tau$ but the projection is exact only in the continuous limit $\tau\to 0$.
Thus, $\tau$ has to be small enough to yield converged results, but as large
as possible to make DMC efficient. The time step dependence is system
dependent, so it has to be evaluated on a case by case basis. In periodic
systems this can be computationally very expensive, because each supercell is
possibly affected by the time step differently. As noted, an improved DMC
algorithm Zen  _et al._ (2016) was recently presented. The new algorithm,
denoted ZSGMA from the authors’ initials, gives better convergence with
respect to $\tau$ than the one proposed by Umrigar, Nightingale and Runge
(UNR) Umrigar  _et al._ (1993) which is implemented as standard in DMC codes.
In the evaluation of $E_{\text{latt}}$ an important practical point is the
influence of the simulation cell size on the time step error. This is shown in
Fig. 4 for the example of the ammonia crystal. Specifically, in Fig. 4 the
dependence of $E_{\text{latt}}$ on $\tau$ is shown for a range of different
unit cells. First, it can be seen that with ZSGMA, $E_{\text{latt}}$ exhibits
almost no dependence on $\tau$ for the range of $\tau$ reported. In contrast,
values of $E_{\text{latt}}$ from UNR show a pronounced and non-linear
dependence on $\tau$. This means that for UNR simulations for small values of
$\tau$ (say $\tau\leq 0.001$ au) are required in order to generate a reliable
$\tau\to 0$ extrapolation. Second, and shown here for the first time, we find
that with ZSGMA the time step error on $E_{\text{latt}}$ is independent of the
size of the unit cell (c.f. the 1x1x1, the 2x2x2 and the 3x3x3 cells in Fig.
4). In contrast with UNR the time step error increases significantly as the
size of the simulation cell is increased. This general behavior can be
rationalized by the fact that ZSGMA is (approximatively) size consistent up to
relatively large values of $\tau$, while UNR is size consistent only in the
limit $\tau\to 0$.

In this work we have verified the time step convergence with the ZSGMA
algorithm for each molecular crystal, as reported in the SI Appendix. It is
the larger time step that ZSGMA facilitates and the insensitivity of the time
step error to the size of the cell that enables the converged DMC calculations
on large crystals reported in this study. All results reported in the work are
obtained with the ZSGMA algorithm and a time step that yields a bias of $<1$
kJ/mol. With the UNR algorithm the same accuracy would have required difficult
extrapolations and a computational cost around two orders of magnitude larger.

Figure 4:  Converged lattice energies from DMC in small and large unit cells.
Lattice energies, $E_{\text{latt}}$, for the ammonia crystal, as obtained from
DMC (with the MPC interaction to reduce FSE) by using different time steps and
cell sizes, from the primitive 1x1x1 cell (comprising 4 molecules) to a 3x3x3
supercell (108 molecules). Blue squares represent the results obtained by
employing the algorithm by Umrigar, Nightingale and Runge (UNR) Umrigar  _et
al._ (1993); red circles correspond to the algorithm by Zen  _et al._ (2016)
(ZSGMA), which consistently yields accurate results. The black dashed line is
the value obtained from experimental sublimation enthalpies by Reilly and
Tkatchenko (2013).

We now describe the scheme employed here to compute accurate values of
$E_{\text{latt}}$ with DMC. We recommend the use of the ZSGMA algorithm Zen
_et al._ (2016) in all DMC calculations, and MPC for all DMC calculations in
periodic systems. The following 5 step procedure can be used to assess the
lattice energy for a given molecular crystal.

  * (i)

Geometries – Obtain geometries for the molecular crystal and the isolated
molecule. Since geometry optimizations of large systems are challenging with
DMC, we recommend the use of DFT and an exchange-correlation functional that
accounts for vdW dispersion forces. If reliable experimental structures are
available, the optimization should be performed at the experimental volume.

  * (ii)

Independent particle finite size error (IPFSE) – Using the structure obtained
in (i), converge the energy per molecule in the crystal,
$E_{\text{crys}}^{\text{DFT},\infty}$, using the functional that is used to
obtain the single particle orbitals for the DMC calculations (we generally use
LDA). Convergence is reached by considering $l$x$m$x$n$ Monkhorst-Pack grids
of increasing size. The difference
$E_{\text{crys}}^{\text{DFT},l\text{x}m\text{x}n}$ \-
$E_{\text{crys}}^{\text{DFT},\infty}$=IPFSE${}^{\text{DFT}}_{l\text{x}m\text{x}n}$
provides a good indication of the independent particle contribution to the FSE
in DMC calculations for a $l$x$m$x$n$ supercell, see SI Appendix.

  * (iii)

Jastrow optimization – Take the smallest supercell that is compatible with the
Jastrow factor (typically the Jastrow factor has cutoffs related to the size
of the simulated cell; we suggest to use supercells with the maximum radius of
a sphere inscribed within the Wigner-Seitz cell $>$5 Å) and optimize the
Jastrow factor of the quantum Monte Carlo wave function by minimizing the
variance (or, alternatively, the variational energy). An optional test of the
reliability of the Jastrow can be performed by calculating the DMC binding
energy in a molecular dimer extracted from the crystal, and comparing it with
a reference value obtained from CCSD(T).

  * (iv)

DMC time step – Check the time step dependence either on the cell used in step
(iii) or on the molecular dimer.

  * (v)

Final DMC calculation of $E_{\text{latt}}$ – Take a supercell from step (ii)
with the estimation IPFSE${}^{\text{DFT}}_{l\text{x}m\text{x}n}$ smaller that
10 kJ/mol. Perform the DMC simulation for this crystal using MPC, and the DMC
calculation with open conditions for the molecule. Calculate $E_{\text{latt}}$
and correct for the independent particle FSE using
IPFSE${}^{\text{DFT}}_{l\text{x}m\text{x}n}$. This yields the final DMC(sc)
result. Optionally, consider larger supercells in order to reduce
IPFSE${}^{\text{DFT}}_{l\text{x}m\text{x}n}$ and the MPC correction.

The threshold IPFSE${}^{\text{DFT}}_{l\text{x}m\text{x}n}$ $<$ 10 kJ/mol (in
step v) is motivated by the target accuracy of $\sim$1 kJ/mol and a $\sim$10 %
reliability of the DFT based IPFSE correction. A more accurate alternative to
evaluate the IPFSE is possible (twist averaging) and discussed in the SI
Appendix.

### Supporting Information (SI)

Supporting Information Appendix provides details of the setup for the DMC, RPA
and RPA+GWSE calculations, a discussion of the finite size errors and
additional DMC results. In addition, the computational cost of DMC is
discussed, as well as of RPA and RPA+GWSE. An extended version of Table 1 is
given, and the evaluation of lattice energies from experimental sublimation
enthalpies is discussed. Geometries of the molecular crystals and reference
molecules used for the DMC, RPA and RPA-GWSE calculations are given.

###### Acknowledgements.

A.Z. and A.M. are supported by the European Research Council under the
European Union’s Seventh Framework Program (FP/2007-2013) / ERC Grant
Agreement number 616121 (HeteroIce project). A.Z. and A.M.’s work is also
sponsored by the Air Force Office of Scientific Research, Air Force Material
Command, USAF, under grant number FA8655-12-1-2099. A.M. is also supported by
the Royal Society through a Royal Society Wolfson Research Merit Award. J.G.B
acknowledges support by the Alexander von Humboldt foundation within the
Feodor-Lynen program. J.K. is supported by the European Union’s Horizon 2020
research and innovation program under the Marie Sklodowska-Curie grant
agreement No 658705. We are also grateful for computational resources to
ARCHER, UKCP consortium (EP/ F036884/1), the London Centre for Nanotechnology,
UCL Research Computing, Oak Ridge Leadership Computing Facility (No. DE-
AC05-00OR22725), and IT4Innovations Centre of Excellence
(CZ.1.05/1.1.00/02.0070 and LM2015070).

## Supporting Information for “Fast and accurate quantum Monte Carlo for
molecular crystals”

As supporting information we provide details about the setup for the DMC
calculations in Section S1; a discussion of the finite size errors in DMC in
Section S2; additional results on the DMC calculations performed for carbon
dioxide in Section S3, ammonia in Section S4, benzene in Section S5,
naphthalene in Section LABEL:sec:dmc_C10H8, anthracene in Section
LABEL:sec:dmc_C14H10, and the three ice polymorphs in Section
LABEL:sec:ice_dmc. A discussion on the DMC scaling and computational cost is
reported in Section LABEL:sec:DMCcost. Furthermore, details on the RPA and
RPA+GWSE calculations are given in Section LABEL:sec:RPA. An extended version
of Table I of the main paper is given in Section LABEL:sec:othermethods. We
describe in Section LABEL:sec:ElattEXP how the lattice energy can be obtained
from experimental evaluations of the sublimation enthalpy, and we discuss the
corresponding uncertainty. Finally, we provide the geometries of the molecular
crystals used for the DMC, RPA and RPA-GWSE calculations in Section
LABEL:sec:geoCrys, as well as the geometries of the reference molecules in
Section LABEL:sec:geoMol.

##  S1  DMC setup

The geometries used for the molecular crystals and for the reference molecule
are reported in Sec. LABEL:sec:geoCrys and LABEL:sec:geoMol, respectively. In
all the reported DMC calculations we used a trial wavefunctions of the Slater-
Jastrow type with a single Slater determinants, the Trail-Needs-Dirac-Fock
pseudopotentials (NTDF) Trail and Needs (2005a, b) and the locality
approximation. Single particle orbitals were obtained using DFT with the LDA
functional and a plane-wave cutoff of 600 Ry, re-expanded in terms of
B-splines with the natural grid spacing $a=\pi/G_{\rm max}$, where $G_{\rm
max}$ is the magnitude of the largest plane wave in the expansion. The Jastrow
factor used in the trial wavefunction of the system included a two-body
electron-electron (e-e) term; two-body electron-nucleus (e-n) terms and three-
body electron-electron-nucleus (e-e-n) terms specific for any atom type. The
variational parameters of the Jastrow have been optimized in order to minimize
the variance in the smallest considered simulated cell for each molecular
crystal, and the same Jastrow has been used for larger supercells and the
isolated molecule. The size of the simulation cell impose some constraints of
the Jastrow variational freedom, in the form of cut-offs in the e-n, e-e and
e-e-n terms. For this reason, the primitive cell could be too small and a
larger simulation cell should be considered in order to do not impose a
suboptimal Jastrow. In particular, the smallest simulation cell that we
considered is the primitive cell for the ice polymorphs, ammonia and benzene;
the $2\times 2\times 2$ for the carbon dioxide and the $1\times 2\times 1$ for
naphthalene and anthracene. The Jastrow parametrization that we used for ice
is the same as already used previously in Ref. Santra  _et al._ , 2011a; for
naphthalene and anthracene we used the Jastrow optimized for the benzene
crystals in the primitive cell. The properties (e.g., variance, the time step)
used for the calculations reported in Table I of the main paper are given in
Table LABEL:tab:crysprop, and an estimate of the computational cost for each
of the molecular crystals with different supercells are reported in Table
LABEL:tab:Supercells.

The chosen setup can surely influence the absolute DMC energy value, but in
this work we are interested in energy differences, so a big error cancellation
is expected and slightly different choices should not yield drastic
differences. For instance, in Fig. S1 we compare the values of the binding
energy of the ammonia dimer, obtained using the locality approximation and the
TNDF pseudopotentials (i.e., the set up used for all the other systems) with
those obtained using the T-move Casula (2006); Casula  _et al._ (2010) and the
recently introduced correlated electron pseudopotentials (CEPP) Trail and
Needs (2013). The latter choice shows a slightly larger time step dependance,
but for a timestep $\tau\leq 0.01$ au the results are in good agreement.

Figure S1:  Binding energy for the ammonia dimer calculated via DMC as a
function of the DMC time steps. Results are obtained using locality
approximation and Trail-Needs-Dirac-Fock pseudopotentials (TNDF) Trail and
Needs (2005a, b) (red circles), or T-move and correlated electron
pseudopotentials (CEPP) Trail and Needs (2013) (blue squares). We report as a
black dashed line the reference CCSD(T) value.

The geometries used for the isolated molecule and the molecular crystal can
also influence the value of $E_{\text{latt}}$, as discussed in the main paper.
In Fig. S2 we report the equation of state (EOS) for ammonia and benzene, as
obtained with DMC calculations using DFT geometries from the optB88-vdW and
the HSE-3c functionals. The fitted parameters of the Murnaghan EOS for the two
sets of geometries are in quite good agreement.

Figure S2:  EOS for ammonia and benzene, evaluated using DMC calculations on the $2\times 2\times 2$ cell, FSE corrected via MPC, and geometries from DFT optimizations with the optB88-vdW and the HSE-3c functionals. Dashed lines are obtained fitting the DMC values with the Murnaghan EOS, where the $B_{0}^{\prime}$ term (the derivative of the bulk modulus with the pressure) is set to a value of 4.  Table S1:  Results from DMC and CCSD(T) at complete basis set (CBS), for molecular dimers. The geometry of the dimers is extracted from the molecular crystals. Zero-point effects are not included.  | Carbon dioxide | Ammonia | Benzene | Naphthalene | Anthracene
---|---|---|---|---|---
| CO2 | NH3 | C6H6 | C10H8 | C14H10
DMC($\tau$=0.1) | -14.6 $\pm$ 0.2 | -11.1 $\pm$ 0.2 |  |  |
DMC($\tau$=0.03) | -6.9 $\pm$ 0.3 | -10.2 $\pm$ 0.1 | -10.6 $\pm$ 0.2 | -19.5 $\pm$ 0.6 | -29.7 $\pm$ 0.8
DMC($\tau$=0.01) | -5.4 $\pm$ 0.2 | -9.9 $\pm$ 0.1 | -9.9 $\pm$ 0.3 | -18.6 $\pm$ 0.5 | -30.0 $\pm$ 0.8
DMC($\tau$=0.003) | -4.6 $\pm$ 0.4 | -10.0 $\pm$ 0.2 |  |  |
DMC($\tau\to 0$) | -4.2 $\pm$ 0.2 |  |  |  |
CCSD(T)a | -4.1 | -9.6 | -9.9 $\pm$ 0.3 | -19.1 $\pm$ 0.4 | -29.8 $\pm$ 0.6
a In all the systems, chemical core is kept fixed in the correlation
treatment, no relativistic are included; CO2 is obtained at the CCSD(T)/CP/CBS
level; NH3 at CCSD(T)-F12/CBS(cc-pVTZ-F12,cc-pVQZ)/CP; C6H6, C10H8 and C14H10
at L-CCSD(T) interaction energy with TZ, QZ extrapolation.

In order to test the reliability of all the choices made, system by system, we
extracted from each molecular crystal a dimer, and we compared the DMC results
from our approach, and for different choices of the time step $\tau$, with the
predictions of CCSD(T). Results are provided in Table S1. From the results in
the dimer, we have an indication that in order to have a target accuracy of
around 1 kJ/mol, we can use a time step of $\tau=0.03$ au for ammonia,
benzene, naphthalene and anthracene, while in carbon dioxide there is a strong
time step dependance, thus we have performed calculations at different values
of $\tau$, down to 0.003 au, and we have extrapolated for $\tau\to 0$.
Additional tests on the DMC time step have also been performed on these
molecular crystals, as reported in Sec. S3–LABEL:sec:dmc_C14H10, confirming
the reliability of the over mentioned choices. In the three ice polymorphs we
have used $\tau=0.003$ au, as discussed in Sec. LABEL:sec:ice_dmc. However, we
first discuss in Sec. S2 details about the FSE corrections.

##  S2  Finite size errors in DMC

In the computation of $E_{\textrm{crys}}$, periodic boundary conditions (PBC)
are employed and the value $E_{\textrm{crys}}^{N_{\textrm{mol}}}$ computed is
obtained from a cell with only $N_{\textrm{mol}}$ molecules. The quantity that
we need is the thermodynamic limit for $E_{\textrm{crys}}^{\infty}$, obtained
for $N_{\textrm{mol}}\to\infty$, that is by using an infinitively large
simulation cell. The difference
$E_{\textrm{crys}}^{\infty}-E_{\textrm{crys}}^{N_{\textrm{mol}}}$ is the
finite size error (FSE). FSEs in DMC are quite tricky, but they have been
studied in detail and are well characterized Drummond  _et al._ (2008);
Holzmann  _et al._ (2016). They essentially arise from an independent particle
contribution (IPFSE) due to quantization of momentum, and a spurious Coulomb
interaction (CIFSE) between periodic images.

Methods employing independent orbitals (even effectively, such as using the
Kohn-Sham orbitals in DFT) often allows to calculate the energy (per electron,
per atom or per primitive cell) by performing a reduction of the problem to
the primitive cell. Many-body approaches, such as QMC, do not allow this
reduction, because correlations can have long ranges, thus there is often the
need to use large simulation cells. Assuming that we are simulating a
supercells having $N$ electrons, with positions $\left\\{{\bf
r}_{1},\ldots,{\bf r}_{N}\right\\}$, and that ${\bf R}_{p}$ are the primitive-
cell lattice vectors, then the Hamiltonian has to satisfy the following
translational symmetry:

$\hat{H}({\bf r}_{1}+{\bf R}_{p},\ldots,{\bf r}_{i}+{\bf R}_{p},\ldots,{\bf
r}_{N}+{\bf R}_{p})=\hat{H}({\bf r}_{1},\ldots,{\bf r}_{i},\ldots,{\bf
r}_{N})$ (1)

which leads to the many-body Bloch condition:

$\Psi_{{\bf k}_{p}}({\bf r}_{1},\ldots,{\bf r}_{N})=W_{{\bf k}_{p}}({\bf
r}_{1},\ldots,{\bf r}_{N})\;{\rm e}^{i\,{\bf k}_{p}\cdot{1\over N}\sum_{i}{\bf
r}_{i}}$ (2)

where $W$ has the same translational symmetry of the Hamiltonian in Eq. 1, and
the Bloch momentum ${{\bf k}_{p}}$ can be restricted to the Brillouin zone of
the primitive cell. On top of that, we need to impose periodic boundary
conditions (PBC) across the simulations cell, yielding the following
additional translational symmetry for the Hamiltonian:

$\hat{H}({\bf r}_{1},\ldots,{\bf r}_{i}+{\bf R}_{s},\ldots,{\bf
r}_{N})=\hat{H}({\bf r}_{1},\ldots,{\bf r}_{i},\ldots,{\bf
r}_{N})\quad\textrm{for each }i=1,\ldots,N$ (3)

where ${\bf R}_{s}$ are the simulation cell lattice vectors. Notice that the
symmetry in Eq. 3 is an artifact of the PBC, which leads to the many-body
Bloch condition:

$\Psi_{{\bf k}_{s}}({\bf r}_{1},\ldots,{\bf r}_{N})=U_{{\bf k}_{s}}({\bf
r}_{1},\ldots,{\bf r}_{N})\;{\rm e}^{i\,{\bf k}_{s}\cdot\sum_{i}{\bf r}_{i}}$
(4)

where $U$ has the periodicity of the simulation cell for each electron, and
the simulation cell Bloch momentum ${\bf k}_{s}$, also called twist vector,
can be restricted to the Brillouin zone for the simulation-cell lattice. The
employment of ${\bf k}_{s}$ vectors other than zero is ofter referred as the
application of twisted boundary conditions Lin  _et al._ (2001), and when the
expectation value is taken as the average over many vectors ${\bf k}_{s}$ in
the Brillouin zone of a gives simulation cell, we talk about twist averaging
boundary condition (TABC) 111This process can become tricky for metals, but in
this work we are only interested to insulators, for which there are no
expected difficulties.. The fluctuations of the QMC energy due to the use of
different twists is almost proportional to the corresponding fluctuations
observed by performing Hartree-Fock or DFT calculations at the same ${\bf
k}_{s}$ point. For instance, this can be seen for the specific case of ammonia
crystal in the lower panel of Fig. S3. Moreover, in insulators it is often
possible to find a special point in the Brillouin zone Rajagopal  _et al._
(1994, 1995); Dagrada  _et al._ (2016), such as the Baldereschi point
Baldereschi (1973), which provides an energy very close to the one obtained by
averaging over many ${\bf k}_{s}$ points in the Brillouin zone. For the
ammonia crystal, this aspect is clear in the lower panel in Fig. S3.

Figure S3:  (Upper left panel) Difference between the lattice energy
$E_{\text{latt}}$ obtained with the a finite simulation cell, with respect to
the extrapolated lattice energy $E_{\text{latt}}^{\infty}$ for the infinite
system, as obtained from a DFT calculation with LDA functional, plotted as a
function of the simulation cell size (where $1\times 1\times 1$ indicated the
primitive cell). We show the convergence for three special points in the
Brillouin zone (notice that the system is a simple cubic with lattice
$a\sim$5.13 Å): the $\Gamma$-point corresponding to ${\bf k}_{s}=(0,0,0)$, the
corner point R corresponding to ${\bf k}_{s}={\pi\over a}(1,1,1)$, and the
Baldereschi point ${\bf k}_{s}={\pi\over 2a}(1,1,1)$. (Upper right panel) Same
as the left panel, but for DMC calculations with MPC. Notice that the $\Gamma$
and R points correspond to twists that make the QMC wave function real (see
Eq. 4), so we have used the fixed-node approximation. On the other hand, the
Baldereschi point makes the QMC wave function complex and we used the fixed-
phase approximation. (Lower panel) scatter plot of the lattice energy obtained
using DMC with MPC versus DFT with LDA, using only the primitive cell and 320
random twists. The special twists corresponding to the $\Gamma$, R and
Baldereschi points are highlighted. All the twists (excluded the $\Gamma$ and
R points) correspond to a complex wave function and thus employ the fixed-
phase approximation. The big gray cross indicates $E_{\text{latt}}^{\infty}$
for the two approaches. The red fitting line has an angular coefficient of
0.860$\pm$0.003.

It has to be noticed that in general the QMC wave function of a periodic
system is a complex function, according to Eq. 4. Indeed, there are a finite
set of points in the Brillouin zone which make $\Psi_{{\bf k}_{s}}$ a real
wave function, such as the $\Gamma$-point ${\bf k}_{s}=(0,0,0)$, or the corner
points, or the center of edge or of face points, but a general $\Psi_{{\bf
k}_{s}}$ yield a complex wave function. Thus, it has to be remembered that
with complex wave functions the fixed-node (FN) approximation Anderson (1976);
Reynolds  _et al._ (1982) cannot be used and has to be substituted by the
fixed-phase (FP) approximation Ortiz  _et al._ (1993). The two approximations
are strictly related but not equivalent. Thus, whenever we use a complex wave
function, in the calculation of the lattice energy
$E_{\text{latt}}={E_{\textrm{crys}}}-E_{\textrm{gas}}$, also
$E_{\textrm{gas}}$ needs to be computed with the FP approximation. In order to
do so, we have employed a big simulation cell with the single molecule into
the cell, such that the interaction between the molecule and the periodic
images is negligible, and performed a DMC/MPC calculation with PBC and at a
twist ${\bf k}_{s}=\epsilon(1,1,1)$, with $\epsilon\sim 10^{-6}$, in order to
use the FP approximation also for the reference energy. In this work we have
employed TABC with complex wave functions only for the cases of water and
ammonia crystals. In these systems the have that, if there is any difference
in $E_{\textrm{gas}}$ due to the employment of FP versus FN approximation, it
is smaller than the statistical accuracy of the results. However, this cannot
be assumed in general.

Earlier we mentioned that DMC can have large FSEs in periodic systems, and
that the error can be split into a CIFSE and a IPFSE. At this point we can
clarify what are the contributions of the CIFSE and IPFSE terms. Given a
finite supercell, the IPFSE is the error that we make by taking only one ${\bf
k}_{s}$ point (for instance, the $\Gamma$-point) instead of performing TABC.
The remaining of the FSE of due to the CIFSE. Fig. S3 gives the impression
that the IPFSE is the largest part FSE. In fact, it is the other way around:
CIFSE is much larger, one to two orders of magnitude, than the IPFSE. The
feeling in Fig. S3 that there is no CIFSE is due to the fact that the DMC
calculations are employing the model periodic Coulomb (MPC) model, which is a
very effective method recover from the CIFSE. If we would have used the
standard Ewald summation, the lattice energies for the primitive cell would
have been severely overestimated, as it can be seen for the ammonia crystal in
Fig. S4.

Figure S4:  $E_{\text{latt}}$ for the ammonia crystal versus the inverse of
the number $N_{\textrm{mol}}$ of molecules in the supercell. The reported DMC
values are obtained via bare Ewald interaction, KZK correction Kwee  _et al._
(2008), CCMH correction Chiesa  _et al._ (2006), and the MPC interaction
Fraser  _et al._ (1996); Williamson  _et al._ (1997); Kent  _et al._ (1999).
Primitive and $2\times 2\times 2$ cells are at the Baldereschi point; $3\times
3\times 3$ cell is at the R-point. Clearly the MPC interaction yields accurate
results for all simulation cells considered.

The source of CIFSE can be explained in different ways (see Refs Drummond  _et
al._ (2008); Holzmann  _et al._ (2016)). An intuitive explanation is that the
sampling of configurations with electrons that are always in the same relative
position in the periodic images of the simulation cell will induce a spurious
Coulomb interaction between images (e.g., dipole-dipole). Large supercells
decorrelate the electron configurations and eliminate the issue, however this
is not a practical solution because the convergence is slow and the
calculation becomes very expensive. There are at least three complementary
ways available to deal with the CIFSE: the employment in the simulation of a
model periodic Coulomb (MPC) potential for the long-range interaction Fraser
_et al._ (1996); Williamson  _et al._ (1997); Kent  _et al._ (1999) which
eliminates the spurious Coulomb interaction between images; or a posteriori
corrections to standard DMC simulations employing Ewald summation, as proposed
by Chiesa  _et al._ (2006) (CCMH) or by Kwee  _et al._ (2008) (KZK). KZK is
computationally the fastest but the least accurate; it employs the usage of a
special DFT functional that mimics the FSEs of QMC. It is of usage for an
initial estimate of the order of magnitude of the FSE in a specific system,
before any DMC calculation. MPC and CCMH are almost equivalent, as it is shown
theoretically in Ref. Drummond  _et al._ , 2008, and computationally they
involve only a little overhead in the DMC simulation. In Fig. S4 we show the
behavior of the three correction schemes and of the standard (Ewald)
simulation as a function of the supercell size for the representative case of
ammonia crystal. In the primitive cell the Ewald estimate is three times as
large as the converged value, and also in the $2\times 2\times 2$ the FSE is
much larger than the chemical accuracy. KZK overcorrects, by around 40%, the
FSE. Both MPC and CCMH provide very effective corrections and the are within 1
kJ/mol already in the $2\times 2\times 2$, whilst in the primitive cell MPC is
slightly more accurate than CCMH. We have reported in the following paragraphs
the results from the three corrections for any of the studied molecular
crystals. In general, MPC appears as the most accurate method.

##  S3 DMC results for Carbon Dioxide

Table S2 reports the lattice energies values obtained using DMC with Ewald,
KZK, CCMH or MPC approaches. The values show a strong dependence on $\tau$,
thus we also performed the extrapolation for $\tau\to 0$. FSEs for Ewald and
KZK are quite large for the $2\times 2\times 2$ simulation cell, and are
larger than 1 kJ/mol also for a $3\times 3\times 3$ simulation cell. On the
other hand, both CCMH and MPC appears already very close for a $2\times
2\times 2$ simulation cell. Concerning the IPFSE, a DFT with LDA calculation
shows that the expected error on $2\times 2\times 2$ simulation cell at the
$\Gamma$-point is $<0.1$ kJ/mol, thus negligible, and it is even smaller for
the $3\times 3\times 3$ simulation cell.

Table S2:  $E_{\text{latt}}$ of Carbon Dioxide, in kJ/mol, from DMC results
using Ewald, KZK, CCMH or MPC. $\sigma$ indicated the associated stochastic
error. Values in boldface have been used for Table I of the main paper.
$2\times 2\times 2$ simulation cell, $\Gamma$-point
---
$\tau$ | Ewald | KZK | CCMH | MPC | $\sigma$
0.100 | -50.5 | -34.0 | -39.4 | -38.9 | 0.3
0.030 | -42.8 | -26.3 | -31.7 | -31.5 | 0.3
0.010 | -40.6 | -24.0 | -29.5 | -29.3 | 0.3
0.003 | -40.5 | -24.0 | -29.4 | -29.2 | 0.4
$\to 0$ | -39.8 | -23.3 | -28.7 | -28.5 | 0.4
$3\times 3\times 3$ simulation cell, $\Gamma$-point
$\tau$ | Ewald | KZK | CCMH | MPC | $\sigma$
0.030 | -34.9 | -30.0 | -31.5 | -31.5 | 0.7
0.003 | -31.9 | -27.0 | -28.5 | -28.5 | 1.2
$\to 0$ | -31.6 | -26.6 | -28.2 | -28.2 | 1.3

##  S4 DMC results for Ammonia

The results for the ammonia crystal are shown in Figs. 1 and 2 of the main
paper, and in Figs. S3 and S4. The results reported in Table I of the main
paper are obtained with DMC/MPC and a $\tau=$0.03 au; DMC(lc) is from a
$3\times 3\times 3$ simulation cell at the R-point (for which the IPFSE is
negligible); DMC(sc) is from the primitive cell at the Baldereschi point.

##  S5 DMC results for Benzene

The results for the benzene crystal are reported in Table S3, using Ewald,
KZK, CCMH or MPC. The IPFSE associated to each simulation cell is also
reported. It has been estimated via DFT/LDA calculation, as the difference
IPFSE${}^{\textrm{cell},{\bf k_{s}}}=E_{\textrm{crys}}^{\textrm{cell},{\bf
k_{s}}}-E_{\textrm{crys}}^{\textrm{BZ}}$ between the DMC/LDA energy per
molecule $E_{\textrm{crys}}^{\textrm{cell},{\bf k_{s}}}$ calculated at the
${\bf k_{s}}$-point of the Brillouin zone for the specified simulation cell
minus the DMC/LDA energy per molecule $E_{\textrm{crys}}^{\textrm{BZ}}$ for
the fully sampled Brillouin zone (which is equivalent to
$E_{\textrm{crys}}^{\infty}$ in DFT). Thus, in order to correct both for the
CIFSE and the IPFSE, we have to subtract from the MPC value the IPFSE
contribution.

Table S3:  $E_{\text{latt}}$ of benzene, in kJ/mol, from DMC with $\tau=$0.03 au and using Ewald, KZK, CCMH or MPC. $\sigma$ indicated the associated stochastic error. The IPFSE has been estimated via a DFT/LDA calculation, and it should be subtracted to any $E_{\text{latt}}$ value in the row in order to correct for that, as reported in the last column. Values in boldface have been used for Table I of the main paper.  Cell | Ewald | KZK | CCMH | MPC | $\sigma$ | IPFSE | DMC-FSE
---|---|---|---|---|---|---|---
$1\times 1\times 1(\Gamma)$ | -143.1 | -18.7 | -52.1 | -52.4 | 0.2 | -1.2 | -51.2$\pm$0.2
$1\times 1\times 1$(R) | -134.7 | -10.5 | -43.6 | -44.3 | 0.2 | +7.0 | -51.3$\pm$0.2
$2\times 2\times 2(\Gamma)$ | -63.7 | -48.2 | -51.3 | -52.1 | 0.4 | 0.0 | -52.1$\pm$0.4
$2\times 2\times 2$(R) | -64.1 | -48.6 | -51.6 | -52.5 | 0.6 | 0.0 | -52.5$\pm$0.6
$3\times 2\times 3$(R) | -58.5 | -51.6 | -53.7 | -54.0 | 4.4 | 0.0 | -54.0$\pm$4.4

In addition to the test with the dimer, in Table S1, we also tested the size-
consistency error. In order to do that, we considered a benzene dimer, with
the two benzene molecules far away, such that their residual interaction is
negligible. The size-consistency error, defined as the residual energy per
molecule (i.e., a half of the energy of the far away dimer minus the the
energy of the benzene molecule) is reported in Fig. S5 as a function of the
DMC time step, both for the UNR and the ZSGMA algorithms. It is clear that the
ZSGMA algorithm is more accurate, and the error is smaller than 0.5 kJ/mol for
$\tau<0.03$ au, while in the UNR algorithm a $\tau<0.003$ au is still too
large. The results reported in Table I of the main paper are obtained with
DMC/MPC and a $\tau=$0.03 au; DMC(lc) is from a $2\times 2\times 2$ simulation
cell at the $\Gamma$-point; DMC(sc) is from the primitive cell at the
$\Gamma$-point.

Figure S5:  Estimation of DMC time step error in $E_{\text{latt}}$ for
benzene, obtained from benzene dimer with the two molecules far away (distance
$>12$ Å).

## References

  * Burke (2012) K. Burke, J. Chem. Phys. 136, 150901 (2012).
  * Curtarolo  _et al._ (2013) S. Curtarolo, G. L. W. Hart, M. B. Nardelli, N. Mingo, S. Sanvito, and O. Levy, Nat. Mat. 12, 191 (2013).
  * Marzari (2016) N. Marzari, Nat. Mat. 15, 381 (2016).
  * Sun  _et al._ (2016) J. Sun, R. C. Remsing, Y. Zhang, Z. Sun, A. Ruzsinszky, H. Peng, Z. Yang, A. Paul, U. Waghmare, X. Wu, M. L. Klein, and J. P. Perdew, Nat. Chem. 8, 831 (2016).
  * Cohen  _et al._ (2008) A. J. Cohen, P. Mori-Sanchez, and W. Yang, Science (New York, NY) 321, 792 (2008).
  * Peverati and Truhlar (2014) R. Peverati and D. G. Truhlar, Phil. Trans. R. Soc. A 372, 20120476 (2014).
  * Medvedev  _et al._ (2017) M. G. Medvedev, I. S. Bushmarinov, J. Sun, J. P. Perdew, and K. A. Lyssenko, Science 355, 49 (2017).
  * Hammes-Schiffer (2017) S. Hammes-Schiffer, Science (New York, NY) 355, 28 (2017).
  * Grimme  _et al._ (2016) S. Grimme, A. Hansen, J. G. Brandenburg, and C. Bannwarth, Chem. Rev. 116, 5105 (2016).
  * Klimeš and Michaelides (2012) J. Klimeš and A. Michaelides, J. Chem. Phys. 137, 120901 (2012).
  * Beran (2016) G. J. O. Beran, Chem. Rev. 116, 5567 (2016).
  * Hermann  _et al._ (2017) J. Hermann, R. A. DiStasio, Jr., and A. Tkatchenko, Chem. Rev. 117, 4714 (2017).
  * Sun  _et al._ (2015) J. Sun, A. Ruzsinszky, and J. P. Perdew, Phys. Rev. Lett. 115, 036402 (2015).
  * Mardirossian and Head-Gordon (2014) N. Mardirossian and M. Head-Gordon, Phys. Chem. Chem. Phys. 16, 9904 (2014).
  * Wang  _et al._ (2017) Y. Wang, X. Jin, H. S. Yu, D. G. Truhlar, and X. He, Proc Natl Acad Sci USA 114, 8487 (2017).
  * Cruz-Cabeza  _et al._ (2015) A. J. Cruz-Cabeza, S. M. Reutzel-Edens, and J. Bernstein, Chem. Soc. Rev. 44, 8619 (2015).
  * A. M. Reilly et al. (2016) A. M. Reilly et al., Acta Cryst. B 72, 439 (2016).
  * Booth  _et al._ (2013) G. H. Booth, A. Grüneis, G. Kresse, and A. Alavi, Nature 493, 365 (2013).
  * Wen  _et al._ (2012) S. Wen, K. Nanda, Y. Huang, and G. J. O. Beran, Phys. Chem. Chem. Phys. 14, 7578 (2012).
  * Bygrave  _et al._ (2012) P. J. Bygrave, N. L. Allan, and F. R. Manby, J. Chem. Phys. 137, 164102 (2012).
  * Yang  _et al._ (2014) J. Yang, W. Hu, D. Usvyat, D. Matthews, M. Schutz, and G. K. L. Chan, Science (New York, NY) 345, 640 (2014).
  * Werner and Schütz (2011) H.-J. Werner and M. Schütz, J. Chem. Phys. 135, 144116 (2011).
  * Schimka  _et al._ (2010) L. Schimka, J. Harl, A. Stroppa, A. Grüneis, M. Marsman, F. Mittendorfer, and G. Kresse, Nat. Mater. 9, 741 (2010).
  * Ren  _et al._ (2011) X. Ren, A. Tkatchenko, P. Rinke, and M. Scheffler, Phys. Rev. Lett. 106, 153003 (2011).
  * Klimeš (2016) J. Klimeš, J. Chem. Phys. 145, 094506 (2016).
  * Foulkes  _et al._ (2001) W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal, Rev. Mod. Phys. 73, 33 (2001).
  * Dubecký  _et al._ (2016) M. Dubecký, L. Mitas, and P. Jurečka, Chem. Rev. 116, 5188 (2016).
  * Zen  _et al._ (2016) A. Zen, S. Sorella, M. J. Gillan, A. Michaelides, and D. Alfè, Phys. Rev. B 93, 241118(R) (2016).
  * Fraser  _et al._ (1996) L. M. Fraser, W. M. C. Foulkes, G. Rajagopal, R. J. Needs, S. D. Kenny, and A. J. Williamson, Phys. Rev. B 53, 1814 (1996).
  * Beran  _et al._ (2016) G. J. O. Beran, J. D. Hartman, and Y. N. Heit, Acc. Chem. Res. 49, 2501 (2016).
  * Reilly and Tkatchenko (2013) A. M. Reilly and A. Tkatchenko, J. Chem. Phys. 139, 024705 (2013).
  * Otero-de-la Roza and Johnson (2012) A. Otero-de-la Roza and E. R. Johnson, J. Chem. Phys. 137, 054103 (2012).
  * Price and Reutzel-Edens (2016) S. L. Price and S. M. Reutzel-Edens, Drug Discov. Today 21, 912 (2016).
  * Pulido  _et al._ (2017) A. Pulido, L. Chen, T. Kaczorowski, D. Holden, M. A. Little, S. Y. Chong, B. J. Slater, D. P. McMahon, B. Bonillo, C. J. Stackhouse, A. Stephenson, C. M. Kane, R. Clowes, T. Hasell, A. I. Cooper, and G. M. Day, Nature 543, 657 (2017).
  * Santra  _et al._ (2011a) B. Santra, J. Klimeš, D. Alfè, A. Tkatchenko, B. Slater, A. Michaelides, R. Car, and M. Scheffler, Phys. Rev. Lett. 107, 185701 (2011a).
  * Gillan  _et al._ (2013) M. J. Gillan, D. Alfè, P. J. Bygrave, C. R. Taylor, and F. R. Manby, J. Chem. Phys. 139, 114101 (2013).
  * Cutini  _et al._ (2016) M. Cutini, B. Civalleri, M. Corno, R. Orlando, J. G. Brandenburg, L. Maschio, and P. Ugliengo, J. Chem. Theory Comput. 12, 3340 (2016).
  * Wen and Beran (2011) S. Wen and G. J. O. Beran, J. Chem. Theory Comput. 7, 3733 (2011).
  * Podeszwa  _et al._ (2008) R. Podeszwa, B. M. Rice, and K. Szalewicz, Phys. Rev. Lett. 101, 115503 (2008).
  * Bludský  _et al._ (2008) O. Bludský, M. Rubeš, and P. Soldán, Phys. Rev. B 77, 092103 (2008).
  * Ringer and Sherrill (2008) A. L. Ringer and C. D. Sherrill, Chem. Eur. J. 14, 2542 (2008).
  * Klimeš  _et al._ (2010) J. Klimeš, D. R. Bowler, and A. Michaelides, J. Phys.: Cond. Mat. 22, 022201 (2010).
  * Needs  _et al._ (2010) R. J. Needs, M. D. Towler, N. D. Drummond, and P. L. Rios, J. Phys.: Condens. Matter 22, 023201 (2010).
  * Trail and Needs (2005a) J. R. Trail and R. J. Needs, J. Chem. Phys. 122, 014112 (2005a).
  * Trail and Needs (2005b) J. R. Trail and R. J. Needs, J. Chem. Phys. 122, 174109 (2005b).
  * Mitas  _et al._ (1991) L. Mitas, E. L. Shirley, and D. M. Ceperley, J. Chem. Phys. 95, 3467 (1991).
  * (47) S. Baroni, A. Dal Corso, S. de Gironcoli, and P. Giannozzi, http://www.pwscf.org.
  * Alfè and Gillan (2004) D. Alfè and M. J. Gillan, Phys. Rev. B 70, 161101 (2004).
  * Trail and Needs (2013) J. R. Trail and R. J. Needs, J. Chem. Phys. 139, 014101 (2013).
  * Lin  _et al._ (2001) C. Lin, F. H. Zong, and D. M. Ceperley, Phys. Rev. E 64, 016702 (2001).
  * Williamson  _et al._ (1997) A. J. Williamson, G. Rajagopal, R. J. Needs, L. M. Fraser, W. M. C. Foulkes, Y. Wang, and M.-Y. Chou, Phys. Rev. B 55, R4851 (1997).
  * Kent  _et al._ (1999) P. R. C. Kent, R. Q. Hood, A. J. Williamson, R. J. Needs, W. M. C. Foulkes, and G. Rajagopal, Phys. Rev. B 59, 1917 (1999).
  * Chiesa  _et al._ (2006) S. Chiesa, D. M. Ceperley, R. M. Martin, and M. Holzmann, Phys. Rev. Lett. 97, 076404 (2006).
  * Kwee  _et al._ (2008) H. Kwee, S. Zhang, and H. Krakauer, Phys. Rev. Lett. 100, 126404 (2008).
  * Umrigar  _et al._ (1993) C. J. Umrigar, M. P. Nightingale, and K. J. Runge, J. Chem. Phys. 99, 2865 (1993).
  * Casula (2006) M. Casula, Physical Review B 74, 161102 (2006).
  * Casula  _et al._ (2010) M. Casula, S. Moroni, S. Sorella, and C. Filippi, J. Chem. Phys. 132, 154113 (2010).
  * Drummond  _et al._ (2008) N. D. Drummond, R. J. Needs, A. Sorouri, and W. M. C. Foulkes, Phys. Rev. B 78, 125106 (2008).
  * Holzmann  _et al._ (2016) M. Holzmann, R. C. Clay, M. A. Morales, N. M. Tubman, D. M. Ceperley, and C. Pierleoni, Phys. Rev. B 94, 035126 (2016).
  * Note (1) This process can become tricky for metals, but in this work we are only interested to insulators, for which there are no expected difficulties.
  * Rajagopal  _et al._ (1994) G. Rajagopal, R. J. Needs, S. Kenny, W. M. C. Foulkes, and A. James, Phys Rev Lett 73, 1959 (1994).
  * Rajagopal  _et al._ (1995) G. Rajagopal, R. J. Needs, A. James, S. D. Kenny, and W. M. C. Foulkes, Phys. Rev. B 51, 10591 (1995).
  * Dagrada  _et al._ (2016) M. Dagrada, S. Karakuzu, V. L. Vildosola, and M. Casula, Phys. Rev. B 94, 245108 (2016).
  * Baldereschi (1973) A. Baldereschi, Phys. Rev. B 7, 5212 (1973).
  * Anderson (1976) J. B. Anderson, J. Chem. Phys. 65, 4121 (1976).
  * Reynolds  _et al._ (1982) P. J. Reynolds, D. M. Ceperley, B. J. Alder, and W. A. Lester, J. Chem. Phys. 77, 5593 (1982).
  * Ortiz  _et al._ (1993) G. Ortiz, D. M. Ceperley, and R. M. Martin, Phys. Rev. Lett. 71, 2777 (1993).
  * Note (2) For instance, the primitive cell of the ammonia crystal has 4 molecules, each having 8 electrons, for a total of 32, but if FSE is large we have to take the $2\times 2\times 2$ supercell, having $4\times 2^{3}=32$ molecules, 256 electrons; if FSE are still non negligible we need the $3\times 3\times 3$ supercell, having 108 molecules, 864 electrons, and so on. The situation is even worse with large molecules, where already a $2\times 2\times 2$ supercell may have over a thousand electrons.
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * Kaltak  _et al._ (2014a) M. Kaltak, J. Klimeš, and G. Kresse, Phys. Rev. B 90, 054115 (2014a).
  * Klimeš  _et al._ (2015) J. Klimeš, M. Kaltak, E. Maggio, and G. Kresse, J. Chem. Phys. 143, 102816 (2015).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Kaltak  _et al._ (2014b) M. Kaltak, J. Klimeš, and G. Kresse, J. Chem. Theo. Comput. 10, 2498 (2014b).
  * Rozzi  _et al._ (2006) C. A. Rozzi, D. Varsano, A. Marini, E. K. U. Gross, and A. Rubio, Phys. Rev. B 73, 205119 (2006).
  * Harl and Kresse (2008) J. Harl and G. Kresse, Phys. Rev. B 77, 045136 (2008).
  * Macher  _et al._ (2014) M. Macher, J. Klimeš, C. Franchini, and G. Kresse, J. Chem. Phys. 140, 084502 (2014).
  * Roux  _et al._ (2008) M. V. Roux, M. Temprado, J. S. Chickos, and Y. Nagano, J. Phys. Chem. Ref. Data 37, 1855 (2008).
  * Acree and Chickos (2010) W. Acree, Jr. and J. S. Chickos, J. Phys. Chem. Ref. Data 39, 043101 (2010).
  * Růžička  _et al._ (2014) K. Růžička, M. Fulem, and C. Červinka, J. Chem. Thermodyn. 68, 40 (2014).
  * Chickos (2003) J. S. Chickos, Netsu Sokutei (2003).
  * Whalley (1984) E. Whalley, J. Chem. Phys. 81, 4087 (1984).
  * (83) Thermodynamics Research Center, Selected Values of Properties of Chemical Compounds., Thermodynamics Research Center, Texas A&M University, College Station, Texas, 1997.
  * Ahlberg  _et al._ (1937) J. E. Ahlberg, E. R. Blanchard, and W. O. Lundberg, J. Chem. Phys. 5, 539 (1937).
  * Woolley (1954) H. W. Woolley, J. Res. NBS 52, 289 (1954).
  * Manzheli (1971) V. G. Manzheli, Phys. Stat. Sol. 44, 39 (1971).
  * Popov  _et al._ (1971) V. A. Popov, V. G. Manzhelii, and M. I. Bagatskii, J. Low Temp. Phys. 5, 4 (1971).
  * Goursot  _et al._ (1970) P. Goursot, H. L. Girdhar, and E. F. Westrum, J. Phys. Chem. 74, 2538 (1970).
  * McCullough  _et al._ (1957) P. McCullough, H. L. Finke, J. F. Messerly, S. S. Todd, T. C. Kincheloe, and G. Waddington, J. Phys. Chem. 61, 1105 (1957).
  * Santra  _et al._ (2013) B. Santra, J. Klimeš, A. Tkatchenko, D. Alfè, B. Slater, A. Michaelides, R. Car, and M. Scheffler, J. Chem. Phys. 139 (2013).
  * Santra  _et al._ (2011b) B. Santra, J. Klimeš, D. Alfè, A. Tkatchenko, B. Slater, A. Michaelides, R. Car, and M. Scheffler, Phys. Rev. Lett. 107, 185701 (2011b).