We have modeled a large sample of infrared starburst galaxies using both the PEGASE v2.0 and STARBURST99 codes to generate the spectral energy distribution (SED) of the young star clusters. PEGASE utilizes the Padova group tracks, while STARBURST99 uses the Geneva group tracks, allowing comparison between the two. We used our MAPPINGS III code to compute photoionization models that include a self-consistent treatment of dust physics and chemical depletion. We use the standard optical diagnostic diagrams as indicators of the hardness of the EUV radiation field in these galaxies. These diagnostic diagrams are most sensitive to the spectral index of the ionizing radiation field in the 1-4 ryd region. We find that warm infrared starburst galaxies contain a relatively hard EUV field in this region. The PEGASE ionizing stellar continuum is harder in the 1-4 ryd range than that of STARBURST99. As the spectrum in this regime is dominated by emission from Wolf-Rayet (W-R) stars, this discrepancy is most likely due to the differences in stellar atmosphere models used for the W-R stars. The PEGASE models use the Clegg & Middlemass planetary nebula nuclei (PNN) atmosphere models for the W-R stars, whereas the STARBURST99 models use the Schmutz, Leitherer, & Gruenwald W-R atmosphere models. We believe that the Schmutz et al. atmospheres are more applicable to the starburst galaxies in our sample; however, they do not produce the hard EUV field in the 1-4 ryd region required by our observations. The inclusion of continuum metal blanketing in the models may be one solution. Supernova remnant (SNR) shock modeling shows that the contribution by mechanical energy from SNRs to the photoionization models is 20%. The models presented here are used to derive a new theoretical classification scheme for starbursts and active galactic nucleus (AGN) galaxies based on the optical diagnostic diagrams.