Location of Artinite (Mg 2 CO 3 (OH) 2 ·3H 2 O) within the MgO-CO 2 -H 2 O System using Ab initio Thermodynamics

,


Introduction
Magnesium hydrates and carbonates within the MgO-CO 2 -H 2 O system have a wide variety of applications including catalysis [1][2][3][4] , the treatment of industrial carcinogenic heavy-metal waste 5 , the immobilisation of nuclear waste 6 , fire retardation 7 , construction 8 , high-magnesian biogenic calcites 9 , and carbon-mineralization reactions for permanent storage of anthropogenic carbon dioxide. [10][11][12][13] Magnesium is also one of the most abundant elements in the Earth's mantle 14 , and Mg-rich minerals are thought to exist in significant quantities on dwarf planets 15 and other extra-terrestrial bodies 16 . Measuring phase stability is essential to predicting which material phase will form under a given set of conditions, and hence for understanding and controlling the structure and resultant properties. This problem has drawn intensive research 17 from both experiment and theory and across many fields of organic 18 and inorganic chemistry [19][20][21][22][23][24][25] . The stability of different phases is conceptually described by phase diagrams, which capture the relationships between equilibrium and metastable phases, the activities of component species, and external conditions including pressure, volume and temperature, in an elegant way that can be used to predict the nature and composition of a material. 26,27 Whereas the major species in the MgO-CO 2 -H 2 O phase diagram such as magnesia and magnesite are relatively well studied, minor phases have received comparatively less attention. One such example is the Artinite phase, which was first observed by Brugnatelli in 1902 28 and is a naturally-occurring metastable phase formed at low temperature. The crystal structure of Artinite has been analysed using X-ray diffraction 29 and the thermodynamic properties have been determined using calorimetry. 30 A few studies have also examined the thermodynamic stability of Artinite at low temperature, [31][32][33] but comparatively little is known on the stability of Artinite with respect to other phases in the MgO-CO 2 -H 2 O system. Given that the materials in this system are found under a wide range of conditions, both natural and manmade, a better understanding of where this phase sits within the MgO-CO 2 -H 2 O phase diagram is potentially important to a range of disciplines. For prediction and generation of phase diagrams, the "brute force" option is to use so-called global exploration methods to sample the potential-energy landscape, but this approach quickly becomes intractable as the number of variables (i.e. phases) increases. [34][35][36][37] When data on experimentally-identified materials and their structures is available, there are alternative, cheaper methodologies (namely, using computer simulation) that can be employed whose capabilities are often overlooked. 21,[38][39][40][41][42] The most accurate and reliable way of calculating phase diagrams, i.e. the current "gold standard" is to employ a consistent and accurate ab initio methodology to all the phases, so that the resulting predictions are independent of experimental thermodynamic data. It is however critical that the choice of methodology is appropriate to the system being studied, i.e. transferable across different material phases, which may display significant differences in chemical bonding, but also feasible given the available computational power. Transferability can be a particular challenge for systems such as MgO-CO 2 -H 2 O that display multiple types of bonding across the phase space. For example, the bonding between Mg and different ligands (e.g. hydroxyl groups, water molecules and carbonate anions) is strongly ionic, whereas the intramolecular bonding within the ligands is predominantly covalent, and some of the layered magnesium hydroxides may also display interlayer van der Waals dispersion forces. In this work, we compare a series of different approaches, all without the inclusion of experimentally derived corrections, to construct first-principles phase diagrams of the MgO-CO 2 -H 2 O in order to study the thermodynamics of the under-explored Artinite phase. To account for the hierarchy of chemical bonding, we compare predictions made using density-functional theory with the popular PBE generalised-gradient approximation functional and three van der Waals corrected functionals used successfully in the literature for these materials. [43][44][45][46] We also assess the impact of using the quasi-harmonic approximation to derive the temperature dependence of the Gibbs free energy. This comprehensive study allows us to compare the predictions made from ab initio thermodynamic models with different approximations, and to assess the accuracy and reliability of the different flavours of DFT for studying the hydrated and carbonated phases of magnesian and other mineral systems.

Methodology
The VASP code [47][48][49] was used to carry out DFT calculations on the seven structures shown in Figure 1. These structures are known, experimentally reported phases within the MgO-CO 2 -H 2 O system, and are discussed in detail in the Supporting Information. The valence wavefunctions were described by a plane-wave basis and the core electrons were described using projected augmentedwave pseudopotentials. 50 As our baseline, we used the Perdew-Burke-Ernzerhof 51 generalised-gradient approximation (PBE GGA) to model electron exchange and correlation. We also considered three different functionals including van der Waals forces, viz. PBE with the DFT-D3 correction (PBE-D3), 52 optB86b-vdW 53,54 and optB88-vdW. 53 These approaches are improvements over the wellestablished PBE-D2, 55 where a dispersion correction is computed from the atomic polarizabilities and added to the PBE total energy. In PBE-D3, the polarizability is modulated by the local environments of the atoms, and for optB86b-vdW and optB88-vdW the correction is implemented as a modification to the exchange correlation functional itself. 56 The plane-wave cut-off energy was set to 500 eV. A 4 4 4 Monkhorst-Pack k-point mesh was used × × to integrate the electronic Brillouin zone for all phases apart from Brucite where a 6 6 4 k-point mesh was used. These settings × × were chosen based on explicit convergence tests on the total energies. During the electronic wavefunction optimisations the total energies were converged to 10 −6 eV. For geometry optimisations the relaxation of the atomic structures was deemed to have converged when the forces acting on the atoms were below 0.01 eV/Å.  The vibrational frequencies and elastic-constant matrices were calculated using the finite-displacement routines implemented in VASP, the latter of which were used to calculate the bulk moduli by averaging the first nine elastic constants (c 11 :c 33 ). [57][58][59] For quasi-harmonic approximation (QHA) calculations, harmonic phonon calculations were performed using the Phonopy code 60 on a set of 10 expansions and contractions around the equilibrium volume (approx. ±5%). The forces were determined in a 2 2 2 × × supercell for all phases apart from Brucite, for which we used a 3 × 3 2 expansion, with VASP as the force calculator. For calculating × the thermal properties, the Brillouin zone was integrated by sampling with the following regular Γ-centred q-point meshes, which were converged individually for each structure: 30 30  The phase stability plots where generated using SurfinPy 2 61,62 . We note here several limitations in the procedure used here, viz.: (1) we do not account for the possibility of the coexistence of several phases, miscibility gaps, and the possible existence of different polymorphs (unless explicitly considered); (2) we similarly do not account for the possibility of unknown phases that may arise at high pressure; and (3) we do not account for the anharmonic effects that may play an important role at elevated temperature. it is not feasible to replicate this using the small unit cell used in our DFT calculations, and we have therefore constructed two models with different, fixed distributions ( Figure 2). We suggest that the most facile means to properly account for the random distribution of these groups would be classical Monte Carlo and energy-minimization techniques based on force fields, but this is beyond the scope of the present study.    Physical properties for the two models including the predicted lattice parameters and elastic constants are presented in Table 1.

Results and discussion
The predictions for both Artinite structures made using all four DFT techniques give generally excellent agreement with available experimental measurements, with the notable exception of a 6 % underestimation of the enthalpy of formation using PBE. This highlights the importance of accurately describing vdW forces. For both Artinite structures the DFT methods including dispersion corrections predict a decrease in the cell volume compared to experiments. We also performed similar calculations on the other phases shown in Figure 1, and structural data and predictions of the physical properties are given in the Supporting Information (Table  S1-S4).

Phase stability plots of Mg-Rich Phases I: Theoretical
Framework. To derive our phase stability plots, we assume that the Mg-rich phases in Figure 1 can be formed through the reaction scheme in Equation 1: The system is in equilibrium when the chemical potentials of the reactants and products in Equation 1 are equal; i.e. when the change in Gibbs free energy is δG T,p = 0. ( Assuming that H 2 O and CO 2 are gaseous species, µ CO2 and µ H2O can be written as The chemical potentials are the partial molar free energies of 0 the reactants or products in their standard states. For an ideal gas, ∆µ x = RT ln(p x /p θ ) where p θ is the standard pressure of 1 bar. As both MgO and the products in Equation 1 are solid phases, we assume that: (5) MgO = 0 MgO and (6) = 0 Substituting Equations 3-6 into Equation 2, we obtain: A corresponds to the partial molar free energy of product A, we can replace the left-hand side of Equation 7 with the Gibbs free energy ( ): At equilibrium δG T,p = 0, and hence: In the case of an ideal gas: (10) We can therefore find the values of ∆µ H2O and ∆µ CO2 (or (p H2O ) y and (p CO2 ) z ) for which the Mg-rich phases are in thermodynamic equilibrium, i.e. when they are more or less stable than MgO. This procedure can be applied to all the solid phases to identify which is the most stable, provided that the free energy is known for Δ 0 each Mg-rich phase. These free energies can be calculated using Equation 11: (11) Δ 0 = ∑Δ 0, products -∑ΔG 0, reactants In this work, we calculate the free energy of each component from ab initio modelling, i.e. we do not rely on experimental thermodynamic data for any of the solid phases. This is important to provide an independent verification of experimental data and to make predictions where suitable data is not available. This approach is particularly suitable for use with structural searches, which could potentially identify previously-unreported phases. We used two approaches for calculating the formation energies of the solid phases, viz. the harmonic approximation (HA) and the quasiharmonic approximation (QHA). In the harmonic model, the free energies are calculated from: (12) ( ) = 0 + phonon ( ) = 0 + ZPE + vib ( ) where U 0 is the calculated internal energy from a DFT calculation, U ZPE is the zero point energy and A vib (T) is the temperaturedependent component of the vibrational Helmholtz free energy: Both sums run over the 3 vibrational modes, where is the number of atoms in the unit cell, with characteristic vibrational temperatures θ i (equivalent to the frequencies in Kelvin): For these harmonic calculations, we make two approximations: (a) we only consider the frequencies at the Brillouin zone centre, and we therefore neglect the acoustic modes; and (b) we assume a fixed volume, allowing the Gibbs and Helmholtz energies to be equated. In the QHA model, the harmonic approximation is modified to include the effect of volume changes due to thermal expansion at finite temperature. This is practically achieved by applying the harmonic approximation to a series of contracted and expanded unit cells about the equilibrium volume V 0 . Using this method, the Gibbs free energy is obtained at a specific temperature T and pressure p as: 65 (16) where U(V) is the lattice energy as a function of volume, A phonon (T;V) is the corresponding phonon contribution to the Helmoltz free energy, and the function is minimised with respect to volume for each value of T and p. This approximation has been shown to be appropriate up to around 2/3 of the melting temperature, beyond which anharmonic effects begin to dominate. 66,67 For gaseous species, the standard free energy varies significantly with temperature, and as plane-wave DFT simulations are designed for condensed-phase systems we use experimental data to determine the temperature-dependent free energy term for the gaseous species: where S expt (T) is specific entropy value for a given T and is obtained from the NIST database. 68 Provided the free energies of solid phases can be calculated, the phase boundaries corresponding to phase transitions can be identified using Equation 17, and we reiterate that the present method does not require any form of experimental corrections for the solid phases. 43,69 3.3 Phase stability plots of Mg-rich phases II: Free Energy at T = 0 K. We first verify our methodology against an experimental phase stability plot composed of the well-characterised Mg-rich phases. Figure 3 (a) shows the phase stability plots of Brucite (Mg(OH) 2 ), Magnesite (MgCO 3 ) and Periclase (MgO) obtained from experiments, while Figure 3 (b) shows the phase stability plot containing all the Mg-rich phases predicted using optB86b-vdW. These diagrams are predicted at T = 0 K, i.e. the A vib terms in Equation 12 are set to zero, which simplifies Equation 17 to: (18) = 0 + ZPE The phase stability plots calculated using the other functionals are included in Supporting Information Figure S1 and S2. The ∆µ H2O and ∆µ CO2 are varied from -2 to +1 eV, which assuming ideal-gas behaviour are equivalent to pressures of 10 −25 bar and 10 27 bar at 298 K. This chemical potential range is in line with previous computational work, 70 and is chosen to illustrate the difference in the free energies of the key phases. For the experimental phase stability plot in Figure 3(a), is Δ 0 calculated using the corrected enthalpies of formation defined in Equation 19, where ∆H expt (298 K → 0K) is the change in enthalpy from 298 K to 0 K. Thermodynamic data was taken from the NIST database, 68 but experimental values are only available for these well studied phases.
The computed phase stability plot in Figure 3(b) agrees relatively well with the experimental phase stability plot in Figure 3(a), with small constant shifts in the values of p H2O and p CO2 . This may indicate that a correction factor can be applied to each technique to better reproduce the experimental data. This first test for a thermodynamic modelling strategy is to evaluate the relative stability of Mg(OH) 2 , MgO and MgCO 3 following the procedure outlined in Equation 18. The triple point (where all three phases co-exist; Table 2) was calculated from the experimental data to be at ∆µ CO2 = −1.18 eV and ∆µ H2O = −0.83 eV, whereas using optB86b we predict the point to occur at ∆µ CO2 = −1.27 eV and ∆µ H2O = −0.66 eV. This is excellent agreement, with differences of only 0.08 eV for ∆µ CO2 and 0.17 eV for ∆µ H2O . Due to the limited literature data on Lansfordite (MgCO 3 ·5H 2 O) and Artinite (Mg 2 CO 3 (OH) 2 ·3H 2 O), we are not able to include these phases on the experimental phase stability plot. However, at p CO2 = p H2O = 1 bar at 298 K (where the dashed white lines cross in Figure  3(b)), the most thermodynamically stable phase is predicted to be Magnesite (MgCO 3 ), whereas at p CO2 = p H2O = 1 bar at 0 K (where the solid white lines cross in Figure 3 (b)) the predicted stable phase is Lansfordite.  As noted above, one of our aims is to evaluate the thermodynamic and structural data neglecting experimental data for the condensed phases, not least as comparison with experiment is then a rigorous check on the reliability. Nevertheless, inclusion of experimentally-derived corrections does represent a valuable tool for demonstrating the viability of DFT to the system concerned as shown in Chaka et al. 43 Comparing the performance of the different DFT methods, we find that optB86b-vdW gives the most accurate representation of the thermodynamics as the triple point between the Brucite (Mg(OH) 2 ), Magnesite (MgCO 3 ) and Periclase (MgO) phases is closest to the experimental one. Comparative data for the other three DFT techniques tested is given in the SI, Figures S1 and S2.

Phase stability plots of Mg-Rich Phases III: Harmonic and
Quasi-Harmonic Free Energies. Accounting for the temperature dependence of the free energy in the phase stability plot, in particularly the entropy term, is important, as the most thermodynamically stable phase is likely to vary with temperature. We choose a temperature range between 273 K and 373 K, i.e. where liquid water is present, which are appropriate conditions for industrial 6 and nuclear waste reprocessing 5 applications, for construction materials, 8,72 and for catalysis. [2][3][4] The vibrational entropy of the Mg-rich phases is included in the calculation of Δ 0 (Equation 17) for ∆µ = -2 to +1 eV, with entropies for the solid phases calculated using Equations 13 to 15 and entropies for gaseous species taken from the NIST database. 68 We showed in the previous sections that optB86b-vdW provides a good description of the experimental phase stability plot of the well-characterised Mg-rich phases, and that the predicted phase stability plots obtained using the different DFT methods are all shifted relative to each other. We therefore consider here only the results obtained with optB86b-vdW, but as before the phase stability plots predicted using the other DFT techniques are available in Figures S3-S5 for comparison.
We organise the results in this section into three parts, showing first the free energies of Mg-rich phases at different p CO2 as a function of temperature at constant p H2O = 1 bar, followed by the phase stability plots for p CO2 as a function of p H2O at a fixed T = 298 K and finally p CO2 as a function of temperature at constant p H2O . Figure 4 compares the formation free energies of the different phases as a function of temperature, computed using optB86b-vdW with the harmonic and quasi-harmonic free energies, over the 273-373 K temperature range and a variety of p CO2 and p H2O . These phase stability plots consider the effect of varying p CO2 and thus provide insight into how this variable impacts the synthesis of the Mg-rich phases, that is, we can identify all the phases appearing at different values of p CO2 , which enables us to resolve very small energy differences between phases. At ∆µ CO2 = −1 eV (low CO 2 pressure, p CO2 = 10 −17 bar) and ∆µ H2O = 0 eV (p H2O = 1 bar) the most thermodynamically stable phase is predicted to be Brucite (Mg(OH) 2 ) (Figure 4a-4b). Using optB86b and the QHA free energies, Artinite shows no phase transitions at low CO 2 pressure, identified by intersections of its free-energy curve with other compounds, over the temperature range modelled. Increasing the temperature stabilises phases that contain less water per Mg 2+ ion, due to the evaporative loss of H 2 O. The harmonic free energies predict that at ∆µ CO2 = ∆µ H2O = 0 eV, corresponding to standard conditions with p CO2 = p H2O = 1 bar, MgCO 3 is the most thermodynamically stable phase at all temperatures (Figure 4c), 43,71 whereas the QHA predicts that Lansfordite becomes the most stable phase at 277 K (Figure 4d) (Figure 4c). These are shifted compared to the 332 and 320 K predicted using the QHA free energies (Figure 4d). We also compared the phase stability plots under atmospheric conditions for the different DFT techniques (Figures S3-S5), which all show very similar features to the phase stability plots computed under standard conditions. Temperatures of the pahse transitions are tabulated in Tables S5-S7. At ∆µ CO2 = 0.6 eV (high CO 2 pressure, p CO2 = 10 10 bar) and ∆µ H2O = 0 eV (p H2O = 1 bar) the most thermodynamically stable phase at all temperatures is predicted to be MgCO 3 (Figure 4e), as with p CO2 = p H2O = 1 bar (Figure 4e-4f). When the formation of MgCO 3 is inhibited, at T < 290 K Lansfordite (MgCO 3 ·5H 2 O) is predicted to be the second most stable phase, with Nesquehonite (MgCO 3 ·3 H 2 O) favoured between 308-316 K, and Mg 5 (CO 3 ) 4 (OH) 2 ·4H 2 O favoured only at T > 316 K. At high p CO3 Artinite has one phase boundary at 364 K (Artinite Mg 2 CO 3 (OH) 2 ·3H 2 O/Lansfordite MgCO 3 ·5H 2 O). Increasing the temperature results in phase transitions associated with the loss of one H 2 O per MgO unit. Again, we see a similar trend using the QHA where the phase transition temperatures are shifted ~5-10 K higher compared to those computed using the harmonic free energies. Figure 5a shows that with p CO2 = p H2O = 1 bar the most thermodynamically stable phase is MgCO 3 . In general, all the DFT techniques show excellent agreement, with small shifts in the predicted phase-transition pressures (see Supporting Information Figure S6). As ∆µ CO2 is increased, phases with a higher carbon content are thermodynamically stabilised, whereas increasing ∆µ H2O stabilises the hydrated phases. As the temperature is raised from 0 to 298 K (Figure 3(b)-5(b)), the QHA predicts that Artinite appears at the triple point between Brucite Mg(OH) 2 , Magnesite MgCO 3 and Lansfordite MgCO 3 ·5H 2 O. We can also compare the predicted pressures of the triple point using the harmonic and QHA free energies. Values for the Artinite triple points obtained from the QHA are given in Table 3. Testing the validity of these phase stability plots using experimental techniques is challenging, as changing the pressure of one reference state independently can be difficult to achieve in practice. One method is to put MgCO 3 under vacuum and decrease the pressure of both CO 2 and H 2 O simultaneously, which would probe a diagonal line on the phase stability plot. The pressure at which MgCO 3 transforms to MgO could then be determined and compared to the modelling results. Alternatively, these phase stability plots could be adjusted to represent temperature as a function of p CO2 , which we consider in the following section.  Table 3: Predicted triple points for the phase stability plots in Figure 5, where three phases coexist, using optB86b and the harmonic approximation (HA) or the quasi-harmonic approximation (QHA).   76,77 As magnesite will be the predominant phase at most CO 2 pressures in the T vs. p CO2 phase stability plot, inhibiting its growth produces a variety of new phases. The appearance of these phases has repercussions for the sequestration of CO 2 and shows that metastable carbonatecontaining phases are more important than the anhydrous magnesium carbonate (magnesite). More generally, this has implications for applications where the evolution over time of phases in the MgO-H 2 O-CO 2 system is an important factor, for example during the storage of industrial 6 or nuclear 5 waste materials, in construction, 8,72 where Mg binders have to retain stability over a long period of time, or during catalytic cycles where the formation of metastable phases may lead to poisoning. [2][3][4] Multiple phase transitions may occur during the process that ultimately leads to MgCO 3 formation. Indeed, Ostwald's rule 78 states that nucleation from solution occurs in steps, where the most thermodynamically unstable phases form first, followed by a step towards the most thermodynamically-stable phase. While this rule does not always hold true, 79 a larger number of phases increases the likelihood that Ostwald's rule applies. 79 Figure 6 compares our calculated phase stability plots as a function of temperature and ∆µ CO2 at a fixed ∆µ H2O = 0 eV (p H2O = 1 bar), computed using the harmonic and quasi-harmonic free energies, and either including all the Mg-rich phases or with the formation of MgCO 3 inhibited. Comparable phase stability plots obtained using the other functionals can be found in Figures S7 and S8.  43 who obtained a transition temperature of ~340 K but countered this by applying an experimental correction to the free energy of MgCO 3 ·5H 2 O. This is possibly an indication of the shortcomings of the description of van der Waals forces in DFT in general, and of the treatment in the DFT+D2 method specifically. Please do not adjust margins Please do not adjust margins There is debate in the literature about the thermodynamic stability of Artinite Mg 2 CO 3 (OH) 2 ·3H 2 O. Our ab initio phase stability plots predict a small stability region for Artinite at low temperature ( Figure  6c), which is only found in some of the experimental results. 74,75 Using the QHA free energies increases the stability of this phase over a broader range of temperatures (Figure 6d). We note however that the complexity of the Artinite structure that arises from the random distribution of the and H 2 O ligands is only partially considered CO 2 -3 in this study, as we use a single model structure. A random distribution would give some additional configurational entropy that would further stabilise Artinite. This is consistent with the work of Schott et al., 74 which suggests that Mg 2 CO 3 (OH) 2 ·3H 2 O is stable at low temperatures and low p CO2 . We also note that there are multiple entropy values for Mg 2 CO 3 (OH) 2 ·3H 2 O in the literature. 31,71,74 This makes it difficult to compare unambiguously our predictions to experiments, but, as our ab initio phase stability plots do not require any experimental corrections for the solid phases, and are therefore not biased towards particular sets of measurements, we can consider them truly predictive. With this in mind our results confirm that Mg 2 CO 3 (OH) 2 ·3H 2 O is likely to remain a rare phase and one deserving of future investigation, but the favourable comparison between our predicted phase stability plots and literature data suggests that our predictions may be accurate enough to identify the synthesis conditions required to form specific Mg-rich phases. However, we highlight a key assumption of our methodology, namely that the vibrational properties of the different phases are purely harmonic. This assumption may lead to deviation from experiment under some conditions, particularly at elevated temperature. However, using the quasi-harmonic approximation to account for thermal expansion, as we have done in the present study, did not significantly change the results despite a ten-fold increase in the computational cost. 81 However, despite these potential shortcomings, a key advantage of this thermodynamic framework is that we do not use any experimental corrections for the solid phases, and, as a result, it can be used unambiguously to compare the DFT functionals and the different techniques to account for entropy. It could also be used with a global-exploration structure search to identify potential unknown phases and place them on the phase stability plot. Similarly, although this work examines the MgO-CO 2 -H 2 O system, the methodology is completely general and could be applied to a wide variety of other systems, subject to choosing an appropriate theoretical description of the compounds under study.

Conclusions
In this work, we have developed an ab initio thermodynamics framework to situate the mineral Artinite (Mg 2 CO 3 (OH) 2 ·3H 2 O) within the phase space of other Mg hydroxides and carbonates. We find that the widely-overlooked Artinite phase is metastable and can be stabilised by inhibiting the formation of Magnesite (MgCO 3 ), thus providing an explanation for the experimental controversy over the conditions under which this rare phase can be observed. 31,71,74 Our methodology has several advantages over alternatives, viz.: (1) it can identify overlooked materials that may be stable under specific conditions of temperature and partial pressures of water and carbon dioxide; (2) it can be used to identity competing and metastable phases, such as Artinite; and (3) it does not rely on the availability and accuracy of solid-phase experimental thermodynamic data. The latter in particular means that the present results can inform the future development of synthetic strategies to target and selectively stabilise specific phases of interest. Our theoretical results also show that the inclusion of van der Waals dispersion corrections to density-functional theory is necessary to accurately predict the thermodynamic properties of the MgO-CO 2 -H 2 O system. Such corrections are particularly important for phases that display large numbers of hydrogen bonds, such as Artinite. On the other hand, applying the quasi-harmonic approximation does not lead to a significant improvement in the accuracy of our predictions over the harmonic approximation. Finally, we note that a potential limitation of our study is that we do not fully (because of the usage of ab initio calculations) treat the disorder in Artinite, which may provide additional stabilisation. Future work on this system should therefore consider applying classical molecular dynamics or Monte Carlo simulations to better take disorder into account.

Conflicts of interest
There are no conflicts to declare.