Encoding biomolecular systems requires a large, inconsistently-sized parameter space that is intractable to train and use when considered as a whole. ML approaches commonly reduce complex systems into their component parts ("tokens"), then infer properties of the system as a whole based on relationships between tokens. This approach is complementary to fragmentation methods in quantum chemistry, which approximate supersystem properties by systematically partitioning that system into numerous fragments, for which it is relatively inexpensive to perform high-quality calculations. This provides a hierarchy of well-defined fragments and a database can be used to train ML model for large systems. Generating high-quality training data for protein-ligand binding is complicated, however, by the requisite size of the models involved. Furthermore, the non-covalent nature of many protein-ligand interactions means that the electronic structure model must be chosen carefully. Fragment-based quantum chemistry leverages distributed computing by means of physics-based approximations, as an alternative to parallelization of conventional algorithms. In this way, O(N p ) computational scaling (where N measures system size and the exponent p depends on the electronic structure model) is reduced to N sub ×O(n p ), where n is a fixed subsystem size that does not grow with N , and N sub is the number of subsystems, which increases with N . This is an attractive approach to parallelization, in part because the storage footprint (i.e., memory and/or disk space), which is often the practical limitation, is reduced to that of the largest subsystem and checkpointing can be organized at the level of individual subsystem calculations. However, the number of fragments can be prohibitive for large systems, and must be culled via some kind of screening algorithm. This is necessary not only to reduce cost but also to forestall precision issues associated with calculations that might involve 10 5 or more individual separate subsystems. To this end, we have recently introduced a new software framework, Fragme∩t, with inherent database management, parallelization, and screening capabilities. It is built upon a generalized many-body expansion (MBE) and interfaced with numerous quantum chemistry codes. In recent applications, Fragme∩t has been used to investigate enzyme thermochemistry in large active-site models, and to perform high-order n-body calculations on water clusters and ionwater clusters. In the present work, we aim to apply fragmentation to protein-ligand interaction energies using enzyme models that include not just nearest-neighbor residues but which afford energetically converged interaction energies. Even at the level of density functional theory (DFT), there have been few studies with converged results for full-protein models of ligand binding. There has been other work applying fragment-based quantum chemistry to calculate protein-ligand interaction energies, mostly using DFT although a few studies using second-order Møller-Plesset perturbation theory in small basis sets have been published. The purpose of the present work is to establish protocols that are robust and reliable, which could eventually be used at better levels of theory. Crucially, we aim to compute interaction energies (∆E int ) that are faithful to a supramolecular calculation performed at the same level of theory (method and basis set), and to use sufficiently large molecular models so that ∆E int is converged with respect to further increases in system size. Our approach is based on the MBE truncated at n-body interactions [MBE(n)] and we examine convergence for n = 2-4 using single-residue fragments, in models containing up to 3,124 atoms. This represents unprecedented size and scope for application of MBE(n).