Efficient molecular dynamics simulation of 2D penta-silicene materials using machine learning potentials

Machine Learning Interatomic Potentials (MLIPs) are a modern computational method that allows achieving near-quantum mechanical accuracy (DFT) while still describing large-scale systems in molecular dynamics (MD) simulations. In this work, we use MLIP from DeepMD package and the classical Tersoff potential for SiC (Tersoff.SiC potential) to fully and accurately describe atomic interactions and apply them to molecular dynamics simulations of penta silicene sheet. The results show that the melting points (T$_g$) temperatures of the system in the canonical NVT and isobaric NPT sets are 632 K and 606 K, while the Tersoff.SiC potential have the high melting points, respectively. In addition, the radial distribution function exhibits characteristic peaks at interatomic distances of 2.275 Å\text{} and 2.375 Å, while the Tersoff.SiC potential only describe distance of 2.375 Å. Furthermore, penta silicene was also simulated using on-the-fly machine learning for 10 ps to evaluate the structural stability of the system. This study investigates the thermodynamic properties of two-dimensional penta silicene sheets with pentagonal structures using a high-precision, cost-effective method, contributing further evidence to support experimental synthesis and opening up potential future applications of this material.

💡 Research Summary

This paper presents a comprehensive study on the use of machine‑learning interatomic potentials (MLIPs) to enable large‑scale molecular dynamics (MD) simulations of two‑dimensional penta‑silicene, a silicon allotrope composed of pentagonal rings. The authors combine a deep‑learning potential generated with the DeepMD‑kit package and a classical Tersoff SiC potential to evaluate structural, electronic, and thermodynamic properties of the material.

First, density‑functional theory (DFT) calculations were performed with VASP (PBE functional, PAW pseudopotentials, 400 eV cutoff, 12 × 12 × 1 k‑point mesh) on a 4 × 4 × 1 supercell containing 96 Si atoms. Ab‑initio molecular dynamics (AIMD) simulations were carried out at temperatures ranging from 100 K to 900 K in 100 K increments, each for 0.5 ps with a 0.5 fs timestep, using a Nose‑Hoover thermostat in the canonical (NVT) ensemble. Only the Γ‑point was sampled to keep the data‑generation cost low. The resulting trajectories provided forces, energies, and virial stresses for training.

The DeepMD‑kit workflow constructs atomic descriptors D_ij that encode inter‑atomic distances and directional components within a cutoff radius of 6 Å (smoothing radius 5.8 Å). A neural network with three hidden layers (25, 50, 100 neurons) and a 240‑neuron output layer was trained. The loss function combined energy, force, and virial terms with weighting factors 1, 0.1, and 0.001, respectively. Training proceeded for 1 × 10⁶ steps, with an exponentially decaying learning rate from 1 × 10⁻³ to 3.5 × 10⁻⁸. The final model achieved root‑mean‑square errors (RMSE) of 6 × 10⁻⁶–2 × 10⁻⁵ eV/atom for energy and 4 × 10⁻³–1 × 10⁻² eV Å⁻¹ for forces, the lowest errors occurring at 300 K and 400 K where the training data were most diverse.

Using the trained MLIP, the authors performed classical MD simulations in LAMMPS on a much larger penta‑silicene sheet (2400 atoms, 111.63 Å × 111.63 Å). Simulations were run in both NVT and isothermal‑isobaric (NPT) ensembles with a 1 fs timestep. The melting (glass transition) temperatures were identified as 632 K (NVT) and 606 K (NPT), considerably lower than the values obtained with the Tersoff‑SiC potential, indicating that the MLIP captures the weaker Si–Si bonding in the pentagonal network more accurately.

Radial distribution function (RDF) analysis revealed two distinct peaks at 2.275 Å and 2.375 Å, corresponding to the two Si–Si bond lengths observed in the optimized geometry (2.28 Å and 2.39 Å). The classical Tersoff potential reproduced only the 2.375 Å peak, missing the shorter bond contribution. Phonon calculations (using PHONOPY) confirmed dynamical stability with only a negligible imaginary mode near Γ, while electronic band‑structure calculations showed an indirect band gap of ~0.24 eV, consistent with previous reports.

To assess the robustness of the potential during a simulation, an on‑the‑fly learning scheme was employed. AIMD runs from 300 K to 700 K (1 fs timestep, 10 ps total) were used to continuously update the MLIP while the system evolved. The structure remained intact up to 700 K, with only minor distortions, demonstrating that the potential can be refined in real time and still preserve structural fidelity over short timescales.

The paper concludes that (i) a relatively small AIMD dataset (0.5 ps per temperature) is sufficient to train a high‑accuracy MLIP for a complex 2D material, (ii) MLIPs outperform traditional empirical potentials in predicting melting points and detailed RDF features, and (iii) on‑the‑fly training provides a practical route for real‑time validation of structural stability. Limitations include the short AIMD sampling window, the use of only the Γ‑point for DFT calculations, and the modest 10 ps on‑the‑fly simulation length, which may not capture long‑time diffusion or defect dynamics. Future work should expand the training set with longer trajectories, multiple k‑points, and diverse defect or strain configurations to generalize the approach to other low‑dimensional pentagonal materials and to enable reliable predictions for experimental synthesis and device applications.