Bugfix: Issue with spin-polarized calculations

* The code had several errors to do with detecting whether the calculation was
  spin-polarized or not. This has been fixed and the detection simplified.

HISTORY.rst

@@ -1,6 +1,10 @@
 =======
 History
 =======
+2025.1.21 -- Bugfix: Issue with spin-polarized calculations
+    * The code had several errors to do with detecting whether the calculation was
+      spin-polarized or not. This has been fixed and the detection simplified.
+
 2024.12.14 -- Bugfix: Issue with initialization in subflowcharts
 
 2024.10.20 -- Added the standard results for drivers

dftbplus_step/base.py

@@ -113,7 +113,13 @@ def exe_config(self):
         return self.parent.exe_config
 
     def band_structure(
-        self, input_path, sym_points, sym_names, Efermi=[0.0, 0.0], DOS=None
+        self,
+        input_path,
+        sym_points,
+        sym_names,
+        Efermi=0.0,
+        DOS=None,
+        spin_polarized=False,
     ):
         """Prepare the graph for the band structure.
 
@@ -123,7 +129,11 @@ def band_structure(
             The path to the band output from DFTB+.
         """
         Band_Structure = self.create_band_structure_data(  # noqa: F841
-            input_path, sym_points, sym_names, Efermi=Efermi
+            input_path,
+            sym_points,
+            sym_names,
+            Efermi=Efermi,
+            spin_polarized=spin_polarized,
         )
 
         figure = cms_plots.band_structure(
@@ -141,7 +151,7 @@ def band_structure(
             figure.dump(wd / "band_structure.html")
 
     def create_band_structure_data(
-        self, input_path, sym_points, sym_names, Efermi=[0.0, 0.0]
+        self, input_path, sym_points, sym_names, Efermi=0.0, spin_polarized=False
     ):
         """Massage the band structure data into a standard form
 
@@ -155,8 +165,6 @@ def create_band_structure_data(
 
         seamm_options = self.parent.global_options
 
-        spin_polarized = len(Efermi) == 2
-
         # Read configuration file for DFTB+
         ini_dir = Path(seamm_options["root"]).expanduser()
         full_config = configparser.ConfigParser()
@@ -250,7 +258,7 @@ def create_band_structure_data(
 
         return BandStructure
 
-    def create_dos_data(self, input_path, Efermi=[0.0]):
+    def create_dos_data(self, input_path, Efermi=0.0, spin_polarized=False):
         """Create the DOS data
 
         Parameters
@@ -260,13 +268,6 @@ def create_dos_data(self, input_path, Efermi=[0.0]):
         Efermi : float
             The Fermi energy in eV
         """
-        spin_polarized = len(Efermi) == 2
-
-        if spin_polarized and Efermi[0] != Efermi[1]:
-            raise NotImplementedError(
-                f"Cannot handle different Fermi energies yet: {Efermi}"
-            )
-
         logger.info("Preparing DOS")
 
         # Total DOS
@@ -444,13 +445,13 @@ def create_dos_data(self, input_path, Efermi=[0.0]):
                 DOS[element] = total.array
 
         # Shift the Fermi level to 0
-        DOS.index -= Efermi[0]
+        DOS.index -= Efermi
 
         DOS.to_csv(Path(self.directory) / "DOS.csv")
 
         return DOS
 
-    def dos(self, input_path, Efermi=[0.0]):
+    def dos(self, input_path, Efermi=0.0, spin_polarized=False):
         """Prepare the graph for the density of states.
 
         Parameters
@@ -459,8 +460,12 @@ def dos(self, input_path, Efermi=[0.0]):
             The path to the band output from DFTB+.
         Efermi : float
             The Fermi energy in eV
+        spin_polarized : bool
+            Whether the calculation is spin polarized
         """
-        DOS = self.create_dos_data(input_path, Efermi=Efermi)
+        DOS = self.create_dos_data(
+            input_path, Efermi=Efermi, spin_polarized=spin_polarized
+        )
 
         figure = cms_plots.dos(DOS, template="line.graph_template")
 
@@ -660,6 +665,11 @@ def parse_results(self, lines):
                         # The Fermi level has one or two values if spin-polarized, but
                         # normally they are the same, so turn into a scalar.
                         property_data[key] = values[0]
+                        property_data["spin polarized"] = len(values) == 2
+                        if len(values) == 2 and values[0] != values[1]:
+                            raise NotImplementedError(
+                                f"Cannot handle different Fermi energies yet: {values}"
+                            )
                     else:
                         property_data[key] = redimension(values, dims)
                 if key not in properties:
@@ -671,7 +681,7 @@ def parse_results(self, lines):
 
         # Create the standard properties needed for energy, gradients, etc.
         property_data["energy"] = property_data["total_energy"]
-        property_data["energy,units"] = "e_H"
+        property_data["energy,units"] = "E_h"
         if "forces" in property_data:
             property_data["gradients"] = [
                 [-v for v in row] for row in property_data["forces"]

dftbplus_step/energy.py

@@ -15,6 +15,7 @@
 from pathlib import Path
 import shutil
 import textwrap
+import traceback
 
 import hsd
 from tabulate import tabulate
@@ -643,11 +644,11 @@ def analyze(self, indent="", data={}, out=[]):
         # Prepare the DOS graph(s)
         if "fermi_level" in data:
             # Efermi = list(Q_(data["fermi_level"], "hartree").to("eV").magnitude)
-            Efermi = [Q_(data["fermi_level"], "hartree").to("eV").magnitude]
+            Efermi = Q_(data["fermi_level"], "hartree").to("eV").magnitude
         else:
-            Efermi = [0.0]
+            Efermi = 0.0
         wd = Path(self.directory)
-        self.dos(wd / "band.out", Efermi=Efermi)
+        self.dos(wd / "band.out", Efermi=Efermi, spin_polarized=data["spin polarized"])
 
         text_lines = []
         # Get charges and spins, etc.
@@ -733,6 +734,8 @@ def analyze(self, indent="", data={}, out=[]):
             try:
                 text += self.make_plots(data)
             except Exception as e:
+                traceback.print_exc()
+                print(f"There was an error making the plots: {str(e)}")
                 text += f"There was an error making the plots: {str(e)}"
         else:
             text += f"Orbital plots not supported for {self.model}"