Commit 98559c84 authored by Lauri Himanen's avatar Lauri Himanen
Browse files

Added structure information for 2D systems in EncyclopediaNormalizer, updated...

Added structure information for 2D systems in EncyclopediaNormalizer, updated SystemNormalizer so that it does not try to fetch symemtry information for non-bulk systems, added improved system type detection from MatID to SystemNormalizer.
parent 30b071d2
Pipeline #67235 failed with stages
in 13 minutes and 32 seconds
......@@ -193,7 +193,11 @@ normalize = NomadConfig(
# cluster. Used in detecting system type.
cluster_threshold=3.1,
# Defines the "bin size" for rounding cell angles for the material hash
angle_rounding=float(10.0) # unit: degree
angle_rounding=float(10.0), # unit: degree
# The threshold for a system to be considered "flat". Used e.g. when
# determining if a 2D structure is purely 2-dimensional to allow extra rigid
# transformations that are improper in 3D but proper in 2D.
flat_dim_threshold=0.1,
)
client = NomadConfig(
......
......@@ -21,8 +21,10 @@ import math
import json
import ase
import ase.data
from ase import Atoms
import numpy as np
from matid import SymmetryAnalyzer
import matid.geometry
from nomad.normalizing.normalizer import Normalizer, s_scc, s_system, s_frame_sequence, r_frame_sequence_to_sampling, s_sampling_method, r_frame_sequence_local_frames, r_scc_to_system
from nomad.metainfo.encyclopedia import Encyclopedia, Material, Calculation
......@@ -389,11 +391,11 @@ class EncyclopediaNormalizer(Normalizer):
# Fill structure related metainfo
struct = None
if system_type == Material.system_type.type.bulk:
struct = StructureBulk()
struct = StructureBulk(self._backend, self.logger)
elif system_type == Material.system_type.type.two_d:
struct = Structure2D()
struct = Structure2D(self._backend, self.logger)
if struct is not None:
struct.fill(self._backend, representative_system)
struct.fill(representative_system)
def normalize(self, logger=None) -> None:
super().normalize(logger)
......@@ -428,17 +430,16 @@ class Structure():
"""A base class that is used for processing structure related information
in the Encylopedia.
"""
def __init__(self, backend, logger):
self.backend = backend
self.logger = logger
def atom_labels(self, material: Material, std_atoms: ase.Atoms) -> None:
material.atom_labels = std_atoms.get_chemical_symbols()
def atom_positions(self, material: Material, std_atoms: ase.Atoms) -> None:
material.atom_positions = std_atoms.get_scaled_positions(wrap=False)
def atomic_density(self, calculation: Calculation, repr_system: ase.Atoms) -> None:
orig_n_atoms = len(repr_system)
orig_volume = repr_system.get_volume() * (1e-10)**3
calculation.atomic_density = float(orig_n_atoms / orig_volume)
@abstractmethod
def bravais_lattice(self, material: Material, section_system: Dict) -> None:
pass
......@@ -474,18 +475,13 @@ class Structure():
def has_free_wyckoff_parameters(self, material: Material, wyckoff_sets: Dict) -> None:
pass
@abstractmethod
def lattice_parameters(self, calculation: Calculation, std_atoms: ase.Atoms) -> None:
pass
def material_hash(self, material: Material, symmetry_analyzer: SymmetryAnalyzer) -> None:
wyckoff_sets = symmetry_analyzer.get_wyckoff_sets_conventional()
space_group_number = symmetry_analyzer.get_space_group_number()
def mass_density(self, calculation: Calculation, repr_system: ase.Atoms) -> None:
orig_volume = repr_system.get_volume() * (1e-10)**3
mass = structure.get_summed_atomic_mass(repr_system.get_atomic_numbers())
calculation.mass_density = float(mass / orig_volume)
@abstractmethod
def material_hash(self, material: Material, section_system: Dict) -> None:
pass
# Create and store hash based on SHA512
norm_hash_string = structure.get_symmetry_string(space_group_number, wyckoff_sets)
material.material_hash = sha512(norm_hash_string.encode('utf-8')).hexdigest()
@abstractmethod
def material_name(self, material: Material, symbols: list, numbers: list) -> None:
......@@ -494,10 +490,6 @@ class Structure():
def number_of_atoms(self, material: Material, std_atoms: ase.Atoms) -> None:
material.number_of_atoms = len(std_atoms)
@abstractmethod
def periodicity(self, material: Material) -> None:
pass
@abstractmethod
def point_group(self, material: Material, section_symmetry: Dict) -> None:
pass
......@@ -510,15 +502,10 @@ class Structure():
class StructureBulk(Structure):
"""Processes structure related metainfo for Encyclopedia bulk structures.
"""
def material_hash(self, material: Material, section_symmetry: Dict) -> None:
# Get symmetry information from the section
space_group_number = section_symmetry["space_group_number"]
section_std_system = section_symmetry["section_std_system"][0]
wyckoff_sets = section_std_system.tmp["wyckoff_sets"]
# Create and store hash based on SHA512
norm_hash_string = structure.get_symmetry_string(space_group_number, wyckoff_sets)
material.material_hash = sha512(norm_hash_string.encode('utf-8')).hexdigest()
def atomic_density(self, calculation: Calculation, repr_system: ase.Atoms) -> None:
orig_n_atoms = len(repr_system)
orig_volume = repr_system.get_volume() * (1e-10)**3
calculation.atomic_density = float(orig_n_atoms / orig_volume)
def bravais_lattice(self, material: Material, section_symmetry: Dict) -> None:
bravais_lattice = section_symmetry["bravais_lattice"]
......@@ -554,6 +541,11 @@ class StructureBulk(Structure):
cell_normalized = std_atoms.get_cell()
calculation.lattice_parameters = structure.get_lattice_parameters(cell_normalized)
def mass_density(self, calculation: Calculation, repr_system: ase.Atoms) -> None:
orig_volume = repr_system.get_volume() * (1e-10)**3
mass = structure.get_summed_atomic_mass(repr_system.get_atomic_numbers())
calculation.mass_density = float(mass / orig_volume)
def material_name(self, material: Material, symbols: list, numbers: list) -> None:
# Systems with one element are named after it
if len(symbols) == 1:
......@@ -698,26 +690,27 @@ class StructureBulk(Structure):
wyckoff_list.append(data)
material.wyckoff_groups = json.dumps(wyckoff_list, sort_keys=True)
def fill(self, backend, representative_system: Dict) -> None:
def fill(self, representative_system) -> None:
# Fetch resources
sec_enc = backend.get_mi2_section(Encyclopedia.m_def)
sec_enc = self.backend.get_mi2_section(Encyclopedia.m_def)
material = sec_enc.material
calculation = sec_enc.calculation
sec_symmetry = representative_system["section_symmetry"][0]
symmetry_analyzer = representative_system["section_symmetry"][0].tmp["symmetry_analyzer"]
std_atoms = sec_symmetry["section_std_system"][0].tmp["std_atoms"] # Temporary value stored by SystemNormalizer
prim_atoms = sec_symmetry["section_primitive_system"][0].tmp["prim_atoms"] # Temporary value stored by SystemNormalizer
repr_atoms = sec_symmetry["section_original_system"][0].tmp["orig_atoms"] # Temporary value stored by SystemNormalizer
wyckoff_sets = sec_symmetry["section_std_system"][0].tmp["wyckoff_sets"]
std_atoms = symmetry_analyzer.get_conventional_system()
prim_atoms = symmetry_analyzer.get_primitive_system()
repr_atoms = representative_system.tmp["representative_atoms"] # Temporary value stored by SystemNormalizer
wyckoff_sets = symmetry_analyzer.get_wyckoff_sets_conventional()
names, counts = structure.get_hill_decomposition(prim_atoms.get_chemical_symbols(), reduced=False)
greatest_common_divisor = reduce(gcd, counts)
reduced_counts = np.array(counts) / greatest_common_divisor
# Fill structural information
self.mass_density(calculation, repr_atoms)
self.material_hash(material, sec_symmetry)
self.material_hash(material, symmetry_analyzer)
self.number_of_atoms(material, std_atoms)
self.atom_labels(material, std_atoms)
self.atom_positions(material, std_atoms)
self.atomic_density(calculation, repr_atoms)
self.bravais_lattice(material, sec_symmetry)
self.cell_normalized(material, std_atoms)
......@@ -738,48 +731,110 @@ class StructureBulk(Structure):
class Structure2D(Structure):
"""Processes structure related metainfo for Encyclopedia 2D structures.
"""
def material_hash(self, material: Material, section_symmetry: Dict) -> None:
pass
def get_symmetry_analyzer(self, original_system: Atoms) -> SymmetryAnalyzer:
# Determine the periodicity by examining vacuum gaps
vacuum_directions = structure.find_vacuum_directions(original_system, threshold=7.0)
periodicity = np.invert(vacuum_directions)
# If two axis are not periodic, return. This only happens if the vacuum
# gap is not aligned with a cell vector.
if sum(periodicity) != 2:
self.logger.warn("Could not detect the periodic dimensions in a 2D system.")
return False
# Center the system in the non-periodic direction, also taking
# periodicity into account. The get_center_of_mass()-function in MatID
# takes into account periodicity and can produce the correct CM unlike
# the similar function in ASE.
pbc_cm = matid.geometry.get_center_of_mass(original_system)
cell_center = 0.5 * np.sum(original_system.get_cell(), axis=0)
translation = cell_center - pbc_cm
translation[periodicity] = 0
symm_system = original_system.copy()
symm_system.translate(translation)
# Set the periodicity according to detected periodicity in order for
# SymmetryAnalyzer to use the symmetry analysis designed for 2D
# systems.
symm_system.set_pbc(periodicity)
symmetry_analyzer = SymmetryAnalyzer(
symm_system,
config.normalize.symmetry_tolerance,
config.normalize.flat_dim_threshold
)
return symmetry_analyzer
def get_non_periodic_index_std(self, symmetry_analyzer: SymmetryAnalyzer) -> int:
# Get the periodicity as detected by classification
pbc = symmetry_analyzer._original_system.get_pbc()
non_periodic_index_orig = np.where(pbc == False)[0][0] # noqa: E712
# The index of the originally non-periodic dimension may not correspond
# to the one in the normalized system, because the normalized system
# may use a different coordinate system. We will have to get this
# transformation_matrix = self.symmetry_dataset["transformation_matrix"]
transformation_matrix = symmetry_analyzer._get_spglib_transformation_matrix()
for i_axis, axis in enumerate(transformation_matrix):
if axis[non_periodic_index_orig] != 0 and \
axis[(non_periodic_index_orig + 1) % 3] == 0.0 and \
axis[(non_periodic_index_orig + 2) % 3] == 0.0:
non_periodic_index_std = i_axis
break
return non_periodic_index_std
def cell_angles_string(self, calculation: Calculation) -> None:
pass
def cell_normalized(self, material: Material, std_atoms: ase.Atoms) -> None:
pass
cell_normalized = std_atoms.get_cell()
cell_normalized *= 1e-10
material.cell_normalized = cell_normalized
def cell_primitive(self, material: Material, prim_atoms: ase.Atoms) -> None:
pass
cell_prim = prim_atoms.get_cell()
cell_prim *= 1e-10
material.cell_primitive = cell_prim
def lattice_parameters(self, calculation: Calculation, std_atoms: ase.Atoms) -> None:
pass
def lattice_parameters(self, calculation: Calculation, symmetry_analyzer: SymmetryAnalyzer, std_atoms: Atoms, non_periodic_index_std: int) -> None:
# Eliminate the parameters in the non-periodic dimension
periodicity = np.array([True, True, True])
periodicity[non_periodic_index_std] = False
full_parameters = structure.get_lattice_parameters(std_atoms.get_cell())
lengths = np.array(full_parameters[:3])
angles = np.array(full_parameters[3:])
def periodicity(self, material: Material) -> None:
pass
a, b = lengths[periodicity] * 1e-10
alpha = angles[non_periodic_index_std]
calculation.lattice_parameters = np.array([a, b, 0.0, alpha, 0.0, 0.0])
def periodicity(self, material: Material, non_periodic_index_std: int) -> None:
periodic_indices = [0, 1, 2]
del periodic_indices[non_periodic_index_std]
material.periodicity = np.array(periodic_indices, dtype=np.int8)
def fill(self, backend, representative_system: Dict) -> None:
def fill(self, representative_system) -> None:
# Fetch resources
sec_enc = backend.get_mi2_section(Encyclopedia.m_def)
sec_enc = self.backend.get_mi2_section(Encyclopedia.m_def)
material = sec_enc.material
calculation = sec_enc.calculation
sec_symmetry = representative_system["section_symmetry"][0]
std_atoms = sec_symmetry["section_std_system"][0].tmp["std_atoms"] # Temporary value stored by SystemNormalizer
prim_atoms = sec_symmetry["section_primitive_system"][0].tmp["prim_atoms"] # Temporary value stored by SystemNormalizer
repr_atoms = sec_symmetry["section_original_system"][0].tmp["orig_atoms"] # Temporary value stored by SystemNormalizer
repr_atoms = representative_system.tmp["representative_atoms"] # Temporary value stored by SystemNormalizer
symmetry_analyzer = self.get_symmetry_analyzer(repr_atoms)
std_atoms = symmetry_analyzer.get_conventional_system()
prim_atoms = symmetry_analyzer.get_primitive_system()
names, counts = structure.get_hill_decomposition(prim_atoms.get_chemical_symbols(), reduced=False)
greatest_common_divisor = reduce(gcd, counts)
reduced_counts = np.array(counts) / greatest_common_divisor
non_periodic_index_std = self.get_non_periodic_index_std(symmetry_analyzer)
# Fill structural information
self.mass_density(calculation, repr_atoms)
self.material_hash(material, sec_symmetry)
# Fill metainfo
self.material_hash(material, symmetry_analyzer)
self.number_of_atoms(material, std_atoms)
self.atom_labels(material, std_atoms)
self.atomic_density(calculation, repr_atoms)
self.atom_positions(material, std_atoms)
self.cell_normalized(material, std_atoms)
self.cell_volume(calculation, std_atoms)
self.cell_primitive(material, prim_atoms)
self.formula(material, names, counts)
self.formula_reduced(material, names, reduced_counts)
self.lattice_parameters(calculation, std_atoms)
self.cell_angles_string(calculation)
self.periodicity(material)
self.lattice_parameters(calculation, symmetry_analyzer, std_atoms, non_periodic_index_std)
self.periodicity(material, non_periodic_index_std)
......@@ -16,6 +16,7 @@ from math import gcd as gcd
from typing import List, Dict, Tuple
from functools import reduce
import numpy as np
from ase import Atoms
from nomad.constants import NUMBER_TO_MASS_MAP_KG
......@@ -208,3 +209,43 @@ def get_formula_string(symbols: List[str], counts: List[int]) -> str:
else:
formula += symbol
return formula
def find_vacuum_directions(system: Atoms, threshold: float) -> np.array:
"""Searches for vacuum gaps that are separating the periodic copies.
Args:
system: The structure to analyze
threshold: Vacuum threshold in angstroms
Returns:
np.ndarray: An array with a boolean for each lattice basis
direction indicating if there is enough vacuum to separate the
copies in that direction.
"""
rel_pos = system.get_scaled_positions()
pbc = system.get_pbc()
# Find the maximum vacuum gap for all basis vectors
gaps = np.empty(3, dtype=bool)
for axis in range(3):
if not pbc[axis]:
gaps[axis] = True
continue
comp = rel_pos[:, axis]
ind = np.sort(comp)
ind_rolled = np.roll(ind, 1, axis=0)
distances = ind - ind_rolled
# The first distance is from first to last, so it needs to be
# wrapped around
distances[0] += 1
# Find maximum gap in cartesian coordinates
max_gap = np.max(distances)
basis = system.get_cell()[axis, :]
max_gap_cartesian = np.linalg.norm(max_gap * basis)
has_vacuum_gap = max_gap_cartesian >= threshold
gaps[axis] = has_vacuum_gap
return gaps
......@@ -242,6 +242,9 @@ class SystemNormalizer(SystemBasedNormalizer):
set_value('configuration_raw_gid', configuration_id)
if is_representative:
# Save the Atoms as a temporary variable
self._backend.add_tmp_value("section_system", "representative_atoms", atoms)
# system type analysis
if atom_positions is not None:
with utils.timer(
......@@ -249,8 +252,10 @@ class SystemNormalizer(SystemBasedNormalizer):
system_size=atoms.get_number_of_atoms()):
self.system_type_analysis(atoms)
system_type = self._backend.get_value("system_type")
# symmetry analysis
if atom_positions is not None and (lattice_vectors is not None or not any(pbc)):
if atom_positions is not None and (lattice_vectors is not None or not any(pbc)) and system_type == "bulk":
with utils.timer(
self.logger, 'symmetry analysis executed',
system_size=atoms.get_number_of_atoms()):
......@@ -294,13 +299,11 @@ class SystemNormalizer(SystemBasedNormalizer):
self._backend.addValue('system_type', system_type)
def symmetry_analysis(self, atoms: ase.Atoms) -> None:
"""Analyze the symmetry of the material being simulated.
"""Analyze the symmetry of the material being simulated. Only performed
for bulk materials.
We feed in the parsed values in section_system to the
the symmetry analyzer. We then use the MatID library
to classify the system as 0D, 1D, 2D or 3D and more specific
when possible. When lattice vectors or simulation cells are not present
we skip this analysis. The analysis results are written to the backend.
We feed in the parsed values in section_system to the the symmetry
analyzer. The analysis results are written to the backend.
Args:
atoms: The atomistic structure to analyze.
......@@ -371,8 +374,6 @@ class SystemNormalizer(SystemBasedNormalizer):
self._backend.addArrayValues('atomic_numbers_std', conv_num)
self._backend.addArrayValues('wyckoff_letters_std', conv_wyckoff)
self._backend.addArrayValues('equivalent_atoms_std', conv_equivalent_atoms)
self._backend.add_tmp_value("section_std_system", "wyckoff_sets", wyckoff_sets)
self._backend.add_tmp_value("section_std_system", "std_atoms", conv_sys)
self._backend.closeSection('section_std_system', std_gid)
prim_gid = self._backend.openSection('section_primitive_system')
......@@ -381,14 +382,12 @@ class SystemNormalizer(SystemBasedNormalizer):
self._backend.addArrayValues('atomic_numbers_primitive', prim_num)
self._backend.addArrayValues('wyckoff_letters_primitive', prim_wyckoff)
self._backend.addArrayValues('equivalent_atoms_primitive', prim_equivalent_atoms)
self._backend.add_tmp_value("section_primitive_system", "prim_atoms", prim_sys)
self._backend.closeSection('section_primitive_system', prim_gid)
orig_gid = self._backend.openSection('section_original_system')
self._backend.addArrayValues('wyckoff_letters_original', orig_wyckoff)
self._backend.addArrayValues('equivalent_atoms_original', orig_equivalent_atoms)
self._backend.closeSection('section_original_system', orig_gid)
self._backend.add_tmp_value("section_original_system", "orig_atoms", atoms)
self._backend.closeSection('section_symmetry', symmetry_gid)
self.springer_classification(atoms, space_group_number) # Springer Normalizer
......
......@@ -50,28 +50,28 @@ def test_system_type(geometry_optimization: Encyclopedia):
def test_bulk_information(geometry_optimization: Encyclopedia):
"""Tests that information for bulk systems is correctly processed."
"""
go = geometry_optimization
assert go.material.system_type == "bulk"
assert go.material.number_of_atoms == 4
assert go.material.atom_labels == ["Na", "Na", "Na", "Na"]
assert go.material.crystal_system == "cubic"
assert go.material.bravais_lattice == "cF"
assert go.material.formula == "Na"
assert go.material.formula_reduced == "Na"
assert go.material.has_free_wyckoff_parameters is False
assert go.material.material_name == "Sodium"
assert go.material.point_group == "m-3m"
assert go.material.cell_normalized is not None
assert go.material.cell_primitive is not None
assert np.array_equal(go.material.periodicity, [0, 1, 2])
assert go.material.wyckoff_groups is not None
assert go.calculation.atomic_density == pytest.approx(4.0e+30, rel=0.000001, abs=None)
assert go.calculation.lattice_parameters is not None
assert go.calculation.cell_angles_string is not None
assert go.calculation.mass_density == 4 * 22.98976928 * 1.6605389e-27 / 1e-30 # Atomic mass in kg / cell volume
assert go.calculation.cell_volume == 1e-30
enc = geometry_optimization
assert enc.material.system_type == "bulk"
assert enc.material.number_of_atoms == 4
assert enc.material.atom_labels == ["Na", "Na", "Na", "Na"]
assert enc.material.atom_positions is not None
assert enc.material.crystal_system == "cubic"
assert enc.material.bravais_lattice == "cF"
assert enc.material.formula == "Na"
assert enc.material.formula_reduced == "Na"
assert enc.material.has_free_wyckoff_parameters is False
assert enc.material.material_name == "Sodium"
assert enc.material.point_group == "m-3m"
assert enc.material.cell_normalized is not None
assert enc.material.cell_primitive is not None
assert np.array_equal(enc.material.periodicity, [0, 1, 2])
assert enc.material.wyckoff_groups is not None
assert enc.calculation.atomic_density == pytest.approx(4.0e+30, rel=0.000001, abs=None)
assert enc.calculation.lattice_parameters is not None
assert enc.calculation.cell_angles_string is not None
assert enc.calculation.mass_density == 4 * 22.98976928 * 1.6605389e-27 / 1e-30 # Atomic mass in kg / cell volume
assert enc.calculation.cell_volume == 1e-30
def test_2d_information(twod: Encyclopedia):
......@@ -79,3 +79,13 @@ def test_2d_information(twod: Encyclopedia):
"""
enc = twod
assert enc.material.system_type == "2D"
assert enc.material.number_of_atoms == 2
assert enc.material.atom_labels == ["C", "C"]
assert enc.material.atom_positions is not None
assert enc.material.cell_normalized is not None
assert enc.material.cell_primitive is not None
assert enc.material.formula == "C2"
assert enc.material.formula_reduced == "C"
assert enc.calculation.lattice_parameters is not None
assert np.array_equal(enc.material.periodicity, [0, 1])
# assert go.calculation.cell_angles_string is not None
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