The phase field crystal equation (PFC) is a popular model for simulating micro-structures in materials science and is very computationally expensive to solve. A highly scalable solver for PFC modeling is presented in this paper. The equation is discredited with a stabilized implicit finite difference method and the time step size is adaptively controlled to obtain physically meaningful solutions. The nonlinear system arising at each time step is solved by using a parallel Newton-Krylov-Schwarz algorithm. In order to achieve good performance, low-order homogeneous boundary conditions are imposed on the sub domain boundary in the Schwarz preconditioner. Experiments are carried out to exploit optimal choices of the preconditioner type, the sub domain solver and the overlap size. Numerical results are provided to show that the solver is scalable to thousands of processor cores.