From 41f196b2b46022c179bdce0dcd44ab338d956551 Mon Sep 17 00:00:00 2001 From: Lorenzo <79980269+bastonero@users.noreply.github.com> Date: Tue, 14 Jan 2025 17:23:23 +0100 Subject: [PATCH] `Hubbard`: add useful utility functions (#996) General routines to perform and handle better Hubbard parameters are added, especially for the inter-site parameters. This functions can be used in a number of use-cases, with special application to their initialization and their robust definition during their self-consistent calculations. In particular, we: * Add the possibility of finding first (intersites) neighbours using nearest neighbour finders via Voronoi algorithm (pymatgen). This is helpful when using the `hp.x` code and having a structure containing onsite Hubbard atoms with different number of neighbours. * Add useful routines to allow for simple Hubbard parameters initialization, by nearest neighbour analysis (pymatgen). Standardization of the lattice and of the atomic positions can also be performed to have all the atoms within the unit cell. The latter might be needed to run the `pw.x` binary, which might complain about inter-sites defined beyond the 3x3x3 supercell. --- docs/source/conf.py | 1 + .../data/hubbard_structure.py | 8 +- src/aiida_quantumespresso/utils/hubbard.py | 392 +++++++++++++++++- tests/fixtures/structures/Fe3O4.cif | 263 ++++++++++++ tests/fixtures/structures/LMT.cif | 30 ++ tests/fixtures/structures/MnCoS.cif | 28 ++ tests/utils/test_hubbard.py | 254 +++++++++++- 7 files changed, 961 insertions(+), 15 deletions(-) create mode 100644 tests/fixtures/structures/Fe3O4.cif create mode 100644 tests/fixtures/structures/LMT.cif create mode 100644 tests/fixtures/structures/MnCoS.cif diff --git a/docs/source/conf.py b/docs/source/conf.py index 02fe239f6..ff86dbd75 100755 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -51,6 +51,7 @@ 'python': ('https://docs.python.org/3.8', None), 'aiida': ('https://aiida.readthedocs.io/projects/aiida-core/en/latest/', None), 'aiida_pseudo': ('http://aiida-pseudo.readthedocs.io/en/latest/', None), + 'pymatgen': ('https://pymatgen.org/', None), } # Settings for the `autoapi.extenstion` automatically generating API docs diff --git a/src/aiida_quantumespresso/data/hubbard_structure.py b/src/aiida_quantumespresso/data/hubbard_structure.py index a5c8d2b63..552b78816 100644 --- a/src/aiida_quantumespresso/data/hubbard_structure.py +++ b/src/aiida_quantumespresso/data/hubbard_structure.py @@ -133,9 +133,11 @@ def append_hubbard_parameter( parameters = HubbardParameters.from_tuple(hp_tuple) hubbard = self.hubbard - if parameters not in hubbard.parameters: - hubbard.parameters.append(parameters) - self.hubbard = hubbard + if parameters in hubbard.parameters: + hubbard.parameters.remove(parameters) + + hubbard.parameters.append(parameters) + self.hubbard = hubbard def pop_hubbard_parameters(self, index: int): """Pop Hubbard parameters in the list. diff --git a/src/aiida_quantumespresso/utils/hubbard.py b/src/aiida_quantumespresso/utils/hubbard.py index 65565191b..2e23d9374 100644 --- a/src/aiida_quantumespresso/utils/hubbard.py +++ b/src/aiida_quantumespresso/utils/hubbard.py @@ -3,19 +3,22 @@ # pylint: disable=no-name-in-module from itertools import product import os -from typing import List, Tuple, Union +from typing import Dict, List, Tuple, Union from aiida.orm import StructureData +import numpy as np from aiida_quantumespresso.common.hubbard import Hubbard from aiida_quantumespresso.data.hubbard_structure import HubbardStructureData __all__ = ( 'HubbardUtils', + 'initialize_hubbard_parameters', 'get_supercell_atomic_index', 'get_index_and_translation', 'get_hubbard_indices', 'is_intersite_hubbard', + 'max_number_of_neighbours', ) QE_TRANSLATIONS = list(list(item) for item in product((-1, 0, 1), repeat=3)) @@ -251,8 +254,6 @@ def get_hubbard_for_supercell(self, supercell: StructureData, thr: float = 1e-3) :returns: a new ``HubbbardStructureData`` with all the mapped Hubbard parameters """ - import numpy as np - uc_pymat = self.hubbard_structure.get_pymatgen_structure() sc_pymat = supercell.get_pymatgen_structure() uc_positions = uc_pymat.cart_coords # positions in Cartesian coordinates @@ -313,6 +314,376 @@ def get_hubbard_for_supercell(self, supercell: StructureData, thr: float = 1e-3) return HubbardStructureData.from_structure(structure=supercell, hubbard=new_hubbard) + def get_interacting_pairs(self) -> Dict[str, List[str]]: + """Return tuple of kind name interaction pairs. + + :returns: dictionary of onsite kinds with a list of V kinds + """ + pairs_dict = dict() + sites = self.hubbard_structure.sites + parameters = self.hubbard_structure.hubbard.parameters + + onsite_indices = [p.atom_index for p in parameters if p.atom_index == p.neighbour_index] + + for index in onsite_indices: + name = sites[index].kind_name + if name not in pairs_dict: + pairs_dict[name] = [] + + for parameter in parameters: + onsite_name = sites[parameter.atom_index].kind_name + neigh_name = sites[parameter.neighbour_index].kind_name + + if onsite_name in pairs_dict and onsite_name != neigh_name: + if not neigh_name in pairs_dict[onsite_name]: + pairs_dict[onsite_name].append(neigh_name) + + return pairs_dict + + def get_pairs_radius( + self, + onsite_index: int, + neighbours_names: List[str], + number_of_neighbours: int, + radius_max: float = 7.0, + thr: float = 1.0e-2, + ) -> Tuple[float, float]: + """Return the minimum and maximum radius of the first neighbours of the onsite site. + + :param onsite_index: index in the structure of the onsite Hubbard atom + :param neighbours_names: kind names of the neighbours + :param number_of_neighbours: number of neighbours coming to select + :param radius_max: maximum radius (in Angstrom) to use for looking for neighbours + :param thr: threshold (in Angstrom) for defining the shells + :return: (radius min +thr, radius max -thr) defining the shells containing only the first neighbours + """ + rmin = 0 + pymat = self.hubbard_structure.get_pymatgen_structure() + _, neigh_indices, _, distances = pymat.get_neighbor_list(sites=[pymat[onsite_index]], r=radius_max) + + sort = np.argsort(distances) + neigh_indices = neigh_indices[sort] + distances = distances[sort] + + count = 0 + for i in range(len(neigh_indices)): # pylint: disable=consider-using-enumerate + index = i + if self.hubbard_structure.sites[neigh_indices[i]].kind_name in neighbours_names: + rmin = max(rmin, distances[i]) + count += 1 + if count == number_of_neighbours: + break + + rmax = distances[index + 1] + if np.abs(rmax - rmin) < thr: + rmax = rmin + 2 * rmin + + return rmin + thr, rmax - thr + + def get_intersites_radius( + self, + nn_finder: str = 'crystal', + nn_inputs: Union[Dict, None] = None, + radius_max: float = 7.0, + thr: float = 1.0e-2, + **_, + ) -> float: + """Return the radius (in Angstrom) for intersites from nearest neighbour finders. + + It peforms a nearest neighbour analysis (via pymatgen modules) to find the first inersite + neighbours for all the onsite atoms. A radius is returned which can be used to + run an ``hp.x`` calculation. Such radius defines a shell including only the first + neighbours of each onsite Hubbard atom. + + :param nn_finder: string defining the nearest neighbour finder; options are: + * `crystal`: use :class:`pymatgen.analysis.local_env.CrystalNN` + * `voronoi`: use :class:`pymatgen.analysis.local_env.VoronoiNN` + :param nn_inputs: inputs for the nearest neighbours finder; when None, standard inputs + are used to find geometric first neighbours (recommended) + :param radius_max: max radius where to look for neighbouring atoms, in Angstrom + :param thr: threshold (in Angstrom) for defining the shells + :return: radius defining the shell containing only the first neighbours + """ + import warnings + + from pymatgen.analysis.local_env import CrystalNN, VoronoiNN + + rmin, rmax = 0.0, radius_max + + if nn_finder not in ['crystal', 'voronoi']: + raise ValueError('`nn_finder` must be either `cyrstal` or `voronoi`') + if nn_inputs is None: + if nn_finder == 'crystal': + nn_inputs = {'distance_cutoffs': None, 'x_diff_weight': 0, 'porous_adjustment': False} + if nn_finder == 'voronoi': + nn_inputs = {'tol': 0.1, 'cutoff': radius_max} + + voronoi = CrystalNN(**nn_inputs) if nn_finder == 'crystal' else VoronoiNN(**nn_inputs) # pylint: disable=unexpected-keyword-arg + + sites = self.hubbard_structure.sites + name_to_specie = {kind.name: kind.symbol for kind in self.hubbard_structure.kinds} + pymat = self.hubbard_structure.get_pymatgen_structure() + pairs = self.get_interacting_pairs() + + for i, site in enumerate(sites): + if site.kind_name in pairs: + neigh_species = voronoi.get_cn_dict(pymat, i) # e.g. {'O': 4, 'S': 2, ...} + number_of_neighs = 0 + + for neigh_name in pairs[site.kind_name]: + specie = name_to_specie[neigh_name] + + if specie in neigh_species: + number_of_neighs += neigh_species[specie] + neigh_species.pop(specie) # avoid 'duplicating' same specie but different (kind) name + + rmin_, rmax_ = self.get_pairs_radius(i, pairs[site.kind_name], number_of_neighs, radius_max, thr) + + rmin = max(rmin_, rmin) # we want the largest to include them all + rmax = min(rmax_, rmax) # we want the smallest to check whether such radius exist + + if rmin > rmax: + warnings.warn('A common radius seems to not exist! Try lowering `thr`.') + + return min(rmin, rmax) + + def get_intersites_list( + self, + nn_finder: str = 'crystal', + nn_inputs: Union[Dict, None] = None, + radius_max: float = 7.0, + **_, + ) -> List[Tuple[int, int, Tuple[int, int, int]]]: + """Return the list of intersites from nearest neighbour finders. + + It peforms a nearest neighbour analysis (via pymatgen modules) to find the first inersite + neighbours for all the onsite atoms. A list is returned with all the nearest neighbours + providing all the information about the couples indices and the associated trasnaltion vector. + Also on-site information is included. + + :param nn_finder: string defining the nearest neighbour finder; options are: + * `crystal`: use :class:`pymatgen.analysis.local_env.CrystalNN` + * `voronoi`: use :class:`pymatgen.analysis.local_env.VoronoiNN` + :param nn_inputs: inputs for the nearest neighbours finder; when None, standard inputs + are used to find geometric first neighbours (recommended) + :param radius_max: max radius where to look for neighbouring atoms, in Angstrom + :param thr: threshold (in Angstrom) for defining the shells + :return: list of lists, each having (atom index, neighbouring index, translation vector) + """ + from pymatgen.analysis.local_env import CrystalNN, VoronoiNN + + if nn_finder not in ['crystal', 'voronoi']: + raise ValueError('`nn_finder` must be either `cyrstal` or `voronoi`') + + if nn_inputs is None: + if nn_finder == 'crystal': + nn_inputs = {'distance_cutoffs': None, 'x_diff_weight': 0, 'porous_adjustment': False} + if nn_finder == 'voronoi': + nn_inputs = {'tol': 0.1, 'cutoff': radius_max} + + voronoi = CrystalNN(**nn_inputs) if nn_finder == 'crystal' else VoronoiNN(**nn_inputs) # pylint: disable=unexpected-keyword-arg + + sites = self.hubbard_structure.sites + name_to_specie = {kind.name: kind.symbol for kind in self.hubbard_structure.kinds} + pymat = self.hubbard_structure.get_pymatgen_structure() + pairs = self.get_interacting_pairs() + neigh_list = [] + + for i, site in enumerate(sites): # pylint: disable=too-many-nested-blocks + if site.kind_name in pairs: + neigh_species = voronoi.get_cn_dict(pymat, i) # e.g. {'O': 4, 'S': 2, ...} + neigh_list.append([i, i, (0, 0, 0)]) + + for neigh_name in pairs[site.kind_name]: + try: # we want an API that allows for specifying neighbours that are not present + specie = name_to_specie[neigh_name] + + count = 0 + if specie in neigh_species: + _, neigh_indices, images, distances = pymat.get_neighbor_list( + sites=[pymat[i]], r=radius_max + ) + + sort = np.argsort(distances) + neigh_indices = neigh_indices[sort] + images = images[sort] + + for index, image in zip(neigh_indices, images): + if pymat[index].specie.name == specie: + neigh_list.append([ + i, + int(index), + tuple(np.array(image, dtype=np.int64).tolist()), + ]) + count += 1 + + if count >= neigh_species[specie]: + break + except KeyError: + pass + + return neigh_list + + def get_max_number_of_neighbours( + self, + nn_finder: str = 'crystal', + nn_inputs: Union[Dict, None] = None, + radius_max: float = 7.0, + **_, + ) -> list: + """Return the maximum number of nearest neighbours, aslo counting the non-interacting ones. + + It peforms a nearest neighbour analysis (via pymatgen modules) to find the first inersite + neighbours for all the onsite atoms. A list is returned with all the nearest neighbours + providing all the information about the couples indices and the associated trasnaltion vector. + Also on-site information is included. + + :param nn_finder: string defining the nearest neighbour finder; options are: + * `crystal`: use :class:`pymatgen.analysis.local_env.CrystalNN` + * `voronoi`: use :class:`pymatgen.analysis.local_env.VoronoiNN` + :param nn_inputs: inputs for the nearest neighbours finder; when None, standard inputs + are used to find geometric first neighbours (recommended) + :param radius_max: max radius where to look for neighbouring atoms, in Angstrom + :param thr: threshold (in Angstrom) for defining the shells + :return: list of lists, each having (atom index, neighbouring index, translation vector) + """ + from pymatgen.analysis.local_env import CrystalNN, VoronoiNN + + if nn_finder not in ['crystal', 'voronoi']: + raise ValueError('`nn_finder` must be either `cyrstal` or `voronoi`') + + if nn_inputs is None: + if nn_finder == 'crystal': + nn_inputs = {'distance_cutoffs': None, 'x_diff_weight': 0, 'porous_adjustment': False} + if nn_finder == 'voronoi': + nn_inputs = {'tol': 0.1, 'cutoff': radius_max} + + voronoi = CrystalNN(**nn_inputs) if nn_finder == 'crystal' else VoronoiNN(**nn_inputs) # pylint: disable=unexpected-keyword-arg + + sites = self.hubbard_structure.sites + pymat = self.hubbard_structure.get_pymatgen_structure() + pairs = self.get_interacting_pairs() + max_num_of_neighs = 0 + + for i, site in enumerate(sites): # pylint: disable=too-many-nested-blocks + if site.kind_name in pairs: + neigh_species = voronoi.get_cn_dict(pymat, i) # e.g. {'O': 4, 'S': 2, ...} + max_num_of_neighs = max(max_num_of_neighs, np.sum(list(neigh_species.values()))) + + return max_num_of_neighs + + +def initialize_hubbard_parameters( + structure: StructureData, + pairs: Dict[str, Tuple[str, float, float, Dict[str, str]]], + nn_finder: str = 'crystal', + nn_inputs: Union[Dict, None] = None, + fold: bool = True, + standardize: bool = False, + radius_max: float = 7.0, + thr: float = 1e-5, + use_kinds: bool = True, + **_, +) -> HubbardStructureData: + """Initialize the on-site and intersite parameters using nearest neighbour finders. + + It peforms a nearest neighbour analysis (via pymatgen modules) to find the first inersite + neighbours for all the onsite atoms. Only the atoms in the `pair`, and that are nearest + neighbours from the analysis, are initialized. + + :param structure: a StructureData instance + :param pairs: dictionary of the kind + {onsite name: [onsite manifold, onsite value, intersites value, {neighbour name: neighbour manifold}], ...} + For example: {'Fe': ['3d', 5.0, 1.0, {'O':'2p', 'O1':'1s', 'Se':'4p'}]} + :param nn_finder: string defining the nearest neighbour finder; options are: + * `crystal`: use :class:`pymatgen.analysis.local_env.CrystalNN` + * `voronoi`: use :class:`pymatgen.analysis.local_env.VoronoiNN` + :param nn_inputs: inputs for the nearest neighbours finder; when None, standard inputs + are used to find geometric first neighbours (recommended) + :param fold: whether to fold in within the cell the atoms + :param standardize: whether to standardize the atoms and the cell via spglib (symmetry analysis) + :param radius_max: max radius where to look for neighbouring atoms, in Angstrom + :param thr: threshold to refold the atoms with crystal coordinates close to 1.0 + :param use_kinds: whether to use kinds for initializing the parameters; when False, it + initializes all the ``Kinds`` matching the given specie + :return: HubbardStructureData with initialized Hubbard parameters + """ + from aiida.tools.data import spglib_tuple_to_structure, structure_to_spglib_tuple + from pymatgen.analysis.local_env import CrystalNN, VoronoiNN + from spglib import standardize_cell + + if nn_finder not in ['crystal', 'voronoi']: + raise ValueError('`nn_finder` must be either `cyrstal` or `voronoi`') + + if nn_inputs is None: + if nn_finder == 'crystal': + nn_inputs = {'distance_cutoffs': None, 'x_diff_weight': 0, 'porous_adjustment': False} + if nn_finder == 'voronoi': + nn_inputs = {'tol': 0.1, 'cutoff': radius_max} + + voronoi = CrystalNN(**nn_inputs) if nn_finder == 'crystal' else VoronoiNN(**nn_inputs) # pylint: disable=unexpected-keyword-arg + + if not standardize and not fold: + hubbard_structure = HubbardStructureData.from_structure(structure=structure) + + if fold: + cell, kind_info, kinds = structure_to_spglib_tuple(structure) + cell = list(cell) + cell[1] = cell[1] % 1.0 + cell[1] = np.where(np.abs(cell[1] - 1.0) < thr, 0.0, cell[1]) + + folded_struccture = spglib_tuple_to_structure(cell, kind_info, kinds) + hubbard_structure = HubbardStructureData.from_structure(structure=folded_struccture) + + if standardize: + cell, kind_info, kinds = structure_to_spglib_tuple(structure) + new_cell = standardize_cell(cell) + new_structure = spglib_tuple_to_structure(new_cell, kind_info, kinds) + hubbard_structure = HubbardStructureData.from_structure(structure=new_structure) + + sites = hubbard_structure.sites + name_to_specie = {kind.name: kind.symbol for kind in hubbard_structure.kinds} + pymat = hubbard_structure.get_pymatgen_structure() + + for i, site in enumerate(sites): # pylint: disable=too-many-nested-blocks + site_name = site.kind_name if use_kinds else name_to_specie[site.kind_name] + if site_name in pairs: + neigh_species = voronoi.get_cn_dict(pymat, i) # e.g. {'O': 4, 'S': 2, ...} + onsite = pairs[site.kind_name] + hubbard_structure.append_hubbard_parameter(i, onsite[0], i, onsite[0], onsite[1]) + + for neigh_name in onsite[3]: + try: # we want an API that allows for specifying neighbours that are not present + specie = name_to_specie[neigh_name] + + count = 0 + if specie in neigh_species: + _, neigh_indices, images, distances = pymat.get_neighbor_list(sites=[pymat[i]], r=radius_max) + + sort = np.argsort(distances) + neigh_indices = neigh_indices[sort] + images = images[sort] + + for index, image in zip(neigh_indices, images): + if pymat[index].specie.name == specie: + hubbard_structure.append_hubbard_parameter( + atom_index=i, + atom_manifold=onsite[0], + neighbour_index=int(index), # otherwise the validator complains + neighbour_manifold=onsite[3][neigh_name], + value=onsite[2], + translation=np.array(image, dtype=np.int64).tolist(), + hubbard_type='V', + ) + count += 1 + + if count >= neigh_species[specie]: + break + except KeyError: + pass + + return hubbard_structure + def get_supercell_atomic_index(index: int, num_sites: int, translation: List[Tuple[int, int, int]]) -> int: """Return the atomic index in 3x3x3 supercell. @@ -351,3 +722,18 @@ def is_intersite_hubbard(hubbard: Hubbard) -> bool: """Return whether `Hubbard` contains intersite interactions (+V).""" couples = [(param.atom_index != param.neighbour_index) for param in hubbard.parameters] return any(couples) + + +def max_number_of_neighbours(intersites_list: List[Tuple[int, int]]) -> int: + """Return the maximum number of neighbours found. + + .. note:: it assumes only one onsite parameter is defined per atom index, + that means `intra` Hubbard interactions are not defined, such as `V Fe 3d Fe 2p 1 1 1.0` + + :param intersites_list: list of lists of shape (atom index, neigbours index) + """ + from collections import Counter + + first_indices = np.array(intersites_list, dtype='object')[:, 0] + counts = Counter(first_indices) + return max(counts.values()) - 1 # assuming there is always 1 and only 1 onsite parameter diff --git a/tests/fixtures/structures/Fe3O4.cif b/tests/fixtures/structures/Fe3O4.cif new file mode 100644 index 000000000..51ac1b00c --- /dev/null +++ b/tests/fixtures/structures/Fe3O4.cif @@ -0,0 +1,263 @@ +#------------------------------------------------------------------------------ +#$Date: 2017-10-13 02:32:00 +0300 (Fri, 13 Oct 2017) $ +#$Revision: 201954 $ +#$URL: file:///home/coder/svn-repositories/cod/cif/1/01/03/1010369.cif $ +#------------------------------------------------------------------------------ +# +# This file is available in the Crystallography Open Database (COD), +# http://www.crystallography.net/ +# +# All data on this site have been placed in the public domain by the +# contributors. +# +data_1010369 +loop_ +_publ_author_name +'Montoro, V' +_publ_section_title +; +Miscibilita fra gli ossidi salini di ferro e di manganese +; +_journal_coden_ASTM GCITA9 +_journal_name_full 'Gazzetta Chimica Italiana' +_journal_page_first 728 +_journal_page_last 733 +_journal_volume 68 +_journal_year 1938 +_chemical_formula_structural 'Fe3 O4' +_chemical_formula_sum 'Fe3 O4' +_chemical_name_systematic 'Iron diiron(III) oxide' +_space_group_IT_number 227 +_symmetry_cell_setting cubic +_symmetry_space_group_name_Hall 'F 4d 2 3 -1d' +_symmetry_space_group_name_H-M 'F d -3 m :1' +_cell_angle_alpha 90 +_cell_angle_beta 90 +_cell_angle_gamma 90 +_cell_formula_units_Z 8 +_cell_length_a 8.384 +_cell_length_b 8.384 +_cell_length_c 8.384 +_cell_volume 589.3 +_cod_original_sg_symbol_H-M 'F d -3 m S' +_cod_database_code 1010369 +loop_ +_symmetry_equiv_pos_as_xyz +x,y,z +y,z,x +z,x,y +x,z,y +y,x,z +z,y,x +x,-y,-z +y,-z,-x +z,-x,-y +x,-z,-y +y,-x,-z +z,-y,-x +-x,y,-z +-y,z,-x +-z,x,-y +-x,z,-y +-y,x,-z +-z,y,-x +-x,-y,z +-y,-z,x +-z,-x,y +-x,-z,y +-y,-x,z +-z,-y,x +1/4-x,1/4-y,1/4-z +1/4-y,1/4-z,1/4-x +1/4-z,1/4-x,1/4-y +1/4-x,1/4-z,1/4-y +1/4-y,1/4-x,1/4-z +1/4-z,1/4-y,1/4-x +1/4-x,1/4+y,1/4+z +1/4-y,1/4+z,1/4+x +1/4-z,1/4+x,1/4+y +1/4-x,1/4+z,1/4+y +1/4-y,1/4+x,1/4+z +1/4-z,1/4+y,1/4+x +1/4+x,1/4-y,1/4+z +1/4+y,1/4-z,1/4+x +1/4+z,1/4-x,1/4+y +1/4+x,1/4-z,1/4+y +1/4+y,1/4-x,1/4+z +1/4+z,1/4-y,1/4+x +1/4+x,1/4+y,1/4-z +1/4+y,1/4+z,1/4-x +1/4+z,1/4+x,1/4-y +1/4+x,1/4+z,1/4-y +1/4+y,1/4+x,1/4-z +1/4+z,1/4+y,1/4-x +x,1/2+y,1/2+z +1/2+x,y,1/2+z +1/2+x,1/2+y,z +y,1/2+z,1/2+x +1/2+y,z,1/2+x +1/2+y,1/2+z,x +z,1/2+x,1/2+y +1/2+z,x,1/2+y +1/2+z,1/2+x,y +x,1/2+z,1/2+y +1/2+x,z,1/2+y +1/2+x,1/2+z,y +y,1/2+x,1/2+z +1/2+y,x,1/2+z +1/2+y,1/2+x,z +z,1/2+y,1/2+x +1/2+z,y,1/2+x +1/2+z,1/2+y,x +x,1/2-y,1/2-z +1/2+x,-y,1/2-z +1/2+x,1/2-y,-z +y,1/2-z,1/2-x +1/2+y,-z,1/2-x +1/2+y,1/2-z,-x +z,1/2-x,1/2-y +1/2+z,-x,1/2-y +1/2+z,1/2-x,-y +x,1/2-z,1/2-y +1/2+x,-z,1/2-y +1/2+x,1/2-z,-y +y,1/2-x,1/2-z +1/2+y,-x,1/2-z +1/2+y,1/2-x,-z +z,1/2-y,1/2-x +1/2+z,-y,1/2-x +1/2+z,1/2-y,-x +-x,1/2+y,1/2-z +1/2-x,y,1/2-z +1/2-x,1/2+y,-z +-y,1/2+z,1/2-x +1/2-y,z,1/2-x +1/2-y,1/2+z,-x +-z,1/2+x,1/2-y +1/2-z,x,1/2-y +1/2-z,1/2+x,-y +-x,1/2+z,1/2-y +1/2-x,z,1/2-y +1/2-x,1/2+z,-y +-y,1/2+x,1/2-z +1/2-y,x,1/2-z +1/2-y,1/2+x,-z +-z,1/2+y,1/2-x +1/2-z,y,1/2-x +1/2-z,1/2+y,-x +-x,1/2-y,1/2+z +1/2-x,-y,1/2+z +1/2-x,1/2-y,z +-y,1/2-z,1/2+x +1/2-y,-z,1/2+x +1/2-y,1/2-z,x +-z,1/2-x,1/2+y +1/2-z,-x,1/2+y +1/2-z,1/2-x,y +-x,1/2-z,1/2+y +1/2-x,-z,1/2+y +1/2-x,1/2-z,y +-y,1/2-x,1/2+z +1/2-y,-x,1/2+z +1/2-y,1/2-x,z +-z,1/2-y,1/2+x +1/2-z,-y,1/2+x +1/2-z,1/2-y,x +1/4-x,3/4-y,3/4-z +3/4-x,1/4-y,3/4-z +3/4-x,3/4-y,1/4-z +1/4-y,3/4-z,3/4-x +3/4-y,1/4-z,3/4-x +3/4-y,3/4-z,1/4-x +1/4-z,3/4-x,3/4-y +3/4-z,1/4-x,3/4-y +3/4-z,3/4-x,1/4-y +1/4-x,3/4-z,3/4-y +3/4-x,1/4-z,3/4-y +3/4-x,3/4-z,1/4-y +1/4-y,3/4-x,3/4-z +3/4-y,1/4-x,3/4-z +3/4-y,3/4-x,1/4-z +1/4-z,3/4-y,3/4-x +3/4-z,1/4-y,3/4-x +3/4-z,3/4-y,1/4-x +1/4-x,3/4+y,3/4+z +3/4-x,1/4+y,3/4+z +3/4-x,3/4+y,1/4+z +1/4-y,3/4+z,3/4+x +3/4-y,1/4+z,3/4+x +3/4-y,3/4+z,1/4+x +1/4-z,3/4+x,3/4+y +3/4-z,1/4+x,3/4+y +3/4-z,3/4+x,1/4+y +1/4-x,3/4+z,3/4+y +3/4-x,1/4+z,3/4+y +3/4-x,3/4+z,1/4+y +1/4-y,3/4+x,3/4+z +3/4-y,1/4+x,3/4+z +3/4-y,3/4+x,1/4+z +1/4-z,3/4+y,3/4+x +3/4-z,1/4+y,3/4+x +3/4-z,3/4+y,1/4+x +1/4+x,3/4-y,3/4+z +3/4+x,1/4-y,3/4+z +3/4+x,3/4-y,1/4+z +1/4+y,3/4-z,3/4+x +3/4+y,1/4-z,3/4+x +3/4+y,3/4-z,1/4+x +1/4+z,3/4-x,3/4+y +3/4+z,1/4-x,3/4+y +3/4+z,3/4-x,1/4+y +1/4+x,3/4-z,3/4+y +3/4+x,1/4-z,3/4+y +3/4+x,3/4-z,1/4+y +1/4+y,3/4-x,3/4+z +3/4+y,1/4-x,3/4+z +3/4+y,3/4-x,1/4+z +1/4+z,3/4-y,3/4+x +3/4+z,1/4-y,3/4+x +3/4+z,3/4-y,1/4+x +1/4+x,3/4+y,3/4-z +3/4+x,1/4+y,3/4-z +3/4+x,3/4+y,1/4-z +1/4+y,3/4+z,3/4-x +3/4+y,1/4+z,3/4-x +3/4+y,3/4+z,1/4-x +1/4+z,3/4+x,3/4-y +3/4+z,1/4+x,3/4-y +3/4+z,3/4+x,1/4-y +1/4+x,3/4+z,3/4-y +3/4+x,1/4+z,3/4-y +3/4+x,3/4+z,1/4-y +1/4+y,3/4+x,3/4-z +3/4+y,1/4+x,3/4-z +3/4+y,3/4+x,1/4-z +1/4+z,3/4+y,3/4-x +3/4+z,1/4+y,3/4-x +3/4+z,3/4+y,1/4-x +loop_ +_atom_site_label +_atom_site_type_symbol +_atom_site_symmetry_multiplicity +_atom_site_Wyckoff_symbol +_atom_site_fract_x +_atom_site_fract_y +_atom_site_fract_z +_atom_site_occupancy +_atom_site_attached_hydrogens +_atom_site_calc_flag +Fe1 Fe2+ 16 d 0.625 0.625 0.625 0.5 0 d +Fe2 Fe3+ 16 d 0.625 0.625 0.625 0.5 0 d +Fe3 Fe3+ 8 a 0. 0. 0. 1. 0 d +O1 O2- 32 e 0.375 0.375 0.375 1. 0 d +loop_ +_atom_type_symbol +_atom_type_oxidation_number +Fe2+ 2.000 +Fe3+ 3.000 +O2- -2.000 +loop_ +_cod_related_entry_id +_cod_related_entry_database +_cod_related_entry_code +1 ChemSpider 4937312 diff --git a/tests/fixtures/structures/LMT.cif b/tests/fixtures/structures/LMT.cif new file mode 100644 index 000000000..3ed2120f0 --- /dev/null +++ b/tests/fixtures/structures/LMT.cif @@ -0,0 +1,30 @@ +# generated using pymatgen +data_LiMnTe2 +_symmetry_space_group_name_H-M 'P 1' +_cell_length_a 4.31097682 +_cell_length_b 4.31097682 +_cell_length_c 5.22251802 +_cell_angle_alpha 90.00000000 +_cell_angle_beta 90.00000000 +_cell_angle_gamma 120.00000000 +_symmetry_Int_Tables_number 1 +_chemical_formula_structural LiMnTe2 +_chemical_formula_sum 'Li1 Mn1 Te2' +_cell_volume 84.05469084 +_cell_formula_units_Z 1 +loop_ + _symmetry_equiv_pos_site_id + _symmetry_equiv_pos_as_xyz + 1 'x, y, z' +loop_ + _atom_site_type_symbol + _atom_site_label + _atom_site_symmetry_multiplicity + _atom_site_fract_x + _atom_site_fract_y + _atom_site_fract_z + _atom_site_occupancy + Mn Mn0 1 0.33333333 0.66666667 0.47975861 1.0 + Te Te1 1 0.33333333 0.66666667 -0.01976285 1.0 + Te Te2 1 0.66666667 0.33333333 0.45389956 1.0 + Li Li3 1 -0.00000000 -0.00000000 0.73880468 1.0 diff --git a/tests/fixtures/structures/MnCoS.cif b/tests/fixtures/structures/MnCoS.cif new file mode 100644 index 000000000..2d2b5c4d5 --- /dev/null +++ b/tests/fixtures/structures/MnCoS.cif @@ -0,0 +1,28 @@ +data_image0 +_chemical_formula_structural MnCoS +_chemical_formula_sum "Mn1 Co1 S1" +_cell_length_a 4.02254 +_cell_length_b 4.02254 +_cell_length_c 4.02254 +_cell_angle_alpha 60 +_cell_angle_beta 60 +_cell_angle_gamma 60 + +_space_group_name_H-M_alt "P 1" +_space_group_IT_number 1 + +loop_ + _space_group_symop_operation_xyz + 'x, y, z' + +loop_ + _atom_site_type_symbol + _atom_site_label + _atom_site_symmetry_multiplicity + _atom_site_fract_x + _atom_site_fract_y + _atom_site_fract_z + _atom_site_occupancy + Mn Mn1 1.0 0.00000 0.00000 0.00000 1.0000 + Co Co1 1.0 0.75000 0.75000 0.75000 1.0000 + S S1 1.0 0.50000 0.50000 0.50000 1.0000 diff --git a/tests/utils/test_hubbard.py b/tests/utils/test_hubbard.py index 8846b6a21..8dd2ed2f4 100644 --- a/tests/utils/test_hubbard.py +++ b/tests/utils/test_hubbard.py @@ -3,10 +3,14 @@ # pylint: disable=redefined-outer-name import os +from aiida.orm import StructureData +from ase.io import read +import numpy as np import pytest from aiida_quantumespresso.common.hubbard import Hubbard -from aiida_quantumespresso.utils.hubbard import HubbardUtils +from aiida_quantumespresso.data.hubbard_structure import HubbardStructureData +from aiida_quantumespresso.utils.hubbard import HubbardUtils, initialize_hubbard_parameters @pytest.fixture @@ -14,10 +18,7 @@ def generate_hubbard_structure(): """Return a `HubbardStructureData` instance.""" def _generate_hubbard_structure(parameters=None, projectors=None, formulation=None): - from aiida.orm import StructureData - - from aiida_quantumespresso.data.hubbard_structure import HubbardStructureData - + """Return a `HubbardStructureData` instance.""" a, b, c, d = 1.40803, 0.81293, 4.68453, 1.62585 # pylint: disable=invalid-name cell = [[a, -b, c], [0.0, d, c], [-a, -b, c]] # pylint: disable=invalid-name positions = [[0, 0, 0], [0, 0, 3.6608], [0, 0, 10.392], [0, 0, 7.0268]] @@ -44,8 +45,8 @@ def _generate_hubbard_structure(parameters=None, projectors=None, formulation=No def generate_hubbard_utils(generate_hubbard_structure): """Return a `HubbardUtils` instance.""" - def _generate_hubbard_utils(): - return HubbardUtils(hubbard_structure=generate_hubbard_structure()) + def _generate_hubbard_utils(**kwargs): + return HubbardUtils(hubbard_structure=generate_hubbard_structure(**kwargs)) return _generate_hubbard_utils @@ -186,7 +187,6 @@ def test_is_to_reorder(generate_hubbard_structure): @pytest.mark.usefixtures('aiida_profile') def test_reorder_supercell_atoms(generate_hubbard_structure): """Test the `reorder_atoms` method with a supercell.""" - from aiida.orm import StructureData parameters = [ (0, '3d', 0, '3d', 5.0, (0, 0, 0), 'U'), (0, '3d', 1, '2p', 5.0, (0, 0, 0), 'U'), @@ -214,7 +214,6 @@ def test_hubbard_for_supercell(generate_hubbard_structure): .. warning:: we only test for the ``U`` case, assuming that if U is transfered to the supercell correctly, any parameter will. """ - from aiida.orm import StructureData parameters = [[3, '1s', 3, '2s', 5.0, [0, 0, 0], 'U']] # Li atom projectors = 'atomic' formulation = 'liechtenstein' @@ -235,3 +234,240 @@ def test_hubbard_for_supercell(generate_hubbard_structure): assert parameters[1] == '1s' assert parameters[3] == '2s' assert hubbard_supercell.sites[parameters[0]].kind_name == 'Li' + + +@pytest.fixture +def get_non_trivial_hubbard_structure(filepath_tests): + """Return a multi-coordination number `HubbardStructureData`.""" + + def _get_non_trivial_hubbard_structure(name=None, use_kinds=False): + """Return a multi-coordination number `HubbardStructureData`.""" + path = os.path.join(filepath_tests, 'fixtures', 'structures') + + if name is None: + atoms = read(os.path.join(path, 'Fe3O4.cif')) + + if use_kinds: + tags = [0] * atoms.get_global_number_of_atoms() + for i in range(10): + tags[i] = 1 + tags[0] = 2 + for i in [26, 27, 28, 30, 31, 32]: + tags[i] = 1 + atoms.set_tags(tags) + + hubbard_structure = HubbardStructureData.from_structure(StructureData(ase=atoms)) + pymat = hubbard_structure.get_pymatgen_structure() + + tag_names = ['Fe'] + if use_kinds: + tag_names = ['Fe', 'Fe1', 'Fe2'] + + for i, site in enumerate(hubbard_structure.sites): + if site.kind_name in tag_names: + hubbard_structure.append_hubbard_parameter(i, '3d', i, '3d', 5.0) + pymat_sites = pymat.get_sites_in_sphere(site.position, r=2.96) + + for pymat_site in pymat_sites: + if pymat_site.position_atol != site.position: + # if pymat_site.specie.name != 'Fe': + translation = np.array(pymat_site.image, dtype=np.int64).tolist() + args = (i, '3d', int(pymat_site.index), '2p', 1.0, translation) + hubbard_structure.append_hubbard_parameter(*args) + + return hubbard_structure + + return _get_non_trivial_hubbard_structure + + +def test_get_interacting_pairs(get_non_trivial_hubbard_structure): + """Test the `HubbardUtils.get_interacting_pairs` method.""" + hubbard_structure = get_non_trivial_hubbard_structure() + hubbard_utils = HubbardUtils(hubbard_structure=hubbard_structure) + + pairs = hubbard_utils.get_interacting_pairs() + + assert 'Fe' in pairs + assert pairs['Fe'] == ['O'] + + # using actual kinds + hubbard_structure = get_non_trivial_hubbard_structure(use_kinds=True) + hubbard_utils = HubbardUtils(hubbard_structure=hubbard_structure) + + pairs = hubbard_utils.get_interacting_pairs() + # print(hubbard_utils.get_hubbard_card()) + + assert pairs == {'Fe': ['O', 'O1'], 'Fe1': ['O', 'O1'], 'Fe2': ['O', 'O1']} + + +def test_get_pairs_radius(get_non_trivial_hubbard_structure): + """Test the `HubbardUtils.get_pairs_radius` method.""" + hubbard_structure = get_non_trivial_hubbard_structure() + hubbard_utils = HubbardUtils(hubbard_structure=hubbard_structure) + + rmin, rmax = hubbard_utils.get_pairs_radius(0, ['O'], 6) + assert rmin < 2.2 < rmax + + rmin, rmax = hubbard_utils.get_pairs_radius(23, ['O'], 4) + assert rmin < 2.2 < rmax + + # using actual kinds + hubbard_structure = get_non_trivial_hubbard_structure(use_kinds=True) + hubbard_utils = HubbardUtils(hubbard_structure=hubbard_structure) + + rmin, rmax = hubbard_utils.get_pairs_radius(0, ['O', 'O1'], 6) + assert rmin < 2.2 < rmax + + +def test_get_intersites_radius(get_non_trivial_hubbard_structure): + """Test the `HubbardUtils.get_intersites_radius` method.""" + hubbard_structure = get_non_trivial_hubbard_structure() + hubbard_utils = HubbardUtils(hubbard_structure=hubbard_structure) + + radius = hubbard_utils.get_intersites_radius() + + assert abs(radius - 2.106) < 1.0e-5 + + pymat = hubbard_structure.get_pymatgen_structure() + for i in range(24): + assert len(pymat.get_neighbors(pymat[i], r=radius)) in [4, 6] + + +def test_get_intersites_radius_five_nn(filepath_tests): + """Test the `HubbardUtils.get_intersites_radius` method against LiMnTe (5 neighbours).""" + path = os.path.join(filepath_tests, 'fixtures', 'structures', 'LMT.cif') + atoms = read(path) + + hubbard_structure = HubbardStructureData.from_structure(StructureData(ase=atoms)) + hubbard_structure.initialize_onsites_hubbard('Mn', '3d') + hubbard_structure.initialize_intersites_hubbard('Mn', '3d', 'Te', '2p') + hubbard_utils = HubbardUtils(hubbard_structure=hubbard_structure) + + radius = hubbard_utils.get_intersites_radius() + + pymat = hubbard_structure.get_pymatgen_structure() + assert len(pymat.get_neighbors(pymat[0], r=radius)) == 5 + + +def test_initialize_hubbard_parameters(get_non_trivial_hubbard_structure): + """Test the `HubbardUtils.initialize_hubbard_parameters` method.""" + structure = get_non_trivial_hubbard_structure() + hubbard_structure = initialize_hubbard_parameters(structure=structure, pairs={'Fe': ['3d', 5.0, 1e-8, {'O': '2p'}]}) + + assert len(hubbard_structure.hubbard.parameters) == 8 * (4 + 1) + 16 * (6 + 1) + + +def test_initialize_hubbard_parameters_five_nn(filepath_tests): + """Test the `HubbardUtils.initialize_hubbard_parameters` method against LiMnTe (5 neighbours).""" + path = os.path.join(filepath_tests, 'fixtures', 'structures', 'LMT.cif') + atoms = read(path) + + structure = StructureData(ase=atoms) + + hubbard_structure = initialize_hubbard_parameters( + structure=structure, pairs={'Mn': ['3d', 5.0, 1e-8, { + 'Te': '4p' + }]} + ) + + assert len(hubbard_structure.hubbard.parameters) == 6 + + +@pytest.mark.parametrize( + ('pairs', 'number_of_parameters'), + ( + ( + { + 'Mn': ['3d', 5.0, 1e-8, { + 'S': '3p' + }], + 'Co': ['3d', 5.0, 1e-8, { + 'S': '3p' + }], + }, + (1 + 6) + (1 + 4) # 2 onsites, 6 + 4 intersites + ), + ( + { + 'Mn': ['3d', 5.0, 1e-8, { + 'S': '3p', + 'Co': '3d' + }], + 'Co': ['3d', 5.0, 1e-8, { + 'S': '3p', + 'Mn': '3d' + }], + }, + (1 + 6 + 4) + (1 + 4 + 4) # 2 onsites, 6 + 4 TM-S, 4 + 4 TM-TM + ), + ) +) +def test_get_intersites_list(filepath_tests, pairs, number_of_parameters): + """Test the `HubbardUtils.get_intersites_list` method against MnCoS.""" + path = os.path.join(filepath_tests, 'fixtures', 'structures', 'MnCoS.cif') + atoms = read(path) + + structure = StructureData(ase=atoms) + hubbard_structure = initialize_hubbard_parameters(structure=structure, pairs=pairs) + + assert len(hubbard_structure.hubbard.parameters) == number_of_parameters + + hubbard_utils = HubbardUtils(hubbard_structure=hubbard_structure) + intersites = hubbard_utils.get_intersites_list() + + init_intersites = np.array(hubbard_structure.hubbard.to_list(), dtype='object')[:, [0, 2, 5]].tolist() + + assert len(intersites) == len(init_intersites) + + for intersite in intersites: + assert intersite in init_intersites + + +def test_get_max_number_of_neighbours(filepath_tests): + """Test the `HubbardUtiles.get_max_number_of_neighbours` method.""" + path = os.path.join(filepath_tests, 'fixtures', 'structures', 'MnCoS.cif') + atoms = read(path) + + pairs = { + 'Mn': ['3d', 5.0, 1e-8, { + 'S': '3p', + }], + 'Co': ['3d', 5.0, 1e-8, { + 'S': '3p', + }], + } + structure = StructureData(ase=atoms) + hubbard_structure = initialize_hubbard_parameters(structure=structure, pairs=pairs) + hubbard_utils = HubbardUtils(hubbard_structure=hubbard_structure) + + assert hubbard_utils.get_max_number_of_neighbours() == 10 + + +def test_max_number_of_neighbours(filepath_tests): + """Test the `max_number_of_neighbours` method.""" + from aiida_quantumespresso.utils.hubbard import max_number_of_neighbours + + array = [[0, 1], [0, 1], [0, 0], [0, 1], [0, 2], [0, 2]] + + assert max_number_of_neighbours(array) == 5 + + path = os.path.join(filepath_tests, 'fixtures', 'structures', 'MnCoS.cif') + atoms = read(path) + + pairs = { + 'Mn': ['3d', 5.0, 1e-8, { + 'S': '3p', + 'Co': '3d' + }], + 'Co': ['3d', 5.0, 1e-8, { + 'S': '3p', + 'Mn': '3d' + }], + } + structure = StructureData(ase=atoms) + hubbard_structure = initialize_hubbard_parameters(structure=structure, pairs=pairs) + + hubbard_utils = HubbardUtils(hubbard_structure=hubbard_structure) + intersites = np.array(hubbard_utils.get_intersites_list(), dtype='object')[:, [0, 1]] + + assert max_number_of_neighbours(intersites) == 10