Following a review, the only safe model to fit and have sigmas always…

… be positive is EQ106 with a positive A value.

chemicals/interface.py

@@ -413,6 +413,14 @@ def Watson_sigma(T, Tc, a1, a2, a3=0.0, a4=0.0, a5=0.0):
     -----
     This expression is also used for enthalpy of vaporization in [1]_.
     The coefficients from NIST TDE for enthalpy of vaporization are kJ/mol.
+    
+    This model is coded to return 0 values at Tr >= 1. It is otherwise not
+    possible to evaluate this expression at Tr = 1, as log(0) is undefined
+    (although the limit shows the expression converges to 0).
+    
+    This equation does not have any term forcing it to become near-zero
+    at the critical point, but it cannot be fit so as to produce negative 
+    values.
 
     Examples
     --------
@@ -427,6 +435,8 @@ def Watson_sigma(T, Tc, a1, a2, a3=0.0, a4=0.0, a5=0.0):
        https://trc.nist.gov/TDE/Help/TDE103b/Eqns-Pure-SurfaceTension/HVPExpansion-SurfaceTension.htm
     '''
     Tr = T/Tc
+    if Tr >= 1.0:
+        return 0.0
     l = log(1.0 - Tr)
     return exp(a1 + l*(a2 + Tr*(a3 + Tr*(a4 + a5*Tr))))
 
@@ -461,6 +471,9 @@ def ISTExpansion(T, Tc, a1, a2, a3=0.0, a4=0.0, a5=0.0):
 
     Notes
     -----
+    This equation hsa a term term forcing it to become zero
+    at the critical point, but it can easily be fit so as to produce negative 
+    values at any reduced temperature.
 
     Examples
     --------
@@ -475,6 +488,8 @@ def ISTExpansion(T, Tc, a1, a2, a3=0.0, a4=0.0, a5=0.0):
        https://trc.nist.gov/TDE/Help/TDE103b/Eqns-Pure-SurfaceTension/ISTExpansion-SurfaceTension.htm
     '''
     tau = 1.0 - T/Tc
+    if tau <= 0.0:
+        return 0.0
     return tau*(a1 + tau*(a2 + tau*(a3 + tau*(a4 + a5*tau))))
 
 def Somayajulu(T, Tc, A, B, C):

tests/test_interface.py

@@ -207,13 +207,29 @@ def test_Watson_sigma():
     sigma = Watson_sigma(T=350.0, Tc=543.836, a1=-3.02417, a2=1.21792, a3=-5.26877e-9, a4=5.62659e-9, a5=-2.27553e-9)
     assert_close(sigma, 0.013834092660564925, rtol=1e-14)
 
+    # Test a couple calls very near to Tc
+    Watson_sigma(T=543.835999999, Tc=543.836, a1=-3.02417, a2=1.21792, a3=-5.26877e-9, a4=5.62659e-9, a5=-2.27553e-9)
+    Watson_sigma(T=543.839999999999, Tc=543.836, a1=-3.02417, a2=1.21792, a3=-5.26877e-9, a4=5.62659e-9, a5=-2.27553e-9)
+
+    assert 0 == Watson_sigma(T=543.836, Tc=543.836, a1=-3.02417, a2=1.21792, a3=-5.26877e-9, a4=5.62659e-9, a5=-2.27553e-9)
+    assert 0 == Watson_sigma(T=600, Tc=543.836, a1=-3.02417, a2=1.21792, a3=-5.26877e-9, a4=5.62659e-9, a5=-2.27553e-9)
+
+    # This call will give very large numbers
+    Watson_sigma(T=543.836*(1-1e-10), Tc=543.836, a1=-3.02417, a2=-1.21792, a3=-5.26877e-9, a4=5.62659e-9, a5=-2.27553e-9)
+
+
 def test_ISTExpansion():
     sigma = ISTExpansion(T=400.0, Tc=776.0, a1=0.037545, a2=0.0363288)
     assert_close(sigma, 0.02672100905515996, rtol=1e-13)
 
     sigma = ISTExpansion(T=400.0, Tc=776.0, a1=0.037545, a2=0.0363288, a3=1e-4, a4=1e-3, a5=1e-4)
     assert_close(sigma, 0.02679017489704166, rtol=1e-15)
 
+    assert 0 == ISTExpansion(T=777.0, Tc=776.0, a1=0.037545, a2=0.0363288)
+
+    # Point near Tc but with really high values
+    ISTExpansion(T=775.999999, Tc=776.0, a1=1e20, a2=0.0363288)
+
 
 def test_sigma_Gharagheizi_1():
     # point from article supporting material exactly matches