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