A method for computing the singlet-triplet energy difference (STED) of molecules with nearly degenerate states is described. It is based on the shifted contour auxiliary field Monte Carlo, implemented with plane waves and pseudopotentials. Two additional critical elements are: a variational multireference approach, for treating non-dynamical correlation and direct correlated sampling of energy differences. Applications to the H-He-H and CH2 systems are given, estimating STED of - 14.3 +/-0.1 and 9.8 +/-0.3 kcal/mol, respectively. The overall accuracy is limited by the type of pseudopotential used, with preference to those based on generalized gradients.