1) It could be calculated in vacuo and kept fixed ({q vac j }); 2) be calculated in vacuo for each simulation step ({q vac j (R QM )}, where (R QM ) indicates the charges are a function of the coordinates of the QM atoms); or 3) calculated in presence of the environment (namely, the solvent) for each simulation step ({q sol j (R QM )}). In the first case, the atomic partial charges of the QM region would play the role of fixed parameters that need to be assigned when building the system topology. In options 2 and 3, the charges are recomputed on-the-fly. In this work, the atomic partial charges will always be dependent on the system geometry, and the (R QM ) notation will be omited.