diff --git a/tests/test_stream.py b/tests/test_stream.py
index 66d63ad3..271999b5 100644
--- a/tests/test_stream.py
+++ b/tests/test_stream.py
@@ -31,24 +31,24 @@ def test_registration_alias():
 
 def test_vlle():
     tmo.settings.set_thermo(['Water', 'Ethanol', 'Octane'], cache=True)
-    s = tmo.Stream(None, Water=1, Ethanol=1, Octane=2, vlle=True, T=350)
+    s = tmo.Stream(None, Water=1, Ethanol=1, Octane=2, vlle=True, T=351)
     assert_allclose(s.mol, [1, 1, 2]) # mass balance
     total = s.F_mol
     xl = s.imol['l'].sum() / total
     xL = s.imol['L'].sum() / total
     xg = s.imol['g'].sum() / total
-    assert_allclose(xl + xL, 0.6773951226773982, atol=0.05) # Convergence
-    assert_allclose(xg, 0.32260487732260174, atol=0.05)
+    assert_allclose(xl + xL, 0.533447771190868, atol=0.05) # Convergence
+    assert_allclose(xg, 0.46655222880913216, atol=0.05)
     # Make sure past equilibrium conditions do not affect result of vlle
-    s = tmo.Stream(None, Water=1, Ethanol=1, Octane=2, T=350)
-    s.vle(T=350, P=101325)
-    s.vlle(T=350, P=101325)
+    s = tmo.Stream(None, Water=1, Ethanol=1, Octane=2, T=351)
+    s.vle(T=351, P=101325)
+    s.vlle(T=351, P=101325)
     assert_allclose(s.mol, [1, 1, 2]) # mass balance
     xl = s.imol['l'].sum() / total
     xL = s.imol['L'].sum() / total
     xg = s.imol['g'].sum() / total
-    assert_allclose(xl + xL, 0.6773951226773982, atol=0.05) # Convergence
-    assert_allclose(xg, 0.32260487732260174, atol=0.05)
+    assert_allclose(xl + xL, 0.533447771190868, atol=0.05) # Convergence
+    assert_allclose(xg, 0.46655222880913216, atol=0.05)
     
     s = tmo.Stream(None, Water=1, Ethanol=1, Octane=2, vlle=True, T=300)
     assert set(s.phases) == set(['l', 'L']) # No gas phase
diff --git a/thermosteam/equilibrium/lle.py b/thermosteam/equilibrium/lle.py
index 4310dd07..4284cfc9 100644
--- a/thermosteam/equilibrium/lle.py
+++ b/thermosteam/equilibrium/lle.py
@@ -60,7 +60,8 @@ def psuedo_equilibrium_inner_loop(Kgammay, z, T, n, f_gamma, gamma_args, phi):
     K = gammax / gammay 
     y = K * x
     y /= y.sum()
-    K = gammax / f_gamma(y, T, *gamma_args)
+    gammay = f_gamma(y, T, *gamma_args)
+    K = gammax / gammay
     Kgammay_new[:n] = K
     Kgammay_new[n:] = gammay
     return Kgammay_new
@@ -157,7 +158,7 @@ class LLE(Equilibrium, phases='lL'):
                  '_K',
                  '_phi'
     )
-    default_method = 'shgo'
+    default_method = 'pseudo equilibrium'
     shgo_options = dict(f_tol=1e-6, minimizer_kwargs=dict(f_tol=1e-6))
     differential_evolution_options = {'seed': 0,
                                       'popsize': 12,
diff --git a/thermosteam/separations.py b/thermosteam/separations.py
index 233148a3..8f069ac3 100644
--- a/thermosteam/separations.py
+++ b/thermosteam/separations.py
@@ -381,8 +381,8 @@ def lle_partition_coefficients(top, bottom):
     >>> IDs, K = tmo.separations.lle_partition_coefficients(s['L'], s['l'])
     >>> IDs
     ('Water', 'Ethanol', 'Octanol')
-    >>> K[2] # Octanol
-    3324.4
+    >>> round(K[2], -1) # Octanol
+    3330.0
 
     """
     IDs = tuple([i.ID for i in bottom.lle_chemicals])