Fix sfd nan at large negative x from log(0) (#347)

sfd(x) computes the Fermi-Dirac entropy -f*log(f) - (1-f)*log(1-f).
For x <= -39, expm1(x) rounds to -1.0 in float64, so f -> 1.0 and
1-f -> 0.0, producing log(0) -> nan and RuntimeWarnings (divide by zero
in log, invalid value in multiply). The physical entropy limit there is
0 (x*log(x) -> 0 as x -> 0).

Mask the f=0/1 boundaries to 0 instead of computing the log. Interior
values are unchanged. This was a real numerical bug, not just noise:
test_band and other tests consumed nan-valued entropy from sfd.

Add a regression test covering x in [-100, 100] asserting finiteness and
the vanishing boundary. Drop the NaN-tolerating filter in
test_sfd_properties since sfd no longer returns nan.

Author: QG-phy

dptb/tests/test_occupy.py

@@ -161,10 +161,29 @@ def test_sfd_properties(self):
         x = np.linspace(-5, 5, 100)
         s = sfd(x)
 
-        # Entropy should be non-negative (handling log singularities)
-        # Filter out NaN values which can occur at boundaries
-        s_valid = s[np.isfinite(s)]
-        assert np.all(s_valid >= 0)
+        # Entropy should be non-negative
+        assert np.all(s >= 0)
+
+        # Maximum at Fermi level (x=0)
+        assert np.isclose(s[50], np.max(s), atol=1e-6)
+
+    def test_sfd_finite_at_large_negative_x(self):
+        """sfd must stay finite for large negative x where expm1 rounds to -1.
+
+        For x <= -39, expm1(x) rounds to -1.0 in float64, so f -> 1 and
+        1-f -> 0. The entropy has physical limit 0 there (x*log(x) -> 0),
+        so the result must be finite 0 rather than log(0) -> nan.
+        """
+        x = np.array([-100, -50, -40, -39, -20, 0, 20, 40, 50, 100])
+        s = sfd(x)
+        assert np.all(np.isfinite(s))
+        # Extreme values and negative round-to-one boundaries vanish exactly.
+        assert np.all(s[[0, 1, 2, 3, 8, 9]] == 0.0)
+        # x=40 is inside the compute range; f is tiny but nonzero, so entropy
+        # is a small positive number rather than 0
+        assert s[7] > 0 and s[7] < 1e-10
+        # interior unchanged: x=0 is the entropy maximum ln(2)
+        assert np.isclose(s[5], np.log(2.0), rtol=1e-12)
 
     def test_smp1_zero(self):
         """Test that MP1 entropy is zero (by design)."""

dptb/utils/occupy.py

@@ -147,7 +147,18 @@ def sfd(x):
     x_valid = x[mask_in_limit]
     f_valid = 1.0 / (np.expm1(x_valid) + 2.0)  # stable ffd
     f1_valid = 1.0 - f_valid
-    out[mask_in_limit] = -f_valid * np.log(f_valid) - f1_valid * np.log(f1_valid)
+    # Guard the f=0/1 boundaries: for sufficiently negative x, expm1(x) rounds
+    # to -1.0 in float64, so f_valid -> 1.0 and f1_valid -> 0.0, which would
+    # produce log(0) -> nan. The physical entropy limit there is 0
+    # (x*log(x) -> 0 as x -> 0), so mask those entries to 0 instead of
+    # computing the log. The interior values are unchanged.
+    interior = (f_valid > 0.0) & (f_valid < 1.0)
+    out_valid = np.zeros_like(f_valid)
+    out_valid[interior] = (
+        -f_valid[interior] * np.log(f_valid[interior])
+        - f1_valid[interior] * np.log(f1_valid[interior])
+    )
+    out[mask_in_limit] = out_valid
     return out
 
 def fmp1(x):