Precise predictions of the radiation environment inside space vehicles, and inside the human body, are essential when planning for long term deep space missions. Since these predictions include complex geometries, as well as the contributions from many different types of radiation, including neutrons, 3-D Monte Carlo codes with precise physics models are needed. In this paper, we present improvements and developments of some physics models used in the general purpose 3-D Monte Carlo code PHITS [1]. The total reaction cross section (sigma(R)) and the decay lifetime of a projectile particle are the first essential quantities in MC calculations, since these determine the mean free path of the transported particles and the probability function according to which a projectile particle will collide within a certain distance in the matter depends on the sigma(R). This will also scale the calculated partial fragmentation cross sections. In this paper we present comparisons of calculated and measured sigma(R) using the Kurotama Hybrid sigma(R), model [2] which is incorporated into PHITS. The prediction of the fragmentation reactions of relativistic heavy ions is also essential for ensuring radiation safety of astronauts. The default model for nuclear-nuclear reactions is JQMD in PHITS. However, JQMD cannot accurately enough describe the nucleon and d, t, He-3 and 4He induced reactions. Therefore the Intra-Nuclear Cascade of Liege (INCL) [3] has been selected as the default model for these reactions. Moreover, it has been realized that the production of light fragments is underestimated by conventional simulation codes based on a combination of intranuclear cascade and statistical decay models. This is because this combination cannot reproduce the high multiplicity events that are responsible for the production of light fragments. To better reproduce high multiplicity events, we have simulated fragmentation cross sections using a combination of JQMD/INCL, statistical multi-fragmentation model (SMM) [4,5] and the generalized evaporation model (GEM). Examples of these simulations will be presented. A new approach to describe neutron spectra of deuteron-induced reactions in the Monte Carlo simulations has also been developed by combining the INCL and the Distorted Wave Born Approximation (DWBA) calculation [6]. We have incorporated this combined method into PHITS and applied it to estimate (d,xn) spectra on light targets at incident energies ranging from 10 to 40 MeV. In this paper, we will show that the double differential cross sections obtained by INCL and DWBA successfully reproduced broad peaks and discrete peaks, respectively.