Neural network molecular dynamics (NNMD) simulations could revolutionize atomistic modeling of materials with quantum-mechanical accuracy at a fraction of computational cost. However, popular NNMD frameworks are generally implemented for a single computing node, and conventional energy-based NN models still suffer from large time-to-solution (T2S), prohibiting the application of NNMD to challenging materials simulations encompassing large spatiotemporal scales. Consequently, no leadership-scale NNMD simulation has thus far been reported. Here, we present a scalable parallel NNMD software (RXMD-NN) based on our scalable reactive molecular dynamics simulation engine named RXMD. RXMD-NN has achieved high scalability up to 786,432 IBM BlueGene/Q cores involving 1.7 billion atoms. Furthermore, we have achieved 4.6-fold reduction of T2S by using a novel network that directly predicts atomic forces from feature vectors. Reduced T2S has for the first time allowed the study of large-scale off-stoichiometry effect in a widely used phase change material, Ge2 Se2 Te5, thereby resolving its “firstsharp diffraction peak mystery”.