MB-PIP: Many-Body Permutational Invariant Polynomials

• MB-PIP is a framework that decomposes total energy into explicit many-body contributions using permutationally invariant polynomials in transformed intersite distances.
• It enforces correct asymptotic behavior by purifying the polynomial basis, ensuring each many-body term vanishes appropriately in dissociation limits.
• The methodology enables efficient gradient evaluations and scalable fitting strategies, making it applicable to water clusters, linear alkanes, and other molecular systems.

Many-body permutationally invariant polynomial (MB-PIP) denotes a class of potential-energy representations in which a system is decomposed into explicit low-order many-body contributions and selected terms are written as permutationally invariant polynomials in transformed intersite distances. In the canonical water formulation associated with MB-pol, the PIPs are the short-range 2-body and 3-body corrections within a many-body expansion that otherwise contains a Partridge–Schwenke 1-body term, physics-based long-range electrostatics, induction, dispersion, and higher-body polarization (Nguyen et al., 2018). Subsequent work extended the same symmetry-adapted polynomial logic to an explicit 4-body water interaction, to atom-based expansions for linear alkanes, and to monomer-centered PIP-descriptor hybrids, so MB-PIP now refers not to a single architecture but to a family of closely related constructions that share body ordering, algebraic permutation symmetry, and explicit polynomial bases (Nandi et al., 2021, Qu et al., 22 Sep 2025, Yu et al., 2024).

1. Formal structure and scope

In the strictest usage, MB-PIP is tied to the many-body expansion (MBE),

$$E_N(1,\dots,N)=\sum_i V^{\mathrm{1B}(i)}+\sum_{i<j}V^{\mathrm{2B}(i,j)}+\sum_{i<j<k}V^{\mathrm{3B}(i,j,k)}+\cdots ,$$

with higher-order terms defined recursively by subtracting all lower-order contributions from the energy of each $n$-mer subsystem. In this formulation, a fitted PIP does not ordinarily learn total energy directly; it learns a specific many-body component, or a short-range residual of that component after a physics-based baseline has been removed (Nguyen et al., 2018).

That definition is exact for monomer-based water models, but the literature also contains broader body-ordered interpretations. In linear alkanes, where there is “no obvious monomer partitioning,” the expansion is written over atoms, so the “monomers are atoms,” all 2-body, 3-body, and 4-body atomic tuples are enumerated, and the total energy is the sum of those contributions (Qu et al., 22 Sep 2025). In monomer-centered MB-PIPNet, by contrast, the total energy is written as a sum of effective monomeric energies predicted from PIP descriptors; this retains a many-body viewpoint but is not an explicit MBE over fitted $n$-body energy surfaces (Yu et al., 2024, Li et al., 15 Feb 2026).

Formulation	Decomposition unit	Representative use
Explicit MBE MB-PIP	Molecules or monomers	Water 2B/3B MB-pol; explicit 4B water
Atom-based MB-PIP	Atoms	Linear alkane 2B/3B/4B expansion
Monomer-centered PIP hybrid	Monomers or fragments	MB-PIPNet for water, CO$_2$, alkanes

A persistent conceptual distinction follows from this table. In explicit MB-PIP, the fitted object is an $n$-body energy term with its own dissociation constraints. In atom-based or monomer-centered descendants, the same PIP machinery is retained, but the definition of “body” and the role of the polynomial basis change.

2. Polynomial construction, symmetry, and asymptotics

The core MB-PIP ansatz is an algebraic expansion

$V_{\mathrm{ML,PIP}}=\sum_l c_l \eta_l,$

where the $\eta_l$ are symmetrized monomials in transformed geometric variables. In water MB-pol, these variables are functions of pair distances involving real atoms and auxiliary lone-pair sites $L_1$ and $L_2$, with elementary forms $\xi_i=e^{-k d_i}$, $n$0, and $n$1 (Nguyen et al., 2018). The resulting polynomial is invariant under permutation of water molecules, permutation of equivalent H atoms within each molecule, and permutation of equivalent lone-pair sites within each molecule.

The same principle reappears in other MB-PIP realizations, but with different local symmetry classes. The atom-based alkane model assigns separate PIP bases to tuple types such as CCC, CCH, CHH, and HHH, with permutation symmetries $n$2 or $n$3 at 3-body order and $n$4, $n$5, or $n$6 at 4-body order (Qu et al., 22 Sep 2025). The important point is that permutation invariance is imposed at basis-construction time, not learned statistically.

Correct asymptotic behavior is equally central. In the water 4-body PES, the basis is purified so that the genuine 4-body interaction vanishes rigorously in all monomer+trimer and dimer+dimer dissociation channels, exactly as required by the recursive definition of $n$7 (Nandi et al., 2021). The 2021 PIP methodology paper generalizes this idea algorithmically: basis purification sorts generated polynomials into those with correct large-distance limits and those without, and compaction then deletes monomials and polynomials that are no longer needed after purification (Houston et al., 2021). The alkane MB-PIP imposes the same physical principle in a simpler statement: any polynomial that does not go to zero when an atom is infinitely far away from the remaining atoms is removed from the basis (Qu et al., 22 Sep 2025).

These asymptotic constraints are not decorative. They define whether a fitted term is genuinely many-body rather than merely permutation-invariant. A polynomial with the wrong dissociation limit can still interpolate training data well, but it no longer represents a valid $n$8-body contribution in the sense of the MBE.

3. Fitting strategies, gradients, and computational realization

Once the nonlinear variables are fixed, MB-PIP fitting is typically a linear regression problem in the coefficients $n$9. In water PIP-MB-pol, the 2-body and 3-body short-range residuals are fitted by ridge-regularized least squares,

$n$0

with $n$1 for 2B and $n$2 for 3B; linear parameters are optimized by singular value decomposition and nonlinear parameters by the simplex algorithm. The reference datasets contain 42,069 dimer configurations and 12,347 trimer configurations, randomly split into training, validation, and test subsets in the ratio $n$3 (Nguyen et al., 2018).

Not all MB-PIP models fit isolated many-body components. The transferable alkane model is trained all at once against total molecular energies with a standard over-determined linear least-squares fit, using 734 undetermined linear coefficients for the low-level 2b/3b/4b model and 383 coefficients for the $n$4-ML correction PES (Qu et al., 22 Sep 2025). This is a substantive architectural difference from water MB-pol-style fitting, where the PIPs are trained on short-range residuals of individual body terms.

Analytical gradients are a major practical advantage of PIPs. The reverse-differentiation framework evaluates the energy as

$n$5

and computes all derivatives with respect to transformed internal variables in one backward sweep through the computational graph, followed by Cartesian chain-rule contraction. On ethanol benchmarks, this yielded energy and force evaluations that were reported as $n$6 to $n$7 times faster than several alternative ML methods while maintaining comparable precision (Houston et al., 2021).

This computational pattern explains a recurrent MB-PIP tradeoff. Training can require careful basis generation, purification, and nonlinear parameter selection, but inference is then a fixed algebraic evaluation with exact gradients. That separation between expensive model construction and cheap deterministic deployment is one of the defining practical features of the method family.

4. Canonical explicit MB-PIP realizations

The best-known explicit MB-PIP is the short-range part of MB-pol for water. The total model contains explicit terms through 3-body order and treats all higher-body contributions via classical polarization. The 2-body term is decomposed into

$n$8

where the long-range baseline contains permanent electrostatics, induction, and damped dispersion; the 3-body term is written as

$n$9

The PIPs enter only through the switched short-range corrections, so the polynomial is explicitly confined to the overlap-driven region associated with “charge transfer and penetration, and Pauli repulsion” (Nguyen et al., 2018).

For the 2-body correction, the final PIP uses 31 primitive variables, 1153 symmetrized monomials, and maximum degree 4. For the 3-body correction, it uses 36 primitive exponential variables, 1163 symmetrized monomials, and maximum degree 4. The switching cutoffs are $_2$0 Å and $_2$1 Å for dimers, and $_2$2 Å and $_2$3 Å for trimers. Test RMSEs for the PIP model are $_2$4 kcal/mol per dimer for 2B and $_2$5 kcal/mol per trimer for 3B. For water clusters $_2$6, $_2$7–6, errors in 2B and 3B contributions are below $_2$8 kcal/mol in all cases, total interaction-energy errors never exceed $_2$9 kcal/mol, and the relative ordering of isomers is preserved (Nguyen et al., 2018).

An important result of that study is that the many-body scaffold, not the specific regressor, is the dominant structural ingredient. PIP-MB-pol, BPNN-MB-pol, and GAP-MB-pol all show very similar accuracy on the same 2B/3B water problem, and the largest residual errors arise in sparsely sampled highly repulsive regions rather than from any decisive weakness of the PIP ansatz itself (Nguyen et al., 2018).

The explicit 4-body water PES extends this logic to a genuine tetramer interaction. It is a 12-atom, maximum-order-3 PIP built from all 66 internuclear distances, using Morse-type variables for the 12 intramolecular coordinates and inverse-distance variables for the 54 intermolecular coordinates. The fit is based on 2119 symmetry-unique CCSD(T)-F12a/haTZ 4-body energies. Under full 84 symmetry, purification leaves only 2 basis functions, so the final model uses reduced $n$0 symmetry plus replication over all $n$1 monomer permutations; this produces 10,737 candidate PIPs before purification and 1649 retained basis functions afterward. The overall RMS fitting error is $n$2, and the resulting 4-body energies of water-hexamer isomers are in better agreement with benchmark CCSD(T) results than those from MB-pol (Nandi et al., 2021).

The 4-body study also shows that explicit higher-body MB-PIP terms can be computationally practical when their asymptotics are correct. A single-core implementation evaluates 100,000 PES calls in $n$3 s, and in a 256-water cluster a simple screen based on the largest O–O distance below $n$4 Å discards $n$5 of tetramers, reducing the estimated 4-body evaluation time from about 1 hour to 14 seconds (Nandi et al., 2021).

5. Extensions to covalent molecules, local body-order models, and PIP-descriptor hybrids

The atom-based MB-PIP for linear alkanes demonstrates that explicit many-body PIPs are not limited to noncovalent monomer expansions. Here the total energy is written as

$n$6

with $n$7 the number of atoms, truncation at 4-body, and physical range cutoffs to control the combinatorial growth of tuples. The model uses Morse variables $n$8, maximum polynomial powers 10, 8, and 6 for 2-body, 3-body, and 4-body terms, respectively, and a total of 734 fitted coefficients. Trained on $n$9, it was previously reported transferable from $V_{\mathrm{ML,PIP}}=\sum_l c_l \eta_l,$0 to $V_{\mathrm{ML,PIP}}=\sum_l c_l \eta_l,$1, and in the $V_{\mathrm{ML,PIP}}=\sum_l c_l \eta_l,$2-ML extension the same framework is used for conformational energy differences from $V_{\mathrm{ML,PIP}}=\sum_l c_l \eta_l,$3 to $V_{\mathrm{ML,PIP}}=\sum_l c_l \eta_l,$4. The $V_{\mathrm{ML,PIP}}=\sum_l c_l \eta_l,$5-correction itself is another MB-PIP with 383 coefficients and fit errors of $V_{\mathrm{ML,PIP}}=\sum_l c_l \eta_l,$6 for $V_{\mathrm{ML,PIP}}=\sum_l c_l \eta_l,$7-B3LYP and $V_{\mathrm{ML,PIP}}=\sum_l c_l \eta_l,$8 for $V_{\mathrm{ML,PIP}}=\sum_l c_l \eta_l,$9-PBE0+MBD; over the tested chain-length range, the corrected B3LYP PES yields a mean absolute error of $\eta_l$0 kcal/mol for the hairpin-minus-linear energy difference (Qu et al., 22 Sep 2025).

Closely related local body-ordered variants use the same polynomial-invariant philosophy but are not strict fragment MBEs. The molecular aPIP force-field framework extends atomic PIPs to multi-element molecules with generalized 2-body, 3-body, and 4-body channels, distance-angle coordinates, smooth cutoffs, QR-based linear fitting, and explicit regularization. In a combined alkane fit trained through hexane, the model retained less than $\eta_l$1 meV/atom error at 300 K for unseen heptane and octane and high-temperature errors of $\eta_l$2 and $\eta_l$3 meV/atom, respectively (Allen et al., 2020). The materials-oriented aPIP framework makes the same point in elemental condensed phases: explicit low-dimensional body-order terms, combined with Tikhonov and integral-function regularization, can improve transferability and remove catastrophic “holes,” and the authors report better transferability than SOAP-GAP in sparse-data regimes for W, Si, and Ti (Oord et al., 2019).

MB-PIPNet moves further from explicit MB-PIP while preserving its descriptor logic. In the 2024 monomer-centered formulation, the total energy is

$\eta_l$4

where each monomer energy depends on a self descriptor $\eta_l$5 and an environmental descriptor

$\eta_l$6

both built from 1-body and 2-body PIPs in Morse-type transformed distances. For the water trimer, this architecture achieved training and test RMSEs of $\eta_l$7 and $\eta_l$8 meV/atom; for a liquid-water model trained on MB-pol data, the reported training and test RMSEs are $\eta_l$9 and $L_1$0 meV/atom (Yu et al., 2024). The 2026 extension to covalent systems applies the same monomeric logic to $L_1$1 and $L_1$2 fragments in $L_1$3, using 445-dimensional PIP descriptors and $L_1$4 neural networks. On 100,000 geometries, the reported timings are 1792 s for DeepMD energies and gradients, 1248 s for MB-PES energies and gradients, and 240 s for MB-PIPNet energies and gradients, while MB-PIPNet remains more accurate than DeepMD and somewhat less accurate than explicit MB-PES on total energies (Li et al., 15 Feb 2026).

These extensions show that the PIP component of MB-PIP can serve three distinct roles: direct representation of explicit many-body terms, local body-ordered basis for transferable force fields, or compact chemically motivated descriptor for a neural network.

6. Accuracy, transferability, recurring misconceptions, and open problems

Several conclusions recur across the literature. First, the decisive factor is often the physically structured decomposition rather than the choice of regression engine alone. In water, PIP, BPNN, and GAP give very similar 2B/3B accuracy inside the same many-body scaffold, and the dominant residual errors are associated with nonuniform training-set coverage in sparse, highly repulsive regions (Nguyen et al., 2018). This directly contradicts the common simplification that MB-PIP performance can be understood by viewing the polynomial only as a standalone regressor.

Second, the literature does not support a blanket claim that PIPs are restricted to very small systems. Explicit 4-body water fitting, transferable 2b/3b/4b alkane models, regularized atomic body-order formulations, and monomer-centered PIP-descriptor networks all demonstrate workable scaling strategies. Those strategies, however, are not interchangeable: explicit high-order MB-PIPs rely on purification and exact asymptotics; local aPIP models rely on cutoffs and low-dimensional body terms; MB-PIPNet gives up the explicit MBE in favor of effective monomer energies (Nandi et al., 2021, Qu et al., 22 Sep 2025, Oord et al., 2019, Yu et al., 2024).

Third, computational efficiency depends strongly on basis engineering. Reverse differentiation, purification, compaction, and pruning have been shown to make PIP energies and forces much faster to evaluate than several alternative ML models while preserving high precision (Houston et al., 2021). At the same time, the raw basis size still grows rapidly with variable count and polynomial order, and extending explicit PIPs to larger or more heterogeneous systems remains substantially more difficult than deploying generic atom-centered or graph-based models (Nguyen et al., 2018, Li et al., 15 Feb 2026).

The current limitations are correspondingly specific. In the alkane $L_1$5-MB-PIP, the correction was fitted to direct electronic energy differences rather than to the discrepancy between the low-level MB-PIP PES and CC energies, and the authors identify this as a main reason the corrected result is not perfect; the residual conformational error remains about $L_1$6 kcal/mol (Qu et al., 22 Sep 2025). In MB-PIPNet, higher-body effects are only implicit, long-range interactions are not separately modeled by a dedicated physically motivated term, and the current alkane implementation modestly overestimates the C–H stretch region near $L_1$7 (Yu et al., 2024, Li et al., 15 Feb 2026). In explicit higher-body PIPs, the need to enforce all relevant dissociation channels becomes increasingly intricate as body order rises (Nandi et al., 2021).

A plausible synthesis is that MB-PIP is best regarded not as one model class but as a design principle: decompose the energy into chemically or physically meaningful pieces, encode the exact permutation symmetry algebraically, enforce the asymptotics appropriate to the chosen body term, and exploit the resulting explicit polynomial structure for stable fitting and exact differentiation. The literature suggests that this principle remains competitive across water clusters, flexible covalent molecules, and condensed-phase systems, but that its most successful realizations are those in which the many-body decomposition, the asymptotic constraints, and the data design are all treated as first-class components of the model rather than as secondary implementation details.