DescriptionReactive molecular dynamics is a powerful simulation method for describing chemical reactions. Here, we introduce a new generalized polarizable reactive force-field (ReaxPQ+) model to significantly improve the accuracy by accommodating the reorganization of surrounding media. The increased computation is accelerated by (1) extended Lagrangian approach to eliminate the speed-limiting charge iteration, (2) shiftcollapse computation of many-body renormalized n-tuples, which provably minimizes data transfer, (3) multithreading with roundrobin data privatization, and (4) data reordering to reduce computation and allow vectorization. The new code achieves (1) weak-scaling parallel efficiency of 0.989 for 131,072 cores, and (2) eight-fold reduction of time-to-solution (T2S) compared with the original code, on an Intel Knights Landing-based computer. The reduced T2S has for the first time allowed purely computational synthesis of atomically-thin transition metal dichalcogenide layers assisted by machine learning to discover a novel synthetic pathway.