We extracted at least 5000 snapshots from the AIMD trajectories, evenly spaced by 2.5 fs. For each snapshot, we computed MLWFs, Berry phase dipole moments and polarizabilities using the CP2K program package. For all CP2K calculations, we employed the PBE exchange-correlation functional, Goedecker-Teter-Hutter (GTH) pseudopotentials and the DZVP-MOLOPT-SR-GTH basis sets. The CP2K calculations of the dipole moments and polarizabilities were conducted only at the Γ-point. Snapshots from the AIMD trajectories with the SC setup were used without modification, while the snapshots from the KP trajectories were extended to an 1 × 1 × 2 supercell.