A theoretical approach is presented which represents the first attempt known to the authors to develop molecular dynamics algorithms capable of modeling phase coexistence between two nonequilibrium steady state phases confined in a closed (E, V, N) system. We deal exclusively with shearing liquids because of their importance in rheology. In the present paper, as in the equilibrium Gibbs ensemble Monte Carlo technique for systems at equilibrium, the coexisting phases have no physical contact but their dynamics are coupled in order to reach mechanical, thermal, and composition balance between bulk regions of the two phases. The thermal balance is maintained by requiring zero net heat flow across a hypothetical boundary. This can be achieved by starting from equilibrium and gradually increasing the strength of the external field (the shear rate) in a quasistatic process. For particle interchanges we invoke the Evans-Baranyai variational principle which is at the very least a good approximation for similar simulated steady state systems far from equilibrium. Results of several model calculations are presented. The Limitations and the implications of the methods are discussed.