where ! is a mixing coe cient introduced to prevent polarization "sloshing", encountered when modelling the crystal. Similar mixing schemes have been implemented for accelerating the convergence. Here, t indexes the SCF step. A more general formulation and introduction to the notation of the perturbative expansion of the electrostatic interactionresulting in the multipole moment model -can be found elsewhere. The inclusion of polarization requires the multiplication of damping functions on the interaction tensors T ij in equation 7. They are derived from considering the Gaussian charge density overlap to avoid the so-called 'polarization catastrophe'. The spatial extent of these damping functions are described by a parameter g -damping length. See the Supporting Information in this work, and in references for more details.