For QM calculations the DFTB+ software package is employed, which computes the electronic structure at the SCC-DFTB level . For this class of systems, we have previously found that a relatively loose set of simulation parameters may be used to obtain satisfactorily converged atomic charges. This helps to keep the cost of the QM part of the simulation to a minimum. To this end we sample the 1st Brillouin zone at the point only and we apply a Fermi temperature of 1×10 -6 K for smearing the occupancies of the electronic states. The convergence threshold of the self-consistent charge optimization is set to 1 × 10 -2 Hartree. In the MD simulation, we run the model under NVT ensemble using GROMACS (version 2018.4) . The Nosé-Hoover thermostat is employed fixing the temperature of 298 K with relaxation times of 1ps. The cutoff distance of short range Lennard-Jones and electrostatic interactions is set to 1.4 nm. The long-range electrostatic interaction is treated by the Particle Mesh Ewald approach with the standard (3D) geometry. The switch function of the Lennard Jones 12-6 potential is onset from 1 nm to smoothly truncate at 1.4 nm. The TIP4P-2005 is used as the water model with the rigid constraint by SETTLE algorithm. The parameters of ions (Li + and Cl -) and carbon atoms are taken from the Madrid-2019 force field and the Amber force-field , respectively. The