Commit 2d5a84b1 authored by Lauri Himanen's avatar Lauri Himanen
Browse files

Added tests for band structure normalization.

parent e1d57e90
Pipeline #73541 canceled with stages
in 1 minute and 49 seconds
......@@ -41,7 +41,7 @@ def get_summed_atomic_mass(atomic_numbers: np.ndarray) -> float:
return mass
def find_match(self, pos: np.array, positions: np.array, eps: float) -> Union[int, None]:
def find_match(pos: np.array, positions: np.array, eps: float) -> Union[int, None]:
"""Attempts to find a position within a larger list of positions.
Args:
......
......@@ -20,7 +20,7 @@ default all constants are reported in SI units.
import numpy as np
import scipy.constants
atomic_mass_constant = scipy.constants.physical_constants["atomic mass constant"]
atomic_mass_constant = scipy.constants.physical_constants["atomic mass constant"][0]
pi = scipy.constants.pi
# List of atomic masses (natural isotope dist.) in order, atomic mass units.
......
......@@ -4337,11 +4337,6 @@ class section_single_configuration_calculation(MSection):
repeats=True,
a_legacy=LegacyDefinition(name='section_energy_van_der_Waals'))
section_k_band_normalized = SubSection(
sub_section=SectionProxy('section_k_band_normalized'),
repeats=True,
a_legacy=LegacyDefinition(name='section_k_band_normalized'))
section_k_band = SubSection(
sub_section=SectionProxy('section_k_band'),
repeats=True,
......
......@@ -35,11 +35,12 @@ There is one ABC for all normalizer:
from typing import List, Any, Iterable, Type
from .dos import DosNormalizer
from .system import SystemNormalizer
from .optimade import OptimadeNormalizer
from .fhiaims import FhiAimsBaseNormalizer
from .dos import DosNormalizer
from .normalizer import Normalizer
from .optimade import OptimadeNormalizer
from .system import SystemNormalizer
from .band_structure import BandStructureNormalizer
from .encyclopedia.encyclopedia import EncyclopediaNormalizer
normalizers: Iterable[Type[Normalizer]] = [
......@@ -47,5 +48,6 @@ normalizers: Iterable[Type[Normalizer]] = [
OptimadeNormalizer,
FhiAimsBaseNormalizer,
DosNormalizer,
EncyclopediaNormalizer
BandStructureNormalizer,
EncyclopediaNormalizer,
]
......@@ -17,7 +17,7 @@ import json
import numpy as np
from typing import List
from nomad.datamodel.metainfo.public import section_k_band, section_single_configuration_calculation, section_band_gap, section_system
from nomad.datamodel.metainfo.public import section_k_band, section_band_gap, section_system
from nomad.normalizing.normalizer import Normalizer
from nomad.constants import pi
from nomad import config, atomutils
......@@ -30,9 +30,6 @@ class BandStructureNormalizer(Normalizer):
- Creates labels for special points within the band path (band_path_labels).
- Determines if the path is a standard one or not (is_standard)
"""
def __init__(self):
pass
def normalize(self, logger=None) -> None:
# Setup logger
if logger is not None:
......@@ -51,14 +48,15 @@ class BandStructureNormalizer(Normalizer):
# In order to resolve the special points and the reciprocal cell,
# we need informatoin about the system.
system = scc.section_single_configuration_calculation_to_system_ref
system = scc.single_configuration_calculation_to_system_ref
for band in scc.section_k_band:
self.add_reciprocal_cell(band, system)
self.add_brillouin_zone(band)
self.add_band_gaps(band, valence_band_maximum)
self.add_path_labels(band, system)
self.add_is_standard(band)
if band.band_structure_kind != "vibrational":
self.add_reciprocal_cell(band, system)
self.add_brillouin_zone(band)
self.add_band_gaps(band, valence_band_maximum)
self.add_path_labels(band, system)
# self.add_is_standard(band)
def add_reciprocal_cell(self, band: section_k_band, system: section_system):
"""A reciprocal cell for this calculation. If the original unit cell is
......@@ -286,7 +284,7 @@ class BandStructureNormalizer(Normalizer):
start_point = segment.band_k_points[0]
end_point = segment.band_k_points[-1]
start_index = atomutils.find_match(start_point, special_k_points, eps)
end_index = atomutils.find_match(end_point, special_points, eps)
end_index = atomutils.find_match(end_point, special_k_points, eps)
if start_index is None:
start_label = None
else:
......@@ -297,101 +295,101 @@ class BandStructureNormalizer(Normalizer):
end_label = special_point_labels[end_index]
segment.band_path_labels = [start_label, end_label]
def add_is_standard(self, band: section_k_band) -> None:
pass
def crystal_structure_from_cell(self, cell: ase.cell.Cell, eps=1e-4):
"""Return the crystal structure as a string calculated from the cell.
"""
cellpar = ase.geometry.cell_to_cellpar(cell=cell)
abc = cellpar[:3]
angles = cellpar[3:] / 180 * pi
a, b, c = abc
alpha, _, gamma = angles
# According to:
# https://www.physics-in-a-nutshell.com/article/6/symmetry-crystal-systems-and-bravais-lattices#the-seven-crystal-systems
# If a=b=c and alpha=beta=gamma=90degrees we have cubic.
if abc.ptp() < eps and abs(angles - pi / 2).max() < eps:
return 'cubic'
elif abc.ptp() < eps and abs(angles - pi / 3).max() < eps:
return 'fcc'
elif abc.ptp() < eps and abs(angles - np.arccos(-1 / 3)).max() < eps:
return 'bcc'
# If a=b!=c, alpha=beta=gamma=90deg, tetragonal.
elif abs(a - b) < eps and abs(angles - pi / 2).max() < eps:
return 'tetragonal'
elif abs(angles - pi / 2).max() < eps:
return 'orthorhombic'
# if a = b != c , alpha = beta = 90deg, gamma = 120deg, hexagonal
elif (abs(a - b) < eps and abs(gamma - pi / 3 * 2) < eps and abs(angles[:2] - pi / 2).max() < eps):
return 'hexagonal'
elif (c >= a and c >= b and alpha < pi / 2 and abs(angles[1:] - pi / 2).max() < eps):
return 'monoclinic'
else:
raise ValueError('Cannot find crystal structure')
def get_special_points(self, cell, eps=1e-4):
"""Return dict of special points.
The definitions are from a paper by Wahyu Setyawana and Stefano
Curtarolo::
http://dx.doi.org/10.1016/j.commatsci.2010.05.010
lattice: str
One of the following: cubic, fcc, bcc, orthorhombic, tetragonal,
hexagonal or monoclinic.
cell: 3x3 ndarray
Unit cell.
eps: float
Tolerance for cell-check.
"""
lattice = self.crystal_structure_from_cell(cell)
cellpar = ase.geometry.cell_to_cellpar(cell=cell)
abc = cellpar[:3]
angles = cellpar[3:] / 180 * pi
a, b, c = abc
alpha, _, gamma = angles
# Check that the unit-cells are as in the Setyawana-Curtarolo paper:
if lattice == 'cubic':
assert abc.ptp() < eps and abs(angles - pi / 2).max() < eps
elif lattice == 'fcc':
assert abc.ptp() < eps and abs(angles - pi / 3).max() < eps
elif lattice == 'bcc':
angle = np.arccos(-1 / 3)
assert abc.ptp() < eps and abs(angles - angle).max() < eps
elif lattice == 'tetragonal':
assert abs(a - b) < eps and abs(angles - pi / 2).max() < eps
elif lattice == 'orthorhombic':
assert abs(angles - pi / 2).max() < eps
elif lattice == 'hexagonal':
assert abs(a - b) < eps
assert abs(gamma - pi / 3 * 2) < eps
assert abs(angles[:2] - pi / 2).max() < eps
elif lattice == 'monoclinic':
assert c >= a and c >= b
assert alpha < pi / 2 and abs(angles[1:] - pi / 2).max() < eps
if lattice != 'monoclinic':
return special_points[lattice]
# Here, we need the cell:
eta = (1 - b * np.cos(alpha) / c) / (2 * np.sin(alpha)**2)
nu = 1 / 2 - eta * c * np.cos(alpha) / b
return {'Γ': [0, 0, 0],
'A': [1 / 2, 1 / 2, 0],
'C': [0, 1 / 2, 1 / 2],
'D': [1 / 2, 0, 1 / 2],
'D1': [1 / 2, 0, -1 / 2],
'E': [1 / 2, 1 / 2, 1 / 2],
'H': [0, eta, 1 - nu],
'H1': [0, 1 - eta, nu],
'H2': [0, eta, -nu],
'M': [1 / 2, eta, 1 - nu],
'M1': [1 / 2, 1 - eta, nu],
'M2': [1 / 2, eta, -nu],
'X': [0, 1 / 2, 0],
'Y': [0, 0, 1 / 2],
'Y1': [0, 0, -1 / 2],
'Z': [1 / 2, 0, 0]}
# def add_is_standard(self, band: section_k_band) -> None:
# pass
# def crystal_structure_from_cell(self, cell: ase.cell.Cell, eps=1e-4):
# """Return the crystal structure as a string calculated from the cell.
# """
# cellpar = ase.geometry.cell_to_cellpar(cell=cell)
# abc = cellpar[:3]
# angles = cellpar[3:] / 180 * pi
# a, b, c = abc
# alpha, _, gamma = angles
# # According to:
# # https://www.physics-in-a-nutshell.com/article/6/symmetry-crystal-systems-and-bravais-lattices#the-seven-crystal-systems
# # If a=b=c and alpha=beta=gamma=90degrees we have cubic.
# if abc.ptp() < eps and abs(angles - pi / 2).max() < eps:
# return 'cubic'
# elif abc.ptp() < eps and abs(angles - pi / 3).max() < eps:
# return 'fcc'
# elif abc.ptp() < eps and abs(angles - np.arccos(-1 / 3)).max() < eps:
# return 'bcc'
# # If a=b!=c, alpha=beta=gamma=90deg, tetragonal.
# elif abs(a - b) < eps and abs(angles - pi / 2).max() < eps:
# return 'tetragonal'
# elif abs(angles - pi / 2).max() < eps:
# return 'orthorhombic'
# # if a = b != c , alpha = beta = 90deg, gamma = 120deg, hexagonal
# elif (abs(a - b) < eps and abs(gamma - pi / 3 * 2) < eps and abs(angles[:2] - pi / 2).max() < eps):
# return 'hexagonal'
# elif (c >= a and c >= b and alpha < pi / 2 and abs(angles[1:] - pi / 2).max() < eps):
# return 'monoclinic'
# else:
# raise ValueError('Cannot find crystal structure')
# def get_special_points(self, cell, eps=1e-4):
# """Return dict of special points.
# The definitions are from a paper by Wahyu Setyawana and Stefano
# Curtarolo::
# http://dx.doi.org/10.1016/j.commatsci.2010.05.010
# lattice: str
# One of the following: cubic, fcc, bcc, orthorhombic, tetragonal,
# hexagonal or monoclinic.
# cell: 3x3 ndarray
# Unit cell.
# eps: float
# Tolerance for cell-check.
# """
# lattice = self.crystal_structure_from_cell(cell)
# cellpar = ase.geometry.cell_to_cellpar(cell=cell)
# abc = cellpar[:3]
# angles = cellpar[3:] / 180 * pi
# a, b, c = abc
# alpha, _, gamma = angles
# # Check that the unit-cells are as in the Setyawana-Curtarolo paper:
# if lattice == 'cubic':
# assert abc.ptp() < eps and abs(angles - pi / 2).max() < eps
# elif lattice == 'fcc':
# assert abc.ptp() < eps and abs(angles - pi / 3).max() < eps
# elif lattice == 'bcc':
# angle = np.arccos(-1 / 3)
# assert abc.ptp() < eps and abs(angles - angle).max() < eps
# elif lattice == 'tetragonal':
# assert abs(a - b) < eps and abs(angles - pi / 2).max() < eps
# elif lattice == 'orthorhombic':
# assert abs(angles - pi / 2).max() < eps
# elif lattice == 'hexagonal':
# assert abs(a - b) < eps
# assert abs(gamma - pi / 3 * 2) < eps
# assert abs(angles[:2] - pi / 2).max() < eps
# elif lattice == 'monoclinic':
# assert c >= a and c >= b
# assert alpha < pi / 2 and abs(angles[1:] - pi / 2).max() < eps
# if lattice != 'monoclinic':
# return special_points[lattice]
# # Here, we need the cell:
# eta = (1 - b * np.cos(alpha) / c) / (2 * np.sin(alpha)**2)
# nu = 1 / 2 - eta * c * np.cos(alpha) / b
# return {'Γ': [0, 0, 0],
# 'A': [1 / 2, 1 / 2, 0],
# 'C': [0, 1 / 2, 1 / 2],
# 'D': [1 / 2, 0, 1 / 2],
# 'D1': [1 / 2, 0, -1 / 2],
# 'E': [1 / 2, 1 / 2, 1 / 2],
# 'H': [0, eta, 1 - nu],
# 'H1': [0, 1 - eta, nu],
# 'H2': [0, eta, -nu],
# 'M': [1 / 2, eta, 1 - nu],
# 'M1': [1 / 2, 1 - eta, nu],
# 'M2': [1 / 2, eta, -nu],
# 'X': [0, 1 / 2, 0],
# 'Y': [0, 0, 1 / 2],
# 'Y1': [0, 0, -1 / 2],
# 'Z': [1 / 2, 0, 0]}
......@@ -223,3 +223,10 @@ def hash_exciting() -> Backend:
@pytest.fixture(scope='session')
def hash_vasp(bands_unpolarized_gap_indirect) -> Backend:
return bands_unpolarized_gap_indirect
@pytest.fixture(scope='session')
def band_path_fcc(bands_unpolarized_gap_indirect) -> Backend:
"""Band structure calculation for a cubic lattice.
"""
return bands_unpolarized_gap_indirect
# Copyright 2018 Markus Scheidgen
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an"AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
import numpy as np
import pytest
from tests.normalizing.conftest import ( # pylint: disable=unused-import
phonon,
bands_unpolarized_no_gap,
bands_polarized_no_gap,
bands_unpolarized_gap_indirect,
bands_polarized_gap_indirect,
band_path_fcc,
)
from pint import UnitRegistry
ureg = UnitRegistry()
def test_band_gaps(bands_unpolarized_no_gap, bands_polarized_no_gap, bands_unpolarized_gap_indirect, bands_polarized_gap_indirect):
"""Tests that band gaps are correctly identified for different cases.
"""
def test_generic(bs, n_channels):
"""Generic tests for band structure data."""
assert bs.brillouin_zone is not None
assert bs.reciprocal_cell.shape == (3, 3)
# Unpolarized, no gaps
bs = bands_unpolarized_no_gap.entry_archive.section_run[0].section_single_configuration_calculation[0].section_k_band[0]
test_generic(bs, n_channels=1)
assert bs.section_band_gap is None
assert bs.section_band_gap_spin_up is None
assert bs.section_band_gap_spin_down is None
# Polarized, no gaps
bs = bands_polarized_no_gap.entry_archive.section_run[0].section_single_configuration_calculation[0].section_k_band[0]
test_generic(bs, n_channels=2)
assert bs.section_band_gap is None
assert bs.section_band_gap_spin_up is None
assert bs.section_band_gap_spin_down is None
# Unpolarized, finite gap, indirect
bs = bands_unpolarized_gap_indirect.entry_archive.section_run[0].section_single_configuration_calculation[0].section_k_band[0]
test_generic(bs, n_channels=1)
gap_ev = (bs.section_band_gap.value * ureg.J).to(ureg.eV).magnitude
assert gap_ev == pytest.approx(0.62, 0.01)
assert bs.section_band_gap.type == "indirect"
assert bs.section_band_gap_spin_up is None
assert bs.section_band_gap_spin_down is None
# Polarized, finite gap, indirect
bs = bands_polarized_gap_indirect.entry_archive.section_run[0].section_single_configuration_calculation[0].section_k_band[0]
test_generic(bs, n_channels=2)
gap = bs.section_band_gap
gap_up = bs.section_band_gap_spin_up
gap_down = bs.section_band_gap_spin_down
gap_ev = (gap.value * ureg.J).to(ureg.eV).magnitude
gap_up_ev = (gap_up.value * ureg.J).to(ureg.eV).magnitude
gap_down_ev = (gap_down.value * ureg.J).to(ureg.eV).magnitude
assert gap_up.type == "indirect"
assert gap_down.type == "indirect"
assert gap_up_ev != gap_down_ev
assert gap_up_ev == gap_ev
assert gap_up_ev == pytest.approx(0.956, 0.01)
assert gap_down_ev == pytest.approx(1.230, 0.01)
def test_band_paths(band_path_fcc):
"""Tests that the path labeling workds for different lattices.
"""
# FCC
scc = band_path_cubic.entry_archive.section_run[0].section_single_configuration_calculation[0]
system = scc.single_configuration_calculation_to_system_ref
space_group_number = system.section_symmetry[0].space_group_number
assert space_group_number == 227
assumed_path = np.array([
["G", "X"],
["X", "W"],
["W", "K"],
["K", "G"],
["G", "L"],
["L", "U"],
["U", "W"],
["W", "L"],
["L", "K"],
["U", "X"],
])
bs = band_path_cubic.entry_archive.section_run[0].section_single_configuration_calculation[0].section_k_band[0]
for i, segment in enumerate(bs.section_k_band_segment):
labels = segment.band_path_labels
assert np.array_equal(labels, assumed_path[i, :])
def test_phonon_band(phonon):
"""Ensures that phonon bands are not touched.
"""
bs = phonon.entry_archive.section_run[0].section_single_configuration_calculation[0].section_k_band[0]
assert bs.is_standard_path is None
assert bs.section_band_gap is None
assert bs.section_band_gap_spin_up is None
assert bs.section_band_gap_spin_down is None
......@@ -29,10 +29,6 @@ from tests.normalizing.conftest import ( # pylint: disable=unused-import
phonon,
two_d,
bulk,
bands_unpolarized_no_gap,
bands_polarized_no_gap,
bands_unpolarized_gap_indirect,
bands_polarized_gap_indirect,
dos_unpolarized_vasp,
dos_polarized_vasp,
hash_exciting,
......@@ -454,61 +450,6 @@ def test_method_gw_metainfo(gw):
assert enc.method.gw_starting_point == "GGA_C_PBE+0.75*GGA_X_PBE+0.25*HF_X"
def test_band_structure(bands_unpolarized_no_gap, bands_polarized_no_gap, bands_unpolarized_gap_indirect, bands_polarized_gap_indirect):
def test_generic(bs, n_channels):
"""Generic tests for band structure data."""
assert bs.brillouin_zone is not None
assert bs.reciprocal_cell.shape == (3, 3)
# Unpolarized, no gaps
enc = bands_unpolarized_no_gap.entry_archive.section_encyclopedia
properties = enc.properties
bs = properties.electronic_band_structure
test_generic(bs, n_channels=1)
assert bs.band_gap is None
assert bs.band_gap_spin_up is None
assert bs.band_gap_spin_down is None
# Polarized, no gaps
enc = bands_polarized_no_gap.entry_archive.section_encyclopedia
properties = enc.properties
bs = properties.electronic_band_structure
test_generic(bs, n_channels=2)
assert bs.band_gap is None
assert bs.band_gap_spin_up is None
assert bs.band_gap_spin_down is None
# Unpolarized, finite gap, indirect
enc = bands_unpolarized_gap_indirect.entry_archive.section_encyclopedia
properties = enc.properties
bs = properties.electronic_band_structure
test_generic(bs, n_channels=1)
gap_ev = (bs.band_gap.value * ureg.J).to(ureg.eV).magnitude
assert gap_ev == pytest.approx(0.62, 0.01)
assert bs.band_gap.type == "indirect"
assert bs.band_gap_spin_up is None
assert bs.band_gap_spin_down is None
# Polarized, finite gap, indirect
enc = bands_polarized_gap_indirect.entry_archive.section_encyclopedia
properties = enc.properties
bs = properties.electronic_band_structure
test_generic(bs, n_channels=2)
gap = bs.band_gap
gap_up = bs.band_gap_spin_up
gap_down = bs.band_gap_spin_down
gap_ev = (gap.value * ureg.J).to(ureg.eV).magnitude
gap_up_ev = (gap_up.value * ureg.J).to(ureg.eV).magnitude
gap_down_ev = (gap_down.value * ureg.J).to(ureg.eV).magnitude
assert gap_up.type == "indirect"
assert gap_down.type == "indirect"
assert gap_up_ev != gap_down_ev
assert gap_up_ev == gap_ev
assert gap_up_ev == pytest.approx(0.956, 0.01)
assert gap_down_ev == pytest.approx(1.230, 0.01)
def test_hashes_exciting(hash_exciting):
"""Tests that the hashes has been successfully created for calculations
from exciting.
......
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