The formation of massive stars remains one of the most significant unsolved problems in astrophysics, with implications for the formation of the elements and the structure and evolution of galaxies. It is these stars, with masses greater than 8-10 solar masses, that eventually explode as supernovae and produce most of the heavy elements in the universe, dominate the energy injection into the interstellar medium of galaxies and by injecting both heavy elements and energy into the surrounding medium, shape the evolution of galaxies. High mass star formation poses a major theoretical challenge: How is it possible to sustain a sufficiently high mass accretion rate into a protostellar core despite the radiationpressure on the accreting envelope? I will first summarize our work on the first 3D simulations of massive star formation that resolve this longstanding problem using our high resolution 3D magneto-radiation-hydrodynamic adaptive mesh refinement code ORION, identifying 2 new mechanisms that efficiently solve the problem of the Eddington barrier to high mass star formation. By including the feedback effects of radiation pressure, protostellar radiation heating and protostellar outflows, I will discuss our work on stellar multiplicity. I will then present the first simulations of the effects of radiation feedback from high mass protostars on the universality of the IMF across a wide range of environments. Finally, I will present the first simulations of the formation of high mass clusters with radiation feedback.