Monte Carlo (MC) radiation transport methods are used for dose calculation as 'gold standard.' However, the method is computationally time-consuming and thus impractical for normal tissue dose reconstructions for the large number of proton therapy patients required for epidemiologic investigations of late health effects. In the present study, we developed a new dose calculation method for the rapid reconstruction of out-of-field neutron dose to patients undergoing pencil beam scanning (PBS) proton therapy. The new dose calculation method is based on neutron dose voxel kernels (DVKs) generated by MC simulations of a proton pencil beam irradiating a water phantom (60 × 60 × 300 cm3), which was conducted using a MC proton therapy simulation code, TOPAS. The DVKs were generated for 19 beam energies (from 70 to 250 MeV with the 10 MeV interval) and three range shifter thicknesses (1, 3, and 5 cm). An in-house program was written in C++ to superimpose the DVKs onto a patient CT images according to proton beam characteristics (energy, position, and direction) available in treatment plans. The DVK dose calculation method was tested by calculating organ/tissue-specific neutron doses of 1- and 5-year-old whole-body computational phantoms where intracranial and craniospinal irradiations were simulated. The DVK-based doses generally showed reasonable agreement with those calculated by direct MC simulations with a detailed PBS model that were previously published, with differences mostly less than 30% and 10% for the intracranial and craniospinal irradiations, respectively. The computation time of the DVK method for one patient ranged from 1 to 30 min on a single CPU core of a personal computer, demonstrating significant improvement over the direct MC dose calculation requiring several days on high-performance computing servers. Our DVK-based dose calculation method will be useful when dosimetry is needed for the large number of patients such as for epidemiologic or clinical research.