For valence electrons, a plane-wave basis set was adopted with an energy cutoff of 450 eV, and the ionic cores were described with the projector augmented-wave (PAW) method. The total energy was converged to an accuracy of 1×10 -6 eV, and a force tolerance of 0.03 eV/Å was used in all structure optimization. The NNPs were trained via multistream extended Kalman filter algorithm using energies and forces. To improve the quality of the NNPs, an expanded dataset was built via active learning using a modified framework based on RuNNerActiveLearn . The active learning procedure was stopped when no extrapolation was found in the runs and all energy (force) differences between the two NNPs are lower than 5 (10) times the RMSE of energy (forces). If the potential is not applicable to the above three systems, the active learning would be launched again, and the temperature of NN-MDs would be increased to increase structural diversity. Further tests regarding temperature stability and melting point are presented in Section S2. The simulation of Neural Network based Molecular Dynamic (NN-MD) were performed using LAMMPS code with the NNP interface from n2p2 , For equilibrium, the simulations were run for 1 ns at the indicated temperature. For annealing, the simulations were run at 700 K for 1ns, then the systems were cooled to 300K in 1 ns, and finally another 1 ns for equilibrium at 300 K. The canonical (NVT) ensemble simulations were modeled with the Nose-Hoover thermostat for surface systems, and the Nose-Hoover thermostat and barostat were employed for isobaric-isothermic (NPT) ensemble simulations of bulk systems .