Hierarchy Configuration Interaction: Combining Seniority Number and Excitation Degree

F ábris Kossoski, 1, a) Yann Damour, [START_REF] Szabo | Modern quantum chemistry[END_REF] and Pierre-Franc ¸ois Loos 1, b) Laboratoire de Chimie et Physique Quantiques (UMR 5626), Université de Toulouse, CNRS, UPS, France

We propose a novel partitioning of the Hilbert space, hierarchy configuration interaction (hCI), where the excitation degree (with respect to a given reference determinant) and the seniority number (i.e., the number of unpaired electrons) are combined in a single hierarchy parameter. The key appealing feature of hCI is that each hierarchy level accounts for all classes of determinants whose number share the same scaling with system size. By surveying the dissociation of multiple molecular systems, we found that the overall performance of hCI usually exceeds or, at least, parallels that of excitation-based CI. For higher orders of hCI and excitation-based CI, the additional computational burden related to orbital optimization usually do not compensate the marginal improvements compared with results obtained with Hartree-Fock orbitals. The exception is orbital-optimized CI with single excitations, a minimally correlated model displaying the qualitatively correct description of single bond breaking, at a very modest computational cost.

Hierarchy configuration interaction (hCI)

Excitation degree e Seniority number s Hierarchy parameter h = e+s/2 In electronic structure theory, configuration interaction (CI) methods allow for a systematic way to obtain approximate and exact solutions of the electronic Hamiltonian, by expanding the wave function as a linear combination of Slater determinants (or configuration state functions). [START_REF] Szabo | Modern quantum chemistry[END_REF][START_REF] Helgaker | Molecular Electronic-Structure Theory[END_REF] At the full CI (FCI) level, the complete Hilbert space is spanned in the wave function expansion, leading to the exact solution for a given one-electron basis set. Except for very small systems, [START_REF] Knowles | [END_REF]4 the FCI limit is unattainable, and in practice the expansion of the CI wave function must be truncated. The question is then how to construct an effective and computationally tractable hierarchy of truncated CI methods that quickly recover the correlation energy, understood as the energy difference between the FCI and the mean-field Hartree-Fock (HF) solutions.

Excitation-based CI is surely the most well-known and popular class of CI methods. In this context, one accounts for all determinants generated by exciting up to e electrons from a given reference, which is usually the HF determinant, but does not have to. In this way, the excitation degree e defines the following sequence of models: CI with single excitations (CIS), CI with single and double excitations (CISD), CI with single, double, and triple excitations (CISDT), and so on. Excitationa) Electronic mail: fkossoski@irsamc.ups-tlse.fr b) Electronic mail: loos@irsamc.ups-tlse.fr based CI manages to quickly recover weak (dynamic) correlation effects, but struggles in strong (static) correlation regimes. It also famously lacks size-consistency which explains issues, for example, when dissociating chemical bonds. Importantly, the number of determinants N det (which is the key parameter governing the computational cost, as discussed later) scales polynomially with the number of basis functions N as N 2e .

Alternatively, seniority-based CI methods (sCI) have been proposed in both nuclear [START_REF] Ring | The Nuclear Many-Body Problem, Theoretical and Mathematical Physics, The Nuclear Many-Body Problem[END_REF] and electronic 6 structure calculations. In short, the seniority number s is the number of unpaired electrons in a given determinant. By truncating at the seniority zero (s = 0) sector (sCI0), one obtains the well-known doublyoccupied CI (DOCI) method, [START_REF] Bytautas | [END_REF][7][8][9] which has been shown to be particularly effective at catching static correlation, while higher sectors tend to contribute progressively less. [START_REF] Bytautas | [END_REF][10][11][12] In addition, sCI0 is size-consistent, a property that is not shared by higher orders of seniority-based CI. However, already at the sCI0 level, N det scales exponentially with N, since excitations of all degrees are included. Therefore, despite the encouraging successes of seniority-based CI methods, their unfavorable computational scaling restricts applications to very small systems. 13 Besides CI, other methods that exploit the concept of seniority number have been pursued. [14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][START_REF] Johnson | Bivariational principle for an antisymmetrized product of nonorthogonal geminals appropriate for strong electron correlation[END_REF][START_REF] Fecteau | Near-exact treatment of seniority-zero ground and excited states with a richardson-gaudin mean-field[END_REF] At this point, we notice the current dichotomy. When targeting static correlation, seniority-based CI methods tend to have a better performance than excitation-based CI, despite their higher computational cost. The latter class of methods, in contrast, are well-suited for recovering dynamic correlation, and only at polynomial cost with system size. Ideally, we aim for a method that captures most of both static and dynamic correlation, with as few determinants as possible. With this goal in mind, we propose a new partitioning of the Hilbert space, named hierarchy CI (hCI). It combines both the excitation degree e and the seniority number s into one single hierarchy parameter

h = e + s/2 2 , (1)
which assumes half-integer values. Here we only consider systems with an even number of electrons, meaning that s takes only even values as well. Figure 1 shows how the Hilbert space is progressively populated in excitation-based CI, senioritybased CI, and our hybrid hCI methods.

We have three key justifications for this new CI hierarchy. The first one is physical. We know that the lower degrees of excitations and lower seniority sectors, when looked at individually, often carry the most important contribution to the FCI expansion. By combining e and s as is Eq. ( 1), we ensure that both directions in the excitation-seniority map (see Fig. 1) are contemplated. Rather than filling the map top-bottom (as in excitation-based CI) or left-right (as in seniority-based CI), the hCI methods fills it diagonally. In this sense, we hope to recover dynamic correlation by moving right in the map (increasing the excitation degree while keeping a low seniority number), at the same time as static correlation, by moving down (increasing the seniority number while keeping a low excitation degree).

The second justification is computational. In the hCI class of methods, each level of theory accommodates additional determinants from different excitation-seniority sectors (each block of same color tone in Fig. 1). The key insight behind hCI is that the number of additional determinants presents the same scaling with respect to N, for all excitation-seniority sectors entering at a given hierarchy h. This justifies the numerator in the definition of h [Eq. (1)].

Finally, the third justification for our hCI method is empirical and closely related to the computational motivation. There are many possible ways to populate the Hilbert space starting from a given reference determinant, and one can in principle formulate any systematic recipe that includes progressively more determinants. Besides a physical or computational perspective, the question of what makes for a good recipe can be framed empirically. Does our hCI class of methods perform better than excitation-based or seniority-based CI, in the sense of recovering most of the correlation energy with the least computational effort?

Hybrid approaches based on both excitation degree and seniority number have been proposed before. 12,[START_REF] Van Raemdonck | [END_REF]32 In these works, the authors established separate maximum values for the excitation and the seniority, and either the union or the intersection between the two sets of determinants have been considered. For the union case, N det grows exponentially with N, while in the intersection approach the Hilbert space is filled rectangle-wise in our excitation-seniority map. In the latter case, the scaling of N det would be dominated by the rightmost bottom block. Bytautas et al. 10 explored a different hybrid scheme combining determinants having a maximum seniority number and those from a complete active space. In comparison to previous approaches, our hybrid hCI scheme has two key advantages. First, it is defined by a single parameter that unifies excitation degree and seniority number [see Eq. ( 1)]. Second and most importantly, each next level includes all classes of determinants whose number share the same scaling with system size, as discussed before, thus preserving the polynomial cost of the method.

Each level of excitation-based CI has a hCI counterpart with the same scaling of N det with respect to N, justifying the denominator in the definition of h [Eq. (1)]. For example, N det = O(N 4 ) in both hCI2 and CISD, whereas N det = O(N 6 ) in hCI3 and CISDT, and so on. From this computational perspective, hCI can be seen as a more natural choice than the traditional excitation-based CI, because if one can afford for, say, CISDT, then one could probably afford hCI3, due to the same scaling of N det . Of course, in practice an integer-h hCI method has more determinants than its excitation-based counterpart (despite the same scaling of N det ), and thus one should first ensure whether including the lower-triangular blocks (going from CISDT to hCI3 in our example) is a better strategy than adding the next column (going from CISDT to CISDTQ). Therefore, here we decided to discuss the results in terms of N det , rather than the formal scaling of N det as a function of N, which could make the comparison somewhat biased toward hCI. It is also interesting to compare the lowest levels of hCI (hCI1) and excitation-based CI (CIS). Since single excitations do not connect with the reference (at least for HF orbitals), CIS provides the same energy as HF. In contrast, the paired double excitations of hCI1 do connect with the reference (and the singles contribute indirectly via the doubles). Therefore, while CIS based on HF orbitals does not improve with respect to the mean-field HF wave function, the hCI1 counterpart already represents a minimally correlated model, with the same and favorable N det = O(N 2 ) scaling. hCI also allows for half-integer values of h, with no equivalent in excitation-based CI. This gives extra flexibility in terms of methodological choice. For a particular application with excitation-based CI, CISD might be too inaccurate, for example, while the improved accuracy of CISDT might be too expensive. hCI2.5 could represent an alternative, being more accurate than hCI2 and less expensive than hCI3.

Our main goal here is to assess the performance of hCI against excitation-based and seniority-based CI. To do so, we have evaluated how fast different observables converge to the FCI limit as a function of N det . In particular, we have calculated the potential energy curves (PECs) for the dissociation of six systems: HF, F 2 , N 2 , ethylene, H 4 , and H 8 , which display a variable number of bond breaking. For the latter two molecules, we have considered linearly arranged with equally spaced hydrogen atoms, and computed PECs along the symmetric dissociation coordinate. For ethylene, we consider the C --C double bond breaking, while freezing the remaining internal coordinates. Its equilibrium geometry was taken from Ref. 33 for electronic structure methods, being often considered when assessing novel methodologies. More precisely, we have evaluated the convergence of four observables: the non-parallelity error (NPE), the distance error, the vibrational frequencies, and the equilibrium geometries. The NPE is defined as the maximum minus the minimum differences between the PECs obtained at a given CI level and the exact FCI result. We define the distance error as the maximum plus the minimum differences between a given PEC and the FCI result. Thus, while the NPE probes the similarity regarding the shape of the PECs, the distance error measures how their overall magnitudes compare. From the PECs, we have also extracted the vibrational frequencies and equilibrium geometries (details can be found in the Supporting Information).

The hCI method was implemented in quantum package via a straightforward adaptation of the configuration interaction using a perturbative selection made iteratively (CIPSI) algorithm, [34][35][36][37] by allowing only for determinants having a given maximum hierarchy h to be selected. The excitationbased CI, seniority-based CI, and FCI calculations presented here were also performed with the CIPSI algorithm implemented in quantum package. 38 It is worth mentioning that the determinant-driven framework of quantum package allows the inclusion of any arbitrary set of determinants. In practice, we consider, for a given CI level, the ground state energy to be converged when the second-order perturbation correction computed in the truncated Hilbert space (which approximately measures the error between the selective and complete calculations) lies below 0.01 mE h . 37 These selected versions of CI require considerably fewer determinants than the formal number of determinants (understood as all those that belong to a given CI level, regardless of their weight or symmetry) of their complete counterparts. Nevertheless, we decided to present the results as functions of the formal number of determinants (see above), which are not related to the particular algorithmic choices of the CIPSI calculations. The ground-state CI energy is obtained with the Davidson iterative algorithm, 39 which in the present implementation of quantum package means that the computation and storage cost us O(N 3/2 det ) and O(N det ), re-spectively. This shows that the determinant-driven algorithm is not optimal in general. However, the selected nature of the CIPSI algorithm drastically reduces the actual number of determinants and therefore calculations are technically feasible.

The CI calculations were performed with both canonical HF orbitals and optimized orbitals. In the latter case, the energy is obtained variationally in the CI space and in the orbital parameter space, hence defining orbital-optimized CI (oo-CI) methods. We employed the algorithm described elsewhere 40 and also implemented in quantum package for optimizing the orbitals within a CI wave function. In order to avoid converging to a saddle point solution, we employed a similar strategy as recently described in Ref. 41. Namely, whenever the eigenvalue of the orbital rotation Hessian is negative and the corresponding gradient component g i lies below a given threshold g 0 , then this gradient component is replaced by g 0 |g i |/g i . Here we took g 0 = 1 µE h , and considered the orbitals to be converged when the maximum orbital rotation gradient lies below 0.1 mE h . While we cannot ensure that the obtained solutions are global minima in the orbital parameter space, we verified that all stationary solutions surveyed here correspond to real minima (rather than maxima or saddle points). All CI calculations were performed with the cc-pVDZ basis set and within the frozen core approximation. For the HF molecule we have also tested basis set effects, by considering the larger cc-pVTZ and cc-pVQZ basis sets.

It is worth mentioning that obtaining smooth PECs for the orbital optimized calculations proved to be far from trivial. First, the orbital optimization was started from the HF orbitals of each geometry. This usually led to discontinuous PECs, meaning that distinct solutions were found by our algorithm. Then, at some geometries that seem to present the lowest lying solution, the optimized orbitals were employed as the guess orbitals for the neighboring geometries, and so on, until a new PEC is obtained. This protocol was repeated until the PEC built from the lowest lying oo-CI solution becomes continuous. We recall that saddle point solutions were purposely avoided in our orbital optimization algorithm. If that was not the case, then even more stationary solutions would have been found.

While the full set of PECs and the corresponding energy differences with respect to FCI are shown in the Supporting Information, in Fig. 2 we present the PECs for F 2 , which display many of the features also observed for the other systems. It already gives a sense of the performance of three classes of CI methods, clearly showing the overall superiority of hCI over excitation-based CI. It further illustrates several important features which will be referenced to in the upcoming discussion. We first discuss the results for HF orbitals. In Fig. 3, we present the NPEs for the six systems studied, and for the three classes of CI methods, as functions of N det . The main result contained in Fig. 3 concerns the overall faster convergence of hCI when compared to excitation-based and seniority-based CI. This is observed for single bond breaking (HF and F 2 ) as well as the more challenging double (ethylene), triple (N 2 ), and quadruple (H 4 ) bond breaking. For H 8 , hCI and excitationbased CI perform similarly. The convergence with respect to N det is slower in the latter, more challenging cases, irrespective of the class of CI methods, as expected. 42,[START_REF] Motta | on the Many-Electron Problem)[END_REF] But more importantly, the superiority of hCI appears to be highlighted in the one-site multiple bond break systems (compare ethylene and N 2 with HF and F 2 in Fig. 3).

For all systems (specially ethylene and N 2 ), hCI2 is better than CISD, two methods where N det scales as N 4 . hCI2.5 is better than CISDT (except for H 8 ), despite its lower computational cost, whereas hCI3 is much better than CISDT, and comparable in accuracy with CISDTQ (again for all systems). Inspection of the PECs (see Fig. 2 for the case of F 2 or the Supporting Information for the other systems) reveals that the lower NPEs observed for hCI stem mostly from the contribution of the dissociation region. This result demonstrates the importance of higher-order excitations with low seniority number in this strong correlation regime, which are accounted for in hCI but not in excitation-based CI (for a given scaling of N det ). These determinants are responsible for alleviating the size-consistency problem when going from excitation-based CI to hCI.

Meanwhile, the first level of seniority-based CI (sCI0, which is the same as DOCI) tends to offer a rather low NPE when compared to the other CI methods with a similar N det (hCI2.5 and CISDT). However, convergence is clearly slower for the next levels (sCI2 and sCI4), whereas excitation-based CI and specially hCI converge faster. Furthermore, seniority-based CI becomes less attractive for larger basis set in view of its exponential scaling. This can be seen in Figs. S2 andS3 of the Supporting Information, which shows that augmenting the basis set leads to a much steeper increase of N det for senioritybased CI.

It is worth mentioning the surprisingly good performance of hCI1 and hCI1.5. For HF, F 2 , and ethylene, they yield lower NPEs than the much more expensive CISDT method, and only slightly higher in the case of N 2 . For the same systems, we also see the NPEs increase from hCI1.5 to hCI2, and decreasing to lower values only at the hCI3 level. (Even than, it is important to remember that the hCI2 results remain overall superior to their excitation-based counterparts.) Both findings are not observed for H 4 and H 8 . It seems that both the relative worsening of hCI2 and the success of hCI1 and hCI1.5 become less apparent as progressively more bonds are being broken (compare, for instance, F 2 , N 2 , and H 8 in Fig. 3). This reflects the fact that higher-order excitations are needed to properly describe multiple bond breaking, and also hints at some cancelation of errors in low-order hCI methods for single bond breaking.

In Fig. S5 of the Supporting Information, we present the distance error, which is also found to decrease faster with hCI. Most of the observations discussed for the NPE also hold for the distance error, with two main differences. The convergence is always monotonic for the latter observable (which is expected from its definition), and the performance of seniority-based CI is much poorer (due to the slow recovery of dynamic correlation).

In Figs. S6 and S7 of the Supporting Information, we present the convergence of the equilibrium geometries and vibrational frequencies, respectively, as functions of N det , for the three classes of CI methods. For the equilibrium geometries, hCI performs slightly better overall than excitation-based CI. A more significant advantage of hCI can be seen for the vibrational frequencies. For both observables, hCI and excitation-based CI largely outperform seniority-based CI. Similarly to what we have observed for the NPEs, the convergence of hCI is also found to be non-monotonic in some cases. This oscillatory behavior is particularly evident for F 2 (notice in Fig. 2 how the potential well goes from shallow in hCI1, to deep in hCI2, and then shallow again in hCI3). It is also noticeable for HF, becoming less apparent for ethylene, virtually absent for N 2 , and showing up again for H 4 and H 8 . Interestingly, equilibrium geometries and vibrational frequencies of HF and F 2 (single bond breaking), are rather accurate when evaluated at the hCI1.5 level, bearing in mind its relatively modest computational cost.

For the HF molecule we have also evaluated how the convergence is affected by increasing the size of the basis set, going from cc-pVDZ to cc-pVTZ and cc-pVQZ (see Figs. S2 and S3 in the Supporting Information). While a larger N det is required to achieve the same level of convergence, as expected, the convergence profiles remain very similar for all basis sets. Vibrational frequency and equilibrium geometry present less oscillations for hCI. We thus believe that the main findings discussed here for the other systems would be equally basis set independent.

Up to this point, all results and discussions have been based on CI calculations with HF orbitals. We recall that senioritybased CI (in contrast to excitation-based CI) is not invariant with respect to orbital rotations within the occupied and virtual subspaces, [START_REF] Bytautas | [END_REF] and for this reason it is customary to optimize the corresponding wave function by performing such rotations. Similarly, hCI wave functions are not invariant under orbital rotations within each subspace. Thus, we decided to further assess the role of orbital optimization (occupied-virtual rotations included) for each class of CI methods. Due to the significantly higher computational cost and numerical difficulties associated with orbital optimization at higher CI levels, such calculations were typically limited up to oo-CISD (for excitation-based), oo-DOCI (for seniority-based), and oo-hCI2 (for hCI). The PECs and convergence of properties as function of N det are shown in the Supporting Information.

Of course, at a given CI level, orbital optimization will lead to lower energies than with HF orbitals. However, even though the energy is lowered (thus improved) at each geometry, such improvement may vary largely along the PEC, which may or may not decrease the NPE. More often than not, the NPEs do decrease upon orbital optimization, though not always. For example, compared with their non-optimized counterparts, oo-hCI1 and oo-hCI1.5 provide somewhat larger NPEs for HF and F 2 , similar NPEs for ethylene, and smaller NPEs for N 2 , H 4 , and H 8 . Following the same trend, oo-CISD presents smaller NPEs than HF-CISD for the multiple bond breaking systems, but very similar ones for the single bond breaking cases. oo-CIS has significantly smaller NPEs than HF-CIS, being comparable to oo-hCI1 for all systems except for H 4 and H 8 , where the latter method performs better. (We will come back to oo-CIS later.) Based on the present oo-CI results, hCI still has the upper hand when compared with excitation-based CI, though by a smaller margin.

Orbital optimization usually reduces the NPE for senioritybased CI (in this case we only considered oo-DOCI) as well. The gain is specially noticeable for H 4 and H 8 (where the orbitals become symmetry-broken 24 ), and much less so for HF, ethylene, and N 2 (where the orbitals remain symmetrypreserved). This is in line with what has been observed before for N 2 . [START_REF] Bytautas | [END_REF] For F 2 , we found that orbital optimization actually increases the NPE (though by a small amount), due to the larger energy lowering in the Franck-Condon region than at dissociation (see Fig. 2). These results suggest that, when bond breaking involves one site, orbital optimization at the DOCI level does not have such an important role, at least in the sense of decreasing the NPE.

Optimizing the orbitals at the CI level also tends to benefit the convergence of vibrational frequencies and equilibrium geometries. The impact is often somewhat larger for hCI than for excitation-based CI, by a small margin. Also, the large oscillations observed in the hCI convergence with HF orbitals (for HF and F 2 ) are significantly suppressed upon orbital optimization.

We come back to the surprisingly good performance of oo-CIS, which is interesting due to its low computational cost. The PECs are compared with those of HF and FCI in Fig. S12 of the Supporting Information. At this level, the orbital rotations provide an optimized reference (different from the HF determinant), from which only single excitations are performed. Since the reference is not the HF determinant, Brillouin's theorem no longer holds, and single excitations actually connect with the reference. Thus, with only single excitations (and a reference that is optimized in the presence of these excitations), one obtains a minimally correlated model. Interestingly, oo-CIS recovers a non-negligible fraction (15%-40%) of the correlation energy around the equilibrium geometries. For all systems, significantly more correlation energy (25%-65% of the total) is recovered at dissociation. In fact, the larger account of cor-relation at dissociation is responsible for the relatively small NPEs encountered at the oo-CIS level. We also found that the NPE drops more significantly (with respect to the HF one) for the single bond breaking cases (HF and F 2 ), followed by the double (ethylene) and triple (N 2 ) bond breaking, then H 4 , and finally H 8 .

The above findings can be understood by looking at the character of the oo-CIS orbitals. At dissociation, the closed-shell reference is actually ionic, with orbitals assuming localized atomic-like characters. The reference has a decreasing weight in the CI expansion as the bond is stretched, becoming virtually zero at dissociation. However, that is the reference one needs to achieve the correct open-shell character of the fragments when the single excitations of oo-CIS are accounted for. Indeed, the most important single excitations promote the electron from the negative to the positive fragment, resulting in two singly open-shell radicals. This is enough to obtain the qualitatively correct description of single bond breaking, hence the relatively low NPEs observed for HF and F 2 . In contrast, the oo-CIS method can only explicitly account for one unpaired electron on each fragment, such that multiple bond breaking become insufficiently described. Nevertheless, double (ethylene) and even triple (N 2 ) bond breaking still appear to be reasonably well-described at the oo-CIS level.

In this Letter, we have proposed a new scheme for truncating the Hilbert space in configuration interaction calculations, named hierarchy CI (hCI). By merging the excitation degree and the seniority number into a single hierarchy parameter h, the hCI method ensures that all classes of determinants sharing the same scaling of N det with the number of basis functions are included in each level of the hierarchy. We evaluated the performance of hCI against excitation-based CI and senioritybased CI, by comparing PECs and derived quantities for six systems, ranging from single to multiple bond breaking.

Our key finding is that the overall performance of hCI either surpasses or equals that of excitation-based CI, in the sense of convergence with respect to N det . The superiority of hCI is more noticeable for the non-parallelity and distance errors, but also observed to a lesser extent for the vibrational frequencies and equilibrium geometries. The comparison to senioritybased CI is less trivial. DOCI (the first level of seniority-based CI) often provides even lower NPEs for a similar N det , but it falls short in describing the other properties investigated here. In addition, if higher accuracy is desired, convergence was found to be faster with hCI (and also excitation-based CI) than seniority-based CI, at least for HF orbitals. Finally, the exponential scaling of seniority-based CI in practice precludes this approach for larger systems and basis sets, while the favorable polynomial scaling and encouraging performance of hCI is an alternative.

We found surprisingly good results for the first level of hCI (hCI1) and the orbital optimized version of CIS (oo-CIS), two methods with very favorable computational scaling. In particular, oo-CIS correctly describes single bond breaking. We hope to report on generalizations to excited states in the future. In contrast, orbital optimization at higher CI levels is not necessarily a recommended strategy, given the overall modest improvement in convergence when compared to re-sults with canonical HF orbitals. One should bear in mind that optimizing the orbitals is always accompanied with wellknown challenges (several solutions, convergence issues, etc) and may imply a significant computational burden (associated with the calculations of the orbital gradient and Hessian, and the many iterations that are often required), specially for larger CI spaces. In this sense, stepping up in the CI hierarchy might be a more straightforward and possibly a cheaper alternative than optimizing the orbitals. One possibility to explore is to first optimize the orbitals at a lower level of CI, and then to employ this set of orbitals at a higher level of CI.

The hCI pathway presented here offers several interesting possibilities to pursue. One could generalize and adapt hCI for excited states 44 and open-shell systems, 45 develop coupledcluster methods based on an analogous excitation-seniority truncation of the excitation operator, [46][47][48] and explore the accuracy of hCI trial wave functions for quantum Monte Carlo simulations. [49][50][51]

FIG. 1 .

 1 FIG.1. Partitioning of the Hilbert space into blocks of specific excitation degree e (with respect to a closed-shell determinant) and seniority number s. This e-s map is truncated differently in excitation-based CI (left), seniority-based CI (right), and hierarchy-based CI (center). The color tones represent the determinants that are included at a given CI level.

FIG. 2 .

 2 FIG.2. Potential energy curves for F 2 , according to HF, FCI, and the three classes of CI methods: seniority-based CI (blue), excitationbased CI (red), and hierarchy-based CI (green), (dashed lines for half-integer h), with HF orbitals (left) and orbitals optimized at a given CI level (right), and with the cc-pVDZ basis set.

8 FIG. 3 .

 83 FIG.3. Non-parallelity errors as function of the number of determinants, for the three classes of CI methods: seniority-based CI (blue), excitation-based CI (red), and our proposed hybrid hCI (green).

ACKNOWLEDGMENTS

This work was performed using HPC resources from CALMIP (Toulouse) under allocation 2021-18005. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No. 863481).

SUPPORTING INFORMATION AVAILABLE

Potential energy curves, energy differences with respect to FCI results, non-parallelity errors, distance errors, vibrational frequencies, and equilibrium geometries, according to the three classes of CI methods (excitation-based CI, seniority-based CI, and hierarchy-based CI), with Hartree-Fock orbitals and with optimized orbitals, computed for HF, F 2 , ethylene, N 2 , H 4 , and H 8 , with the cc-pVDZ basis set, and also for HF, with cc-pVDZ, cc-pVTZ, cc-pVQZ basis sets. Equilibrium geometry of ethylene, and details about the fitting of potential energy curves.


Michele Ruggeri

Kyle Reeves

Tzu-Yao Hsu

Guillaume Jeanmairet

Mathieu Salanne

Carlo Pierleoni

Multi-scale simulation of the adsorption of lithium ion on graphite surface: from Quantum Monte Carlo to Molecular Density Functional Theory

The structure of the double-layer formed at the surface of carbon electrodes is governed by the interactions between the electrode and the electrolyte species. However, carbon is notoriously difficult to simulate accurately, even with well-established methods such as electronic Density Functional Theory and Molecular Dynamics.

Here we focus on the important case of a lithium ion in contact with the surface of graphite, and we perform a series of reference Quantum Monte Carlo calculations that allow us to benchmark various electronic Density Functional Theory functionals. We then fit an accurate carbon-lithium pair potential, which is used in molecular Density Functional Theory calculations to determine the free energy of the adsorption of the ion on the surface in the presence of water. The adsorption profile in solution differs markedly from the gas phase results, which emphasize the role of the solvent on the properties of the double-layer.

I. INTRODUCTION

Carbon materials play a very important role in many chemistry fields, ranging from catalysis to electrochemistry. In energy storage applications, carbon is used as an electrode in the form of graphite for Li-ion batteries, [START_REF] Armand | Building better batteries[END_REF] of hard carbon for Na-ion batteries, 2 and of nanoporous materials for supercapacitors. [START_REF] Simon | Capacitive energy storage in nanostructured carbon-electrolyte systems[END_REF] In all these examples, the carbon materials are put in contact with a charged electrolyte, and the interfacial structure and dynamics play a crucial role in the operation of the devices. [START_REF] Salanne | Efficient Storage Mechanisms for Building Better Supercapacitors[END_REF][START_REF] Striolo | The carbon-water interface: Modeling challenges and opportunities for the water-energy nexus[END_REF] In order to optimize the performance, molecular simulations can play an important role, [START_REF] Lahrar | Carbon-carbon supercapacitors: Beyond the average pore size or how electrolyte confinement and inaccessible pores affect the capacitance[END_REF] but they need to accurately account for the interactions between the ions and the carbon surface.

Over the years, a large number of electronic density functional theory (eDFT) and molecular dynamics (MD) simulations were devoted to the study of electrochemical interfaces between carbon and ions or liquid electrolytes. Although they provide qualitatively similar results, some discrepancies may be observed when analyzing quantitative properties, such as adsorption energies, preferential "binding" distances or adsorption profiles. For example, as reported by Valencia et al. [START_REF] Valencia | Lithium adsorption on graphite from density functional theory calculations[END_REF] , in the case of bare lithium adsorption on graphite surface, eDFT calculations always display binding energies larger than 1 eV with a preference for the center of the C 6 hexagonal rings (hollow sites) with respect to adsorption above C atoms (top sites) or C-C bonds (bridge sites). [START_REF] Valencia | Lithium adsorption on graphite from density functional theory calculations[END_REF] The lithium atom loses a valence electron which is entirely transferred to the carbon a) Electronic mail: mathieu.salanne@sorbonne-universite.fr b) Electronic mail: carlo.pierleoni@aquila.infn.it

surface. Yet, the adsorption energy may differ by more than 0.5 eV when changing the exchange-correlation functional. [START_REF] Rytkönen | Density functional study of alkali-metal atoms and monolayers on graphite (0001)[END_REF] In particular, the inclusion of dispersion effect further increases the binding energy of the lithium ion. [START_REF] Fan | Adsorption of single li and the formation of small li clusters on graphene for the anode of lithium-ion batteries[END_REF] Similar problems arise in classical MD simulations due to the choice of different interaction potentials. On the one hand, the choice of the electrostatic model for the electrode atoms, which can be treated using fixed-charges or the constant potential method, [START_REF] Scalfi | Molecular simulation of electrode-solution interfaces[END_REF] will impact the localization of the charge induced by the ions on the carbon. On the other hand, several force fields may be used for the short-range repulsion and dispersion interactions.

Important discrepancies were already reported in the case of the adsorption of water molecules on carbon surfaces. [START_REF] Striolo | The carbon-water interface: Modeling challenges and opportunities for the water-energy nexus[END_REF] In two studies on carbon nanotubes [START_REF] Al-Hamdani | How strongly do hydrogen and water molecules stick to carbon nanomaterials?[END_REF] and other carbon nanostructures [START_REF] Brandenburg | Interaction between water and carbon nanostructures: How good are current density functional approximations?[END_REF] , Michaelides and co-workers compared the results of eDFT calculations using a large set of exchange-correlation functionals with diffusion Monte Carlo (DMC) calculations. DMC is a Quantum Monte Carlo (QMC) method, which explicitly accounts for electronic correlation and exchange. Although for fermions the fixed node error prevent DMC from being exact, it is more accurate than eDFT and belongs to the class of variational methods that can be systematically improved with enough effort. This method is much more costly than conventional eDFT in terms of computational time, but it is expected to capture with a high accuracy the intermolecular interactions. It can thus be viewed as a good reference for benchmarking purposes. Based on these results [START_REF] Al-Hamdani | How strongly do hydrogen and water molecules stick to carbon nanomaterials?[END_REF][START_REF] Brandenburg | Interaction between water and carbon nanostructures: How good are current density functional approximations?[END_REF] , it appears that the choice of the exchangecorrelation (XC) approximation in eDFT is crucial but also that it is difficult to predict the accuracy since the adsorption energy results from a subtle balance of the interactions, especially at medium range.

In this work we follow a similar approach for the adsorption of lithium on graphite. The choice of the system was made based on its relevance for energy applications and because the lithium ions has a small number of valence electrons, thus allowing the cost of the reference DMC calculations. We computed the adsorption profiles on hollow, top and bridge sites by varying the distance between the lithium and the carbon material. DMC results are then compared with several XC approximations.

Yet, knowing the adsorption energy of a single ion may not be enough. Indeed, the solvent molecules within the electrolyte will also interact strongly both with the surface and with the ions, which can result in large variations of the adsorption properties under realistic conditions. [START_REF] Merlet | Influence of solvation on the structural and capacitive properties of electrical double layer capacitors[END_REF] The main quantity to be determined is then the free energy of adsorption profile, which should be obtained by sampling the whole liquid degrees of freedom. This quantity is not accessible to electronic structure methods due to their large computational cost, and it is necessary to resort to classical mechanics-based methods instead. Classical MD is usually the method of choice for such purposes, but it suffers from high inefficiency for systems studied under large dilution conditions. In such cases, it may conveniently be replaced by molecular DFT (MDFT). MDFT describes a liquid by its density field, and a functional of this density is introduced to account for entropic, external and solvent-solvent interactions. The external contribution arises from the species in contact with the liquid, which are the graphite electrode and the lithium ion in the present study. The equilibrium density and its related free energy are computed through a variational procedure, which reduces the computational cost by several orders of magnitude w.r.t atomistic MD. Recently, we have extended the MDFT method to account for constant potential electrodes, and we showed it was providing the correct structure and thermodynamics for liquid water in contact with graphite electrodes. [START_REF] Jeanmairet | Study of a water-graphene capacitor with molecular density functional theory[END_REF] In the second step of the present work, we build upon these developments to study the adsorption of lithium on graphite, in the presence of liquid water. The lithiumcarbon interaction is parameterized using the electronic DFT calculations benchmarked on DMC. We show that the energy minimum observed in vacuum vanishes in the presence of solvent, and that the lithium ion does not show anymore a preferred binding distance close to the surface due to the presence of water molecules.

II. COMPUTATIONAL METHODS

A. Quantum Monte Carlo QMC simulations are a set of stochastic computational methods for the evaluation of observables of quantum systems. The fundamental idea behind these techniques is that expectation values of physical quantities of quantum systems can be written as

A = dR Ψ * (R) ÂΨ(R) (1)
where  is a generic observable, Ψ(R) is the (normalized) wave function for the quantum system under analysis , and R represents the 3N -dimensional electronic coordinate. The integral in Equation 1 can then be evaluated using a Monte Carlo sampling [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF] . A Markov chain Monte Carlo simulation consists in sampling a collection of configurations of the system used to estimate configurational integrals such as the one in Equation 1.

In order to obtain meaningful information on the physical system it is important to use a proper wave function. Since the exact ground state for physical systems of interest is unknown it is necessary to guess the form of the ground state wave functions. A large variety of QMC methods exists, implementing different strategies. In this study two QMC techniques were used: Variational Monte Carlo (VMC) and DMC, using the fixed-node approximation.

In VMC we approximate the ground state of the system by defining a trial wave function Ψ T (p; R), which is then plugged in Equation 1; the resulting integral is then computed using a generalised Metropolis sampling. The trial wave function has a functional form that is chosen by taking into account the physical features of the system, and depends on a set of variational parameters p, which are optimized to obtain an approximation of the exact ground state. To optimize the trial wave function we use the variational principle in quantum mechanics: for any given Hamiltonian the ground state has by definition the lowest possible energy. What is done in practice is computing the energy using the trial wave function

E p = dR Ψ * T (p; R) ĤΨ T (p; R) = = dR |Ψ T (p; R)| 2 E loc (p; R) (2)
where we used the local energy E loc (p; R) defined as

E loc (p; R) = ĤΨ T (p; R) Ψ T (p; R) ; (3)
By varying the parameters p we can minimize the energy E p , thus finding an approximation of the exact ground state within the chosen class the trial functions. Another property that is used in wave function optimization is the zero-variance principle: for a given Hamiltonian, the energy variance of each energy eigenstate is zero. In order to use both the variational and the zero-variance principles the quantity that is optimized in practice is a linear combination of the energy and its variance. While VMC allows estimates of ground state properties starting from a relatively simple trial wave function, VMC results are typically not realistically accurate. In order to improve these results projection methods have to be used. Projection methods are a class of QMC methods that use imaginary time evolution to filter out excited states contribution from the trial wave function, allowing accurate computation of ground state properties. While in principle projection techniques are exact, the fermionic sign problem in practice prevents the exact evaluation of quantum observables; to avoid the sign problem the fixednode approximation is used. With this approximation the results will still be approximate, but they will remain variational, and significantly more accurate than the ones obtained with VMC (i.e. the computed energies and energy variance will be lower).

Several different projection QMC methods exist. Among those, one of the most commonly used in electronic structure calculations is DMC. If we consider a wave function Ψ(R), a Hamiltonian with a potential energy term U (R) and we define the imaginary time as τ = it in atomic unit we have

ĤΨ(R) = - 1 2 ∆Ψ(R) + U (R)Ψ(R) = ∂Ψ(R) ∂τ ; (4)
the basic idea at the foundation of DMC is that there is a strong analogy between the imaginary time Schrödinger equation and a diffusion equation (coming from the kinetic part of Ĥ) with an additional branching term (from U (R)). It is thus possible to simulate the imaginary time propagation by using the process of diffusion and branching of classical random walkers, initially distributed according to a trial wave function (obtained via the optimization procedure described above). The practical details of the implementation of DMC can be found in literature [START_REF] Foulkes | Quantum monte carlo simulations of solids[END_REF][START_REF] Umrigar | A diffusion monte carlo algorithm with very small time-step errors[END_REF] . The important things to stress here is that accuracy of fixed-node DMC result depends on the accuracy of the nodal surface of the trial wave function, and that the results are variational with respect to the nodal location.

All QMC computations in the present work were performed using the QMCpack software [START_REF] Kim | QMCPACK: an open source ab initio quantum Monte Carlo package for the electronic structure of atoms, molecules and solids[END_REF] . We used QMC calculations to determine the adsorption energy of Li atom adsorbed on a graphite substrate as a function of the separation between the Li atom from the surface. If z is the distance between the Li atom and the carbon surface the adsorption energy profile E ads (z) is defined as

E ads (z) = E Li+C (z) -(E Li + E C ) (5)
where E Li+C (z) is the energy of a system made of a graphite substrate with a Li atom at a distance z and E Li and E C are the energy of the isolated atom and graphite respectively. The latter is modeled using two graphene layers made of 50 C atoms each, with an AB stacking (which corresponds to a 5×5×1 supercell). The energy profile was computed for three different setups: with the Li atom lying above the centre of a C hexagon (hollow site), above a C atom (top site) and above a C-C bond (bridge site). All simulations were done using a simulation box with cell parameters (in Å)   12.336 0.000 0.000 6.168 10.683 0.000 0.000 0.000 30.000

 

The height of the simulation box was selected after systematically testing the convergence of the total energy as a function of the amount of vacuum between periodic images in the z direction of the graphite bi-layer using eDFT (Section S1 of the Supporting Information) using the PBE functional. Dispersion corrections [START_REF] Grimme | A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu[END_REF] were then used to relax the carbon bilayer structure for the selected box height of 30 Å, yielding a distance between the two planes of 3.47 Å which was used for the QMC computations. This value is slightly larger than the one provided in previous QMC studies (3.43 Å), [START_REF] Mostaani | Quantum monte carlo calculation of the binding energy of bilayer graphene[END_REF][START_REF] Shin | Nature of interlayer binding and stacking of sp-sp 2 hybridized carbon layers: A quantum monte carlo study[END_REF] but this should not impact the computed adsorption energies. Previous investigations systematically studied convergence with the size of the supercell and the number of twists for single and a bilayer systems [START_REF] Mostaani | Quantum monte carlo calculation of the binding energy of bilayer graphene[END_REF][START_REF] Shin | Nature of interlayer binding and stacking of sp-sp 2 hybridized carbon layers: A quantum monte carlo study[END_REF][START_REF] Ganesh | Binding and diffusion of lithium in graphite: Quantum monte carlo benchmarks and validation of van der waals density functional methods[END_REF] . Here for computational reasons, we assumed the mentioned 5x5x1 supercell and we used a 4x4x1 twist grid (corresponding to eight nonequivalent twists) [START_REF] Lin | Twist-averaged boundary conditions in continuum quantum monte carlo algorithms[END_REF] based on the convergence of the corresponding DFT energy. Although convergence of the DFT total energy with twist grid is reached only for a 8x8x1 k-points grid, the residual effect is less than 1mHa. Moreover in the property we are aiming, namely the absorption energy of the Lithium atom, size effects should largely cancel since the same simulation cell and twist grid were used in all the simulations, including the ones of the isolated atom and substrate. In all VMC and DMC simulations we used trial wave function with a Slater-Jastrow form

Ψ T (p; R) = J(p; R)D(R) (6)
where J(p; R) is a Jastrow term describing electronic correlation, with one and two body terms, and D(R) is a Slater determinant, ensuring the correct fermionic antisymmetry. The single particle orbitals used in the Slater determinant were evaluated using DFT with a PBE [START_REF] Perdew | Generalized gradient approximation made simple[END_REF] functional. The orbital calculations were performed using the Quantum ESPRESSO software [START_REF] Giannozzi | Quantum espresso: a modular and open-source software project for quantum simulations of materials[END_REF][START_REF] Giannozzi | Advanced capabilities for materials modelling with quantum espresso[END_REF] . In the DFT calculations a plane wave basis set was used, with a cutoff at 150 Ry, using norm conserving pseudopotentials for both the Li and C atoms. In QMC simulations the Burkatzki-Dolg-Filippi set of pseudopotentials was used [START_REF] Burkatzki | Energy-consistent pseudopotentials for quantum monte carlo calculations[END_REF] . Wave function optimization is performed using the Linear method [START_REF] Toulouse | Optimization of quantum monte carlo wave functions by energy minimization[END_REF] , and iterated until the optimized energy converges. Only the Jastrow part of the trial wave function is optimized. More information on the Jastrow part of the trial wave function can be found in section S2 of the Supporting Information. We report DMC results obtained using an imaginary time step of τ = 0.01 Ha -1 and a population of 12000 to 16000 random walkers. Previous QMC simulations of a similar system [START_REF] Ganesh | Binding and diffusion of lithium in graphite: Quantum monte carlo benchmarks and validation of van der waals density functional methods[END_REF] shown that using these parameters leads to energies converged to a few mHa per atom. We also note that we are performing QMC simulations for a system whose size was deemed large enough to give an accurate description of Li in proximity of carbon. Size consistent, non local T-moves [START_REF] Casula | Size-consistent variational approaches to nonlocal pseudopotentials: Standard and lattice regularized diffusion monte carlo methods revisited[END_REF] are used, as well as the ZSGMA branching scheme [START_REF] Zen | Boosting the accuracy and speed of quantum monte carlo: Size consistency and time step[END_REF] .

B. Electronic Density functional theory

Electronic DFT calculations were performed using the Quantum ESPRESSO electronic structure code. [START_REF] Giannozzi | Quantum espresso: a modular and open-source software project for quantum simulations of materials[END_REF][START_REF] Giannozzi | Advanced capabilities for materials modelling with quantum espresso[END_REF] To be consistent with the QMC result, an identical simulation cell was considered, consisting of one hundred carbon atoms divided amongst two graphite layers with AB stacking (which corresponds to a 5×5×1 supercell). A kinetic energy cutoff of 40 Ry was used.

We compared the adsorption energy profiles for two series of XC functionals. In the first series we used the LDA, 31 PBE 24 and BLYP [START_REF] Becke | Density-functional exchange-energy approximation with correct asymptotic behavior[END_REF][START_REF] Lee | Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density[END_REF] which neglect the London dispersion interaction. Then we included the latter using either the D2 correction parameterized by Grimme, [START_REF] Grimme | Semiempirical GGA-type density functional constructed with a long-range dispersion correction[END_REF] or through the use of the vdW-DF-C09 functional [START_REF] Thonhauser | Spin signature of nonlocal correlation binding in metal-organic frameworks[END_REF][START_REF] Thonhauser | Van der waals density functional: Selfconsistent potential and the nature of the van der waals bondvan der waals density functional: Self-consistent potential and the nature of the van der waals bond[END_REF][START_REF] Berland | Van Der Waals Forces in Density Functional Theory: A Review of the VdW-DF MethodVan Der Waals Forces in Density Functional Theory: A Review of the VdW-DF Method[END_REF] implemented in the Libxc library. [START_REF] Lehtola | Recent developments in libxc -a comprehensive library of functionals for density functional theory[END_REF] Rappe-Rabe-Kaxiras-Joannopoulos ultrasoft (rrkjus) pseudopotentials [START_REF] Rappe | Optimized pseudopotentials[END_REF] were used for both carbon and lithium atoms.

Several uniform Monkhorst-Pack grids of 1×1×1, 2×2×1, 3×3×1 and 5×5×1 k-points were tested for a single Li distance of 2.4 Å from the graphite surface (see Supporting Information Section S3). The difference between the total energy of the 1×1×1 and 2×2×1 grids is roughly 26 meV and thus the 1×1×1 grid was used for the calculations throughout this work. Additionally, an extended system consisting in four carbon layers instead of two was simulated to check the effect of the number of graphite layers; almost no difference was observed for the absorption energy profile as shown in Section 4 of the Supporting Information.

Like for the QMC calculations, the lithium atom was systematically above the three adsorption sites. The energies were converged at each step to an accuracy of 1 × 10 -6 Ry. To align the various curves, the noninteracting systems (i.e. graphite and lithium atom separately) were computed for each XC functional, and the binding energy was obtained according to Equation 5.

C. Molecular Density Functional Theory

Solvation free energies were computed with MDFT [START_REF] Jeanmairet | Study of a water-graphene capacitor with molecular density functional theory[END_REF][START_REF] Jeanmairet | Molecular Density Functional Theory of Water[END_REF][START_REF] Ding | Efficient molecular density functional theory using generalized spherical harmonics expansions[END_REF] while the polarisability of the graphite sheets was handled using fluctuating Gaussian charges method [START_REF] Jeanmairet | Study of a water-graphene capacitor with molecular density functional theory[END_REF][START_REF] Siepmann | Influence of Surface-Topology and Electrostatic Potential on Water Electrode Systems[END_REF][START_REF] Scalfi | Charge fluctuations from molecular simulations in the constant-potential ensemble[END_REF] . MDFT is a flavor of classical density functional theory (cDFT) developed to study the solvation properties of molecular solutes into molecular solvents such as water or acetonitrile. The solvent is described by its density field ρ(r, ω) which measures the average number per unit volume of molecules with an orientation ω at a given position r. The solute acts a perturbation through an external potential V ext (r, ω) causing the solvent to deviate from the homogeneous bulk fluid.

According to the cDFT principles [START_REF] Evans | The nature of the liquid-napour interface and other topics in the statistical mechanics of non-uniform, classical fluids[END_REF][START_REF] Hansen | Theory of Simple Liquids[END_REF] , there exists a unique functional, F , of the solvent density, ρ, that is equal to the solvation free energy at its minimum which is reached for the equilibrium solvent density. To find an expression for the functional, a common practice is to start by splitting it into the following sum

F [ρ(r, ω)] = F id [ρ(r, ω)] + F exc [ρ(r, ω)] + ρ(r, ω)V ext (r, ω)drdω. (7)
In Equation 7, the first term of the rhs is called ideal and corresponds to the entropic contribution of a noninteracting fluid with the same density. The second term is due to solvent-solvent interaction and is often called the excess term while the last term is due to solute-solvent interaction and thus called the external term.

Exact expressions exist for the ideal and external functionals that can be computed numerically. The excess part, however, requires approximations. It can be expressed as an infinite Taylor expansion around the homogeneous bulk solvent density ρ b ,

F excess [ρ] = - k B T 2 ∆ρ(r, ω)c (2) (r -r , ω, ω ) ∆ρ(r , ω )drdωdr dω + F B [ρ]. (8)
In Equation 8, ∆ρ(r, ω) = ρ(r, ω) -ρ b , k B is the Boltzmann constant, T is the temperature and c [START_REF] Dou | Hard carbons for sodium-ion batteries: Structure, analysis, sustainability, and electrochemistry[END_REF] is the two-body direct correlation function of bulk solvent. F B is so-called bridge functional that collects all the terms higher than quadratic, involving many-body direct correlation functions of the bulk solvent. A common way to approximate the excess functional is to ignore the bridge functional, i.e. F B = 0, resulting to the "HNC" functional because it is equivalent to using the hypernetted chain (HNC) closure for the solute-solvent correlations in the molecular Ornstein-Zernike equation [START_REF] Ding | Efficient molecular density functional theory using generalized spherical harmonics expansions[END_REF] . In this work, we use a very simple bridge functional [START_REF] Borgis | Simple Parameter-Free Bridge Functionals for Molecular Density Functional Theory. Application to Hydrophobic Solvation[END_REF][START_REF] Borgis | Accurate prediction of hydration free energies and solvation structures using molecular density functional theory with a simple bridge functional[END_REF] based on weighted density approximation that is known to correct well for the dramatic pressure overestimation of the HNC approximation. Water is modeled with the SPC/E force field while the external potential V ext is created by the graphite electrodes and the lithium ion whose interactions are described as the sum of Lennard-Jones and electrostatic interactions. The Lennard-Jones and charge parameters of the 4 types of atoms involved are collected in Table I. For the carbon atoms, the choice of the force field of Werder et al. [START_REF] Werder | On the Water-Carbon Interaction for Use in Molecular Dynamics Simulations of Graphite and Carbon Nanotubes[END_REF] was made based on a previous QMD study, in which it was shown to provide a good estimate of the energy of adsorption of a water molecule on a graphene surface. [START_REF] Al-Hamdani | How strongly do hydrogen and water molecules stick to carbon nanomaterials?[END_REF] The electrode charges are represented with Gaussian charge distributions with a width of 0.40 Å. These partial charges are calculated by enforcing a uniform potential within the whole carbon electrode, with an overall electroneutrality constraint (hence forcing the total charge on the carbon to be equal to -1). [START_REF] Siepmann | Influence of Surface-Topology and Electrostatic Potential on Water Electrode Systems[END_REF] . Electrode charges are optimized self consistently with the functional minimization through an iterative scheme [START_REF] Jeanmairet | Study of a water-graphene capacitor with molecular density functional theory[END_REF] . In the first step, the functional of Equation 7is minimized with no charges on the lithium and the carbon atoms. Then, carbon charges are optimized in the presence of the inhomogeneous water charge density and of the lithium cation. The functional is minimized again but in the presence of lithium charge and the previously determined electrode charges. The process is repeated until it converges, with a convergence criterion of 5 × 10 -4 on the relative change in solvation free energy between two consecutive steps. MDFT calculations were performed using an in-house Fortran code and electrode charges are optimized using the constant potential molecular dynamics software MetalWalls [START_REF] Marin-Laflèche | Metalwalls: A classical molecular dynamics software dedicated to the simulation of electrochemical systems[END_REF] . We use a 24.672 × 21.366 × 40 Å 3 simulation box (the unit cell used in QMC is replicated twice in x and y directions) with 74 × 64 × 120 grid nodes and an angular grid of 196 orientations per grid node. We run calculations for a distance z between the electrode plane and the lithium varying between z = 1.0 Å and z = 10 Å with an increment of 0.2 Å between 1 Å and 6 Å and of 0.5 Å otherwise.

III. RESULTS AND DISCUSSION

A. Benchmark of the DFT functionals with QMC

We first discuss the QMC results for the three adsorption sites. QMC total energies and adsorption energies are reported in Table II.

From Table II, it appears clearly that the adsorption energies are larger for the top and bridge sites than for the hollow site for distances of 1.5 and 2.0 Å in good agreement with previous eDFT results from the literature and from the current study. At a larger distance of 3.0 Å the three sites display similar energies, which shows that the difference between them has a short-range character. Due to the high computational cost of DMC, further lithium-carbon distances were only considered for the hollow site. We obtained a binding energy E b of -1.08 eV for a lithium-surface distance of 1.8 Å. We also observe a somewhat peculiar behavior since the adsorption energy is rather similar for distances of 3 and 4 Å. By analyzing the corresponding electronic densities as shown on Figure 1A, we observe that this correspond to the region in which the electron transfer occurs. For distances lower than 3 Å the energies correspond to the adsorption of a lithium ion on a polarized carbon surface while for distances greater than 4 Å the system corresponds to neutral lithium atom and carbon material. The adsorption profiles obtained for the hollow site for the LDA, BLYP, PBE, BLYP+D2, PBE+D2 and dW-DF-C09 functionals are compared with the DMC benchmark on Figure 1B. The LDA results in a strong overbinding, which is expected. The comparison with the other functionals is more surprising. The QMC results lie between the BLYP and the BLYP+D2, while all the other functionals predict too low adsorption energies. However, as noted by Valencia et al., the binding energy (defined as the minimum of the adsorption energy profile) should be approximately given by 7

E b ≈ E b (Li + ) -(IP [Li] -W f [Graphite]) (9)
where E b (Li + ) is the binding energy of the lithium ion, IP its ionization potential and W f the work function of graphite. The observed variation may therefore be due to different values for IP and W f from the various functionals.

By visualizing the electronic density of the system using various DFT functionals (Supplementary Figure S6), we also observe that they yield similar results as the QMC calculations up to 3 Å while some discrepancies are observed at a larger distance of 4 Å. In the present work, we

z = 4 Å z = 3 Å z = 1.8 Å 0 1 2 3 4 5 6 7 8
Li-C distance (Å) are mostly interested in developing an accurate potential for the adsorption of the lithium ion. Consequently, we performed a second comparison of the various functionals in which the energy at z = 3 Å is subtracted. The results are shown on Figure 1C. The discrepancy between the various functionals is somewhat lower. The best agreement is now obtained with BLYP, PBE and the vdW-DF-C09 functionals, while the others predicts overbinding. For the two former approximations, the good agreement may be fortuitous since they do not account for the dispersion effects. However, this points towards an overestimation of the dispersion effects when using the D2 correction. Indeed, in the case of the lithium ion, only two semi-core electrons take part in the interaction, which should result in a very weak dispersion term. The vdW-DF-C09 functional, which accounts for these effect explicitly and not through a parameterised term, seems to better catch the interactions.

B. Fitting the carbon-lithium potential

In order to incorporate solvent effects, it is necessary to develop accurate classical interaction potentials. It is not possible to fit them directly on the QMC calculations due to the limited number of data. Instead, we pick the most accurate functional, vdW-DF-C09, and calculate the Ligraphite binding energy for a large number of distances. The intermolecular interaction should in principle account for four different effects: electrostatics, polarization, shortrange repulsion and dispersion. In our electrostatic model, the two former effects are explicitly introduced through the use of a +1 point charge on the lithium and of the calculation of partial (Gaussian) charges on the carbon atoms. These partial charges are calculated for each lithium-carbon distance using the same methodology as described in II C.

Concerning the short-range repulsion and the dispersion effects, the two main potentials used in the literature are the Lennard-Jones and the Born-Huggins-Mayer (BHM) ones. However, it appears that the electrostatic interaction was sufficient to account for the attractive part of the binding energy, as shown on Figure 2. The fitted potential should therefore add very few, if no contribution for the dispersion interaction, which agrees with the previous observation on the use of dispersion-corrected functionals. A well-known drawback of the Lennard-Jones potential is that it is not possible to fit the short-range repulsion and the dispersion term separately since they both involve the same parameters. We have therefore chosen a BHM potential instead, which analytical form is:

V BHM (r) = A exp (-br) - C 6 r 6 (10)
where A, b are the parameters describing the intensity and the range of the repulsion interaction, while C 6 is the dipole-dipole dispersion interactions. In principle higher order terms could be included for dispersion, but as discussed above this term is almost negligible in the case of the lithium ion. The fitted potential reproduces with a very high accuracy the vdW-DF-C09 for the three types of adsorption sites as shown on Figure 2. It also matches very well with the BLYP (not shown on the figure) as well as the QMC results (on which it was not fitted) adsorption curves. The corresponding parameters are A = 91.17, b = 2.518 and C 6 = 1.107 (all numbers are given in atomic units).

C. Adsorption of the lithium ion in the presence of water

The fitted potential can directly be used in any classical molecular simulation, such as MD. Since we focus here on the adsorption free energy of the lithium ion on the carbon surface we prefer to use MDFT which is a computationally more efficient alternative. The solvation free energy of a single system can be computed within a few minutes on a single CPU while it would require tens of CPU hours with MD. The free energy profile obtained for the three adsorption sites in the presence of liquid water are shown on Figure 3A. The profiles are very different from the gas phase results. The minimum at ≈ 1.8 Å completely disappears and is replaced by a strongly repulsive wall. This shows that there is no preferential adsorption of lithium on the graphite surface in aqueous phase. This result is in qualitative agreement with a recent MD study on the adsorption of hydrated ions on graphene, [START_REF] Loche | Breakdown of linear dielectric theory for the interaction between hydrated ions and graphene[END_REF] which provided a repulsive free energy profile over the whole range of considered distances. Yet, the latter study did not include any Coulombic interaction between the ion and the carbon surface, which is the main driving force for adsorption in the gas phase as discussed above. It is thus interesting to examine the various contributions to the total free energy, which are provided on Figure 3B. We observe that the electrostatic attraction between Li and C is in fact counterbalanced by the solvation free energy. The latter contains the electrostatic interactions of the water molecules with both the lithium and the graphite surface, which results in strong screening effect. The extent of this screening effect was studied in details in a recent study focused on gold surfaces: [START_REF] Pireddu | A molecular perspective on induced charges on a metallic surface[END_REF] The presence of the water molecules strongly impacts the polarization of the surface. Consequently, the total free energy is almost equal to the BHM contribution over the whole range of distances, except between 2 and 5 Å where the solvation free energy overcomes the ion-surface Coulombic interaction, resulting in a more repulsive potential.

The effect of the solvent can be further analyzed by plotting the density profiles for various distances between the ion and the surface (Figure 4). At z = 8 Å, two regions with larger densities emerge, corresponding to the surface adsorbed water molecules at a distance of 3 Å from the surface 14 on the one hand and to the lithium ion first solvation shell on the other hand. At a distance of 5 Å the solvation shell starts to overlap with the adsorbed layer at close contact to the electrode, which results in small depletion zones in the latter. These depletion zones remain observable at smaller distances, but the impact on the free energy becomes negligible w.r.t the large shortrange repulsion between the carbon and the lithium.

IV. CONCLUSION

Two main results were obtained in this work. Firstly, we provide an accurate set of benchmark energies for the adsorption of lithium on a graphite surface, using an accurate QMC approach. This data can be used in future works to test the accuracy of XC functionals for such important problems in surface science, or to parameterize new force fields between the two species.

Here we obtained a very good agreement using a BHM potential, but different approaches may be proposed in the future. Secondly, we used the parameterized force field in order to compute the free energy profile for the adsorption of a lithium ion at the graphite surface in the presence of water as a solvent. MDFT was preferred over classical MD for this calculation due to its much lower computational cost. We showed that the low energy minimum obtained in the gas phase completely vanishes, resulting in an overall repulsive profile over the whole range of studied lithium-carbon distances. This is due to the screening of the attractive Coulombic term by the water molecules on the one hand, and on the other hand to the interferences between the densities corresponding to the first adsorbed layer on carbon and to the first solvation shell of the cation.

The proposed multi-scale approach may be extended to a large variety of systems in the future. In particular, it would be interesting to study how adsorption varies within the alkali and alkali-earth cationic series. The study of larger organic ions would be very useful for the scope supercapacitors devices, but this would require additional work to account accurately for the flexibility of the ion. Finally, the study of more polarizable electrode materials such as gold or platinum would be relevant in the context of catalysis. Changing metal would lead to very different adsorption properties for the solvent, so that one can expect large differences in the adsorption free-energy profiles of ions at electrodes.

FIG. 1 :FIG. 2 :

 12 FIG.1: A) Electron density around a Li atom adsorbed on a graphite substrate, computed via quantum Monte Carlo, for different atom-substrate distances. The densities were obtained by computing the overall electronic density of a system with a Li atom adsorbed on graphite, and subtracting the density of the isolated substrate, in absence of the adsorbed atom. All shown isosurfaces correspond to a density of 6 • 10 -4 electrons/Å 3 . B) Comparison of the adsorption energies obtained with various DFT functionals and DMC for the adsorption of the lithium on the hollow site of graphite. C) Same as B) but substracting the adsorption energy at a distance of 3 Å.

FIG. 3 :

 3 FIG.3: A) Adsorption free energy for a lithium ion on the graphite surface in the presence of water, computed using MDFT, for the three adsorption sites. The energy variation in the absence of water is also shown for comparison. B) Contributions to the total free energy for the hollow adsorption site in the MDFT calculation.

FIG. 4 :

 4 FIG.4: Projection of the solvent densities computed using MDFT for various lithium-carbon distances.

TABLE I :

 I Force-field parameters used in the MDFT simulations. Mixed parameters are computed using the Lorentz-Berthelot rules (except for the C-Li interaction, which does not affect the MDFT results).

TABLE II :

 II QMC energies for the adsorption of lithium on graphite, for the three different sites for several lithium-carbon distances reported in the first column (in Å). Second column are VMC results while third column are DMC results. The last column are the adsorption energies computing with DMC using equation 5.

ACKNOWLEDGMENTS

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 771294). This work was granted access to the HPC resources of CINES under Allocation A0100910463 made by GENCI and benefited from the 'Grand Challenge Jean Zay' program. We acknowledge support from EoCoE, a project funded by the European Union Contracts No. H2020-EINFRA-2015-1-676629 and H2020-INFRAEDI-2018-824158.

DATA AVAILABILITY

The input files used in this study for the QMC, DFT and classical simulations are openly available in the repository https://gitlab.com/ampere2, together with the reference energies for the QMC calculations.