Monte Carlo simulation has become an effective method for solving different physical problems, which is especially advantageous in the particle transport problems with complex geometric models. Charged particle transport simulation programs based on Monte Carlo method are widely used in medical therapy, radiation protection, accelerator design and other fields. Although the domestic neutron and photon transport simulation programs are well developed, there are few programs being able to do charged particle transport simulation. In addition, the energy range for programs based on the intranuclear cascade model and evaporative fission model is usually high, above several tens of MeV. For medium and low energy proton transport, the accuracy of nuclear reaction simulations based on nuclear databases is better than that of models. The biggest problem with Monte Carlo methods for simulating charged particle transport calculations is that the number of collisions of charged particles is too frequent. Therefore, this paper was based on the condensed history method to simulate the proton transport. In this study, with the application of key physical model of charged particle transport, such as the continuous slowing down approximation theory, multiple scattering theory and energy straggling model, a framework of heavy charged particles transport simulation program was established. Based on nuclear database of ENDF/B-Ⅶ.0 proton-nucleus interactions, the proton transport simulation function was implemented on the NPTS (neutron-photon transport simulation program). The extended NPTS program supported the simulation of proton transport problems in the energy range from 1 keV to 150 MeV and multi-coupling transport calculations of protons, neutrons, photons and electrons. A comparison of the stopping power of protons, deuterons and alpha particles in common materials with SRIM was carried out by the NPTS program, which proved the correctness of NPTS in calculating the stopping power of charged particles. The numerical results of neutron yield, average energy and angle for the thick target model agree well with GEANT4, MCNP6 and experimental results. In order to fully verify the program, a proton treatment chamber model with larger spatial dimensions and more complex geometry was used for the study. The neutron flux distribution inside the treatment chamber obtained by NPTS was basically consistent with GEANT4. The accuracy of the NPTS simulations of proton transport and multiparticle coupled transport was verified by several test models. Moreover, the proton energy and continuous slowing down approximation range curves obtained from the NPTS calculations can be used to quickly estimate the target thickness at maximum neutron yield, providing a favourable reference for BNCT (boron neutron capture therapy) target design. This paper effectively extends the application of the NPTS program in accelerator target design, shielding analysis, and radiation protection. The research effectively improves the application value of NPTS program, as well as enriching the research method of medium and low energy proton transport programs in domestic.