Commit cb4f01ca authored by Lauri Himanen's avatar Lauri Himanen
Browse files

Fixed issues with DOS normalization and added extended description to the...

Fixed issues with DOS normalization and added extended description to the metainfo description about the normalization procedure.
parent 42f21688
Pipeline #78948 passed with stages
in 23 minutes and 24 seconds
...@@ -239,8 +239,8 @@ normalize = NomadConfig( ...@@ -239,8 +239,8 @@ normalize = NomadConfig(
# The threshold for point equality in k-space. Unit: 1/m. # The threshold for point equality in k-space. Unit: 1/m.
k_space_precision=150e6, k_space_precision=150e6,
# The energy threshold for how much a band can be on top or below the fermi # 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 # level in order to detect a gap. Unit: Joule.
band_structure_energy_tolerance=300 * 1.38064852E-23, band_structure_energy_tolerance=1.6022e-20, # 0.1 eV
springer_db_path=os.path.join( springer_db_path=os.path.join(
os.path.dirname(os.path.abspath(__file__)), os.path.dirname(os.path.abspath(__file__)),
'normalizing/data/springer.msg' 'normalizing/data/springer.msg'
......
...@@ -1264,10 +1264,17 @@ class section_dos(MSection): ...@@ -1264,10 +1264,17 @@ class section_dos(MSection):
shape=['number_of_dos_values'], shape=['number_of_dos_values'],
unit='joule', unit='joule',
description=''' description='''
Array containing the set of discrete energy values with respect to the top of the Array containing the set of discrete energy values with respect to the
valence band for the density (electronic-energy) of states (DOS). This is the highest occupied energy level. This is the total DOS, see
total DOS, see atom_projected_dos_energies and species_projected_dos_energies for atom_projected_dos_energies and species_projected_dos_energies for
partial density of states. 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')) a_legacy=LegacyDefinition(name='dos_energies_normalized'))
......
...@@ -47,7 +47,8 @@ class DosNormalizer(Normalizer): ...@@ -47,7 +47,8 @@ class DosNormalizer(Normalizer):
dos_values = dos.dos_values dos_values = dos.dos_values
dos_energies = dos.dos_energies dos_energies = dos.dos_energies
energy_reference_fermi = scc.energy_reference_fermi 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 continue
# Normalize DOS values to be 1/J/atom/m^3 # Normalize DOS values to be 1/J/atom/m^3
...@@ -68,75 +69,90 @@ class DosNormalizer(Normalizer): ...@@ -68,75 +69,90 @@ class DosNormalizer(Normalizer):
dos_values_normalized = dos_values / (number_of_atoms * unit_cell_volume) dos_values_normalized = dos_values / (number_of_atoms * unit_cell_volume)
# Normalize energies so that they are normalized to HOMO. # Normalize energies so that they are normalized to HOMO.
i_channel = 0 dos_energies_normalized = None
fermi_energy = energy_reference_fermi[i_channel]
fermi_idx = (np.abs(dos_energies - fermi_energy)).argmin() # Primarily normalize to the HOMO level reported by the parser.
energy_threshold = config.normalize.band_structure_energy_tolerance # The first channel is used.
value_threshold = 1e-8 # The DOS value that is considered to be zero energy_reference = None
homo_found = False if energy_reference_highest_occupied is not None:
zero_found = False energy_reference = energy_reference_highest_occupied[0]
energy_reference = fermi_energy # If HOMO is not reported, search for it in the DOS values
else:
# Walk through the energies in descencing direction to see if a i_channel = 0
# gap is nearby (see energy_threshold). If gap found, continue fermi_energy = energy_reference_fermi[i_channel]
# until HOMO found fermi_idx = (np.abs(dos_energies - fermi_energy)).argmin()
idx = fermi_idx energy_threshold = config.normalize.band_structure_energy_tolerance
while True: value_threshold = 1e-8 # The DOS value that is considered to be zero
try: homo_found = False
value = dos_values_normalized[i_channel, idx] zero_found = False
energy_distance = fermi_energy - dos_energies[idx] energy_reference = fermi_energy
except IndexError:
break # First check that the closest dos energy to fermi_energy is not too
if energy_distance.magnitude > energy_threshold and not zero_found: # far away. If it is very far away, the normalization may be very
break # inaccurate and we do not report it.
if value <= value_threshold: fermi_energy_closest = dos_energies[fermi_idx]
zero_found = True distance = np.abs(fermi_energy_closest - fermi_energy)
if zero_found and value > value_threshold: if distance.magnitude <= energy_threshold:
energy_reference = dos_energies[idx + 1]
homo_found = True # Walk through the energies in descencing direction to see if a
break # gap is nearby (see energy_threshold). If gap found, continue
idx -= 1 # until HOMO found
# If gap was not found in descending direction, check the idx = fermi_idx
# ascending direction for a nearby (see energy_threshold) HOMO while True:
# value try:
if not homo_found: value = dos_values_normalized[i_channel, idx]
idx = fermi_idx + 1 energy_distance = fermi_energy_closest - dos_energies[idx]
while True: except IndexError:
try: break
value = dos_values_normalized[i_channel, idx] if energy_distance.magnitude > energy_threshold and not zero_found:
energy_distance = dos_energies[idx] - fermi_energy break
except IndexError: if value <= value_threshold:
break zero_found = True
if energy_distance.magnitude > energy_threshold: if zero_found and value > value_threshold:
break energy_reference = dos_energies[idx + 1]
if value <= value_threshold: homo_found = True
energy_reference = dos_energies[idx] break
break idx -= 1
idx += 1 # If gap was not found in descending direction, check the
# ascending direction for a nearby (see energy_threshold) HOMO
dos_energies_normalized = dos_energies - energy_reference # value
if not homo_found:
# Data for DOS fingerprint idx = fermi_idx + 1
dos_fingerprint = None while True:
try: try:
dos_fingerprint = DOSFingerprint().calculate( value = dos_values_normalized[i_channel, idx]
dos_energies_normalized.magnitude, energy_distance = dos_energies[idx] - fermi_energy
dos_values_normalized, except IndexError:
n_atoms=number_of_atoms break
) if energy_distance.magnitude > energy_threshold:
except Exception as e: break
self.logger.error('could not generate dos fingerprint', exc_info=e) if value <= value_threshold:
energy_reference = dos_energies[idx]
break
idx += 1
# Add quantities to NOMAD's Metainfo # Add quantities to NOMAD's Metainfo
scc_url = '/section_run/0/section_single_configuration_calculation/%d/section_dos/0' % scc.m_parent_index scc_url = '/section_run/0/section_single_configuration_calculation/%d/section_dos/0' % scc.m_parent_index
self._backend.openContext(scc_url) self._backend.openContext(scc_url)
dos.dos_values_normalized = dos_values_normalized dos.dos_values_normalized = dos_values_normalized
dos.dos_energies_normalized = dos_energies_normalized
if dos_fingerprint is not None: # Data for DOS fingerprint
sec_dos_fingerprint = dos.m_create(section_dos_fingerprint) if energy_reference is not None:
sec_dos_fingerprint.bins = dos_fingerprint.bins dos_energies_normalized = dos_energies - energy_reference
sec_dos_fingerprint.indices = dos_fingerprint.indices dos.dos_energies_normalized = dos_energies_normalized
sec_dos_fingerprint.stepsize = dos_fingerprint.stepsize try:
sec_dos_fingerprint.grid_id = dos_fingerprint.grid_id dos_fingerprint = DOSFingerprint().calculate(
sec_dos_fingerprint.filling_factor = dos_fingerprint.filling_factor 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) self._backend.closeContext(scc_url)
...@@ -40,19 +40,19 @@ def test_fingerprint(dos_si_vasp): ...@@ -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): # def test_dos_energies(dos_si_vasp: Backend, dos_si_exciting: Backend, dos_si_fhiaims: Backend):
# """For debugging. # """For debugging.
# """ # """
# x_exciting = dos_si_exciting.get_value('dos_energies_normalized', 0) # x_exciting = dos_si_exciting.get_value('dos_energies_normalized', 0)
# y_exciting = dos_si_exciting.get_value('dos_values_normalized', 0) # y_exciting = dos_si_exciting.get_value('dos_values_normalized', 0)
# x_vasp = dos_si_vasp.get_value('dos_energies_normalized', 0) # x_vasp = dos_si_vasp.get_value('dos_energies_normalized', 0)
# y_vasp = dos_si_vasp.get_value('dos_values_normalized', 0) # y_vasp = dos_si_vasp.get_value('dos_values_normalized', 0)
# x_fhiaims = dos_si_fhiaims.get_value('dos_energies_normalized', 0) # x_fhiaims = dos_si_fhiaims.get_value('dos_energies_normalized', 0)
# y_fhiaims = dos_si_fhiaims.get_value('dos_values_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_vasp, y_vasp[0], label="VASP", marker=".")
# mpl.plot(x_exciting, y_exciting[0], label="exciting") # mpl.plot(x_exciting, y_exciting[0], label="exciting", marker=".")
# mpl.plot(x_fhiaims, y_fhiaims[0], label="FHI-aims") # mpl.plot(x_fhiaims, y_fhiaims[0], label="FHI-aims", marker=".")
# mpl.legend() # mpl.legend()
# mpl.show() # mpl.show()
def test_dos_magnitude(dos_si_vasp: Backend, dos_si_exciting: Backend, dos_si_fhiaims: Backend): def test_dos_magnitude(dos_si_vasp: Backend, dos_si_exciting: Backend, dos_si_fhiaims: Backend):
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment