Reliable calculations of nuclear matrix elements are a prerequisite for the determination of the effective neutrino mass and other particle physics parameters from neutrinoless double beta decay. Here, the operator expansion method is improved by including Coulomb, tenser and central interactions simultaneously. Furthermore, the formalism of the OEM is extended to those matrix elements necessary to extract the right-handed parameters [lambda] and [eta] from O nu beta beta decay. OEM includes the dependence of the nuclear matrix elements on the intermediate states implicitly and can therefore be understood as a step beyond the closure approximation. Numerical studies are carried out for the isotope Ge-76 combining the OEM expressions with ground-state wave functions calculated within a proton-neutron quasiparticle Random Phase Approximation (pn-QRPA) model. The influence and relative importance of central, tenser and Coulomb interactions is investigated. Within the OEM, contributions from the Coulomb force are found to be negligible in O nu beta beta decay, while the tenser force leads to a moderate change of the results, of the order of (10-30) %, giving a better agreement between sets of calculations which employ different NN-interactions. Generally, results of the OEM + QRPA calculation are similar to previous calculations of O nu beta beta decay matrix elements, indicating that O nu beta beta decay is not sensitive to model approximations and might therefore be more accurately calculated than the strongly suppressed 2 nu beta beta decay matrix elements.