The neutron production from alpha particles in galactic cosmic rays (GCR) in the lunar subsurface has not been estimated with reliable precision despite its importance for lunar nuclear spectroscopy and space dosimetry. Here, we report our estimation of neutron production from GCR nuclei (protons and alpha particles) with the Particle and Heavy Ion Transport code System (PHITS), which includes several heavy ion interaction models. PHITS simulations of the equilibrium neutron density profiles in the lunar subsurface are compared with experimental data obtained in the Apollo 17 Lunar Neutron Probe Experiment. Our calculations successfully reproduced the data within an experimental error of 15%. Our estimation of neutron production from GCR nuclei, estimated by scaling that from protons by a factor of 1.27, is in good agreement within an error of 1% with the calculations using two different alpha particle interaction models in PHITS during a period of average activity of the solar cycle. However, we show that the factor depends on the incident GCR spectrum model used in the simulation. Therefore, we conclude that the use of heavy ion interaction models is important for estimating neutron production in the lunar subsurface.