<p>Niobium diselenide (NbSe<sub>2</sub>) has garnered significant attention due to the coexistence of superconductivity and charge density waves (CDWs) down to the monolayer limit. However, realistic modeling of CDWs—capturing effects such as layer number, twist angle, and strain—remains challenging due to the high computational cost of first-principles methods. Here, we develop a physically informed workflow for training machine-learning interatomic potentials (MLIPs) based on the E(3)-equivariant Allegro architecture, tailored to capture the subtle structural and dynamical signatures of CDWs in mono- and bilayer NbSe<sub>2</sub>. We find that while CDW lattice distortions are relatively easy to learn, modeling vibrational properties remains more challenging. It requires targeted dataset design and careful hyperparameter tuning, pushing the boundaries and testing the extensibility of current MLIP frameworks. Our MLIPs enable reliable simulations of commensurate and incommensurate CDW phases, including their sensitivity to dimensionality and stacking, as well as CDW dynamics, phonons, and transition temperatures estimated via the stochastic self-consistent harmonic approximation. This work opens new possibilities for studying and tuning CDWs in NbSe<sub>2</sub> and other two-dimensional systems, with implications for electron-phonon coupling, superconductivity, and advanced materials design.</p>
