Recently, we introduced a novel computer simulation technique to determine the optimal set of parameters of an interaction potential for a simple monatomic liquid. This technique was used to obtain interaction potentials of the Lennard-Jones form that accurately describe argon over its entire liquid phase at a fixed pressure [S. D. Bembenek and B. M. Rice, Mol. Phys. 97, 1085 (1999)]. Here, we extend this technique to a homonuclear diatomic molecular system, liquid oxygen. This technique was first applied to a system in which the oxygen molecules were treated as point masses interacting through a modified Lennard-Jones potential. Simulations using the resulting optimal set of potential parameters of this system predict densities that are within 0.25% of experiment over the entire liquid range of oxygen at a fixed pressure. However, the errors in the internal energy and the enthalpy were as large as 9.8%. The technique was then used to determine the optimal parameters for a system of harmonic molecules, in which each molecule has two interaction sites centered at the atomic nuclei. The intermolecular interaction is the sum of all site-site interactions described by a modified Lennard-Jones potential. Simulations using these parameters reproduce experimental densities with an error no greater than 0.80%. The predictions of the internal energy and enthalpy differ from experiment by no more than 3.0% for temperatures below 90 K; predictions at 90 K differ from experiment by no more than 4.11%. These results seem to suggest that our method for determining parameters for an interaction potential is also applicable to simple molecular systems.