The QM-MM interactions can be divided into two parts, namely bonded and nonbonded interactions, where the former only appear when there are covalent bonds crossing the QM-MM boundary. The bonded interactions are described classically using the parameters of the same force field that is used in the MM subsystem. Open valencies in the QM subsystem are saturated by adding boundary atoms that are described by a monovalent pseudopotential and thus contribute a single electron to the QM subsystem. The nonbonded interactions consist of electrostatic and van der Waals interactions. The latter are described through the MM force field. Since the bonded and van der Waals interactions are purely classical, we delegate these to the external MM engine. The electrostatic interactions, which include interactions between the MM charges (or multipoles) and the QM nuclei and electrons, are computed by MiMiC. The current implementation in MiMiC supports an electrostatic embedding scheme in which the electrostatic potential from the MM charges (or multipoles) is added to the KS-DFT external potential. This involves the electronic density of the QM subsystem, which is mapped on a grid, typically of the order of 10 6 points, and thus requires special attention, in particular with respect to computational cost and parallel scalability. Evaluating the full electrostatic part of the nonbonded energy exactly is not only very computationally expensive but also exhibits poor scalability with the size of the system, which can be of the order of 10 5 atoms or larger, since we have to compute the full integral over the electronic grid for each of the MM charges (see below). To overcome these limitations, we have previously formulated and implemented a generalized version of the electrostatic coupling scheme by Laio et al., which was shown to yield the same accuracy as an exact electrostatic treatment but at a substantially reduced cost. Within the generalized electrostatic coupling scheme, the MM atoms are grouped into two domains according to their distance from the QM subsystem. The MM atoms that are within a specified distance from the QM subsystem (the short-range domain) interact directly with the QM electron density and nuclei, whereas the remaining MM atoms (the long-range domain) interact through a multipole expansion of the electrostatic potential from the QM electrons and nuclei. In multi-index notation, as introduced in Ref. 1, the electrostatic QM-MM interaction energy can thus be written as