Charges were in general converged with a threshold of 0.01 e. Geometries were optimised with dl-find 39 using convergence thresholds of 3•10 -6 a.u. and 1.35•10 -3 a.u. for energy and maximum gradient component, respectively. Strain components (in a.u.) were scaled by a factor of 265 to target roughly a convergence threshold of 3.40 • 10 -6 a.u. (0.1 GPa). For DFT calculations, atomic and point charge gradients from the underlying QM theory could be used to analytically calcu- late AC-QM/MM gradients. For the LCCSD(T0) calculations, however, AC-QM/MM gradients and stresses had to be calculated numerically making use of Equations ( ) and ( ). All DFT calculations were performed using ORCA3, 57 the coupled cluster calculations were performed using MOLPRO. NPA charges were calculated with GENNBO 5.9, Chelpg charges with Gaussian 09. IV. RESULTS