Initial reaction mechanisms of lignin pyrolysis were studied by large-scale ReaxFF molecular dynamics simulations (ReaxFF MD) facilitated by the first GPU-enabled code (GMD-Reax) and the unique reaction analysis tool (VARxMD). Simulations were performed over wide temperature ranges both for heat up at 300-2100 K and for NVT at 500-2100 K with a large lignin model, which contained 15920 atoms and was constructed based on Adler's softwood lignin model. By utilizing the relatively continuous observation for pyrolysate evolution in slow heat up simulations, three stages for lignin pyrolysis are proposed by pyrolysate fractions. The underlying mechanisms for the three stages are revealed by analyzing the species structure evolution and the reactions of linkages, aryl units, propyl chains, and methoxy substituents. Stage I is characterized with the complete decomposition of source lignin molecules at low temperatures dominated by breaking of alpha-O-4 and beta-O-4 linkages. The temperature in stage II is relatively high where cracking of all the linkages occurs, accompanied by conversion of propyl chains and methoxy substituents. Stage III mapping to high temperature shows the formation of heavy pyrolysates by recombination reactions of five-, six-, or seven-membered aliphatic rings. The heterocyclic oxygen-containing rings are revealed as important intermediates for the aryl monomer ring opening into aliphatic rings of five-membered, seven-membered, or even larger. The pathways for small molecule formation observed in this work are broadly in agreement with the literature. This work demonstrates a new methodology for investigating the overall behaviors and the underlying complex mechanisms of lignin pyrolysis.