diff --git a/pyiron_atomistics/sphinx/base.py b/pyiron_atomistics/sphinx/base.py
index dc74ef31f..930c36f4c 100644
--- a/pyiron_atomistics/sphinx/base.py
+++ b/pyiron_atomistics/sphinx/base.py
@@ -373,6 +373,7 @@ def load_default_input(self):
         self.input.KpointCoords = [0.5, 0.5, 0.5]
         self.input.KpointFolding = [4, 4, 4]
         self.input.EmptyStates = "auto"
+        self.input.MethfesselPaxton = 1
         self.input.Sigma = 0.2
         self.input.Xcorr = "PBE"
         self.input.VaspPot = False
@@ -522,10 +523,14 @@ def load_hamilton_group(self):
         """
         self.input.sphinx.PAWHamiltonian.setdefault(
             "nEmptyStates", self.input["EmptyStates"]
-            )
+        )
         self.input.sphinx.PAWHamiltonian.setdefault(
             "ekt", self.input["Sigma"]
-            )
+        )
+        for k in ['MethfesselPaxton', 'FermiDirac']:
+            if k in self.input.list_nodes():
+                self.input.sphinx.PAWHamiltonian.setdefault(k, self.input[k])
+                break
         self.input.sphinx.PAWHamiltonian.setdefault("xc", self.input["Xcorr"])
         self.input.sphinx.PAWHamiltonian["spinPolarized"] = self._spin_enabled
 
@@ -965,20 +970,27 @@ def set_mixing_parameters(
         + set_mixing_parameters.__doc__
     )
 
-    def set_occupancy_smearing(self, smearing=None, width=None):
+    def set_occupancy_smearing(self, smearing=None, width=None, order=1):
         """
         Set how the finite temperature smearing is applied in
         determining partial occupancies
 
         Args:
-            smearing (str): Type of smearing (only fermi is
-                            implemented anything else will be ignored)
+            smearing (str): Type of smearing, 'FermiDirac' or 'MethfesselPaxton'
             width (float): Smearing width (eV) (default: 0.2)
-        """
-        if smearing is not None and not isinstance(smearing, str):
-            raise ValueError(
-                "Smearing must be a string"
-            )
+            order (int): Smearing order
+        """
+        if smearing is not None:
+            if not isinstance(smearing, str):
+                raise ValueError("Smearing must be a string")
+            if smearing.lower().startswith('meth'):
+                self.input.MethfesselPaxton = order
+                if 'FermiDirac' in self.input.list_nodes():
+                    del self.input['FermiDirac']
+            elif smearing.lower().startswith('fermi'):
+                self.input.FermiDirac = order
+                if 'MethfesselPaxton' in self.input.list_nodes():
+                    del self.input['MethfesselPaxton']
         if width is not None and width < 0:
             raise ValueError("Smearing value must be a float >= 0")
         if width is not None:
diff --git a/tests/sphinx/test_base.py b/tests/sphinx/test_base.py
index 56782b84c..939dfadbe 100644
--- a/tests/sphinx/test_base.py
+++ b/tests/sphinx/test_base.py
@@ -168,6 +168,7 @@ def test_write_input(self):
             'PAWHamiltonian {\n',
             '\tnEmptyStates = 6;\n',
             '\tekt = 0.2;\n',
+            '\tMethfesselPaxton = 1;\n',
             '\txc = PBE;\n',
             '\tspinPolarized;\n',
             '}\n',
@@ -386,6 +387,9 @@ def test_set_occupancy_smearing(self):
             ValueError, self.sphinx_band_structure.set_occupancy_smearing, "fermi", -0.1
         )
         self.sphinx_band_structure.set_occupancy_smearing("fermi", 0.1)
+        self.assertTrue('FermiDirac' in self.sphinx_band_structure.input)
+        self.sphinx_band_structure.set_occupancy_smearing("methfessel", 0.1)
+        self.assertTrue('MethfesselPaxton' in self.sphinx_band_structure.input)
 
     def test_load_default_groups(self):
         backup  = self.sphinx_band_structure.structure.copy()