Suspected Typo Fix in pymatgen.io.vasp.optics (#2989)

* Small fix for a suspected typo in handling `ismear`

* typo

* add tests for io.vasp.optics.(delta_func|step_func)

---------

Co-authored-by: Janosh Riebesell <[email protected]>

Author: kavanase

pymatgen/io/vasp/optics.py

@@ -48,9 +48,9 @@ class DielectricFunctionCalculator(MSONable):
         - Perform symmetry operations (this is not implemented here)
         - Calculate the real part
 
-    Currently, this Calculator only works for ``ISYM=0`` calculations since we cannot gauranttee that our
-    externally defined symmetry operations are the same as VASP's.  This can be fixed by printing the
-    symmetry operators into the vasprun.xml file.  If this happens in future versions of VASP,
+    Currently, this Calculator only works for ``ISYM=0`` calculations since we cannot guarantee that our
+    externally defined symmetry operations are the same as VASP's. This can be fixed by printing the
+    symmetry operators into the vasprun.xml file. If this happens in future versions of VASP,
     we can dramatically speed up the calculations here by considering only the irreducible kpoints.
     """
 
@@ -279,7 +279,7 @@ def delta_func(x, ismear):
         raise ValueError("Delta function not implemented for ismear < -1")
     if ismear == -1:
         return step_func(x, -1) * (1 - step_func(x, -1))
-    if ismear < 0:
+    if ismear == 0:
         return np.exp(-(x * x)) / np.sqrt(np.pi)
     return delta_methfessel_paxton(x, ismear)
 
@@ -290,7 +290,7 @@ def step_func(x, ismear):
         raise ValueError("Delta function not implemented for ismear < -1")
     if ismear == -1:
         return 1 / (1.0 + np.exp(-x))
-    if ismear < 0:
+    if ismear == 0:
         return 0.5 + 0.5 * scipy.special.erf(x)
     return step_methfessel_paxton(x, ismear)
 

pymatgen/io/vasp/tests/test_optics.py

@@ -1,8 +1,10 @@
 from __future__ import annotations
 
 import numpy as np
+import pytest
+import scipy.special
 
-from pymatgen.io.vasp.optics import DielectricFunctionCalculator
+from pymatgen.io.vasp.optics import DielectricFunctionCalculator, delta_func, delta_methfessel_paxton, step_func
 from pymatgen.io.vasp.outputs import Vasprun
 from pymatgen.util.testing import PymatgenTest
 
@@ -18,6 +20,10 @@ def test_optics(self):
         dfc = DielectricFunctionCalculator.from_directory(eps_data_path)
         egrid, eps = dfc.get_epsilon(0, 0)
 
+        assert egrid[0] == 0
+        assert egrid[-1] == 59.3802
+        assert len(egrid) == len(eps) == 3000
+
         _, eps_real_ref, eps_imag_ref = vrun.dielectric
         eps_real_ref = np.array(eps_real_ref)[:, 0]
         eps_imag_ref = np.array(eps_imag_ref)[:, 0]
@@ -40,3 +46,39 @@ def test_optics(self):
         mask = np.ones_like(dfc.cder, dtype=float)
         x_val, y_val, text = dfc.plot_weighted_transition_data(0, 0, mask=mask, min_val=0.001)
         assert len(x_val) == len(y_val) == len(text)
+
+
+def test_delta_func():
+    x = np.array([0, 1, 2, 3, 4, 5])
+
+    # ismear < -1
+    with pytest.raises(ValueError, match="Delta function not implemented for ismear < -1"):
+        delta_func(x, -2)
+
+    # ismear == -1
+    assert np.all(delta_func(x, -1) == step_func(x, -1) * (1 - step_func(x, -1)))
+
+    # ismear == 0
+    assert np.all(delta_func(x, 0) == np.exp(-(x * x)) / np.sqrt(np.pi))
+
+    # ismear > 0
+    for ismear in [1, 2, 3]:
+        assert np.all(delta_func(x, ismear) == delta_methfessel_paxton(x, ismear))
+
+
+def test_step_func():
+    # array of positive values
+    x = np.array([1, 2, 3, 4, 5])
+    assert np.allclose(step_func(x, -1), 1 / (1.0 + np.exp(-x)))
+
+    # array of negative values
+    x = np.array([-1, -2, -3, -4, -5])
+    assert np.allclose(step_func(x, -1), 1 / (1.0 + np.exp(-x)))
+
+    # array that includes zero
+    x = np.array([-1, 0, 1])
+    assert np.allclose(step_func(x, -1), 1 / (1.0 + np.exp(-x)))
+
+    # ismear == 0
+    x = np.array([1, 2, 3, 4, 5])
+    assert np.allclose(step_func(x, 0), 0.5 + 0.5 * scipy.special.erf(x))