where 𝑆 is an operator that symmetrizes the monomials 𝑦 "# 9 𝑦 "& : 𝑦 #& ; (the symmetrized terms are invariant to all permutations among triplets of NPs) and 𝐶 9:; are linear fitting parameters representing the coefficients of the symmetrized terms with 𝑎 + 𝑏 + 𝑐 = 𝑚. Training Data. To generate the training data, and later validate the developed potential, we carried out MD simulations of polymer-grafted NPs in a polymer matrix treated using a CG model (Fig. ). In this model, the polymer grafts and the matrix polymer were treated as flexible bead-chains of lengths 𝐿 = = 10 and 𝐿 > = 20 beads; the beads represent short segments of polymer chains and are all of size 𝜎 and mass 𝑚. The grafts and matrix polymers are assumed to be fully miscible with each other. Therefore, the graft-matrix intersegment interactions, and the equally miscible matrix-matrix and graft-graft intersegment interactions, were all treated using a potential that accounts for both excluded volume and attractive interactions, with parameters 𝜎 and 𝜀 specifying the size and attraction strength of the segments. The NP cores were treated as rigid spheres of radius 𝑅 ?@ = 3𝜎 and mass of 216𝑚 that interact with each other using a combined attractive and excluded-volume potential of attraction strength 𝜀 ?@ = 30𝜀 (unless otherwise stated) and interact neutrally with all polymer segments using an excluded-volume potential. Bare and grafted NPs with grafting densities of Γ = = 0.15 to 0.4 chains/𝜎 ! spanning the mushroom to weak-brush grafting regimes were explored. All systems were investigated at a constant temperature of 𝑇 = 𝜀/𝑘 * and at a constant polymer segment density of 0.85 beads/𝜎 % producing melt-like conditions. Based on these grafting densities, graft-matrix miscibility, and relative length of matrix and grafted chains, the polymer-grafted NPs are expected to be in the "wetting regime", where the matrix chains can penetrate (wet) the polymer grafts. This CG model has been widely used for studying polymer-NP systems, as it captures the most essential physics of NP-polymer systems, including many-body effects, while keeping the computational costs reasonable. The 𝑊 ! and ∆𝑊 % training data was obtained from free energies (PMFs) of three polymergrafted NPs (labelled NP1, NP2, and NP3) in a polymer matrix computed as a function of interparticle distances using the blue moon ensemble method. These free energies were calculated relative to those of a reference system in which the three NPs are located far from each other where they exhibit no interactions. The 𝑊 ! (𝑑 "# ) potential is simply the separationdistance dependent PMF between a pair of interacting NPs that are isolated from other NPs. To obtain this PMF-𝑊 ! (𝑑 4! ) in our case as we use NP1 and NP2 to compute it-we held the core of NP1 fixed during the simulation while the core of NP2 was brought closer to NP1 in a stepwise manner along a straight-line path passing through its center (red dashed line, Fig. ). 𝑊 ! was obtained as a function of 𝑑 4! by integrating the ensemble-average normal force experienced by NP2 on a finely spaced grid along the path. To ensure that the PMF was not affected by NP3, this NP also was held fixed at a location far from both NP1 and NP2 during the simulation.