To perform spin-orbit coupling calculations on atoms and molecules, good zeroth-order wavefunctions are necessary. Here, we present the software development of the Monte Carlo Configuration Interaction (MCCI) method, to enable calculation of such properties, where MCCI iteratively constructs a multireference wavefunction using a stochastic procedure. In this initial work, we aim to establish the efficacy of this technique in predicting the splitting of otherwise degenerate energy levels on a range of atoms and small diatomic molecules. It is hoped that this work will subsequently act as a gateway toward using this method to investigate singlet-triplet interactions in larger multireference molecules. We show that MCCI can generate very good results using highly compact wavefunctions compared to other techniques, with no prior knowledge of important orbitals. Higher-order relativistic effects are neglected and spin-orbit coupling effects are incorporated using first-order degenerate perturbation theory with the Breit-Pauli Hamiltonian and effective nuclear charges in the one-electron operator. Results are obtained and presented for B, C, O, F, Si, S, and Cl atoms and OH, CN, NO, and C2 diatomic radicals including spin-orbit coupling constants and the relative splitting of the lowest energy degenerate state for each species. Convergence of MCCI to the full configuration interaction result is demonstrated on the multireference problem of stretched OH. We also present results from the singlet-triplet interaction between the math formula and both the math formula and math formula states of the O2 molecule.