Periodic boundary conditions were applied to the purely classical interactions, employing a 15 Å (14 Å for GB1) cutoff on the short-range nonbonded interactions, and using the smooth particle mesh Ewald (SPME) method for the long-range electrostatic interactions. The electrostatic QM-MM interactions were computed within the minimum image convention without periodic repetition. Isolated system conditions for the QM subsystem were applied using the Martyna-Tuckerman 64 method to solve Poisson's equation. The QM cell sizes and PW cutoffs were determined based on convergence analyses shown in the SI (see Figs. ) and taking into consideration the requirements of the Martyna-Tuckerman method. Thus, we chose a cell size of 30 a.u. for the Wat(aq) system and 40 a.u. for the other three systems, and PW cutoffs of 100, 80, 85, and 85 Ry for Wat(aq), Ace(aq), BuOH(ace), and GB1, respectively. The DIIS method was used in the self-consistent field (SCF) optimization of the wavefunction using a convergence threshold of 1 • 10 -5 a.u. with respect to the largest component of the electronic gradient. The initial guess of the wave function in each MD step is based on an extrapolation from the five previous steps to accelerate convergence and to improve energy conservation. To investigate the effect of the long-range coupling on the energy, forces, and dipole moment of the QM subsystem for different short-range cutoffs and multipole orders, we performed an analysis on 10 structures that were extracted from the last half of the simulations (20,000 steps) at regular intervals, i.e. every 2,000th step (∼ 0.97 ps). In the case of GB1, we calculated the Hirshfeld partial charge of the zinc ion instead of the dipole moment. The SCF convergence threshold was set to 1.0 • 10 -6 a.u. in these single-point calculations.