We use the Grand-reaction Monte Carlo method, 105,106 which allows us to simulate chemical equilibrium within and co-existence between the supernatant and coacervate phases in a broad range of salt concentrations, pH values and in the presence of small molecular or ionic solutes. The Grand-reaction method guarantees the equality of chemical potentials of the partitioning species between the phases, equality of osmotic pressures between the phases, and it also provides an option to assess its self-consistency. Similar to the preceding studies, our model assumes that there is no polymer in the supernatant phase. This approximation works well far from the critical point whereas it fails in the critical region, which defines the range of applicability of our approach. In the current work, we use only conventional trial Monte Carlo insertion moves, which can become inefficient when simulating the partitioning of bulky solutes. This deficiency can be suppressed by coupling our method to the continuous fractional component ansatz 107,108 or to thermodynamic integration, which we plan to do in the future.