Purpose: This study investigates the influence of several Monte Carlo radiation transport codes and nuclear models on the simulation of secondary neutron spectra and its impact on calculating and measuring the neutron doses in proton therapy. Materials and methods: Three different multi-purpose Monte Carlo radiation transport codes (FLUKA, MCNPX, Geant4) were used together with different available nuclear models, to calculate secondary neutron energy spectra at various points inside a water tank phantom with PMMA walls using a 10 × 10 cm2 rectangular, monoenergetic proton beam (110 MeV, 150 MeV, 180 MeV, 210 MeV). Using Kerma approximation secondary neutron doses were calculated applying fluence-to-dose equivalent conversion coefficients in water. Moreover, the impact of varying spectra for electrochemically etched CR39 detector calibration was analyzed for different codes and models. Results: In distal positions beyond the Bragg peak, results show largest variations between the codes, which was up to 53% for the high energy neutron fluence at 16 cm from the Bragg peak of the 110 MeV proton beam. In lateral positions, the variation between the codes is smaller and for the total neutron fluence within 20%. Variation in the nuclear models in MCNPX was only visible for the proton beam energies of 180 and 210 MeV and modeling the high energy neutron fluence which reached up to 23% for 210 MeV at 11 cm lateral from the beam axis. Impact on neutron dose equivalent was limited for the different models used (<8%) while it was pronounced for the different codes (45% at 16 cm from the Bragg peak of the 110 MeV proton beam). CR39 calibration factors in lateral positions were on average varying 10% between codes and 5% between nuclear models. Conclusions: This study demonstrated a large impact on the neutron fluence spectra calculated by different codes while the impact of different models in MCNPX proved to be less prominent for the neutron modeling in proton therapy.