Scalable generalized screening for high-order terms in the many-body expansion: Algorithm, open-source implementation, and demonstration

The many-body expansion lies at the heart of numerous fragment-based methods that are intended to sidestep the nonlinear scaling of ab initio quantum chemistry, making electronic structure calculations feasible in large systems. In principle, inclusion of higher-order n -body terms ought to improve the accuracy in a controllable way, but unfavorable combinatorics often defeats this in practice and applications with n ≥ 4 are rare. Here, we outline an algorithm to overcome this combinatorial bottleneck, based on a bottom-up approach to energy-based screening. This is implemented within a new open-source software application (“Fragme ∩ t”), which is integrated with a lightweight semi-empirical method that is used to cull subsystems, attenuating the combinatorial growth of higher-order terms in the graph that is used to manage the calculations. This facilitates applications of unprecedented size, and we report four-body calculations in ( H 2 O ) 64 clusters that afford relative energies within 0.1 kcal/mol/monomer of the supersystem result using less than 10% of the unique subsystems. We also report n -body calculations in ( H 2 O ) 20 clusters up to n = 8, at which point the expansion terminates naturally due to screening. These are the largest n -body calculations reported to date using ab initio electronic structure theory, and they confirm that high-order n -body terms are mostly artifacts of basis-set superposition


I. INTRODUCTION
Quantum-chemical fragmentation approaches replace an intractable electronic structure calculation with a large number of relatively small subsystem calculations. 12][3][4] The simplest form of the MBE expresses the total energy as n-body corrections to the sum of fragment energies: where is the two-body interaction energy, is a three-body correction, etc. Higher-order terms are expressible in a natural way. 5Truncation of Eq. ( 1) at some finite order, n, affords an approximation that we will call MBE(n).3][34][35][36][37] Finally, MBE(n) forms the basis of the fragment molecular orbital method. 38or a system with N fragments, MBE(n) generates a combinatorial number of n-body subsystems, Assuming that N grows linearly with system size, the asymptotic cost of MBE(n) then grows as O(N n ), absent some kind of screening The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp4][45][46][47][48][49][50][51][52][53][54][55][56][57][58] Very few studies have examined higher-order terms for N ≥ 20 fragments. 43,560][61][62][63] By means of these corrections, or else simply by using very large basis sets to eliminate BSSE, it has been determined that MBE(n) for neat water clusters converges at n = 4, 40,42,[56][57][58] or possibly n = 5. 48Higher-order n-body interactions persist in ion-water clusters 53,54 and in some response properties, 51,52 whereas MBE(n) energies for proteins appear to converge at n = 3, using single-residue fragments. 64These observations suggest that the vast majority of n-body interactions are insignificant so long as BSSE is taken into account.
Distance-based thresholding has often been considered as a means to tame the combinatorics of MBE(n), but for highorder expansions this effectively limits the subsystems to nearestneighbors only.This may omit cooperative dipolar interactions that can sometimes cause energetically-important terms to persist at longer length scales. 39Within a distance-based screening approach, these terms must either be omitted (at the expense of accuracy), or else the distance cutoff must be relaxed, which is detrimental to efficiency.Distance-based screening must furthermore be done carefully (by means of smooth cutoffs) in order to avoid introducing discontinuities in potential energy surfaces. 40,65Ab initio molecular dynamics simulations are quite sensitive to small inconsistencies between the energy and forces, which can lead to catastrophic failure to conserve energy. 66n previous work, 41 we introduced an energy-screened MBE(3) method in which a polarizable force field was used to screen the three-body terms.The target level of ab initio theory was applied to trimer IJK only if where ΔE low IJK indicates the three-body correction evaluated using a low-level method (force field), and the cutoff τ 3 was set to 0.25 kJ/mol following some testing.This approach proves to be far more efficient than distance-based screening, if the cutoffs are set (in either case) so as to achieve ∼1 kcal/mol fidelity with respect to a supersystem calculation at the same level of theory.Four-body terms were omitted in our previous work, 41 with the effect that 1 kcal/mol fidelity requires an additional supersystem correction at a low level of theory, in order to capture four-body polarization that is neglected by MBE(3).
Inclusion of higher-order terms may obviate the need for such a correction, thus the present work is an effort to implement energy screening in a general way at all orders of MBE(n).We describe a new "bottom-up" algorithm, by means of which an arbitrary screening procedure may be extended from low-order to high-order nbody subsystems, for arbitrary n, without the need to define bespoke screening algorithms for different n-body levels.The bottom-up algorithm is demonstrated in applications using MBE(n) up to order n = 8, making use of a new software application in which MBE(n) is tightly integrated with the semi-empirical extended tight-binding (xTB) model. 67,68This replaces the polarizable force field that was used previously as the low-level screening method, and represents a much more flexible and general means to screen the low-order subsystems.

II. METHODS
For a given subsystem IJK ⋅ ⋅ ⋅, energy corrections ΔE low IJK⋅⋅⋅ are computed using a low-level screening method.Here, we use the semi-empirical GFN2-xTB model 67,68 for that purpose.If then the corresponding n-body correction ΔEIJK⋅ ⋅ ⋅ is excluded in Eq. ( 1), and the parameter τn serves as a controllable threshold for screening the n-body interactions.In previous work, 41 we used the effective fragment potential 69 as the low-level screening method, but we have transitioned to GFN2-xTB due to its greater input flexibility and user-friendly software interface.

A. Extensible, generalizable many-body expansion
The generalized (G)MBE of Richard and Herbert [1][2][3][4] extends the simple MBE of Eq. ( 1) to handle overlapping fragments, meaning those that may share some atoms in common.Many other fragmentation schemes can be viewed as special cases of (or approximations to) the GMBE, [1][2][3] recovering the conventional MBE in Eq. ( 1) when the fragments are disjoint.The fragments used in the present work will be disjoint (specifically, single H 2 O molecules in water clusters), yet the set-theoretical formulation of the GMBE is still useful.
Here, we present a variation on the GMBE that can be updated iteratively to include new fragments.We regard each fragment FA as a set of atoms and a fragmentation scheme is composed of a unique set of fragments S = {FA, FB, . . ., FN}.
Starting with a scheme where x represents the scheme's current state in some sequence of additions, the corresponding approximation to the supersystem energy is where Ei is the energy of fragment Fi.Coefficients Ci,x arise from the principle of inclusion-exclusion (PIE) in state x, which prevents double-counting. 2,70,71When a new fragment Fα is added, the scheme Sx is updated according to The Journal of Chemical Physics

ARTICLE pubs.aip.org/aip/jcp
The corresponding energy update is where Here, Ci,x is the current coefficient for fragment Fi and Ei∩α is the energy of the fragment formed from the intersection Fi ∩ Fα.For brevity, we refer to fragments using their subscripts (A ∶= FA) and unions of fragments will be denoted using AB ∶= FA ∪ FB.

B. Fragment taxonomy
A fragmentation scheme consists of three types of fragments: (1) primary fragments are small and user-defined, and may overlap; (2) primitive fragments are the largest non-overlapping fragments derived from the intersections of primary fragments; and (3) auxiliary fragments include the primary fragments and all unions thereof.This taxonomy is illustrated in Fig. 1.
Relationships between fragments can be represented using a directed acyclic graph (DAG), wherein each unique fragment is a node.Edges are drawn from fragments subsystems (parents) to larger fragments constructed from them (children).We refer to this structure as a PIE-DAG, although it is also known in set theory as a Hasse diagram. 72,73PIE-DAGs are constructed using the following heuristics: 1.Each node represents a fragment (a set of atoms).

FIG. 2. PIE-DAG for system with disjoint primary fragments {A, B, C}, meaning that
2. A child node is the union of its parents.3. A node cannot have direct parents that are subsets of one another.
An example of a PIE-DAG for system with three primary fragments is illustrated in Fig. 2. Note that A ∩ B ∩ C = Ø in this example.
In practice, a fragmentation scheme can be bootstrapped starting from a single auxiliary fragment (S 0 = {FA}), followed by sequential additions as in Eq. (11).A screening procedure can be implemented to manage this process.In the present work, the PIE-DAG data structure is used to explore two possible strategies for screening auxiliary fragments, which are described in Sec.IV.

C. Test systems
In Sec.IV, we will evaluate the efficacy of two different energy-based screening algorithms.For demonstration purposes, the screening criterion in Eq. ( 6) will only be used for n ≥ 3, with all one-and two-body terms evaluated explicitly.(This also allows us to test screening approximations for the higher-order terms in a manner that is uncontaminated by any approximations introduced at the pairwise level.)As test systems, we will use several sets of (H 2 O)N clusters with N = 20, 32, and 64, some of which were taken from previous literature. 74,75As in previous work, 41,64 an ONIOM-style 1,76 supersystem correction is used to correct for long-range electrostatic and polarization interactions.
Calculations on (H 2 O) 20 and (H 2 O) 32 clusters were performed at the ωB97X-V/def2-TZVPD level with a supersystem correction at the level of Hartree-Fock (HF)/def2-TZVPD.This combination of methods is not a realistic use case for MBE(n) approximations, since both the fragment method and the supersystem correction exhibit the same formal scaling with respect to system size, but this choice facilitates high throughput and is useful for evaluating different algorithms and thresholds.For a more realistic application of fragmentation, we performed calculations on (H 2 O) 64 at the level of second-order Møller-Plesset perturbation theory (MP2), using the jun-cc-pVDZ basis set and a resolution of identity (RI) approximation, in conjunction with a supersystem correction at the HF/jun-cc-pVDZ level.

III. SOFTWARE
The extensible GMBE algorithm described in Sec.II has been implemented in Fragme∩t, a new Python-based application designed for fragmentation calculations. 772][43]59,60,78,79 The new Python application has already been used to obtain converged thermochemical quantities for enzyme-catalyzed reactions, 64 though only non-covalent fragmentation is considered in the present work.
An overview of the fragmentation process is shown in Fig. 3.The Fragme∩t code loads molecular systems, generates primary fragments, and then iteratively generates and screens auxiliary fragments to build the PIE-DAG that was described in Sec.II A. Fragme∩t also includes code capable interfacing with Q-Chem, 80 PySCF, 81 NWChem, 82 ORCA, 83 GFN-xTB, 67,68 and MOPAC. 84mong these, codes that rely on text input and output files are called as subprocesses and the output files are parsed, whereas codes such as GFN-xTB and PySCF that can be called as libraries are interfaced directly, to bypass overhead from system calls.Parallel calculations are supported through a prototype HTTP-based architecture capable of running ∼60 simultaneous calculations.Calculations reported here were performed by interfacing to Q-Chem for the density functional theory and RIMP2 calculations, and to GFN2-xTB for energy screening.
In developing Fragme∩t, our goal is to provide a framework for rapid implementation and testing of new and existing fragmentation methods.In part, we hope this will facilitate community validation of fragmentation methods, by making it relatively easy to perform side-by-side comparisons of different fragmentation methods.(Few such comparisons exist in the literature, 2,3,43,71,85 and several of those that do were performed using Fragme∩t's precursor C++ code. 2,3,43) To this end, Fragme∩t allows plugins to customize each step of a fragment-based calculation.Python was selected in order to lower the barrier to entry for new developers.Data are saved to a SQLite database, which makes Fragme∩t easy to deploy and allows data to be exchanged between a high-performance computing environment and a local installation, for ease of analysis.The code is presently under rapid development but a pre-release version is available on GitLab. 77

IV. ALGORITHMS AND TESTS
We next describe implementation and testing of generalized energy-based screening.Screening has directionality within the PIE-DAG paradigm, and we will describe both a top-down and a bottom-up algorithm to implement energy-based screening.Tests presented in this section will demonstrate the superior efficacy of the bottom-up approach.

A. Top-down screening
Top-down screening starts with the high-order nodes on the PIE-DAG and progressively extends this to lower-order nodes.This is illustrated at the MBE(3) level in Fig. 4(a), for an example consisting of non-overlapping primary fragments S primary = {A, B, C, D}.The screening procedure is first applied to the trimers ABC, ABD, ACD, and BCD.Terms that meet the screening criteria serve as initial nodes added to the PIE-DAG, and in this example we assume (for illustrative purposes) that only ABC, ACD, and BCD meet the screening criteria, not ABD.However, the fact that ABD is negligible does not guarantee that any of AB, AD, or BD is negligible, so each of these (and all other dimers) need to be tested.(Again for the sake of illustration, we assume in Fig. 4 that AD and BD are negligible but that AB is not.)Lastly, all one-body (primary) terms are added to the graph, ensuring that it contains all atoms.
As mentioned in Sec.II C, we do not actually screen the twobody interactions in this work, and instead evaluate all of the one-FIG.4. Schematic view of PIE-DAG fragmentation algorithms based nonoverlapping primary fragments S primary = {A, B, C, D}, traversed in the direction of the yellow arrows at left.(a) Top-down fragmentation method where all high-order fragments are screened first, followed by successively lower-order fragments.For the sake of example, fragments in red are assumed to be excluded based on screening criteria.Black arrows indicate parent/child relationships.(b) Bottom-up fragmentation where low-order terms are screened and added to the PIE-DAG first (again with fragments in red excluded due to screening), followed by successively higher-order terms.If a high-order term does not have all of its parents, it is not considered for addition, as indicated in gray.(c) Bottom-up screening using the parentage parameter M, with M = 1 in this illustration.A high-order term may be considered for addition if it is missing M or fewer parents.4(a) is intended as a schematic to illustrate that within the top-down approach, any sparsity that may be obtained via screening at the n-body level does not propagate downward, meaning that all of the (n − 1)-body terms must be screened.This is true for higher-order interactions as well, beyond the n = 3 case that is illustrated in Fig. 4(a).Figure 5 illustrates the results of top-down energy screening applied to a test set of (H 2 O) 20 clusters at the MBE(3) and MBE(4) levels.The top-down method requires a screening procedure at each n-body level, for which we use GFN2-xTB with a three-body threshold τ 3 ranging from τ 3 = 0 (meaning no screening) up to 0.15 kcal/mol.For simplicity, the four-body threshold (τ 4 ) is set to τ 4 = ατ 3 , for α = 0.05, 0.10, or 0.20.Larger values of τn indicate more aggressive screening, resulting in fewer subsystems.Throughout this work, error is defined as where Esupersys is a conventional supramolecular calculation at the "target" level of theory that is used for the MBE(n) terms.From Fig. 5(a) we find that MBE(3) maintains sub-kcal/mol accuracy provided that τ 3 ≤ 0.1 kcal/mol, whereas for larger values there is a rapid increase in the error as the screening becomes too aggressive.This loss of accuracy for τ 3 > 0.1 kcal/mol is substantially mitigated by moving to a four-body expansion.This is true for all values of τ 4 that we test, although the most restrictive value (τ 4 = 0.2τ 3 ) does exhibit absolute errors exceeding 1 kcal/mol when τ 3 > 0.1 kcal/mol Setting τ 4 = 0.1τ 3 maintains the absolute error within 1 kcal/mol at the expense of increased fragment counts.The latter settings retain most of the trimers, which are generated from the four-body overlaps, and this masks the influence of τ 3 .For top-down screening to be effective, both τ 3 and τ 4 must be tuned to strike a balance between accuracy and efficiency, where efficiency is measured by the fragment counts that appear in the lower panels of Fig. 5.

B. Bottom-up screening
PIE-DAGs can also be constructed upwards from the lowestorder fragments, as illustrated in Fig. 4(b), eliminating higher-order nodes based on whether low-order ones are significant.Somewhat similar procedures have sometimes been used to eliminate higherorder contributions to potential surfaces, 86 which can be couched in a MBE-type formalism, 87 although in that context the expansion is seldom extended to high order.Here, we will develop a general algorithm and show that it can be carried to arbitrary (self-terminating) order in realistic test cases.
The bottom-up process begins by including all primary fragments in the graph, whereupon dimers are added if they meet the screening criteria.In the hypothetical example of Fig. 4(b), dimers AD and BD are determined to be negligible, as they were in the corresponding top-down example of Fig. 4(a).Trimers and higherorder fragments are added subsequently, with the constraint that in order for a new node (representing a non-negligible subsystem) to be added, all of its parents must exist already within the PIE-DAG, meaning that each parent must be non-negligible according to the screening criteria.Under these conditions, the nodes ABD, ACD, and BCD are not considered for addition to the scheme in Fig. 4(b), meaning that they are not even tested with the low-level screening method.The node ABC (depicted in yellow) is considered, and ΔE low ABC is evaluated; this trimer will be included in the scheme if and only if it satisfies the three-body screening criteria, meaning Eq. ( 6) in the present implementation.
Inclusion of a particular trimer IJK depends upon the presence of its parents: IJ, IK, and JK.In the algorithm that we have just described, which is illustrated in Fig. 4(b), each of these parents must be present in order for IJK to be included in the PIE-DAG.We call this the M = 0 algorithm because no missing parents are  4) with M = 0, 1, or 2. (The screening algorithm is fully specified by τ 3 and M so there is no τ 4 in the bottom-up approach.)Lower panels (e)-(h) show the mean fragment counts in the form of stacked plots.
permitted.Relaxing this constraint, to allow a single missing parent, affords what we call the M = 1 algorithm and this is illustrated in Fig. 4(c).Relative to the M = 0 algorithm, use of M = 1 means considering ACD and BCD, which are thus evaluated using the lowlevel screening method and included in the PIE-DAG if they meet the screening criteria.Use of M = 1 may provide better accuracy in cases where one fragment significantly polarizes two others, such as an ion located between two water molecules, where the corresponding water dimer might not have been strongly interacting (due to the separation between the water molecules) when considered in the absence of the intervening ion.
Figure 6 shows the bottom-up algorithm applied to the same set of (H 2 O) 20 clusters that were used to test the top-down method.Energy screening is applied only to the three-body terms, retaining FIG. 7. Bottom-up energy screening on a collection of 10 (H 2 O) 32 clusters, at the ωB97X-V/def2-TZVPD level with a HF/def2-TZVPD correction.As in Figs. 5 and 6, the traces in the upper panel show the error for each cluster isomer as a function of τ 3 , for (a) MBE( 3 all one-and two-body terms (as in the top-down example).In this case, four-body terms are accepted or rejected based strictly on parentage (M = 0, 1, or 2), rather than using a four-body threshold τ 4 .
For all three values of M, the number of fragments decreases as the screening on the lower-level fragments becomes stricter (larger τ 3 ).A sharp drop in accuracy for MBE(3) is again observed for τ 3 > 0.1 kcal/mol [Fig.6(a)], which remains evident for MBE(4) using either M = 0 [Fig.6(b)] or M = 1 [Fig.6(c)].This indicates that addition of four-body terms does not obfuscate the screening that is performed at the three-body level.However, the choice M = 2 [Fig.6(d)] does appear to interfere with the three-body screening, leading to a significant growth in the number of subsystems [Fig.6(h)].For M = 2, errors do remain tightly clustered even for τ 3 > 0.1 kcal/mol, so in some sense this approach is able to inoculate the method against a poor choice of τ 3 , yet that security comes at significant cost.Setting M = 1 and τ 3 = 0.05 kcal/mol provides a good balance between accuracy and efficiency.
To evaluate the efficacy of bottom-up energy screening using larger systems, a similar analysis is presented for (H 2 O) 32 clusters in Fig. 7 and for (H 2 O) 64 clusters in Fig. 8.In order to make a direct comparison to results presented above, we use the same levels of theory for (H 2 O) 32 as for the (H 2 O) 20 clusters.However, calculations for (H 2 O) 64 are performed at a level of theory that constitutes a more realistic use case for MBE(n), with subsystem calculations evaluated at the RIMP2/jun-cc-pVDZ level of theory in conjunction with a HF/jun-cc-pVDZ supersystem correction.Error is then measured relative to a supersystem calculation at the RIMP2/jun-cc-pVDZ level.As demonstrated in previous energy-screened MBE(3)   calculations, 41 this results in dramatic speedups without large-scale parallelization, despite the need for a supersystem HF calculation, while maintaining 1 kcal/mol fidelity for relative energies.These (H 2 O) 64 calculations represent the largest application of MBE(4) in electronic structure theory of which we are aware.
As cluster size increases a systematic error emerges, in both the (H 2 O) 32 and (H 2 O) 64 calculations, that is larger than the (relatively constant) error in amongst relative energies of different clusters of the same size.To demonstrate that this systematic error is due to BSSE, we performed MBE(3) and MBE(4) calculations for a single (H 2 O) 64 isomer using the full cluster basis set, meaning that n-body subsystem calculations are carried out using the basis set of the (H 2 O) 64 supersystem.Errors in these cluster-basis calculations are plotted in Fig. 8(c) alongside the corresponding errors obtained using the conventional subsystem basis.Use of the cluster basis reduces the error below 0.025 kcal/mol/monomer for both MBE(3) and MBE(4), whereas the use of the subsystem basis afford errors that grow with τ 3 , even when that threshold is set quite conservatively.Although the use of a supersystem basis set may not be practical in large-scale calculations, in future work this can replaced with many-body counterpoise corrections. 59

V. DISCUSSION
Both top-down and bottom-up screening can achieve good fidelity with respect to the corresponding supersystem calculation at the same level of theory.However, top-down fragmentation has a significant drawback in that it requires a well-defined screening protocol for each n-body order, and lower orders in the expansion do The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcpnot benefit directly from the screening performed at higher orders when it comes to constructing the graph of non-negligible subsystem nodes.Even for a screening protocol as simple as Eq. ( 6), where n-body terms are excluded based on an energy threshold τn, the topdown approach is impractical for large systems when n > 3, because every high-order n-body subsystem must be considered individually.
To put this into perspective, and to contrast the efficiency of the bottom-up approach, Fig. 9 presents average timing data to construct the PIE-DAG for (H 2 O) 20 , for the calculations reported in Figs. 5 and 6.Cached GFN2-xTB screening data were used to eliminate variations in run times, such that the only cost associated with energy screening is that of a hash lookup for the appropriate energy.When τ 3 is small, both algorithms exhibit a relatively large fragmentation cost because the PIE-DAG is quite dense; however, top-down screening has a high fixed cost even when the screening is more aggressive (larger τ 3 ).Bottom-up screening is more efficient even when τ 3 is small because it does not attempt to screen fourth-order and higher fragments whose parentage is incomplete.The present implementation can perform a single ΔE low IJKL calculation (using GFN2-xTB) in 0.02 s, and for top-down screening this equates to 4 h just to generate and screen the 679 120 unique subsystems required to apply MBE(4) to (H 2 O) 64 .Bottom-up screening can fragment the same system in less than 5 min.
This has significant implications for counterpoise correction schemes that can be used to address BSSE but which have historically been underutilized in fragment-based quantum chemistry, perhaps due to cost.The data in Fig. 8(c) demonstrate that error increases as the number of high-order subsystems is reduced by more aggressive screening, which suggests that most of the high-order terms in MBE(n) contribute nothing more than incomplete error cancellation.Counterpoise corrections for MBE(n) and its generalizations have been developed 40,59,61-63,88-90 but are not yet implemented in our code.These methods introduce an additional combinatorial scaling factor but aggressive energy screening may render these corrections tractable.

FIG. 9.
Mean fragmentation times for (H 2 O) 20 clusters using bottom-up energy screening (in blue) vs top-down screening (in orange).In both cases, cached data from GFN2-xTB were used for screening, so these timings simply reflect the lookup cost and the cost to construct the PIE-DAG.The top-down method requires data for every tetramer whereas the bottom-up approach does not.The steep cost increase for small values of τ 3 results from construction of a dense PIE-DAG.were performed on laptop with a 2.6 GHz dual-core Intel Core i5 processor and Gb of memory.
Because the bottom-up algorithm propagates the screening procedure from low orders to high orders, it allows the energyscreened MBE(n) to be extended to arbitrarily high orders.Note that each n-body node in the PIE-DAG has n parents, each of which is a (n − 1)-body subsystem.Since higher-order subsystems have more parents than lower-order ones, fewer terms need to be considered as n increases, even prior to evaluating ΔE low IJKL⋅⋅⋅ .This natural attenuation is demonstrated in Fig. 10 using (H 2 O) 20 clusters.For each cluster, we ran the bottom-up algorithm up to n = 20, with τ 3 = 0.05 kcal/mol and M = 1.For 10 out of the 12 clusters in the data set, we found that the procedure terminated naturally (due to a lack of parents) at n = 7, and for the remaining two clusters it terminated at n = 8.
Figure 10(b) shows the number of unique subsystems generated at each n-body order, excluding terms where Ci,x = 0 in Eq. ( 12).An interesting feature of these data is that the number of three-, four-, and five-body terms actually decreases as n increases.In essence, the combination of higher-order n-body terms with energy screening serves to localize inter-fragment interactions into larger domains (characterized by larger n), which are capable of accounting for all inter-domain interactions.The most dramatic examples of this effect are the fused-cube isomers No. 1 and No. 20 (which are provided in the supplementary material), for which the addition of a single eight-body fragment to a MBE (7) scheme reduces the number of subsystems by ≈200 for either cluster.
It is commonly asserted that high-order n-body terms are energetically insignificant.In particular, for homogeneous water clusters it has been assumed that the n-body corrections are negligible beyond n = 4, 39,40,42,[56][57][58] or possibly n = 5. 29,48 (Higher-order terms undoubtedly persist in ion-water clusters. 53,54) A careful analysis of this assertion in large (N ≫ 20) water clusters has been hindered by intractable combinatorics, but should now be feasible and we hope to report on this in due course.For now, we note that that terms beyond n = 5 in (H 2 O) 20 clusters converge to a constant, nonzero error [Fig.10(a)].We ascribe this constant offset to the combined The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcpeffects of BSSE and perhaps some residual error due to the supersystem correction.This analysis suggests that for N = 20 monomers, MBE(n) does indeed converge at the n = 5 level at worst, and perhaps sooner, if higher-quality basis sets were to be employed.Others have suggested that five-body terms in neat water clusters are BSSE artifacts that disappear in the complete-basis limit. 56The present results certainly do not preclude that possibility, which we plan to explore in future work.

VI. CONCLUSIONS
We have described a new, open-source implementation of the straightforward MBE(n) method without any kind of electrostatic embedding.The asceticism is by design, as we aim to circumvent complexities arising from self-consistent charge embedding including overly complicated analytic gradients, 66 instead using high-order n-body expansions and low-level supersystem corrections to capture cooperative polarization.Numerical problems arising from unfavorable combinatorics for n > 3, which we have previously documented in large-scale MBE(n) calculations, 5,43 are avoided here by means of energy screening.Stability in the presence of diffuse basis functions is achieved by iterating the n-body self-consistent field calculations to convergence.
Using a new "bottom-up" screening algorithm based on the semi-empirical GFN2-xTB method, 67,68 we are able to achieve subkcal/mol accuracy for water clusters when MBE(n) is paired with an ONIOM-style supersystem correction.This constitutes a predictive level of accuracy by means of which MBE(n) calculations can reliably rank-order the energies of different hydrogen-bonding networks, which has proven to be a challenging problem for methods that use single-H 2 O fragments. 1 Numerically-complete four-body results have been presented for (H 2 O) 64 , which represents a reasonable size for the periodic cell in a liquid water simulation.We believe these to be the largest complete application of MBE(4) to date.
Unlike the top-down screening method that we reported previously, 41 which requires energy screening for every individual n-body system, the bottom-up approach propagates screening information upwards.As a result, the screening is much more aggressive (and thus the algorithm more performant) without loss of accuracy, and we find that the number of n-body subsystems naturally attenuates as n increases.As such, MBE(n) can be used without specifying a maximum order n, yet subsystem counts remain manageable.For (H 2 O) 20 , terms beyond n = 5 appear to contribute only to BSSE.In future work, this will be investigated further using larger basis sets and many-body counterpoise corrections. 59he software reported here 77 makes MBE(n) into a robust tool for high-fidelity fragment-based quantum chemistry.Our opensource implementation is easy to interface with various quantum chemistry programs and can be used to prototype new fragmentation strategies.We anticipate that this implementation will lower the entry barrier to using fragmentation methods, and furthermore provide a degree of community validation that has heretofore been lacking in fragment-based quantum chemistry.

SUPPLEMENTARY MATERIAL
Coordinates for all test systems used here.

FIG. 5 .FIG. 6 .
FIG. 5. Top-down energy screening applied to a set of 12 different (H 2 O) 20 clusters at the ωB97X-V/def2-TZVPD level using a HF/def2-TZVPD supersystem correction.Traces in the upper panels show the error for each individual cluster as a function of τ 3 , with the mean error in bold, for (a) MBE(3) and (b)-(d) MBE(4) with various choices for τ 4 .The blue shaded region delineates one standard deviation around the mean and the gray shaded region indicates the target error of ±1 kcal/mol.Lower panels (e)-(h) show mean fragment counts at each order in the n-body expansion, in the form of stacked plots.
FIG. 7.Bottom-up energy screening on a collection of 10 (H 2 O) 32 clusters, at the ωB97X-V/def2-TZVPD level with a HF/def2-TZVPD correction.As in Figs.5 and 6, the traces in the upper panel show the error for each cluster isomer as a function of τ 3 , for (a) MBE(3) vs (b)-(c) MBE(4) with either M = 0 or 1. Lower panels (d)-(f) show the mean fragment counts.

FIG. 8 .
FIG. 8. Bottom-up energy screening for a set of 11 (H 2 O) 64 clusters, at the RIMP2/jun-cc-pVDZ level with a HF/jun-cc-pVDZ supersystem correction.Traces in the upper panel show the errors for (a) and (b) MBE(4) with M 1, both as a function of τ 3 and reported in per-monomer units.(c) Errors for both MBE(3) and MBE(4) for one particular snapshot, in either the usual subsystem basis set (in blue) vs the full cluster basis set (in orange).(d)-(e) Mean fragment counts as a function of τ 3 .

FIG. 10 .
FIG. 10.Bottom-up MBE(n) calculations using the (H 2 O) 20 test set, setting τ 3 = 0.05 kcal/mol and M = 1.(a) Errors for all 12 cluster isomers as a function of n.(b) Mean fragment counts for each subsystem size.(c) Close-up view of the mean fragment counts for n = 7 and 8.No additional fragments are generated for n ≥ 9.
Fragmentation workflow as implemented in Fragme∩t.Once the supersystem's geometry has been loaded and split into primary fragments, it is passed to Fragme∩t's main loop, which iteratively proposes new auxiliary fragments.These may or may not be added to the PIE-DAG, depending on the outcome of the screening procedure.New fragments and PIE metadata are stored in a SQLite database.Electronic structure calculations are run in parallel and are also stored in the database.

org/aip/jcp and
two-body terms at the target level of theory.(This means one fewer parameter to test in the present work, and two-body screening can easily be incorporated later.)The graph in Fig.