The search for neutrinoless double beta decay of nuclei is believed to be one of the most promising means to search for new physics. Observation of this very rare nuclear process, which violates Lepton Number conservation, would imply the neutrino sector has a Majorana mass component and may also provide an explanation for the universe matter-antimatter asymmetry of the universe. In the case where a heavy intermediate particle is exchanged in this process, QCD contributions from short range interactions become relevant and the calculation of matrix elements with four-quark operators becomes necessary. In these proceedings we will discuss our current progress in the calculation of these four-quark operators from LQCD.