Merge branch 'bs-dos' into 'v0.8.3'

Improved DOS normalization

See merge request !136

nomad/config.py

@@ -239,8 +239,8 @@ normalize = NomadConfig(
     # The threshold for point equality in k-space. Unit: 1/m.
     k_space_precision=150e6,
     # The energy threshold for how much a band can be on top or below the fermi
-    # level in order to detect a gap. k_B x T at room temperature. Unit: Joule
-    band_structure_energy_tolerance=300 * 1.38064852E-23,
+    # level in order to detect a gap. Unit: Joule.
+    band_structure_energy_tolerance=1.6022e-20,  # 0.1 eV
     springer_db_path=os.path.join(
         os.path.dirname(os.path.abspath(__file__)),
         'normalizing/data/springer.msg'

nomad/datamodel/metainfo/public.py

@@ -1264,10 +1264,17 @@ class section_dos(MSection):
         shape=['number_of_dos_values'],
         unit='joule',
         description='''
-        Array containing the set of discrete energy values with respect to the top of the
-        valence band for the density (electronic-energy) of states (DOS). This is the
-        total DOS, see atom_projected_dos_energies and species_projected_dos_energies for
+        Array containing the set of discrete energy values with respect to the
+        highest occupied energy level. This is the total DOS, see
+        atom_projected_dos_energies and species_projected_dos_energies for
         partial density of states.
+
+        If not available through energy_reference_highest_occupied, the highest
+        occupied energy level is detected by searching for a non-zero DOS value
+        below (or nearby) the reported energy_reference_fermi. In case the
+        highest occupied energy level cannot be detected accurately, the
+        normalized values are not reported. For calculations with multiple
+        spin-channels, the normalization is determined by the first channel.
         ''',
         a_legacy=LegacyDefinition(name='dos_energies_normalized'))
 

nomad/normalizing/dos.py

@@ -47,7 +47,8 @@ class DosNormalizer(Normalizer):
                 dos_values = dos.dos_values
                 dos_energies = dos.dos_energies
                 energy_reference_fermi = scc.energy_reference_fermi
-                if dos_energies is None or dos_values is None or energy_reference_fermi is None:
+                energy_reference_highest_occupied = scc.energy_reference_highest_occupied
+                if dos_energies is None or dos_values is None or (energy_reference_fermi is None and energy_reference_highest_occupied is None):
                     continue
 
                 # Normalize DOS values to be 1/J/atom/m^3
@@ -68,75 +69,90 @@ class DosNormalizer(Normalizer):
                 dos_values_normalized = dos_values / (number_of_atoms * unit_cell_volume)
 
                 # Normalize energies so that they are normalized to HOMO.
-                i_channel = 0
-                fermi_energy = energy_reference_fermi[i_channel]
-                fermi_idx = (np.abs(dos_energies - fermi_energy)).argmin()
-                energy_threshold = config.normalize.band_structure_energy_tolerance
-                value_threshold = 1e-8  # The DOS value that is considered to be zero
-                homo_found = False
-                zero_found = False
-                energy_reference = fermi_energy
-
-                # Walk through the energies in descencing direction to see if a
-                # gap is nearby (see energy_threshold). If gap found, continue
-                # until HOMO found
-                idx = fermi_idx
-                while True:
-                    try:
-                        value = dos_values_normalized[i_channel, idx]
-                        energy_distance = fermi_energy - dos_energies[idx]
-                    except IndexError:
-                        break
-                    if energy_distance.magnitude > energy_threshold and not zero_found:
-                        break
-                    if value <= value_threshold:
-                        zero_found = True
-                    if zero_found and value > value_threshold:
-                        energy_reference = dos_energies[idx + 1]
-                        homo_found = True
-                        break
-                    idx -= 1
-                # If gap was not found in descending direction, check the
-                # ascending direction for a nearby (see energy_threshold) HOMO
-                # value
-                if not homo_found:
-                    idx = fermi_idx + 1
-                    while True:
-                        try:
-                            value = dos_values_normalized[i_channel, idx]
-                            energy_distance = dos_energies[idx] - fermi_energy
-                        except IndexError:
-                            break
-                        if energy_distance.magnitude > energy_threshold:
-                            break
-                        if value <= value_threshold:
-                            energy_reference = dos_energies[idx]
-                            break
-                        idx += 1
-
-                dos_energies_normalized = dos_energies - energy_reference
-
-                # Data for DOS fingerprint
-                dos_fingerprint = None
-                try:
-                    dos_fingerprint = DOSFingerprint().calculate(
-                        dos_energies_normalized.magnitude,
-                        dos_values_normalized,
-                        n_atoms=number_of_atoms
-                    )
-                except Exception as e:
-                    self.logger.error('could not generate dos fingerprint', exc_info=e)
+                dos_energies_normalized = None
+
+                # Primarily normalize to the HOMO level reported by the parser.
+                # The first channel is used.
+                energy_reference = None
+                if energy_reference_highest_occupied is not None:
+                    energy_reference = energy_reference_highest_occupied[0]
+                # If HOMO is not reported, search for it in the DOS values
+                else:
+                    i_channel = 0
+                    fermi_energy = energy_reference_fermi[i_channel]
+                    fermi_idx = (np.abs(dos_energies - fermi_energy)).argmin()
+                    energy_threshold = config.normalize.band_structure_energy_tolerance
+                    value_threshold = 1e-8  # The DOS value that is considered to be zero
+                    homo_found = False
+                    zero_found = False
+                    energy_reference = fermi_energy
+
+                    # First check that the closest dos energy to fermi_energy is not too
+                    # far away. If it is very far away, the normalization may be very
+                    # inaccurate and we do not report it.
+                    fermi_energy_closest = dos_energies[fermi_idx]
+                    distance = np.abs(fermi_energy_closest - fermi_energy)
+                    if distance.magnitude <= energy_threshold:
+
+                        # Walk through the energies in descencing direction to see if a
+                        # gap is nearby (see energy_threshold). If gap found, continue
+                        # until HOMO found
+                        idx = fermi_idx
+                        while True:
+                            try:
+                                value = dos_values_normalized[i_channel, idx]
+                                energy_distance = fermi_energy_closest - dos_energies[idx]
+                            except IndexError:
+                                break
+                            if energy_distance.magnitude > energy_threshold and not zero_found:
+                                break
+                            if value <= value_threshold:
+                                zero_found = True
+                            if zero_found and value > value_threshold:
+                                energy_reference = dos_energies[idx + 1]
+                                homo_found = True
+                                break
+                            idx -= 1
+                        # If gap was not found in descending direction, check the
+                        # ascending direction for a nearby (see energy_threshold) HOMO
+                        # value
+                        if not homo_found:
+                            idx = fermi_idx + 1
+                            while True:
+                                try:
+                                    value = dos_values_normalized[i_channel, idx]
+                                    energy_distance = dos_energies[idx] - fermi_energy
+                                except IndexError:
+                                    break
+                                if energy_distance.magnitude > energy_threshold:
+                                    break
+                                if value <= value_threshold:
+                                    energy_reference = dos_energies[idx]
+                                    break
+                                idx += 1
 
                 # Add quantities to NOMAD's Metainfo
                 scc_url = '/section_run/0/section_single_configuration_calculation/%d/section_dos/0' % scc.m_parent_index
                 self._backend.openContext(scc_url)
                 dos.dos_values_normalized = dos_values_normalized
-                dos.dos_energies_normalized = dos_energies_normalized
-                if dos_fingerprint is not None:
-                    sec_dos_fingerprint = dos.m_create(section_dos_fingerprint)
-                    sec_dos_fingerprint.bins = dos_fingerprint.bins
-                    sec_dos_fingerprint.indices = dos_fingerprint.indices
-                    sec_dos_fingerprint.stepsize = dos_fingerprint.stepsize
-                    sec_dos_fingerprint.grid_id = dos_fingerprint.grid_id
-                    sec_dos_fingerprint.filling_factor = dos_fingerprint.filling_factor
+
+                # Data for DOS fingerprint
+                if energy_reference is not None:
+                    dos_energies_normalized = dos_energies - energy_reference
+                    dos.dos_energies_normalized = dos_energies_normalized
+                    try:
+                        dos_fingerprint = DOSFingerprint().calculate(
+                            dos_energies_normalized.magnitude,
+                            dos_values_normalized,
+                            n_atoms=number_of_atoms
+                        )
+                    except Exception as e:
+                        self.logger.error('could not generate dos fingerprint', exc_info=e)
+                    else:
+                        sec_dos_fingerprint = dos.m_create(section_dos_fingerprint)
+                        sec_dos_fingerprint.bins = dos_fingerprint.bins
+                        sec_dos_fingerprint.indices = dos_fingerprint.indices
+                        sec_dos_fingerprint.stepsize = dos_fingerprint.stepsize
+                        sec_dos_fingerprint.grid_id = dos_fingerprint.grid_id
+                        sec_dos_fingerprint.filling_factor = dos_fingerprint.filling_factor
                 self._backend.closeContext(scc_url)

tests/normalizing/test_dos.py

@@ -40,19 +40,19 @@ def test_fingerprint(dos_si_vasp):
 
 
 # def test_dos_energies(dos_si_vasp: Backend, dos_si_exciting: Backend, dos_si_fhiaims: Backend):
-#  """For debugging.
-#  """
-#  x_exciting = dos_si_exciting.get_value('dos_energies_normalized', 0)
-#  y_exciting = dos_si_exciting.get_value('dos_values_normalized', 0)
-#  x_vasp = dos_si_vasp.get_value('dos_energies_normalized', 0)
-#  y_vasp = dos_si_vasp.get_value('dos_values_normalized', 0)
-#  x_fhiaims = dos_si_fhiaims.get_value('dos_energies_normalized', 0)
-#  y_fhiaims = dos_si_fhiaims.get_value('dos_values_normalized', 0)
-#  mpl.plot(x_vasp, y_vasp[0], label="VASP")
-#  mpl.plot(x_exciting, y_exciting[0], label="exciting")
-#  mpl.plot(x_fhiaims, y_fhiaims[0], label="FHI-aims")
-#  mpl.legend()
-#  mpl.show()
+#    """For debugging.
+#    """
+#    x_exciting = dos_si_exciting.get_value('dos_energies_normalized', 0)
+#    y_exciting = dos_si_exciting.get_value('dos_values_normalized', 0)
+#    x_vasp = dos_si_vasp.get_value('dos_energies_normalized', 0)
+#    y_vasp = dos_si_vasp.get_value('dos_values_normalized', 0)
+#    x_fhiaims = dos_si_fhiaims.get_value('dos_energies_normalized', 0)
+#    y_fhiaims = dos_si_fhiaims.get_value('dos_values_normalized', 0)
+#    mpl.plot(x_vasp, y_vasp[0], label="VASP", marker=".")
+#    mpl.plot(x_exciting, y_exciting[0], label="exciting", marker=".")
+#    mpl.plot(x_fhiaims, y_fhiaims[0], label="FHI-aims", marker=".")
+#    mpl.legend()
+#    mpl.show()
 
 
 def test_dos_magnitude(dos_si_vasp: Backend, dos_si_exciting: Backend, dos_si_fhiaims: Backend):