diff --git a/src/pyEQL/engines.py b/src/pyEQL/engines.py
index 0087072f..e79ceb6a 100644
--- a/src/pyEQL/engines.py
+++ b/src/pyEQL/engines.py
@@ -1021,7 +1021,14 @@ def get_activity_coefficient(self, solution: "solution.Solution", solute: str) -
         return ureg.Quantity(act, "dimensionless")
 
     def get_osmotic_coefficient(self, solution: "solution.Solution") -> ureg.Quantity:
+        if self.ppsol is not None:
+            self.ppsol.forget()
+        self._setup_ppsol(solution)
+
         osmotic = self.ppsol.get_osmotic_coefficient()
+        if osmotic == 0:
+            # PHREEQC returns 0 when it assumes a unit osmotic coefficient, so we need to catch that case and return 1 instead.
+            osmotic = 1
         return ureg.Quantity(osmotic, "dimensionless")
 
     def get_solute_volume(self, solution: "solution.Solution") -> ureg.Quantity:
diff --git a/tests/test_phreeqc.py b/tests/test_engine_phreeqc.py
similarity index 97%
rename from tests/test_phreeqc.py
rename to tests/test_engine_phreeqc.py
index 31828603..22334b3c 100644
--- a/tests/test_phreeqc.py
+++ b/tests/test_engine_phreeqc.py
@@ -3,7 +3,7 @@
 =================================================
 
 This file contains tests for the volume and concentration-related methods
-used by pyEQL's Solution class
+used by pyEQL's Solution class using the `phreeqc` engine (phreeqpython)
 """
 
 import logging
@@ -128,6 +128,11 @@ def test_init_engines():
     s.engine._destroy_ppsol()
     assert s.engine.ppsol is None
 
+def test_osmotic_pressure():
+    s1 = Solution([["Na+", "1 mol/L"], ["SO4-2", "0.5 mol/L"]], engine="phreeqc")
+    s2 = Solution([["Na+", "1 mol/L"], ["SO4-2", "0.5 mol/L"]], engine="ideal")
+    assert np.isclose(s1.osmotic_pressure.to("MPa").magnitude, s2.osmotic_pressure.to("MPa").magnitude), f"PHREEQC = {s1.osmotic_pressure.to('MPa').magnitude} MPa, Ideal = {s2.osmotic_pressure.to('MPa').magnitude} MPa"
+
 
 def test_conductivity(s1):
     # even an empty solution should have some conductivity
diff --git a/tests/test_pyeql_phreeqc.py b/tests/test_engine_phreeqc2026.py
similarity index 96%
rename from tests/test_pyeql_phreeqc.py
rename to tests/test_engine_phreeqc2026.py
index 20ccccd9..ef9be7f7 100644
--- a/tests/test_pyeql_phreeqc.py
+++ b/tests/test_engine_phreeqc2026.py
@@ -3,7 +3,7 @@
 =================================================
 
 This file contains tests for the volume and concentration-related methods
-used by pyEQL's Solution class
+used by pyEQL's Solution class using the `phreeqc2026` engine.
 """
 
 import logging
@@ -122,13 +122,19 @@ def test_init_engines():
     s = Solution([["Na+", "4 mol/L"], ["Cl-", "4 mol/L"]], engine="phreeqc2026")
     assert isinstance(s.engine, Phreeqc2026EOS)
     assert s.get_activity_coefficient("Na+").magnitude * s.get_activity_coefficient("Cl-").magnitude < 1
-    assert s.get_osmotic_coefficient().magnitude == 0
+    # PHREEQC itself returns a value of 0 for this, but we need to use a value of 1 (reprseneting ideal solution)
+    # in order to ensure correct osmotic pressure output.
+    assert s.get_osmotic_coefficient().magnitude == 1
     # with pytest.warns(match="Solute Mg+2 not found"):
     assert s.get_activity_coefficient("Mg+2").magnitude == 1
     assert s.get_activity("Mg+2").magnitude == 0
     s.engine._destroy_ppsol()
     assert s.engine.ppsol is None
 
+def test_osmotic_pressure():
+    s1 = Solution([["Na+", "1 mol/L"], ["SO4-2", "0.5 mol/L"]], engine="phreeqc2026")
+    s2 = Solution([["Na+", "1 mol/L"], ["SO4-2", "0.5 mol/L"]], engine="ideal")
+    assert np.isclose(s1.osmotic_pressure.to("MPa").magnitude, s2.osmotic_pressure.to("MPa").magnitude), f"PHREEQC2026 = {s1.osmotic_pressure.to('MPa').magnitude} MPa, Ideal = {s2.osmotic_pressure.to('MPa').magnitude} MPa"
 
 def test_conductivity(s1):
     # even an empty solution should have some conductivity