{'fraction_non_alphanumeric': 0.03443870466033439, 'fraction_numerical': 0.024384824751430317, 'mean_word_length': 4.6598553913863565, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Van der Waals interactions in the ground state of Mg(BH4)2 from density functional theory

In order to resolve an outstanding discrepancy between experiment and theory regarding the ground-state structure of Mg(BH4)2, we examine the importance of long-range dispersive interactions on the compound's thermodynamic stability. Careful treatment of the correlation effects within a recently developed nonlocal van der Waals density functional (vdW-DF) leads to a good agreement with experiment, favoring the {\alpha}-Mg(BH4)2 phase (P6122) and a closely related Mn(BH4)2-prototype phase (P3112) over a large set of polymorphs at low temperatures. Our study demonstrates the need to go beyond (semi)local density functional approximations for a reliable description of crystalline high-valent metal borohydrides.

In order to resolve an outstanding discrepancy between experiment and theory regarding the ground-state structure of Mg(BH4)2, we examine the importance of long-range dispersive interactions on the compound's thermodynamic stability. Careful treatment of the correlation effects within a recently developed nonlocal van der Waals density functional (vdW-DF) leads to a good agreement with experiment, favoring the α-Mg(BH4)2 phase (P6122) and a closely related Mn(BH4)2-prototype phase (P3112) over a large set of polymorphs at low temperatures. Our study demonstrates the need to go beyond (semi)local density functional approximations for a reliable description of crystalline high-valent metal borohydrides. Light-weight metal borohydrides remain in the spotlight of hydrogen-storage research due to their high gravimetric hydrogen content and abundance of the constituent elements. 1,2 Li, Mg, and Ca tetraborohydrides (M (BH 4 ) n ) hold over 12 wt% of H, but their excessive thermodynamic stability presently renders these materials unsuitable for practical reversible hydrogen-storage solutions. 1 In an effort to adjust the formation enthalpies and the decomposition temperatures, several groups have suggested destabilization routes via reactions with other hydrides (see e.g. Refs. [3][4][5] or synthesis of mixed metal tetraborohydrides (e.g. Li-Cu) 6 . Detailed knowledge of materials' ground states becomes essential, as the possibility of new compound formation depends on small free energy differences for the phases involved.
Metal tetraborohydrides are held together primarily by a combination of covalent (B-H) and ionic ([M ] nδ+ -[BH 4 ] δ− ) interactions. Packing of the fairly rigid [BH 4 ] δ− units and [M ] nδ+ ions into crystalline structures depends strongly on the valency and size of the metal atom, with exceptionally complex configurations occurring for the medium-sized divalent magnesium. Based on the experimental data, the ground states below and above T = 453 K, namely α and β, have been assigned space groups P6 1 22 (330 atoms/u.c.) 7,8 and Fddd (704 atoms/u.c.) 9 , respectively. Strikingly, numerous density functional theory (DFT) studies 7,10-15 have converged on a completely unrelated F222 structure as the most stable low-T polymorph. 12,14,15 Recently, it has been pointed out 16 that the lack of a proper description of dispersion forces in the standard DFT approach may be a factor favoring low-density F222 structure over the experimental α phase.
In this study we take a key step towards resolving the existing discrepancy: we demonstrate that the weak dispersive interactions, important in other ionic (KCl and KBr) 17 and covalent-ionic (Mg(OH) 2 and Ca(OH) 2 ) 18 systems, indeed play a critical role in defining the Mg(BH 4 ) 2 ground state. We show that the commonly used Ceperley-Alder 19 or Perdew, Burke, and Ernzerhof (PBE) 20 functionals artificially stabilize configurations with unusually high and low density, respectively. Inclusion of the dispersive contributions via a nonlocal van der Waals density functional (vdW-DF) 21,22 or as a semiempirical PBE-D* correction 23-25 changes the relative stability of considered structures and favors the layered motif 26 occurring in related α-Mg(BH 4 ) 2 and Mn(BH 4 ) 2 prototypes.
The careful re-examination of Mg(BH 4 ) 2 has become possible due to a recent surge of attempts to include vdW interactions in the DFT. 27,28 The truly nonlocal correlation functional in the vdW-DF approach has shown good transferability for a range of systems reaching from simple dimers 29 and physisorbed molecules 30 to DNA 31 and drug design 32 . An efficient FFT formulation of vdW-DF 33 has allowed us to calculate the T = 0 K relative stability for Mg(BH 4 ) 2 structures of unprecedented size (up to 330 atom/u.c.). Calculation of the Gibbs energy vibrational contributions at finite temperatures is much more computationally demanding. In order to examine the relative stability as a function of temperature we have employed a semi-empirical method developed by Grimme 24,25 , that includes the long-range contributions via damped pairwise f dmp (R)C 6 R −6 terms at a negligible cost compared to standard DFT calculations. The method has gained popularity providing an improved description of molecular systems. However, Civalleri et al. 18,23 demonstrated the need to adjust the parameterization (from PBE-D to PBE-D*) for the application of the method to crystalline solids and we observe a better transferrability of the modified set for the challenging case of Mg(BH 4 ) 2 .
We describe the considered library of structure types and the simulation settings in Sec. II, present our systematic comparison of the performance of five DFTbased methods for Mg(BH 4 ) 2 at T = 0 K in Sec. III, demonstrate the effect of vibrational entropy on the poly-morphs' relative stability at finite temperatures in Sec. IV, and conclude in Sec. V.

II. SIMULATION SETUP
Library of structure types. In addition to a large pool of previously considered candidate structures, we include potentially relevant A(BC 4 ) 2 prototypes found in the Inorganic Crystal Structure Database (ICSD) 34 ; the full list of 36 structures is given in Table I. All energies are referenced to the experimental low-T α phase, which has been recently argued to have a higher symmetry (P6 1 22) 7,8 than originally thought (P6 1 ) 9,35 . Our simulations confirm that, while the two structures are nearly degenerate in energy, only P6 1 22 is dynamically stable (P6 1 has multiple imaginary phonon modes reaching 50i cm −1 ). The structure of the high-T β-phase, observed experimentally to be stable between 453 K and 613 K 9, 35 , remains an open question: the first powder diffraction orthorhombic solution with the Fddd symmetry, based on the positions of the B and Mg atoms, has six imaginary phonon modes at Γ reaching 65i cm −1 . Notable Mg(BH 4 ) 2 candidates identified in previous DFT studies include the trigonal P3m1 structures 10 as well as the unusual lowdensity F222 phase derived from I4m2 12,14 . Figures 1a-c illustrate the diversity of the morphologies as different arrangements of the Mg 2δ+ ions and BH δ− 4 units result in layered, hollow-framework-like, or fairly uniform denselypacked structures.
Simulation settings. The total energy calculations are carried out in the generalized gradient approximation (GGA) with the PBE 20 exchange-correlation (xc) functional and the local density functional approximation (LDA) 19 as implemented in VASP 36 . We employ projector augmented-wave pseudopotentials (PAW) 37 in which the semicore states are treated as valence. An energy cutoff of 500 eV and dense Monkhorst-Pack k-point meshes (∼ 0.03Å −1 in each direction in the Brillouin zone) 38 are applied. The vdW-DF, PBE, PBE-D, and PBE-D* calculations are performed with a modified version of PWscf 39 using ultrasoft pseudopotentials. All structures were fully relaxed with the threshold of 10 −5 [Ry] for energy convergence and 3·10 −1 [Ry a.u. −1 ] for residual force; the residual stress was typically under 1 kbar. Benchmark tests given in Appendix A suggest that for a given xc functional the errors arising from other factors (convergence criteria, choice of pseudopotential, ect.) are rather small allowing us to resolve polymorphs ∼ 0.005 eV/BH 4 apart. Phonon spectra and the Gibbs energy corrections for selected structures are calculated with PBE and PBE-D* with a finite displacement method as implemented in PHON 40 . Forces for phonon calculations were obtained from VASP (PBE) or QMPOT 41,42 linked to VASP (PBE-D*) (for more details see Sec. IV).
The necessity to treat Mg(BH 4 ) 2 beyond the standard DFT approximations becomes evident when the PBE functional shows virtually no interlayer binding for the exemplary layered P3m1 phase (Fig. 1d). The lack of a noticeable covalent or electrostatic interaction between the layers is a result of the particular packing of the BH 4 tetrahedra that makes the interlayer interface consist of two parallel sheets of hydrogen atoms already engaged in covalent B-H bonds. Introduction of The size of the extra binding is substantial considering that the reported energy difference between the experimental (α) and theoretical (F222) ground states within PBE is only 0.024 eV/BH 4 15 . The LDA, known to mimic dispersive interactions (e.g., in graphite 43 ), also gives a 0.058 eV/BH 4 binding, but -as we show below -should not be used as a substitute for a properly constructed nonlocal functional.
We compare the performance of the five different methods, PBE, LDA, PBE-D, PBE-D*, and vdW-DF, by plotting the relative stability and compound density for selected structures in Fig. 2 (the corresponding values for all 36 structures are listed in Table I). A number of unexpected features emerge as we examine the total energy data moving from left to right in Fig. 2b. At the PBE level, we find the highest stability of the least compact F222 phase, in agreement with Ref. 13 . However, in the LDA, arguably better suited for simulating Mg(BH 4 ) 2 based on our test in Fig. 1, there is a dramatic change in the ordering of the polymorphs' enthalpies with the ground state now being the most compact I4 1 /acd structure. It is evident that, compared to PBE, the LDA favors higher packing: it lowers the relative enthalpy of the two densest P3m1 and I4 1 /acd polymorphs by 0.134 and 0.120 eV/BH 4 , respectively, and increases the relative enthalpy of the sparsest F222 by 0.083 eV/BH 4 .  Table I. leads to a significant reordering of the relative stabilities obtained with PBE and promotes the I4 1 /acd structure. In fact, there is more agreement between the LDA vs PBE-D rather than the PBE vs PBE-D sets. The P3m1 polymorph gets an even larger 0.217 eV/BH 4 gain in relative stability dropping below α by 0.023 eV/BH 4 , which may look surprising given the similar size of the interlayer binding for this structure in LDA and PBE-D (Fig. 1).

Inclusion of the dispersive contribution in PBE-D also
To help explain this result we plot pairwise interaction strengths (Fig. 3g) and histograms of nearest neighbors (NNs) for relevant structures (Fig. 3a-f). It has been discussed 23,27 that introduction of the dispersive interactions 'on top' of DFT requires a smooth cut-off of the attractive C ij /R −6 ij terms just below the typical vdW distances to avoid their unphysical contributions at short R ij . Panels f and g in Fig. 3 show that a number of intralayer distances in P3m1 (marked with stars) happen to fall in the critical range, around the corresponding R vdW i + R vdW j values fitted to molecular datasets, which can give rise to artificially strong binding.
The suggested adjustment of the vdW parameters within PBE-D* 23 helps reduce the vdW overbinding for the dense P3m1 and I4 1 /acd structures and for the first time we observe the α-phase to be the most stable polymorph among the considered diverse set of candidates, albeit by a small margin. The classical pairwise representation of the vdW corrections in the Grimme approach (Eq. 1) makes it possible to examine the individual vdW contributions. According to Fig. 3h Fig. 4 shows that the BH 4 -BH 4 interaction still dominates the long-range binding, which -not unexpectedly -scales almost linearly with the crystal density. The second largest contribution, Mg-BH 4 , deviates from linear dependence due to a variable coordination of Mg ions. The classical attractive terms cause up to ∼ 20% volume reduction compared to structures optimized with PBE.
Finally, the vdW-DF method correctly predicts α to be the low-T ground state and separates it from all other unrelated structures by a considerable 0.040 eV/BH 4 ; the near degeneracy of α and P3 1 12 (Mn(BH 4 ) 2 prototype) can be traced to the close structural relationship between the two layered structures discussed in detail in Ref. 26 . The method places the hypothetical F222 structure 0.076 eV/BH 4 above α and differentiates between the two highdensity 2D P3m1 and 3D I4 1 /acd structures by improving the stability of just the latter. This suggests that the overestimation of the intralayer binding in the difficult P3m1 case might still not be fully corrected within PBE-D*. Comparison of crystal densities for PBE-D* and vdW-DF in Fig. 2d shows that the latter results in less compact structures; the effect of the nonlocal correlations on the volume change is not easy to quantify because the method is based on the exchange and correlation functionals taken from revPBE 44 and LDA, respectively. The calculated vdW-DF value of ρ = 0.758 gcm −3 for the α phase is 3% below the experimental value which is consistent with a typical 1% overestimation of bond lengths for GGA-based methods.

IV. INCLUSION OF THE VIBRATIONAL CONTRIBUTIONS TO G(T )
Inclusion of finite temperature contributions due to the vibrational entropy is orders of magnitude more expensive than the calculation of the total energy at zero temperature. We estimate the vibrational corrections using two DFT-based methods, PBE-D and PBE-D*, which is computationally feasible and allows us to show directly the effect of the vdW interactions. Hessian matrices were calculated in a numerical fashion from analitycal forces and atomic displacements for sufficiently large supercells (at least 88 atoms). 45 According to our extensive tests for Mg(BH 4 ) 2 , the 0.00133Å displacements and the 10 −8 eV SCF energy convergence criterion ensured minimal errors (∼ 10 cm −1 ) for the calculated frequencies coming from the unharmonic effects and numerical factors. qmesh for the dynamical matrices was typically two times denser than the k-mesh for the energy calculations in each direction. The analitycal forces at the PBE level were found with VASP for fully relaxed polymorph structures. The analytical forces at the PBE-D* level were calculated using QMPOT linked with VASP. In this case, since QMPOT does not currently offer stress calculations, the structures were optimized as follows: the unit cells were first fully relaxed with PWscf and then the ionic positions were re-optimized with a combination of QMPOT and VASP. As described in Appendix A, our careful tests for selected structures showed insignificant variation in the equilibrium cell parameters in the two settings. Once the phonon densities of states n(ω) were found, the vibrational contributions G phon (T ) were included via 40  Calculated relative Gibbs energies G(T ) for selected low-energy structures are shown in Fig. 2a and 2c. We find that the zero point energy alone can change the relative stability by over 0.04 eV/BH 4 . Addition of dispersive interactions leads to about 0.01 eV/BH 4 changes in the relative stability at T = 800 K and can affect the ordering of the polymorphs as happens for the P3 1 12 structure. The related α and P3 1 12 phases remain nearly degenerate in the whole T range and cannot be unambiguously resolved within our vdW-DF or PBE-D* simulations. The search for the high-T ground state may be simplified by our observation of little variation in the en-thalpy difference between α and β for all the five methods in Fig. 2: due to the close structural relationship between the two phases, shown in Fig. 3a and 3c, there must be a cancellation of errors. Hence, despite the remaining noticeable difference between the PBE-D* and vdW-DF relative stability for some polymorphs the former could be an appropriate choice for identification of viable high-T candidates and estimation of the phonon contributions to G(T ).

V. CONCLUSIONS
The presented results illustrate that identification of ground states depends not only on the exhaustive sampling of possible structures and compositions but also on the accuracy of the chosen simulation method. We have demonstrated that inclusion of the dispersive interactions via nonlocal correlation functional within vdW-DF leads to a good agreement with experiment for the Mg(BH 4 ) 2 low-T ground state. The related wellcharacterized Ca(BH 4 ) 2 system should be an interesting test case for the vdW-DF and PBE-D* methods since standard DFT approximations seem to already reproduce the experiment 46 . Proper treatment of the nonlocal correlation effects is expected to be most critical in highervalent metal borohydrides that can form weakly interacting molecular-type complexes 16 . Until such large crystalline systems can be handled with accurate quantum Monte Carlo-level or quantum chemistry methods, the continued improvement of DFT-based methods 47 may provide an attractive alternative for description of the dispersive interactions therein. In any case, it appears to be good practice to systematically test proposed ground states with a range of DFT flavours as this can help spot and avoid potential artifacts of DFT approximations.

VI. ACKNOWLEDGMENTS
A.B. and A.N.K. acknowledge the support of the EPSRC through the DTI Technology Program EP/D062098/1 and the CAF EP/G004072/1. Calculations were performed in the Gdansk Supercomputing Center, on the DEAC cluster of Wake Forest University, and in the Oxford Supercomputing Centre. A.B. thanks Bernhard Seiser for helpful discussions.