A neural network (NN) - based machine learned interatomic potential (MLIP) was developed for simulating the shock response of boron carbide (B­4C). The training data set consisted of atomic species and coordinates of four different polymorphs, their associated total energy, forces, and virial stresses for a broad range of atomic configurations (e.g., various strain levels in tension, compression, and shear loading as well as atomic configurations in different material phases). This data was generated with ab initio methods to maximize accuracy. The MLIP was validated using a three-stage material-agnostic procedure for stability (e.g., energy, force, virial stresses) static property prediction (e.g., elastic constants, equation of state, lattice constants), and dynamic property prediction (e.g., melting point, phase transition and high-pressure response). Shock simulations were then performed at a range of impact velocities to determine the Hugoniot elastic limit (HEL). The evolution of amorphous bands in different polymorphs at different orientation was systematically investigated. It is shown that amorphization causes local increase in the volume and hence causes significant compressive stress within the surrounding crystalline region thus revealing the mechanism for the lack of cracking which has been observed in numerous microscopic observations.