Accurate quantum-chemical fragmentation calculations for ion–water clusters with the density-based many-body expansion

The many-body expansion (MBE) provides an attractive fragmentation method for the efficient quantum-chemical treatment of molecular clusters. However, its convergence with the many-body order is generally slow for molecular clusters that exhibit large intermolecular polarization effects. Ion–water clusters are thus a particularly challenging test case for quantum-chemical fragmentation methods based on the MBE. Here, we assess the accuracy of both the conventional, energy-based MBE and the recently developed density-based MBE [Schmitt-Monreal and Jacob, Int. J. Quantum Chem. 120 , e26228 (2020)] for ion–water clusters. As test cases, we consider hydrated Ca 2+ , F − , OH − , and H 3 O + , and compare both total interaction energies and the relative interaction energies of different structural isomers. We show that an embedded density-based two-body expansion yields highly accurate results compared to supermolecular calculations. Already at the two-body level, the density-based MBE clearly outperforms a conventional, energy-based embedded three-body expansion. We compare different embedding schemes and find that a relaxed frozen-density embedding potential yields the most accurate results. This opens the door to accurate and efficient quantum-chemical calculations for large ion–water clusters as well as condensed-phase systems.


Introduction
Water plays an essential role in chemistry and is a prerequisite for life as we know it [1].
Its exceptional properties have puzzled scientists for generations, and many fundamental questions concerning the properties of water are still unresolved (see, e.g., Refs.2-5 and references therein).Many chemical and most biological processes crucially depend on water as a solvent [6], and solvent effects intricately determine both the structure [7][8][9] and spectroscopic properties [10] of solvated molecular systems.In turn, solutes change the structure and properties of liquid water by influencing water-water interactions [11][12][13].
For studying these interactions at the molecular level, clusters of water as well as ionwater clusters present ideal model systems [14,15].In the past decades, such clusters have been investigated both spectroscopically (see, e.g., Refs.[16][17][18][19][20] as well as computationally (for reviews, see, e.g., Ref. [21,22]).Despite increasing computational resources as well as methodological advances, the accurate, yet sufficiently efficient quantum-chemical treatment of large molecular clusters remains a challenging task.
Fragmentation methods [23][24][25][26][27], which partition such clusters into their molecular constituents, offer an attractive approach for overcoming this bottleneck.The most straightforward realization of such a fragmentation method is the well-known many-body expansion (MBE) (see, e.g., Refs.[28][29][30], which approaches the total interaction energy of a molecular cluster as a sum of n-mer (dimer, trimer, tetramer etc.) interaction energy contributions.If the MBE can be truncated at a sufficiently low order, it provides a low-scaling and highly parallelizable computational approach for the accurate quantumchemical treatment of molecular clusters.Generally, a truncation at second order (i.e., a two-body expansion) would be highly desirable.In addition, the MBE provides the basis for the development of accurate force-field models, in particular for water [31][32][33][34] and ion-water systems [35][36][37].
Ion-water clusters present a particularly challenging test case for the MBE [14,21,33,[38][39][40].There have been numerous computational studies of ion-water clusters that employ a MBE (for recent examples, see Refs.35,36,[39][40][41].Because of the large polarization effects, three-body and four-body effects generally play an important role and the conventional MBE cannot be truncated at low order.To some degree, this can be alleviated with more general fragmentation schemes based on overlapping fragments [42,43]. Recently, starting from the observation that a MBE of the electron density converges faster than a MBE for the total energy, our research group developed a density-based MBE (db-MBE) [44].Instead of expanding the total density, the db-MBE is based on an expansion of the total electron density of a molecular cluster in terms of electron densities of embedded monomers, dimers, trimers etc.This expansion can then be used to obtain a density-based correction to the total energy.We could show that such a db-MBE can capture many-body polarization effects in water clusters already with a two-body expansion and that accurate total and relative energies can be obtained with a densitybased two-body expansion [45].Here, we set out to explore the db-MBE for ion-water clusters, which present an even more challenging test case.We will particularly consider the importance of self-consistent embedding for the fragment calculations, which can be expected to be more important in ion-water clusters.

Computational Methodology 2.1 (Embedded) Many-Body Expansion
To treat a cluster consisting of N molecular fragments with the help of the MBE, instead of doing a quantum-chemical calculation for the full system, calculations are performed for the N monomers, the N (N − 1)/2 dimers that can be formed from them, and possibly also for trimers, tetramers, etc.A property of interest P is then calculated for each of the monomers, P I , for each dimer, P IJ , for each trimer P (3) IJK and so on.Here, the superscript indicates the size of the considered subsystem and the subscript lists the monomers included in it.
The corresponding property P tot of the full cluster is then approximated as [29,30], where ∆P (2) IJK is a trimer interaction contribution, and ∆P • is an n-mer interaction contribution.For explicit expressions for the n-mer interaction contributions we refer to, e.g., Ref. 46.The above expansion is truncated at order n.For such a truncated MBE, the number of subsystem calculations scales as O(N n ).Usually, one aims at truncating the MBE at low order (ideally, already at second order), while the exact, supermolecular result is recovered by construction for n = N .
In the simplest case, the calculations of P (1) IJ , P IJK , . . .are performed for the isolated active subsystem without considering the remaining fragments (isolated MBE ).To accelerate the convergence of the MBE, the effect of the remaining fragments on the active monomer, dimer, trimer etc. can be approximated by means of a suitable embedding potential [47][48][49][50] (embedded MBE ).Here, we consider different local embedding potentials , where the superscript indicates the fragments included in the active subsystem.Note that the inclusion of an embedding potential does not change the scaling of the computational effort of the MBE.
First, we consider a purely electrostatic point-charge (PC) embedding potential, where q k are partial charges placed at the positions of the nuclei R k .Here, we will adopt the point charges provided by suitable force-field models.
Second, we consider a more sophisticated embedding potential depending on the electron densities ρ K (r ′ ) of the fragments in the environment, i.e., the embedding potential of frozen-density embedding (FDE) theory [51,52], with This embedding potential includes the full electrostatic nuclear and electronic Coulomb potentials of the monomers in the environment as well as non-classical contributions due to the exchange-correlation and kinetic energy, which are evaluated using approximate functionals.For obtaining the electron densities of the fragments in the environment, the simplest variant is the use of those obtained in calculations for the isolated monomers.Since for ion-water clusters large polarization effects can be expected, we also iteratively updated all embedded monomer densities using freeze-andthaw (ft) cycles.While the calculation of the relaxed monomer densities comes with some computational effort, this only needs to be done at the level of the monomers, i.e., the scaling of the computational effort remains unchanged.
We will refer to these different variants of the MBE by the following acronyms in the following: 'iso' is used for the isolated MBE, 'PC' refers to the embedded MBE using point-charge embedding, 'FDE' indicates that an embedded MBE with a frozen-density embedding potential based on the unrelaxed frozen densities of the isolated monomers has been used, whereas 'FDE-ft' denotes the use of a frozen-density embedding potential based on relaxed monomer densities optimized self-consistently in freeze-and-thaw cycles.

Energy-Based Many-Body Expansion
Conventionally, the MBE is performed for the total energy of the system, which is approximated as We will refer to this expansion as energy-based many-body expansion (eb-MBE), and to its truncation at order n as eb-MBE(n).
Here, the energies E (1) IJ , E IJK , . . .are obtained from single-point calculations for the monomers, dimers, trimers, etc.In case of an embedded MBE, it is important to ensure that these energies refer to the active subsystem only and that they do not include the interaction with the environment, i.e., the interaction energy or similar contributions are not included in the active subsystem energies [44].

Density-Based Many-Body Expansion
Similarly, the total electron density can be approximated using the MBE as, This expansion forms the basis for the density-based many-body expansion [44], in which the total energy is approximated in the spirit of an ONIOM-style [53,54] composite scheme as Here, is the well known total energy functional of Kohn-Sham density-functional theory (DFT) containing the noninteracting kinetic energy T s [ρ], the electron-nuclei attraction energy tot denotes the n-body expansion of the Kohn-Sham DFT total energy according to Eq. ( 1).Thus, the quantum-chemical method used in the eb-MBE constitutes the high-level method, while the energy functional E tot [ρ] defines the low-level method.
The density-based energy correction E (n) db-corr in Eq. ( 7) is a functional of the monomer, dimer, trimer etc. densities that is given by with the n-body nonadditive kinetic and exchange-correlation energy functionals, and where s , and xc are the n-body expansions according to Eq. ( 1) of the corresponding energy contributions, which are each functionals of the monomer, dimer, trimer etc. electron densities.
The first three terms in the above density-based energy correction E (n) db-corr correct the electrostatic interaction energies such that they correspond to the full n-body electron density.Note that for the nuclear-nuclear attraction energy, a correction only appears at first order.The nonadditive kinetic and exchange-correlation energies account for non-classical interactions.In the spirit of subsystem-DFT [52], these contributions are evaluated using approximate, semi-local density functionals.For further details on the evaluation of the different terms in Eq. ( 9), we refer to Refs.44, 45.We note that the calculation of the density-based energy correction is only done once after the individual subsystem calculations required in the MBE.Therefore, the scaling of the computational effort with cluster size is the same for the eb-MBE(n) and the db-MBE(n), even though in the latter case the need for calculating the subsystem electron densities increases the prefactor [44,45].

Computational Details
All quantum-chemical calculations have been performed using DFT with the Amsterdam Density-Functional (ADF) engine of the Amsterdam Modeling Suite (AMS) [55,56] with the PBE0 hybrid functional [57][58][59] and a triple-zeta plus polarization (TZP) basis set of Slater-type orbitals [60].A Becke integration grid of "good" accuracy [61] and the numerical accuracy setting "good" have been used throughout.All total energies have been obtained with ADF's total energy implementation [62].
For the point-charge embedded MBEs we used TIP3P [63] point charges (q H = +0.417and q O = −0.834)for the environment water molecules.For the Ca 2+ and F − , point charges were set to +2 and −1, respectively.For the hydroxide ion, we selected partial charges of q H = +0.183and q O = −1.183[64], and for the H 3 O + ion we chose q H = +0.524and q O = −0.571[65].For the frozen-density embedded MBEs we employed ADF's implementation of FDE [66] with the PW91k nonadditive kinetic-energy functional [67] and the PBE [57] functional for the nonadditive xc contributions.In the case of the FDE-ft embedded MBEs, three freeze-and-thaw-cycles were performed using monomer subsystems, and the resulting relaxed monomer electron densities were then used to assemble the frozen density in the subsequent calculations of the monomer, dimer, etc. subsystems.
The eb-MBE and db-MBE have been implemented in the PyADF scripting framework [68,69], which was used for all calculations presented in this work.PyADF includes a PyEmbed module, which provides a stand-alone implementation of the subsystem DFT embedding potential and interaction energy terms.For the nonadditive xc and kinetic energy contributions, it makes use of the XCFun library [70,71].Here, we employed the PW91k kinetic-energy functional [67] and the PBE [57] xc functional for the nonadditive energy contributions.The electrostatic interaction energies are evaluated by PyEmbed using numerical integration as described in Ref. 44.The most recent version of PyADF, which can be used for reproducing all calculations presented here, is available at [69].
All figures have been generated using Matplotlib [72] and Seaborn [73].A data set containing coordinates of all considered molecular structure, PyADF input scripts for executing the eb-MBE and db-MBE calculations, raw data from all eb-MBE and db-MBE calculations as well as Jupyter notebooks for data analysis and for generating all tables and figures contained in this article are available at Ref. 74.

Results and Discussion
To assess the accuracy of the db-MBE for ion-water clusters, we will consider different test cases in the following.First, we investigate the absolute error in the total interaction energies of Ca 2+ -water clusters of increasing size (see Section 3.1).Second, we evaluate the accuracy of the relative energies of different structural isomers of selected F − -water clusters (see Section 3.2), OH − -water clusters (see Section 3.3), and H 3 O + -water clusters (see Section 3.4).In each case, the ion (Ca 2+ , F − , OH − , or H 3 O + ) is treated as a separate, charged fragment in the MBE.For the H 3 O + -water clusters, the additional proton is assigned to the closest water molecule to form the H 3 O + ion.
For all clusters, we perform calculations using the eb-MBE of second and third order [eb-MBE(2) and eb-MBE(3)] and using the db-MBE of first, second, and third order [db-MBE(1), db-MBE(2), db-MBE(3)].As reference, a supermolecular calculation of the full cluster is used.In all cases, we calculate interaction energies, E int = E tot − I E (1) I , i.e., the difference between the total energy of the cluster and the energies of the monomers.This also holds for relative interaction energies, which do not include the contribution that is due to changes in the monomer geometries (see the discussion in Ref. 45).As explained above, we compare different embedding schemes in the MBE calculations, namely using the isolated subsystems (no embedding, labelled "iso"), point-charge embedding (labelled "PC"), and frozen-density embedding both with the unrelaxed frozen densities of the isolated environment molecules (labelled "FDE") and with a relaxed frozen environment density updated in freeze-and-thaw cycles (labelled "FDE-ft").
For comparison, we recall some of our earlier results from Ref. 45, in which we benchmarked the accuracy of the db-MBE for neat water clusters.For the total interaction energies of water clusters of increasing size (up to 55 water molecules), we found that the absolute error of the MBEs increases linearly with the number of water molecules.
For the isolated eb-MBE(2), it amounts to about 10 kJ/mol per water molecule (for B3LYP/DZP).This error is reduced to ca. 6 kJ/mol when using point-charge or frozendensity embedding.The db-MBE(1) provides an accuracy that is comparable to the eb-MBE(2), and with the second-order db-MBE(2) the error per water molecule is reduced to ca. 2 kJ/mol, i.e., well below the threshold of chemical accuracy (1 kcal/mol) for the interaction energy per fragment.
For the relative energies of selected isomers of water clusters of different size, we found that the eb-MBE(2) and even the eb-MBE(3) fail to reproduce the ordering of the isomers and/or their energy differences, even if point-charge or frozen-density embedding is used.
In contrast, the embedded db-MBE(2) generally reproduces both the energy ordering of the considered isomers and their energy differences correctly.Across the test set of water cluster isomers investigated in Ref. 45, the mean and the maximum errors in the relative interaction energies amount to only 0.13 kJ/mol and 0.32 kJ/mol per fragment, respectively.
While we found that the use of a suitable embedding potential generally leads to a clear improvement compared to the corresponding isolated MBEs, for water clusters we noticed that the difference between point-charge and frozen-density embedding is rather small.Therefore, the use of freeze-and-thaw cycles (FDE-ft) was not considered previously and is investigated here for the first time.

Ca
As a first test case, we examine hydrated calcium ions in Ca 2+ -water clusters containing three to twenty water molecules.The structures of these clusters have been taken from Ref. 75.This system has been used recently by Veccham et al. to assess the accuracy of an embedded many-body expansion [76].Up to eight water molecules, these are directly coordinated to the calcium ion, while for larger clusters, a second solvation shell starts to build.
The errors in the total interaction energy per water molecule for the different MBEs in combination with different embedding schemes are plotted in Fig. 1 as a function of cluster size.The corresponding mean and maximum errors across all considered cluster sizes are listed in Table I.
For the isolated eb-MBE(2), the error per water molecule increases up to eight water molecules, where it reaches its maximum of 53 kJ/mol.For larger clusters, it decreases and drops to 23 kJ/mol for the largest considered cluster with 20 water molecules.Over all clusters, the mean error per water molecule amounts to 33 kJ/mol.These results are very similar to those reported in Ref. [76] for the same molecular structures, but with The MBEs have been performed using calculations of the isolated subsystems ("iso"), using point-charge embedding ("PC"), as well as using frozen-density embedding with the unrelaxed frozen densities of the isolated fragments ("FDE") and with the relaxed frozen density obtained using freeze-and-thaw cycles ("FDE-ft").For the isolated db-MBEs, the errors are reduced compared to the isolated eb-MBE(2), but remain rather large.Thus, while for water clusters already the isolated db-MBE (2) could reduce the error per fragment below the threshold of chemical accuracy [45], for ion-water clusters the use of a suitable embedding potential is essential.As soon as some embedding potential is included, the accuracy of the db-MBE improves.For the db-MBE(1), which only requires calculations for the monomers, the mean error amounts to ca.Overall, we find that for Ca 2+ -water clusters, the accuracy of the embedded db-MBE (2) previously observed for neat water clusters [45] is retained and the mean error in the total interaction energies per water molecule remains below the threshold of chemical accuracy.With relaxed FDE-ft embedding, it is significantly below this threshold.For the Ca 2+ -water clusters the proper inclusion of an embedding potential is crucial, and the best results are obtained with the parameter-free, relaxed FDE-ft embedding potential.
However, we note that for the clusters considered here, excellent results can already be obtained with the eb-MBE(2), as long as suitable point-charge embedding or a relaxed FDE-ft embedding potential is included.

F − (H 2 O) 10
As a second test case, we consider clusters of a fluoride anion with ten water molecules, Here, we compare ten structural isomers, whose structures have been taken from Ref. [42].These structures had been generated using molecular dynamics simulations and were used earlier as a test case for an MBE with overlapping fragments.The labels used for the different isomers correspond to those used previously in Ref. [42].
Here, we consider the relative interaction energies of the different isomers, where the lowest-energy isomer (isomer 03) is used as reference.We compare the relative energies calculated with the different MBEs to those obtained from supermolecular calculations in Fig. 2. In this figure, the horizontal dashed lines indicate the supermolecular Included are results from both the isolated MBEs ("iso") and the embedded MBEs using pointcharge embedding ("PC") as well as using frozen-density embedding with the unrelaxed frozen densities of the isolated fragments ("FDE") and a relaxed frozen density obtained using freeze-and-thaw cycles ("FDE-ft").Some data points might lie outside the shown energy range.
reference values, where the line colors refer to the respective isomers.These supermolecular relative energies span a range of roughly 45 kJ/mol.The values calculated with the different MBEs are shown by the colored symbols.In the case of perfect agreement, these symbols would be on top of the horizontal line of the same color.The mean and maximum errors in both the total interaction energies and the relative interaction energies across all ten structural isomers are listed in Table II.For the relative interaction energies, the lowest-energy isomer is not included because its error is zero by definition.
Note that in contrast to the previous section, these errors refer to the full clusters and are not normalized to the number of water molecules.
First, we notice that neither the eb-MBE(2) nor the eb-MBE(3) are able to provide the correct energetic ordering of the isomers.While the energy-based three-body expansion Without the inclusion of an embedding potential, the isolated db-MBE does not provide useful results for the relative energies.However, if an embedding potential is included, already the one-body db-MBE(1) leads to a significant improvement compared to the eb-MBE(2) and eb-MBE(3).It groups the isomers correctly into a lower energy group (isomers 01, 02, 03, 04, 06, and 10) and a higher energy group (isomers 05, 07, 08, and 09), even though the ordering within these groups is not fully correct.This is remarkable, given that the db-MBE(1) uses only calculations for the embedded monomers.
With the embedded db-MBE(2), the relative energies are further improved.The energetic ordering of the isomers is mostly correct, except for some close-lying isomers (e.g., isomers 06 and 10, isomers 05 and 07) and some remaining deviations in the energy differences between isomers.For the db-MBE(2) relative energies, all three embedding schemes perform comparably well.The most accurate results are found for the relaxed FDEft embedding potential with with a mean error of 0.6 kJ/mol and a maximum error of 1.4 kJ/mol in the relative energies.The db-MBE(3) further improves the energy differences between isomers slightly.With FDE-ft embedding, the mean and maximum errors in relative energies amount to 0.7 kJ/mol and 1.1 kJ/mol.
Altogether, already the embedded db-MBE(2) clearly improves compared to the eb-MBE(2) and eb-MBE(3), which both do not provide sufficient accuracy for the relative interaction energies.The db-MBE( 2) is able to reproduce the ordering of the isomers as well as their relative energies accurately.Previously, such an accuracy for this test case was only reached with an eb-MBE based on larger, overlapping fragments [42].It is noteworthy that with FDE-ft embedding, excellent agreement is also reached for the total interaction energies, for which the mean and maximum errors are as low as 0.6 kJ/mol and 1.3 kJ/mol, respectively.

OH
As another test case, we consider clusters of a hydroxide anion with three, four, and five water molecules.For each cluster size, different low-energy structural isomers are considered.The structures of these clusters have been taken from Ref. [35], in which Egan and Paesani analyzed many-body effects in these systems using high-level coupledcluster calculations.They found that in OH − (H 2 O) 5 , three-body effects amount to up to 60 kJ/mol, four-body effects to up to 8 kJ/mol, and five-body effects still contribute up to 2 kJ/mol.The lowest-energy clusters of each size were also used as test case in Ref. [76].
At the DFT level, errors of ca.140 kJ/mol were found for OH − (H 2 O) 5 with an isolated eb-MBE(2), which could be reduced to ca.21 kJ/mol by combining the eb-MBE(2) with a sophisticated quantum-embedding scheme.
The relative interaction energies of the different isomers for each cluster size that are calculated with the eb-MBEs and the db-MBEs are compared to the supermolecular reference values in Fig. 3.The mean and maximum errors in the total and in the relative interaction energies are given in Table III.These refer to all considered clusters (the lowestenergy isomer is excluded for the relative interaction energies) and are not normalized to the cluster size.
As for the previous test cases, it is obvious from Fig. 3 that the use of a suitable embedding potential is essential for ion-water clusters.For the isolated eb-MBEs and db- Included are results from both the isolated MBEs ("iso") and the embedded MBEs using point-charge embedding ("PC") as well as using frozendensity embedding with the unrelaxed frozen densities of the isolated fragments ("FDE") and a relaxed frozen density obtained using freeze-and-thaw cycles ("FDE-ft").Some data points might lie outside the shown energy range.However, the embedded eb-MBE( 2) is not able to reproduce the correct ordering of the isomers at all.For the most accurate FDE-ft embedding, the mean and maximum errors in the relative interaction energies still amount to 8 kJ/mol and 17 kJ/mol, respectively.
The embedded eb-MBE(3) improves and results in the correct ordering of the isomers in most cases, but it still shows significant errors in the energy differences, with mean and maximum errors in the relative energies of 1.9 kJ/mol and 3.3 kJ/mol, respectively, for FDE-ft embedding.
The embedded db-MBE(1) performs worse than for the fluoride-water clusters in the previous section, and largely underestimates the energy differences between the structural isomers.In contrast, the embedded db-MBE(2) results in the correct energetic ordering of the structural isomers and yields rather accurate energy differences.With FDE-ft embedding, the mean and maximum errors in the relative energies are reduced to 0.46 kJ/mol and 1.3 kJ/mol, respectively.The embedded db-MBE(3) slightly improves upon this and reduces these errors to 0.34 kJ/mol and 1.0 kJ/mol, respectively.
Again, we note that the excellent accuracy of the db-MBE(2) and db-MBE(3) does not only hold for the relative interaction energies, but also the corresponding total interaction energies turn out to be highly accurate, with maximum errors of only 2.30 kJ/mol and 0.95 kJ/mol, respectively, for the density-based two-body and three-body expansions.

H
As last test case, we examine protonated water clusters with four to six water molecules.
Again, for each cluster size different structural isomers are considered.The structures have been taken from Ref. [36].For applying the MBE, these clusters are partitioned into one H 3 O + cation and N water molecules.We note that some isomers contain a Zundel ion, H 5 O + 2 , as a substructure, which is nevertheless partitioned into an H 3 O + cation and a water molecule.Previously, it was found that with an isolated eb-MBE, three-body contributions amount to up to 40 kJ/mol, four-body contributions can be as large as 5.5 kJ/mol, whereas five-body contributions become negligible [36].
The relative interaction energies obtained with eb-MBEs and db-MBEs in combination with different embedding schemes are visualized in Fig. 4 and compared to the supermolecular reference data.Mean and maximum errors in the total and relative interaction energies are summarized in Table IV.
Overall, the results confirm the picture obtained in the previous sections.The eb-MBE( 2) is not able to reproduce the energetic ordering of the structural isomers correctly.In contrast, the db-MBE(2) yields highly accurate results when combined with a suitable embedding scheme.For the most accurate relaxed FDE-ft embedding, the mean and maximum errors in the relative energies amount to only 1.9 kJ/mol and 2.8 kJ/mol, respectively.A similar error is reached for the total interaction energies with mean and maximum errors of 1.5 kJ/mol and 2.7 kJ/mol, respectively.These errors of the db-MBE(2) are slightly lower than those found for the more expensive eb-MBE(3).Included are results from both the isolated MBEs ("iso") and the embedded MBEs using point-charge embedding ("PC") as well as using frozendensity embedding with the unrelaxed frozen densities of the isolated fragments ("FDE") and a relaxed frozen density obtained using freeze-and-thaw cycles ("FDE-ft").Some data points might lie outside the shown energy range.

Conclusions
In this work, we have evaluated the accuracy of the eb-MBE and of the db-MBE recently developed in our research group, for ion-water clusters.With the conventional eb-MBE, the errors in the total and relative interaction energies of these systems are significantly larger than for neat water clusters, mainly because of the large polarization effects that are induced by the ion.Thus, ion-water clusters present a very challenging test case for the db-MBE.To avoid the steep scaling of the computational effort with cluster size that plagues higher-order MBEs, our goal is to achieve sufficient accuracy already with a two-body expansion.
Our results show that the two-body db-MBE( 2) is able to provide highly accurate relative and total interaction energies for all test cases considered here.The remaining errors are comparable to those found earlier for neat water clusters [45].Thus, the density-based many-body expansion is able to capture the large polarization effects in these systems already at the level of a two-body expansion.With the conventional eb-MBE, a threebody or even a four-body expansion is required to achieve a similar accuracy.
The db-MBE achieves this by performing the MBE not in terms of the energy, but by expanding the electron density and by applying a density-functional to evaluate a densitybased correction to the total energy of the eb-MBE.This MBE of the electron density accounts for the interactions among the polarized electron densities that only appear at higher orders in the conventional eb-MBE [44].Because at any given order, the eb-MBE and the db-MBE rely on the same quantum-chemical calculations for subsystems, the scaling of their computational effort with cluster size remains comparable.Therefore, accounting for many-body polarization effects already in the db-MBE(2) instead of requiring a higher-order eb-MBE reduces the overall computational effort for large systems and thus pushes the limit of the system sizes that can be treated accurately and efficiently.
Further strategies for reducing the computational cost that have been established for the eb-MBE, such as distance-based screening of higher-order contributions [43,47], are also applicable for the db-MBE and will be explored in our future work.
For ion-water clusters, the use of a suitable embedding potential within the MBE subsystem calculations out to be essential.Here, we compared the use of a point-charge embedding potential (with parametrized partial charges) as well as frozen-density embedding with the unrelaxed monomer frozen densities (FDE) and with monomer densities relaxed in freeze-and-thaw cycles (FDE-ft).While we notice that for the db-MBE the resulting energies are less sensitive to the choice of the embedding potential than for the eb-MBE, we still find that in all cases, FDE-ft embedding leads to the most accurate results.Thus, we recommend to use FDE-ft embedding, which does not rely on any parametrization, in future applications of the db-MBE to systems with large polarization effects.
While we only performed DFT calculations in the present study, the db-MBE can be combined with other quantum-chemical methods, in particular with accurate wavefunctionbased methods.As the number of monomer and dimer calculations required for the db-MBE(2) scales only quadratically with the number of fragments while retaining excellent accuracy, it might allow us to overcome the scaling bottleneck of accurate wavefunctionbased quantum chemistry.This will enable the efficient, yet accurate treatment of large ion-water clusters and possibly also of condensed-phase systems with the db-MBE (2).
for executing the eb-MBE and db-MBE calculations, csv files containing the raw data from all eb-MBE and db-MBE calculations as well as Jupyter notebooks for data analysis and for generating all tables and figures contained in this article, are available at Zenodo at https://doi.org/10.5281/zenodo.7113146.
isolated eb-MBE(3), the trend remains similar, but the errors are reduced by a factor of 2-3.Still, the mean and maximum errors remain substantial and amount to 10 kJ/mol and 21 kJ/mol per water molecule, respectively.The accuracy of the eb-MBE(2) improves when it is combined with an embedding potential in the subsystem calculations.With point-charge embedding (using TIP3P charges for the water molecules), the mean and maximum errors drop to only 1.4 kJ/mol and 4.1 kJ/mol per water molecule, respectively.These errors are already below the threshold of chemical accuracy for the interaction energy per water molecule.With an unrelaxed FDE embedding potential, the errors are larger and the mean and maximum errors amount to 19 kJ/mol and 27 kJ/mol per water molecule, respectively.These errors are comparable to those found in Ref.[76] for a point-charge embedded eb-MBE(2) using point charges determined from the electrostatic potential of isolated water molecules.In contrast, TIP3P point charges already account for the mutual polarization of water molecules in the con-densed phase.The fact that the accuracy of electrostatically embedded eb-MBEs shows a strong dependence on the choice of embedding charges is well established in the literature[76][77][78].Accounting for the mutual polarization of the molecules with the cluster by using a FDE-ft embedding potential based on a relaxed frozen density that was obtained in freeze-and-thaw cycles reduces the mean and maximum errors of the eb-MBE(2) to 1.4 kJ/mol and 2.8 kJ/mol per water molecule, which is slightly lower than for TIP3P point-charge embedding.Importantly, in contrast to point-charge embedding, FDE-ft does not rely on any parametrization.When going from the eb-MBE(2) to the eb-MBE(3), the errors are generally reduced.An exception is point-charge embedding, which indicates that the very low error for the pointcharge embedded eb-MBE(2) might be partly accidental.For FDE-ft embedding, the agreement of the eb-MBE(2) with the supermolecular reference becomes almost perfect, with a mean error of only 0.34 kJ/mol per water molecule.

Figure 2 :
Figure 2: Interaction energies (PBE0/TZP) of ten isomers of F − (H 2 O) 10 , relative to the lowest-energy isomer, as calculated with the energy-based MBE [eb-MBE(n), n = 2, 3] as well as the density-based MBE [db-MBE(n), n = 1, 2, 3].The horizontal dashed lines indicate the corresponding reference values from supermolecular calculations.Included are results from both the isolated MBEs ("iso") and the embedded MBEs using pointcharge embedding ("PC") as well as using frozen-density embedding with the unrelaxed frozen densities of the isolated fragments ("FDE") and a relaxed frozen density obtained using freeze-and-thaw cycles ("FDE-ft").Some data points might lie outside the shown energy range.

Figure 3 :
Figure 3: Interaction energies (PBE0/TZP) of different isomers of OH − (H 2 O) N (N = 3, 4, 5), relative to the respective lowest-energy isomer, as calculated with the energybased MBE [eb-MBE(n), n = 2, 3] as well as the density-based MBE [db-MBE(n), n = 1, 2, 3].The horizontal dashed lines indicate the corresponding reference values from supermolecular calculations.Included are results from both the isolated MBEs ("iso") and the embedded MBEs using point-charge embedding ("PC") as well as using frozendensity embedding with the unrelaxed frozen densities of the isolated fragments ("FDE") and a relaxed frozen density obtained using freeze-and-thaw cycles ("FDE-ft").Some data points might lie outside the shown energy range.

Figure 4 :
Figure 4: Interaction energies (PBE0/TZP) of different isomers of H 3 O + (H 2 O) N (N = 3, 4, 5), relative to the respective lowest-energy isomer, as calculated with the energybased MBE [eb-MBE(n), n = 2, 3] as well as the density-based MBE [db-MBE(n), n = 1, 2, 3].The horizontal dashed lines indicate the corresponding reference values from supermolecular calculations.Included are results from both the isolated MBEs ("iso") and the embedded MBEs using point-charge embedding ("PC") as well as using frozendensity embedding with the unrelaxed frozen densities of the isolated fragments ("FDE") and a relaxed frozen density obtained using freeze-and-thaw cycles ("FDE-ft").Some data points might lie outside the shown energy range.
the Coulomb energy J[ρ], the exchange-correlation (xc) energy E xc [ρ], and the nuclear repulsion energy E NN .Finally, E
a different xc functional and basis set.When going to a three-body expansion with the mean and maximum errors of 1.2 kJ/mol and 2.5 kJ/mol per water molecule, respectively, for point-charge embedding, and of 0.23 kJ/mol and 0.44 kJ/mol per water molecule, respectively, for relaxed FDE-ft embedding.In all cases, this is a substantial improvement both compared to the corresponding eb-MBE(3) and db-MBE(2).
8-11 kJ/mol per water molecule.With the db-MBE(2), the mean errors drop Finally, the embedded three-body db-MBE(3) reduces the errors even further, and re-sults in

Table II :
Mean and maximum error (in kJ/mol) in the total and in the relative interaction energy (PBE0/TZP) of ten isomers of F − (H 2 O) 10 for the eb-MBE(n) (n = 2, 3) and the db-MBE(n = 1, 2, 3) in combination with different embedding schemes.

Table III :
Mean and maximum error (in kJ/mol) in the total and in the relative interaction energy (PBE0/TZP) of different isomers of OH − (H 2 O) MBEs, rather large errors are found.The embedded MBEs generally show much better performance.When comparing the mean and maximum errors in TableIII, we notice that point-charge embedding and unrelaxed FDE show comparable errors, with generally slightly better accuracy for unrelaxed FDE.In all cases, the most accurate results are obtained with the relaxed FDE-ft embedding potential, for which the mean and maximum errors are lower than for PC and FDE by at least a factor of two.

Table IV :
Mean and maximum error (in kJ/mol) in the total and in the relative interaction energy (PBE0/TZP) of different isomers of H 3 O + (H 2 O)