"""
Utility code and functions for symmetry analysis of structures and defects.
"""
import contextlib
import math
import os
import warnings
from functools import lru_cache
from itertools import permutations
from typing import Optional, Union
import numpy as np
import pandas as pd
from pymatgen.analysis.defects.core import DefectType
from pymatgen.analysis.structure_matcher import ElementComparator, StructureMatcher
from pymatgen.core.operations import SymmOp
from pymatgen.core.structure import Lattice, PeriodicSite, Structure
from pymatgen.entries.computed_entries import ComputedStructureEntry
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer, SymmetryUndeterminedError
from pymatgen.transformations.standard_transformations import SupercellTransformation
from sympy import Eq, simplify, solve, symbols
from doped.core import DefectEntry
from doped.utils.parsing import (
_get_bulk_supercell,
_get_defect_supercell,
_get_defect_supercell_bulk_site_coords,
_get_unrelaxed_defect_structure,
get_site_mapping_indices,
)
[docs]
@lru_cache(maxsize=int(1e5))
def cached_simplify(eq):
"""
Cached simplification function for ``sympy`` equations, for efficiency.
"""
return simplify(eq)
[docs]
@lru_cache(maxsize=int(1e5))
def cached_solve(equation, variable):
"""
Cached solve function for ``sympy`` equations, for efficiency.
"""
return solve(equation, variable)
def _set_spglib_warnings_env_var():
"""
Set the SPGLIB environment variable to suppress spglib warnings.
"""
os.environ["SPGLIB_WARNING"] = "OFF"
def _check_spglib_version():
"""
Check the versions of spglib and its C libraries, and raise a warning if
the correct installation instructions have not been followed.
"""
import spglib
python_version = spglib.__version__
c_version = spglib.spg_get_version_full()
if python_version != c_version:
warnings.warn(
f"Your spglib Python version (spglib.__version__ = {python_version}) does not match its C "
f"library version (spglib.spg_get_version_full() = {c_version}). This can lead to unnecessary "
f"spglib warning messages, but can be avoided by upgrading spglib with `pip install --upgrade "
f"spglib`."
# No longer required as of spglib v2.5:
# f"- First uninstalling spglib with both `conda uninstall spglib` and `pip uninstall spglib` "
# f"(to ensure no duplicate installations).\n"
# f"- Then, install spglib with `conda install -c conda-forge spglib` or "
# f"`pip install git+https://github.com/spglib/spglib "
# f"--config-settings=cmake.define.SPGLIB_SHARED_LIBS=OFF` as detailed in the doped "
# f"installation instructions: https://doped.readthedocs.io/en/latest/Installation.html"
)
_check_spglib_version()
_set_spglib_warnings_env_var()
def _round_floats(obj, places: int = 5):
"""
Recursively round floats in a dictionary to ``places`` decimal places,
using the ``_custom_round`` function.
"""
if isinstance(obj, float):
return _custom_round(obj, places) + 0.0
if isinstance(obj, dict):
return {k: _round_floats(v, places) for k, v in obj.items()}
if isinstance(obj, (list, tuple)):
return [_round_floats(x, places) for x in obj]
if isinstance(obj, np.ndarray):
return _vectorized_custom_round(obj, places) + 0.0
if isinstance(obj, pd.DataFrame): # if dataframe, convert to dict and round floats
return pd.DataFrame(_round_floats(obj.to_dict(), places))
return obj
def _custom_round(number: float, decimals: int = 3):
"""
Custom rounding function that rounds numbers to a specified number of
decimals, if that rounded number is within 0.15*10^(-decimals) of the
original number, else rounds to [decimals+1] decimals.
Primarily because float rounding with pymatgen/numpy
can give cell coordinates of 0.5001 instead of 0.5
etc, but also can have coordinates of e.g. 0.6125
that should not be rounded to 0.613.
Args:
number (float): The number to round
decimals (int):
The number of decimals to round to (default: 3)
Returns:
float: The rounded number
"""
rounded_number = round(number, decimals)
if abs(rounded_number - number) < 0.15 * float(10) ** (-decimals):
return rounded_number
return round(number, decimals + 1)
_vectorized_custom_round = np.vectorize(_custom_round)
def _get_num_places_for_dist_precision(
structure: Union[Structure, Lattice], dist_precision: float = 0.001
):
"""
Given a structure or lattice, get the number of decimal places that we need
to keep / can round to for _fractional coordinates_ (``frac_coords``), to
maintain a distance precision of ``dist_precision`` in Å.
Intended for use with the ``_round_floats()`` function, to achieve cleanly
formatted structure outputs while ensuring no significant rounding errors are
introduced in site positions (e.g. for very large supercells, small differences
in fraction coordinates become significant).
Args:
structure:
The input structure or lattice.
dist_precision:
The desired distance precision in Å (default: 0.001).
Returns:
int:
The number of decimal places to keep for fractional
coordinates to maintain the desired distance precision.
"""
lattice = structure if isinstance(structure, Lattice) else structure.lattice
frac_precision = dist_precision / max(lattice.abc)
# get corresponding number of decimal places for this precision:
return -1 * min(math.floor(math.log(frac_precision, 10)), 0)
def _round_struct_coords(structure: Structure, dist_precision: float = 0.001, to_unit_cell=False):
"""
Convenience method to round the lattice parameters and fractional
coordinates of a structure to a given distance precision, for cleanly
formatted structure outputs.
Does not apply this operation in-place!
Args:
structure:
The input structure.
dist_precision:
The desired distance precision in Å (default: 0.001).
to_unit_cell:
Whether to round the fractional coordinates to the unit cell
(default: False).
Returns:
Structure:
The structure with rounded lattice parameters and fractional
coordinates.
"""
rounded_struct = structure.copy()
req_places = _get_num_places_for_dist_precision(rounded_struct, dist_precision)
frac_coords = _round_floats(rounded_struct.frac_coords, places=req_places)
lattice = Lattice(_round_floats(rounded_struct.lattice.matrix, places=req_places))
for idx in range(len(rounded_struct)):
orig_site = structure[idx]
rounded_struct._sites[idx] = PeriodicSite(
orig_site.species,
frac_coords[idx],
lattice,
properties=orig_site.properties,
label=orig_site.label,
skip_checks=True,
to_unit_cell=to_unit_cell,
)
return rounded_struct
def _frac_coords_sort_func(coords):
"""
Sorting function to apply on an iterable of fractional coordinates, where
entries are sorted by the number of x, y, z that are (almost) equal (i.e.
between 0 and 3), then by the magnitude of x+y+z, then by the magnitudes of
x, y and z.
"""
if coords is None:
return (1e10, 1e10, 1e10, 1e10, 1e10)
coords_for_sorting = _vectorized_custom_round(
np.mod(_vectorized_custom_round(coords), 1)
) # to unit cell
num_equals = sum(
np.isclose(coords_for_sorting[i], coords_for_sorting[j], atol=1e-3)
for i in range(len(coords_for_sorting))
for j in range(i + 1, len(coords_for_sorting))
)
magnitude = _custom_round(np.linalg.norm(coords_for_sorting))
return (-num_equals, magnitude, *np.abs(coords_for_sorting))
[docs]
def get_sga(struct, symprec=0.01):
"""
Get a ``SpacegroupAnalyzer`` object of the input structure, dynamically
adjusting symprec if needs be.
"""
sga = None
for trial_symprec in [symprec, 0.1, 0.001, 1, 0.0001]:
# if symmetry determination fails, increase symprec first, then decrease, then criss-cross
with contextlib.suppress(SymmetryUndeterminedError):
sga = SpacegroupAnalyzer(struct, symprec=trial_symprec)
if sga:
return sga
import spglib
raise SymmetryUndeterminedError(
f"Could not determine symmetry of input structure! Got spglib error: {spglib.get_error_message()}"
) # well shiiii...
[docs]
def apply_symm_op_to_site(
symm_op: SymmOp,
site: PeriodicSite,
fractional: bool = False,
rotate_lattice: Union[Lattice, bool] = True,
) -> PeriodicSite:
"""
Apply the given symmetry operation to the input site (**not in place**) and
return the new site.
By default, also rotates the lattice accordingly. If you want to apply
the symmetry operation but keep the same lattice definition, set
``rotate_lattice=False``.
Args:
symm_op: ``pymatgen`` ``SymmOp`` object.
site: ``pymatgen`` ``PeriodicSite`` object.
fractional:
If the ``SymmOp`` is in fractional or Cartesian (default)
coordinates (i.e. to apply to ``site.frac_coords`` or
``site.coords``). Default: False
rotate_lattice:
Either a ``pymatgen`` ``Lattice`` object (to use as the new
lattice basis of the transformed site, which can be provided
to reduce computation time when looping) or ``True/False``.
If ``True`` (default), the ``SymmOp`` rotation matrix will be
applied to the input site lattice, or if ``False``, the
original lattice will be retained.
Returns:
``pymatgen`` ``PeriodicSite`` object with the symmetry operation
applied
"""
if isinstance(rotate_lattice, Lattice):
rotated_lattice = rotate_lattice
else:
if rotate_lattice:
if fractional:
rotated_lattice = Lattice(np.dot(symm_op.rotation_matrix, site.lattice.matrix))
else:
rotated_lattice = Lattice(
[symm_op.apply_rotation_only(row) for row in site.lattice.matrix]
)
else:
rotated_lattice = site.lattice
if fractional: # operate in **original** lattice, then convert to new lattice
frac_coords = symm_op.operate(site.frac_coords)
new_coords = site.lattice.get_cartesian_coords(frac_coords)
else:
new_coords = symm_op.operate(site.coords)
return PeriodicSite(
site.species,
new_coords,
rotated_lattice,
coords_are_cartesian=True,
properties=site.properties,
skip_checks=True,
label=site.label,
)
[docs]
def apply_symm_op_to_struct(
symm_op: SymmOp, struct: Structure, fractional: bool = False, rotate_lattice: bool = True
) -> Structure:
"""
Apply a symmetry operation to a structure and return the new structure.
This differs from pymatgen's ``apply_operation`` method in that it
**does not apply the operation in place as well (i.e. does not modify
the input structure)**, which avoids the use of unnecessary and slow
``Structure.copy()`` calls, making the structure manipulation / symmetry
analysis functions more efficient.
Also fixes an issue when applying fractional symmetry operations.
By default, also rotates the lattice accordingly. If you want to apply
the symmetry operation to the sites but keep the same lattice definition,
set ``rotate_lattice=False``.
Args:
symm_op: ``pymatgen`` ``SymmOp`` object.
struct: ``pymatgen`` ``Structure`` object.
fractional:
If the ``SymmOp`` is in fractional or Cartesian (default)
coordinates (i.e. to apply to ``site.frac_coords`` or
``site.coords``). Default: False
rotate_lattice:
If the lattice of the input structure should be rotated
according to the symmetry operation. Default: True.
Returns:
``pymatgen`` ``Structure`` object with the symmetry operation
applied.
"""
# using modified version of ``pymatgen``\'s ``apply_operation`` method:
if rotate_lattice:
if not fractional:
rotated_lattice = Lattice([symm_op.apply_rotation_only(row) for row in struct._lattice.matrix])
else:
rotated_lattice = Lattice(np.dot(symm_op.rotation_matrix, struct._lattice.matrix))
else:
rotated_lattice = struct._lattice
# note could also use ``SymmOp.operate_multi`` for speedup if ever necessary, but requires some more
# accounting of species ordering etc, and this isn't an efficiency bottleneck currently
return Structure.from_sites(
[
apply_symm_op_to_site(symm_op, site, fractional=fractional, rotate_lattice=rotated_lattice)
for site in struct
]
)
[docs]
def summed_rms_dist(struct_a: Structure, struct_b: Structure) -> float:
"""
Get the summed root-mean-square (RMS) distance between the sites of two
structures, in Å.
Note that this assumes the lattices of the two structures are equal!
Args:
struct_a: ``pymatgen`` ``Structure`` object.
struct_b: ``pymatgen`` ``Structure`` object.
Returns:
float:
The summed RMS distance between the sites of the two structures,
in Å.
"""
# orders of magnitude faster than StructureMatcher.get_rms_dist() from pymatgen
# (though this assumes lattices are equal)
# set threshold to a large number to avoid possible site-matching warnings
return sum(get_site_mapping_indices(struct_a, struct_b, threshold=1e10, dists_only=True))
def _get_all_equiv_sites(frac_coords, struct, symm_ops=None, symprec=0.01, dist_tol=0.01):
"""
Get all equivalent sites of the input fractional coordinates in struct.
"""
if symm_ops is None:
sga = get_sga(struct, symprec=symprec)
symm_ops = sga.get_symmetry_operations() # fractional symm_ops by default
dummy_site = PeriodicSite("X", frac_coords, struct.lattice)
x_sites = []
for symm_op in symm_ops:
x_site = apply_symm_op_to_site(
symm_op, dummy_site, fractional=True, rotate_lattice=struct.lattice
).to_unit_cell() # enforce same lattice, as we just want transformed frac coords here
# if distance is >dist_tol for all other sites in x_sites, add x_site to x_sites:
if (
not x_sites
or struct.lattice.get_all_distances(
np.array([site.frac_coords for site in x_sites]), x_site.frac_coords
).min()
> dist_tol
):
x_sites.append(x_site)
return x_sites
def _get_symm_dataset_of_struc_with_all_equiv_sites(
frac_coords, struct, symm_ops=None, symprec=0.01, dist_tol=0.01
):
unique_sites = _get_all_equiv_sites(frac_coords, struct, symm_ops, dist_tol=dist_tol)
struct_with_all_X = _get_struct_with_all_X(struct, unique_sites)
sga_with_all_X = get_sga(struct_with_all_X, symprec=symprec)
return sga_with_all_X.get_symmetry_dataset(), unique_sites
def _get_struct_with_all_X(struct, unique_sites):
"""
Add all sites in unique_sites to a ``copy`` of ``struct``, and return this
new ``Structure``.
"""
struct_with_all_X = struct.copy()
struct_with_all_X.sites += unique_sites
return struct_with_all_X
def _get_equiv_frac_coords_in_primitive(
frac_coords, supercell, primitive, symm_ops=None, equiv_coords=True, dist_tol=0.01
):
"""
Get an equivalent fractional coordinates of ``frac_coords`` from a
supercell, in the primitive cell.
Also returns a list of equivalent fractional coords in the primitive cell
if ``equiv_coords`` is ``True``.
"""
unique_sites = _get_all_equiv_sites(frac_coords, supercell, symm_ops, dist_tol=dist_tol)
supercell_with_all_X = _get_struct_with_all_X(supercell, unique_sites)
prim_with_all_X = get_primitive_structure(supercell_with_all_X, ignored_species=["X"])
# ensure matched to primitive structure:
rotated_struct, matrix = _rotate_and_get_supercell_matrix(prim_with_all_X, primitive)
if rotated_struct is None:
warnings.warn("Could not find a mapping between the primitive and supercell structures!")
return None, None if equiv_coords else None
primitive_with_all_X = rotated_struct * matrix
sm = StructureMatcher(primitive_cell=False, ignored_species=["X"], comparator=ElementComparator())
s2_like_s1 = sm.get_s2_like_s1(primitive, primitive_with_all_X)
s2_really_like_s1 = Structure.from_sites(
[ # sometimes this get_s2_like_s1 doesn't work properly due to different (but equivalent) lattice
PeriodicSite( # vectors (e.g. a=(010) instead of (100) etc.), so do this to be sure
site.specie,
site.frac_coords,
primitive.lattice,
to_unit_cell=True,
)
for site in s2_like_s1.sites
]
)
prim_coord_list = [
_vectorized_custom_round(np.mod(_vectorized_custom_round(site.frac_coords), 1))
for site in s2_really_like_s1.sites
if site.specie.symbol == "X"
]
if equiv_coords:
return ( # sort with _frac_coords_sort_func
sorted(prim_coord_list, key=_frac_coords_sort_func)[0],
prim_coord_list,
)
return sorted(prim_coord_list, key=_frac_coords_sort_func)[0]
def _rotate_and_get_supercell_matrix(prim_struct, target_struct):
"""
Rotates the input prim_struct to match the target_struct orientation, and
returns the supercell matrix to convert from the rotated prim_struct to the
target_struct.
Returns ``(None, None)`` if no mapping is found.
"""
possible_mappings = list(prim_struct.lattice.find_all_mappings(target_struct.lattice))
if not possible_mappings:
return None, None
mapping = next(
iter( # get possible mappings, then sort by R*S, S, R, then return first
sorted(
possible_mappings,
key=lambda x: (
_lattice_matrix_sort_func(np.dot(x[1].T, x[2])),
_lattice_matrix_sort_func(x[2]),
_lattice_matrix_sort_func(x[1]),
),
)
)
)
rotation_matrix = mapping[1]
if np.allclose(rotation_matrix, -1 * np.eye(3)):
# pymatgen sometimes gives a rotation matrix of -1 * identity matrix, which is
# equivalent to no rotation. Just use the identity matrix instead.
rotation_matrix = np.eye(3)
supercell_matrix = -1 * mapping[2]
else:
supercell_matrix = mapping[2]
rotation_symm_op = SymmOp.from_rotation_and_translation(
rotation_matrix=rotation_matrix.T
) # Transpose = inverse of rotation matrices (orthogonal matrices), better numerical stability
output_prim_struct = apply_symm_op_to_struct(rotation_symm_op, prim_struct, rotate_lattice=True)
return _round_struct_coords(output_prim_struct), supercell_matrix
[docs]
def translate_structure(
structure: Structure, vector: np.ndarray, frac_coords: bool = True, to_unit_cell: bool = True
) -> Structure:
"""
Translate a structure and its sites by a given vector (**not in place**).
Args:
structure: ``pymatgen`` ``Structure`` object.
vector: Translation vector, fractional or Cartesian.
frac_coords:
Whether the input vector is in fractional coordinates.
(Default: True)
to_unit_cell:
Whether to translate the sites to the unit cell.
(Default: True)
Returns:
``pymatgen`` ``Structure`` object with translated sites.
"""
translated_structure = structure.copy()
return translated_structure.translate_sites(
indices=list(range(len(translated_structure))),
vector=vector,
to_unit_cell=to_unit_cell,
frac_coords=frac_coords,
)
def _get_supercell_matrix_and_possibly_redefine_prim(
prim_struct, target_struct, sga: Optional[SpacegroupAnalyzer] = None, symprec=0.01
):
"""
Determines the supercell transformation matrix to convert from the
primitive structure to the target structure.
The supercell matrix is defined to be T in ``T*P = S`` where P and S
are the primitive and supercell lattice matrices respectively.
Equivalently, multiplying ``prim_struct * T`` will give the target_struct.
In ``pymatgen``, this requires the output transformation matrix to be
integer.
First tries to determine a simple (integer) transformation matrix with no
basis set rotation required. If that fails, then defaults to using
``_rotate_and_get_supercell_matrix``.
Searches over various possible primitive cell definitions from spglib.
Args:
prim_struct: ``pymatgen`` ``Structure`` object of the primitive cell.
target_struct: ``pymatgen`` ``Structure`` object of the target cell.
sga:
``SpacegroupAnalyzer`` object of the primitive cell. If None,
will be computed from prim_struct.
symprec:
Symmetry precision for ``SpacegroupAnalyzer``, if being generated.
Returns:
prim_struct:
Primitive structure, possibly rotated/redefined.
supercell_matrix:
Supercell transformation matrix to convert from the
primitive structure to the target structure.
"""
def _get_supercell_matrix_and_possibly_rotate_prim(prim_struct, target_struct):
try:
# supercell transform matrix is T in `T*P = S` (P = prim, S = super), so `T = S*P^-1`:
transformation_matrix = np.dot(
target_struct.lattice.matrix, np.linalg.inv(prim_struct.lattice.matrix)
)
if not np.allclose(np.rint(transformation_matrix), transformation_matrix, atol=1e-3):
raise ValueError # if non-integer transformation matrix
return prim_struct, np.rint(transformation_matrix)
except ValueError: # if non-integer transformation matrix
attempt_prim_struct, attempt_transformation_matrix = _rotate_and_get_supercell_matrix(
prim_struct, target_struct
)
if attempt_prim_struct: # otherwise failed, stick with original T matrix
prim_struct = attempt_prim_struct
transformation_matrix = attempt_transformation_matrix
if np.allclose(np.rint(transformation_matrix), transformation_matrix, atol=1e-3):
return prim_struct, np.rint(transformation_matrix)
return prim_struct, transformation_matrix
sga = sga or get_sga(prim_struct, symprec=symprec)
rms_dists_w_candidate_prim_structs_and_T_matrices = []
# Could also apply possible origin shifts to other structs (refined, find_primitive) as well,
# if we find any structures for which this code still fails
candidate_prim_structs = [
*_get_candidate_prim_structs(prim_struct, symprec=symprec),
*_get_candidate_prim_structs(target_struct, symprec=symprec),
]
for possible_prim_struct in candidate_prim_structs:
new_prim_struct, transformation_matrix = _get_supercell_matrix_and_possibly_rotate_prim(
possible_prim_struct, target_struct
)
if not np.allclose(
np.rint(transformation_matrix), transformation_matrix, atol=1e-3
) or not np.allclose(
(new_prim_struct * transformation_matrix).lattice.matrix,
target_struct.lattice.matrix,
atol=1e-3,
):
# not integer or doesn't exactly match bulk supercell, so bad transformation matrix, skip
continue
new_prim_struct = Structure.from_sites([site.to_unit_cell() for site in new_prim_struct])
rms_dist_to_target = summed_rms_dist(
Structure.from_sites(
[
site.to_unit_cell()
for site in (new_prim_struct * transformation_matrix).get_sorted_structure()
]
),
target_struct,
)
rms_dists_w_candidate_prim_structs_and_T_matrices.append(
(rms_dist_to_target, new_prim_struct, transformation_matrix)
)
closest_match = sorted( # sort to get ideal primitive cell definition
rms_dists_w_candidate_prim_structs_and_T_matrices,
key=lambda x: (
round(x[0], 3),
_lattice_matrix_sort_func(x[1].lattice.matrix),
_lattice_matrix_sort_func(x[2]),
_struct_sort_func(x[1]),
),
)[0]
if closest_match[0] > 0.1: # no perfect match has been found. Warn user and return the closest:
warnings.warn(
f"Found the transformation matrix from the primitive cell lattice to the supplied supercell, "
f"but could not determine the transformation to directly match the atomic coordinates ("
f"infinite possible symmetry-equivalent coordinate definitions). Closest match has RMS "
f"distance of {closest_match[0]:.3f} Å.\n"
f"The bulk and defect supercells generated will be equivalent to the input supercell, "
f"but with a different choice of atomic coordinates (e.g. [0.1, 0.1, 0.1] instead of [0.9, "
f"0.9, 0.9]). You should make sure to do the bulk supercell calculation with this "
f"doped-generated supercell (DefectsGenerator.bulk_supercell, which is output to the `Bulk` "
f"folders with the file generation functions), so that the coordinates match those of the "
f"defect supercells (this matters when computing finite-size corrections)."
)
# Note: Could always just get the transformation of the generated supercell to the input supercell, and
# then apply this transformation to each generated bulk/defect supercell at the end of defect
# generation, but means that self.primitive_structure * self.supercell_matrix is no longer
# guaranteed to match self.bulk_supercell... Not the biggest deal though
# Likely way more work than worth
return closest_match[1:]
def _get_candidate_prim_structs(structure, **kwargs):
sga = get_sga(structure, **kwargs)
pmg_prim_struct = structure.get_primitive_structure()
candidate_prim_structs = (
[structure, pmg_prim_struct] if len(structure) == len(pmg_prim_struct) else [pmg_prim_struct]
)
for _i in range(4):
struct = sga.get_primitive_standard_structure()
candidate_prim_structs.append(struct)
spglib_dataset = sga.get_symmetry_dataset()
if not np.allclose(spglib_dataset.origin_shift, 0):
candidate_prim_structs.append(translate_structure(struct, spglib_dataset.origin_shift))
sga = get_sga(struct, sga._symprec) # use same symprec
candidate_prim_structs.append(sga.find_primitive())
for candidate_conv_struct in [sga.get_refined_structure(), sga.get_conventional_standard_structure()]:
if len(candidate_conv_struct) == len(pmg_prim_struct):
# only also try conventional if equivalent to the primitive cell
candidate_prim_structs = [candidate_conv_struct, *candidate_prim_structs]
return candidate_prim_structs
[docs]
def get_wyckoff(frac_coords, struct, symm_ops: Optional[list] = None, equiv_sites=False, symprec=0.01):
"""
Get the Wyckoff label of the input fractional coordinates in the input
structure. If the symmetry operations of the structure have already been
computed, these can be input as a list to speed up the calculation.
Args:
frac_coords:
Fractional coordinates of the site to get the Wyckoff label of.
struct:
pymatgen Structure object for which frac_coords corresponds to.
symm_ops:
List of pymatgen SymmOps of the structure. If None (default),
will recompute these from the input struct.
equiv_sites:
If True, also returns a list of equivalent sites in struct.
symprec:
Symmetry precision for SpacegroupAnalyzer.
"""
symm_dataset, unique_sites = _get_symm_dataset_of_struc_with_all_equiv_sites(
frac_coords, struct, symm_ops, symprec=symprec
)
conv_cell_factor = len(symm_dataset.std_positions) / len(symm_dataset.wyckoffs)
multiplicity = int(conv_cell_factor * len(unique_sites))
wyckoff_label = f"{multiplicity}{symm_dataset.wyckoffs[-1]}"
return (wyckoff_label, unique_sites) if equiv_sites else wyckoff_label
def _struct_sort_func(struct: Union[Structure, np.ndarray]) -> tuple:
"""
Sort by the lattice matrix sorting function, then by (minus) the number of
high-symmetry coordinates (x=y=z, then 2 equal coordinates), then by the
sum of all fractional coordinates, then by the magnitudes of high-symmetry
coordinates (x=y=z, then 2 equal coordinates), then by the summed magnitude
of all x coordinates, then y coordinates, then z coordinates.
Args:
struct:
``pymatgen`` ``Structure`` object, or an array of fractional
coordinates of sites in the structure (in which case the lattice
matrix metric is skipped).
Returns:
tuple: Tuple of sorting criteria values.
"""
if isinstance(struct, Structure):
struct_for_sorting = _round_struct_coords(struct, to_unit_cell=True)
lattice_metric = _lattice_matrix_sort_func(struct_for_sorting.lattice.matrix)
frac_coords = struct_for_sorting.frac_coords
else:
lattice_metric = (False,)
frac_coords = struct
# get summed magnitudes of x=y=z coords:
xyz_matching_coords = frac_coords[ # Find the coordinates where x = y = z:
(frac_coords[:, 0] == frac_coords[:, 1]) & (frac_coords[:, 1] == frac_coords[:, 2])
]
xyz_sum_magnitudes = np.sum(np.linalg.norm(xyz_matching_coords, axis=1))
# get summed magnitudes of x=y / y=z / x=z coords:
xy_matching_coords = frac_coords[
(frac_coords[:, 0] == frac_coords[:, 1])
| (frac_coords[:, 1] == frac_coords[:, 2])
| (frac_coords[:, 0] == frac_coords[:, 2])
]
xy_sum_magnitudes = np.sum(np.linalg.norm(xy_matching_coords, axis=1))
return (
*lattice_metric,
-len(xyz_matching_coords),
-len(xy_matching_coords),
round(np.sum(frac_coords), 2),
round(xyz_sum_magnitudes, 2),
round(xy_sum_magnitudes, 2),
round(np.sum(frac_coords[:, 0]), 2),
round(np.sum(frac_coords[:, 1]), 2),
round(np.sum(frac_coords[:, 2]), 2),
)
def _lattice_matrix_sort_func(lattice_matrix: np.ndarray) -> tuple:
"""
Sorting function to apply on an iterable of lattice matrices.
Matrices are sorted by:
- lattice_matrix is diagonal
- matrix symmetry (around diagonal)
- maximum sum of diagonal element magnitudes.
- minimum number of negative elements
- maximum number of x, y, z that are equal
- maximum number of abs(x), abs(y), abs(z) that are equal
- a, b, c magnitudes (favouring c >= b >= a)
Args:
lattice_matrix (np.ndarray): Lattice matrix to sort.
Returns:
tuple: Tuple of sorting criteria values
"""
is_diagonal = np.allclose(lattice_matrix, np.diag(np.diag(lattice_matrix)))
symmetric = np.allclose(lattice_matrix, lattice_matrix.T)
num_negs = np.sum(lattice_matrix < 0)
diag_sum = np.round(np.sum(np.abs(np.diag(lattice_matrix))), 1)
flat_matrix = lattice_matrix.flatten()
num_equals = sum(
flat_matrix[i] == flat_matrix[j]
for i in range(len(flat_matrix))
for j in range(i, len(flat_matrix))
) # double (square) counting, but doesn't matter (sorting behaviour the same)
num_abs_equals = sum(
abs(flat_matrix[i]) == abs(flat_matrix[j])
for i in range(len(flat_matrix))
for j in range(i, len(flat_matrix))
) # double (square) counting, but doesn't matter (sorting behaviour the same)
a, b, c = np.linalg.norm(lattice_matrix, axis=1)
return (
not is_diagonal,
not symmetric,
-diag_sum,
num_negs,
-num_equals,
-num_abs_equals,
-round(c, 2),
-round(b, 2),
-round(a, 2),
)
[docs]
def get_clean_structure(
structure: Structure, return_T: bool = False, dist_precision: float = 0.001, niggli_reduce: bool = True
):
"""
Get a 'clean' version of the input `structure` by searching over equivalent
cells, and finding the most optimal according to
``_lattice_matrix_sort_func`` (most symmetric, with mostly positive
diagonals and c >= b >= a).
Args:
structure (Structure): Structure object.
return_T (bool):
Whether to return the transformation matrix from the original
structure lattice to the new structure lattice (T * Orig = New).
(Default = False)
dist_precision (float):
The desired distance precision in Å for rounding of lattice
parameters and fractional coordinates. (Default: 0.001)
niggli_reduce (bool):
Whether to Niggli reduce the lattice before searching for the
optimal lattice matrix. If this is set to ``False``, we also
skip the search for the best positive determinant lattice matrix.
(Default: True)
"""
lattice = structure.lattice
if np.all(lattice.matrix <= 0):
lattice = Lattice(lattice.matrix * -1)
possible_lattice_matrices = [
lattice.matrix,
]
for _ in range(4):
lattice = lattice.get_niggli_reduced_lattice() if niggli_reduce else lattice
# want to maximise the number of non-negative diagonals, and also have a positive determinant
# can multiply two rows by -1 to get a positive determinant:
possible_lattice_matrices.append(lattice.matrix)
for i in range(3):
for j in range(i + 1, 3):
new_lattice_matrix = lattice.matrix.copy()
new_lattice_matrix[i] = new_lattice_matrix[i] * -1
new_lattice_matrix[j] = new_lattice_matrix[j] * -1
possible_lattice_matrices.append(new_lattice_matrix)
possible_lattice_matrices.sort(key=_lattice_matrix_sort_func)
new_lattice_matrix = possible_lattice_matrices[0]
if np.all(new_lattice_matrix <= 0):
new_lattice_matrix = new_lattice_matrix * -1
new_structure = Structure(
new_lattice_matrix,
structure.species_and_occu,
structure.cart_coords, # type: ignore
coords_are_cartesian=True,
to_unit_cell=True,
site_properties=structure.site_properties,
labels=structure.labels,
charge=structure._charge,
)
new_structure = _round_struct_coords(new_structure, dist_precision=dist_precision, to_unit_cell=True)
# sort structure to match a desired, deterministic format:
new_structure = new_structure.get_sorted_structure(
key=lambda x: (
x.species.average_electroneg,
x.species_string,
_frac_coords_sort_func(x.frac_coords),
)
)
if niggli_reduce:
new_structure = _get_best_pos_det_structure(new_structure) # ensure positive determinant
if return_T:
# T * Orig = New; T = New * Orig^-1; Orig = T^-1 * New
transformation_matrix = np.matmul(
new_structure.lattice.matrix, np.linalg.inv(structure.lattice.matrix)
)
if not np.allclose(transformation_matrix, np.rint(transformation_matrix), atol=1e-5):
raise ValueError(
"Transformation matrix for clean/reduced structure could not be found! If you are seeing "
"this bug, please notify the `doped` developers"
)
return (new_structure, np.rint(transformation_matrix))
return new_structure
def _get_best_pos_det_structure(structure: Structure):
"""
If the input structure has a negative determinant (corresponding to a left-
hand coordinate system), then find the best possible re-definition of the
lattice vectors which gives a positive determinant, according to
``_struct_sort_func``.
This is to avoid an apparent VASP bug with negative triple products of the
lattice vectors -- not sure if this is only in old versions?
"""
if np.linalg.det(structure.lattice.matrix) < 0:
swap_combo_score_dict = {}
for swap_combo in permutations([0, 1, 2]):
candidate_structure = swap_axes(structure, swap_combo)
if np.linalg.det(candidate_structure.lattice.matrix) > 0:
swap_combo_score_dict[swap_combo] = _struct_sort_func(candidate_structure)
best_swap_combo = min(swap_combo_score_dict, key=lambda x: swap_combo_score_dict[x])
structure = swap_axes(structure, best_swap_combo)
return structure
[docs]
def get_primitive_structure(
structure: Structure,
ignored_species: Optional[list] = None,
clean: bool = True,
return_all: bool = False,
**kwargs,
):
"""
Get a consistent/deterministic primitive structure from a ``pymatgen``
``Structure``.
For some materials (e.g. zinc blende), there are multiple equivalent
primitive cells (e.g. Cd (0,0,0) & Te (0.25,0.25,0.25); Cd (0,0,0) &
Te (0.75,0.75,0.75) for F-43m CdTe), so for reproducibility and in
line with most structure conventions/definitions, take the one with
the cleanest lattice and structure definition, according to
``struct_sort_func``.
If ``ignored_species`` is set, then the sorting function used to determine
the ideal primitive structure will ignore sites with species in
``ignored_species``.
Args:
structure:
``Structure`` object to get the corresponding primitive structure of.
ignored_species:
List of species to ignore when determining the ideal primitive
structure. (Default: None)
clean:
Whether to return a 'clean' version of the primitive structure,
with the lattice matrix in a standardised form. (Default: True)
return_all:
Whether to return all possible primitive structures tested, sorted
by the sorting function. (Default: False)
**kwargs:
Additional keyword arguments to pass to the ``get_sga`` function
(e.g. ``symprec`` etc).
"""
candidate_prim_structs = _get_candidate_prim_structs(structure, **kwargs)
if ignored_species is not None:
pruned_possible_prim_structs = [
Structure.from_sites([site for site in struct if site.specie.symbol not in ignored_species])
for struct in candidate_prim_structs
]
else:
pruned_possible_prim_structs = candidate_prim_structs
# sort and return indices:
sorted_indices = sorted(
range(len(pruned_possible_prim_structs)),
key=lambda i: _struct_sort_func(pruned_possible_prim_structs[i]),
)
prim_structs = [
_get_best_pos_det_structure(_round_struct_coords(candidate_prim_structs[i], to_unit_cell=True))
for i in sorted_indices
]
if clean:
prim_structs = [get_clean_structure(struct) for struct in prim_structs]
return prim_structs if return_all else _get_best_pos_det_structure(prim_structs[0])
[docs]
def get_spglib_conv_structure(sga):
"""
Get a consistent/deterministic conventional structure from a
``SpacegroupAnalyzer`` object. Also returns the corresponding
``SpacegroupAnalyzer`` (for getting Wyckoff symbols corresponding to this
conventional structure definition).
For some materials (e.g. zinc blende), there are multiple equivalent
primitive/conventional cells, so for reproducibility and in line with most
structure conventions/definitions, take the one with the lowest summed norm
of the fractional coordinates of the sites (i.e. favour Cd (0,0,0) and Te
(0.25,0.25,0.25) over Cd (0,0,0) and Te (0.75,0.75,0.75) for F-43m CdTe;
SGN 216).
"""
possible_conv_structs_and_sgas = []
for _i in range(4):
struct = sga.get_conventional_standard_structure()
possible_conv_structs_and_sgas.append((struct, sga))
sga = get_sga(sga.get_primitive_standard_structure(), symprec=sga._symprec)
possible_conv_structs_and_sgas = sorted(
possible_conv_structs_and_sgas, key=lambda x: _struct_sort_func(x[0])
)
return (
_round_struct_coords(possible_conv_structs_and_sgas[0][0], to_unit_cell=True),
possible_conv_structs_and_sgas[0][1],
)
[docs]
def get_BCS_conventional_structure(structure, pbar=None, return_wyckoff_dict=False):
"""
Get the conventional crystal structure of the input structure, according to
the Bilbao Crystallographic Server (BCS) definition. Also returns the
transformation matrix from the ``spglib`` (``SpaceGroupAnalyzer``)
conventional structure definition to the BCS definition.
Args:
structure (Structure): pymatgen Structure object for this to
get the corresponding BCS conventional crystal structure
pbar (ProgressBar): tqdm progress bar object, to update progress.
return_wyckoff_dict (bool): whether to return the Wyckoff label
dict ({Wyckoff label: coordinates})
number.
Returns:
pymatgen Structure object and ``spglib`` -> BCS conv cell transformation matrix.
"""
struc_wout_oxi = structure.copy()
struc_wout_oxi.remove_oxidation_states()
sga = get_sga(struc_wout_oxi)
conventional_structure, conv_sga = get_spglib_conv_structure(sga)
wyckoff_label_dict = get_wyckoff_dict_from_sgn(conv_sga.get_space_group_number())
# determine cell orientation for Wyckoff site determination (needs to match the Bilbao
# Crystallographic Server's convention, which can differ from spglib (pymatgen) in some cases)
sga_wyckoffs = conv_sga.get_symmetrized_structure().wyckoff_symbols
for trial_lattice_vec_swap_array in [ # 3C2 -> 6 possible combinations
# ordered according to frequency of occurrence in the Materials Project
[0, 1, 2], # abc, ~95% of cases
[0, 2, 1], # acb
[2, 1, 0], # cba
[1, 0, 2], # bac
[2, 0, 1], # cab
[1, 2, 0], # bca
None, # no perfect match, default to original orientation
]:
if trial_lattice_vec_swap_array is None:
lattice_vec_swap_array = [0, 1, 2]
break
reoriented_conv_structure = swap_axes(conventional_structure, trial_lattice_vec_swap_array)
if _compare_wyckoffs(
sga_wyckoffs,
reoriented_conv_structure,
wyckoff_label_dict,
):
lattice_vec_swap_array = trial_lattice_vec_swap_array
break
if pbar is not None:
pbar.update(1 / 6 * 10) # 45 up to 55% of progress bar in DefectsGenerator. This part can
# take a little while for low-symmetry structures
bcs_conv_structure = get_clean_structure(
swap_axes(conventional_structure, lattice_vec_swap_array), niggli_reduce=False
)
if return_wyckoff_dict:
return bcs_conv_structure, lattice_vec_swap_array, wyckoff_label_dict
return bcs_conv_structure, lattice_vec_swap_array
[docs]
def get_conv_cell_site(defect_entry):
"""
Gets an equivalent site of the defect entry in the conventional structure
of the host material. If the ``conventional_structure`` attribute is not
defined for defect_entry, then it is generated using ``SpacegroupAnalyzer``
and then reoriented to match the Bilbao Crystallographic Server's
conventional structure definition.
Args:
defect_entry: ``DefectEntry`` object.
"""
bulk_prim_structure = defect_entry.defect.structure.copy()
bulk_prim_structure.remove_oxidation_states() # adding oxidation states adds the
# # deprecated 'properties' attribute with -> {"spin": None}, giving a deprecation warning
prim_struct_with_X = bulk_prim_structure.copy()
prim_struct_with_X.append("X", defect_entry.defect.site.frac_coords, coords_are_cartesian=False)
sga = get_sga(bulk_prim_structure)
# convert to match sga primitive structure first:
sm = StructureMatcher(primitive_cell=False, ignored_species=["X"], comparator=ElementComparator())
sga_prim_struct = sga.get_primitive_standard_structure()
s2_like_s1 = sm.get_s2_like_s1(sga_prim_struct, prim_struct_with_X)
if not s2_like_s1:
warnings.warn(
"The transformation from the DefectEntry primitive cell to the spglib primitive "
"cell could not be determined, and so the corresponding conventional cell site "
"cannot be identified."
)
return None
s2_really_like_s1 = Structure.from_sites(
[ # sometimes this get_s2_like_s1 doesn't work properly due to different (but equivalent) lattice
PeriodicSite( # vectors (e.g. a=(010) instead of (100) etc.), so do this to be sure
site.specie,
site.frac_coords,
sga_prim_struct.lattice,
to_unit_cell=True,
)
for site in s2_like_s1.sites
]
)
conv_struct_with_X = s2_really_like_s1 * np.linalg.inv(
sga.get_conventional_to_primitive_transformation_matrix()
)
# convert to match defect_entry conventional structure definition
s2_like_s1 = sm.get_s2_like_s1(defect_entry.conventional_structure, conv_struct_with_X)
s2_really_like_s1 = Structure.from_sites(
[ # sometimes this get_s2_like_s1 doesn't work properly due to different (but equivalent) lattice
PeriodicSite( # vectors (e.g. a=(010) instead of (100) etc.), so do this to be sure
site.specie,
site.frac_coords,
defect_entry.conventional_structure.lattice,
to_unit_cell=True,
)
for site in s2_like_s1.sites
]
)
conv_cell_site = next(site for site in s2_really_like_s1.sites if site.specie.symbol == "X")
# site choice doesn't matter so much here, as we later get the equivalent coordinates using the
# Wyckoff dict and choose the conventional site based on that anyway (in the DefectsGenerator
# initialisation)
conv_cell_site.to_unit_cell()
conv_cell_site.frac_coords = _vectorized_custom_round(conv_cell_site.frac_coords)
return conv_cell_site
[docs]
def swap_axes(structure, axes):
"""
Swap axes of the given structure.
The new order of the axes is given by the axes parameter. For example,
``axes=(2, 1, 0)`` will swap the first and third axes.
"""
transformation_matrix = [[0, 0, 0], [0, 0, 0], [0, 0, 0]]
for i, axis in enumerate(axes):
transformation_matrix[i][axis] = 1
transformation = SupercellTransformation(transformation_matrix)
return transformation.apply_transformation(structure)
[docs]
def get_wyckoff_dict_from_sgn(sgn):
"""
Get dictionary of {Wyckoff label: coordinates} for a given space group
number.
The database used here for Wyckoff analysis (``wyckpos.dat``) was obtained
from code written by JaeHwan Shim @schinavro (ORCID: 0000-0001-7575-4788)
(https://gitlab.com/ase/ase/-/merge_requests/1035) based on the tabulated
datasets in https://github.com/xtalopt/randSpg (also found at
https://github.com/spglib/spglib/blob/develop/database/Wyckoff.csv).
By default, doped uses the Wyckoff functionality of ``spglib`` (along with
symmetry operations in pymatgen) when possible however.
"""
datafile = _get_wyckoff_datafile()
with open(datafile, encoding="utf-8") as f:
wyckoff = _read_wyckoff_datafile(sgn, f)
wyckoff_label_coords_dict = {}
def _coord_string_to_array(coord_string):
# Split string into substrings, parse each as a sympy expression,
# then convert to list of sympy expressions
return [cached_simplify(x.replace("2x", "2*x")) for x in coord_string.split(",")]
for element in wyckoff["letters"]:
label = wyckoff[element]["multiplicity"] + element # e.g. 4d
wyckoff_coords = [_coord_string_to_array(coords) for coords in wyckoff[element]["coordinates"]]
wyckoff_label_coords_dict[label] = wyckoff_coords
equivalent_sites = [
_coord_string_to_array(coords) for coords in wyckoff.get("equivalent_sites", [])
]
new_coords = [] # new list for equivalent coordinates
for coord_array in wyckoff_coords:
for equivalent_site in equivalent_sites:
# add coord_array and equivalent_site element-wise
equiv_coord_array = coord_array.copy()
equiv_coord_array = equiv_coord_array + np.array(equivalent_site)
new_coords.append(equiv_coord_array)
# add new_coords to wyckoff_label_coords:
wyckoff_label_coords_dict[label].extend(new_coords)
return wyckoff_label_coords_dict
[docs]
def get_wyckoff_label_and_equiv_coord_list(
defect_entry=None, conv_cell_site=None, sgn=None, wyckoff_dict=None
):
"""
Return the Wyckoff label and list of equivalent fractional coordinates
within the conventional cell for the input defect_entry or conv_cell_site
(whichever is provided, defaults to defect_entry if both), given a
dictionary of Wyckoff labels and coordinates (``wyckoff_dict``).
If ``wyckoff_dict`` is not provided, it is generated from the spacegroup
number (sgn) using ``get_wyckoff_dict_from_sgn(sgn)``. If ``sgn`` is not
provided, it is obtained from the bulk structure of the ``defect_entry`` if
provided.
"""
if wyckoff_dict is None:
if sgn is None:
if defect_entry is None:
raise ValueError(
"If inputting `conv_cell_site` and not `defect_entry`, either `sgn` or `wyckoff_dict` "
"must be provided."
)
# get sgn from primitive unit cell of bulk structure:
sgn = get_sga(defect_entry.defect.structure).get_space_group_number()
wyckoff_dict = get_wyckoff_dict_from_sgn(sgn)
def _compare_arrays(coord_list, coord_array):
"""
Compare a list of arrays of sympy expressions (``coord_list``) with an
array of coordinates (``coord_array``).
Returns the matching array from the list.
"""
x, y, z = symbols("x y z")
variable_dicts = [{}] # list of dicts for x,y,z
for sympy_array in coord_list:
match, variable_dict = evaluate_expression_and_update_dict(
sympy_array, coord_array, variable_dicts
)
if match:
# return coord list with sympy expressions subbed with variable_dict:
return [
np.array(
[
np.mod(float(cached_simplify(sympy_expr).subs(variable_dict)), 1)
for sympy_expr in sympy_array
]
)
for sympy_array in coord_list
]
return None # No match found
# get match of coords in wyckoff_label_coords to defect site coords:
def find_closest_match(defect_site, wyckoff_label_coords_dict):
for label, coord_list in wyckoff_label_coords_dict.items():
subbed_coord_list = _compare_arrays(coord_list, np.array(defect_site.frac_coords))
if subbed_coord_list is not None:
# convert coords in subbed_coord_list to unit cell, by rounding to 5 decimal places and
# then modding by 1:
subbed_coord_list = [
_vectorized_custom_round(np.mod(_vectorized_custom_round(coord_array, 5), 1))
for coord_array in subbed_coord_list
]
return label, subbed_coord_list
return None # No match found
def evaluate_expression(sympy_expr, coord, variable_dict):
equation = Eq(sympy_expr, coord)
variable = next(iter(sympy_expr.free_symbols))
variable_dict[variable] = cached_solve(equation, variable)[0]
return cached_simplify(sympy_expr).subs(variable_dict)
def add_new_variable_dict(
sympy_expr_prepend, sympy_expr, coord, current_variable_dict, variable_dicts
):
new_sympy_expr = cached_simplify(sympy_expr_prepend + str(sympy_expr))
new_dict = current_variable_dict.copy()
evaluate_expression(new_sympy_expr, coord, new_dict) # solve for new variable
if new_dict not in variable_dicts:
variable_dicts.append(new_dict)
def evaluate_expression_and_update_dict(sympy_array, coord_array, variable_dicts):
temp_dict = {}
match = False
for variable_dict in variable_dicts:
temp_dict = variable_dict.copy()
match = True
# sort zipped arrays by number of variables in sympy expression:
coord_array, sympy_array = zip(
*sorted(zip(coord_array, sympy_array), key=lambda x: len(x[1].free_symbols))
)
for coord, sympy_expr in zip(coord_array, sympy_array):
# Evaluate the expression with the current variable_dict
expr_value = cached_simplify(sympy_expr).subs(temp_dict)
# If the expression cannot be evaluated to a float
# it means that there is a new variable in the expression
try:
expr_value = np.mod(float(expr_value), 1) # wrap to 0-1 (i.e. to unit cell)
except TypeError:
# Assign the expression the value of the corresponding coordinate, and solve
# for the new variable
# first, special cases with two possible solutions due to PBC:
if sympy_expr == cached_simplify("-2*x"):
add_new_variable_dict("1+", sympy_expr, coord, temp_dict, variable_dicts)
elif sympy_expr == cached_simplify("2*x"):
add_new_variable_dict("-1+", sympy_expr, coord, temp_dict, variable_dicts)
expr_value = evaluate_expression(
sympy_expr, coord, temp_dict
) # solve for new variable
# Check if the evaluated expression matches the corresponding coordinate
if not np.isclose(
np.mod(float(coord), 1), # wrap to 0-1 (i.e. to unit cell)
np.mod(float(expr_value), 1),
atol=0.003,
) and not np.isclose(
np.mod(float(coord), 1) - 1, # wrap to 0-1 (i.e. to unit cell)
np.mod(float(expr_value), 1),
atol=0.003,
):
match = False
break
if match:
break
return match, temp_dict
if defect_entry is not None:
defect_entry.defect.site.to_unit_cell() # ensure wrapped to unit cell
# convert defect site to conventional unit cell for Wyckoff label matching:
conv_cell_site = get_conv_cell_site(defect_entry)
return find_closest_match(conv_cell_site, wyckoff_dict)
def _compare_wyckoffs(wyckoff_symbols, conv_struct, wyckoff_dict):
"""
Compare the Wyckoff labels of a conventional structure to a list of Wyckoff
labels.
"""
def _multiply_wyckoffs(wyckoff_labels, n=2):
return [str(n * int(wyckoff[:-1])) + wyckoff[-1] for wyckoff in wyckoff_labels]
wyckoff_symbol_lists = [_multiply_wyckoffs(wyckoff_symbols, n=n) for n in range(1, 5)] # up to 4x
doped_wyckoffs = []
for site in conv_struct:
wyckoff_label, equiv_coords = get_wyckoff_label_and_equiv_coord_list(
conv_cell_site=site, wyckoff_dict=wyckoff_dict
)
if all(
# allow for sga conventional cell (and thus wyckoffs) being a multiple of BCS conventional cell
wyckoff_label not in wyckoff_symbol_list
for wyckoff_symbol_list in wyckoff_symbol_lists
) and all(
# allow for BCS conv cell (and thus wyckoffs) being a multiple of sga conv cell (allow it fam)
multiplied_wyckoff_symbol not in wyckoff_symbols
for multiplied_wyckoff_symbol in [
_multiply_wyckoffs([wyckoff_label], n=n)[0] for n in range(1, 5) # up to 4x
]
):
return False # break on first non-match
doped_wyckoffs.append(wyckoff_label)
return any(
# allow for sga conventional cell (and thus wyckoffs) being a multiple of BCS conventional cell
set(i) == set(doped_wyckoffs)
for i in wyckoff_symbol_lists
) or any(
set(i) == set(wyckoff_symbols)
for i in [
# allow for BCS conv cell (and thus wyckoffs) being a multiple of sga conv cell (allow it fam)
_multiply_wyckoffs(doped_wyckoffs, n=n)
for n in range(1, 5) # up to 4x
]
) # False if no complete match, True otherwise
def _read_wyckoff_datafile(spacegroup, f, setting=None):
"""
Read the ``wyckpos.dat`` file of specific spacegroup and returns a
dictionary with this information.
"""
if isinstance(spacegroup, int):
pass
elif isinstance(spacegroup, str):
spacegroup = " ".join(spacegroup.strip().split())
else:
raise ValueError("`spacegroup` must be of type int or str")
line = _skip_to_spacegroup(f, spacegroup, setting)
wyckoff_dict = {"letters": [], "multiplicity": [], "number_of_letters": 0}
line_list = line.split()
if line_list[0].isdigit():
wyckoff_dict["spacegroup"] = int(line_list[0])
else:
spacegroup, wyckoff_dict["setting"] = line_list[0].split("-")
wyckoff_dict["spacegroup"] = int(spacegroup)
if len(line.split()) > 1:
eq_sites = line.split("(")[1:]
wyckoff_dict["equivalent_sites"] = ([eq[:-1] for eq in eq_sites])[1:]
wyckoff_dict["equivalent_sites"][-1] = wyckoff_dict["equivalent_sites"][-1][:-1]
while True:
line = f.readline()
if line == "\n":
break
letter, multiplicity = line.split()[:2]
coordinates_raw = line.split()[-1].split("(")[1:]
site_symmetry = "".join(line.split()[2:-1])
wyckoff_dict["letters"].append(letter)
wyckoff_dict["number_of_letters"] += 1
wyckoff_dict["multiplicity"].append(int(multiplicity))
coordinates = [coord[:-1] for coord in coordinates_raw]
wyckoff_dict[letter] = {
"multiplicity": multiplicity,
"site_symmetry": site_symmetry,
"coordinates": coordinates,
}
return wyckoff_dict
def _get_wyckoff_datafile():
"""
Return default path to Wyckoff datafile.
"""
return os.path.join(os.path.dirname(__file__), "wyckpos.dat")
def _skip_to_spacegroup(f, spacegroup, setting=None):
"""
Read lines from f until a blank line is encountered.
"""
name = str(spacegroup) if setting is None else f"{spacegroup!s}-{setting}"
while True:
line = f.readline()
if not line:
raise ValueError(
f"Invalid spacegroup {spacegroup} with setting: {setting}. Not found in the Wyckoff "
f"database!"
)
if line.startswith(name):
break
return line
[docs]
def point_symmetry_from_defect(defect, symm_ops=None, symprec=0.01):
"""
Get the defect site point symmetry from a Defect object.
Note that this is intended only to be used for unrelaxed,
as-generated ``Defect`` objects (rather than parsed defects).
Args:
defect (Defect): ``Defect`` object.
symm_ops (list):
List of symmetry operations of ``defect.structure``, to avoid
re-calculating. Default is None (recalculates).
symprec (float):
Symmetry tolerance for ``spglib``. Default is 0.01.
Returns:
str: Defect point symmetry.
"""
symm_dataset, _unique_sites = _get_symm_dataset_of_struc_with_all_equiv_sites(
defect.site.frac_coords, defect.structure, symm_ops=symm_ops, symprec=symprec
)
spglib_point_group_symbol = schoenflies_from_hermann(symm_dataset.site_symmetry_symbols[-1])
if spglib_point_group_symbol is not None:
return spglib_point_group_symbol
# symm_ops approach failed, just use diagonal defect supercell approach:
defect_diagonal_supercell = defect.get_supercell_structure(
sc_mat=np.array([[2, 0, 0], [0, 2, 0], [0, 0, 2]]),
dummy_species="X",
) # create defect supercell, which is a diagonal expansion of the unit cell so that the defect
# periodic image retains the unit cell symmetry, in order not to affect the point group symmetry
sga = get_sga(defect_diagonal_supercell, symprec=symprec)
return schoenflies_from_hermann(sga.get_point_group_symbol())
[docs]
def point_symmetry_from_defect_entry(
defect_entry: DefectEntry,
symm_ops: Optional[list] = None,
symprec: Optional[float] = None,
relaxed: bool = True,
verbose: bool = True,
return_periodicity_breaking: bool = False,
):
r"""
Get the defect site point symmetry from a ``DefectEntry`` object.
Note: If ``relaxed = True`` (default), then this tries to use the
defect_entry.defect_supercell to determine the site symmetry. This will
thus give the `relaxed` defect point symmetry if this is a DefectEntry
created from parsed defect calculations. However, it should be noted
that this is not guaranteed to work in all cases; namely for non-diagonal
supercell expansions, or sometimes for non-scalar supercell expansion
matrices (e.g. a 2x1x2 expansion)(particularly with high-symmetry materials)
which can mess up the periodicity of the cell. doped tries to automatically
check if this is the case, and will warn you if so.
This can also be checked by using this function on your doped `generated` defects:
.. code-block:: python
from doped.generation import get_defect_name_from_entry
for defect_name, defect_entry in defect_gen.items():
print(defect_name, get_defect_name_from_entry(defect_entry, relaxed=False),
get_defect_name_from_entry(defect_entry), "\n")
And if the point symmetries match in each case, then using this function on your
parsed `relaxed` ``DefectEntry`` objects should correctly determine the final relaxed
defect symmetry - otherwise periodicity-breaking prevents this.
If periodicity-breaking prevents auto-symmetry determination, you can manually
determine the relaxed defect and bulk-site point symmetries, and/or orientational
degeneracy, from visualising the structures (e.g. using VESTA)(can use
``get_orientational_degeneracy`` to obtain the corresponding orientational
degeneracy factor for given defect/bulk-site point symmetries) and setting the
corresponding values in the
``calculation_metadata['relaxed point symmetry']/['bulk site symmetry']`` and/or
``degeneracy_factors['orientational degeneracy']`` attributes.
Note that the bulk-site point symmetry corresponds to that of ``DefectEntry.defect``,
or equivalently ``calculation_metadata["bulk_site"]/["unrelaxed_defect_structure"]``,
which for vacancies/substitutions is the symmetry of the corresponding bulk site,
while for interstitials it is the point symmetry of the `final relaxed` interstitial
site when placed in the (unrelaxed) bulk structure.
The degeneracy factor is used in the calculation of defect/carrier concentrations
and Fermi level behaviour (see e.g. https://doi.org/10.1039/D2FD00043A &
https://doi.org/10.1039/D3CS00432E).
Args:
defect_entry (DefectEntry): ``DefectEntry`` object.
symm_ops (list):
List of symmetry operations of either the defect_entry.bulk_supercell
structure (if ``relaxed=False``) or defect_entry.defect_supercell (if
``relaxed=True``), to avoid re-calculating. Default is None (recalculates).
symprec (float):
Symmetry tolerance for ``spglib``. Default is 0.01 for unrelaxed structures,
0.1 for relaxed (to account for residual structural noise, matching that
used by the ``Materials Project``). You may want to adjust for your system
(e.g. if there are very slight octahedral distortions etc.).
relaxed (bool):
If False, determines the site symmetry using the defect site `in the
unrelaxed bulk supercell` (i.e. the bulk site symmetry), otherwise
tries to determine the point symmetry of the relaxed defect in the
defect supercell. Default is True.
verbose (bool):
If True, prints a warning if the supercell is detected to break
the crystal periodicity (and hence not be able to return a reliable
`relaxed` point symmetry). Default is True.
return_periodicity_breaking (bool):
If True, also returns a boolean specifying if the supercell has been
detected to break the crystal periodicity (and hence not be able to
return a reliable `relaxed` point symmetry) or not. Mainly for
internal ``doped`` usage. Default is False.
Returns:
str: Defect point symmetry (and if ``return_periodicity_breaking = True``,
a boolean specifying if the supercell has been detected to break the crystal
periodicity).
"""
if symprec is None:
symprec = 0.1 if relaxed else 0.01 # relaxed structures likely have structural noise
# May need to adjust symprec (e.g. for Ag2Se, symprec of 0.2 is too large as we have very
# slight distortions present in the unrelaxed material).
# from spglib docs: For atomic positions, roughly speaking, two position vectors x and x' in
# Cartesian coordinates are considered to be the same if |x' - x| < symprec. The angle distortion
# between basis vectors is converted to a length and compared with this distance tolerance.
# we _could_ do bulk-bond length dependent symprec, which seems like it would be physically
# reasonable (basically being a proxy accounting for larger structural/positional noise for same force
# noise in DFT supercell calcs), but from testing this didn't seem to really improve accuracy in
# general (e.g. for Sb2O5 split-interstitial seemed like >0.1 best, while <0.12 required for SrTiO3
# despite smaller bond length)
if not relaxed and defect_entry.defect.defect_type != DefectType.Interstitial:
# then easy, can just be taken from symmetry dataset of defect structure
symm_dataset = get_sga(defect_entry.defect.structure, symprec=symprec).get_symmetry_dataset()
return schoenflies_from_hermann(
symm_dataset.site_symmetry_symbols[defect_entry.defect.defect_site_index]
)
supercell = _get_defect_supercell(defect_entry) if relaxed else _get_bulk_supercell(defect_entry)
if symm_ops is None:
symm_ops = get_sga(supercell).get_symmetry_operations()
# For relaxed = True, often only works for relaxed defect structures if it is a scalar matrix
# supercell expansion of the primitive/conventional cell (otherwise can mess up the periodicity).
# Sometimes works even if not a scalar matrix for _low-symmetry_ systems (counter-intuitively),
# because then the breakdown in periodicity affects the defect site symmetry less. This can be
# checked by seeing if the site symmetry of the defect site in the unrelaxed structure is correctly
# guessed (in which case the supercell site symmetry is not messing up the symmetry detection),
# as shown in the docstrings. Here we test using the 'unrelaxed_defect_structure' in the DefectEntry
# calculation_metadata if present, which, if it gives the same result as relaxed=False, means that
# for this defect at least, there is no periodicity-breaking which is affecting the symmetry
# determination.
matching = True
if relaxed:
if unrelaxed_defect_structure := _get_unrelaxed_defect_structure(defect_entry):
matching = _check_relaxed_defect_symmetry_determination(
defect_entry,
unrelaxed_defect_structure=unrelaxed_defect_structure,
symprec=symprec,
verbose=verbose,
)
else:
warnings.warn(
"`relaxed` was set to True (i.e. get _relaxed_ defect symmetry), but the "
"`calculation_metadata`/`bulk_entry.structure` attributes are not set for `DefectEntry`, "
"suggesting that this DefectEntry was not parsed from calculations using doped. This "
"means doped cannot automatically check if the supercell shape is breaking the cell "
"periodicity here or not (see docstring) - the point symmetry groups are not guaranteed "
"to be correct here!"
)
_failed = False
spglib_point_group_symbol = None
if relaxed:
with contextlib.suppress(Exception):
spglib_point_group_symbol = schoenflies_from_hermann(
get_sga(supercell, symprec=symprec).get_point_group_symbol()
)
if not relaxed or spglib_point_group_symbol is None:
defect_supercell_bulk_site_coords = _get_defect_supercell_bulk_site_coords(
defect_entry, relaxed=relaxed
)
if defect_supercell_bulk_site_coords is not None:
try:
symm_dataset, _unique_sites = _get_symm_dataset_of_struc_with_all_equiv_sites(
defect_supercell_bulk_site_coords,
supercell,
symm_ops=symm_ops, # defect symm_ops needed for relaxed=True, bulk for relaxed=False
symprec=symprec,
dist_tol=symprec,
)
# Note:
# This code works to get the site symmetry of a defect site in a periodicity-breaking
# supercell, but only when the defect has not yet been relaxed. Still has the issue that
# once we have relaxation around the defect site in a periodicity-breaking supercell,
# then the (local) point symmetry cannot be easily determined as the whole supercell
# symmetry is broken.
# Future work could try a local structure analysis to determine the local point symmetry to
# counteract this.
# unique_sites = _get_all_equiv_sites( # defect site but bulk supercell & symm_ops
# site.frac_coords, bulk_supercell, symm_ops
# )
# sga_with_all_X = _get_sga_with_all_X( # defect unique sites but bulk supercell
# bulk_supercell, unique_sites, symprec=symprec
# )
# symm_dataset = sga_with_all_X.get_symmetry_dataset()
except AttributeError:
_failed = True
if defect_supercell_bulk_site_coords is None or _failed:
point_group = point_symmetry_from_defect(
defect_entry.defect, symm_ops=symm_ops, symprec=symprec
)
# possibly pymatgen DefectEntry object without defect_supercell_site set
if relaxed:
warnings.warn(
"Symmetry determination failed with the standard approach (likely due to this being a "
"DefectEntry which has not been generated/parsed with doped?). Thus the _relaxed_ "
"point group symmetry cannot be reliably automatically determined."
)
return (point_group, not matching) if return_periodicity_breaking else point_group
return point_group
if not relaxed:
# `site_symmetry_symbols[-1]` should be used (within this equiv sites approach) for unrelaxed
# defects (rather than `pointgroup`), as the site symmetry can be lower than the crystal point
# group, but not vice versa; so when populating all equivalent sites (of the defect site,
# in the bulk supercell) the overall point group is retained and is not necessarily the defect
# site symmetry. e.g. consider populating all equivalent sites of a C1 interstitial site in a
# structure (such as CdTe), then the overall point group is still the bulk point group,
# but the site symmetry is in fact C1.
# This issue is avoided for relaxed defect supercells as we take the symm_ops of our reduced
# symmetry cell rather than that of the bulk (so no chance of spurious symmetry upgrade from
# equivalent sites), and hence the max point symmetry is the point symmetry of the defect
spglib_point_group_symbol = schoenflies_from_hermann(symm_dataset.site_symmetry_symbols[-1])
# Note that, if the supercell is non-periodicity-breaking, then the site symmetry can be simply
# determined using the point group of the unrelaxed defect structure:
# unrelaxed_defect_supercell = defect_entry.calculation_metadata.get(
# "unrelaxed_defect_structure", defect_supercell
# )
# return schoenflies_from_hermann(
# get_sga(unrelaxed_defect_supercell, symprec).get_symmetry_dataset().pointgroup,
# )
# But current approach works for all cases with unrelaxed defect structures
else:
# For relaxed defects the "defect supercell site" is not necessarily the true centre of mass of
# the defect (e.g. for split-interstitials, split-vacancies, swapped vacancies etc),
# so use 'pointgroup' output (in this case the reduced symmetry avoids the symmetry-upgrade
# possibility with the equivalent sites, as when relaxed=False)
spglib_point_group_symbol = schoenflies_from_hermann(symm_dataset.pointgroup)
# This also works (at least for non-periodicity-breaking supercells) for relaxed defects in
# most cases, but is slightly less robust (more sensitive to ``symprec`` choice) than the
# approach above:
# schoenflies_from_hermann(
# get_sga(
# defect_supercell, symprec=symprec
# ).get_symmetry_dataset().pointgroup
# )
if spglib_point_group_symbol is not None:
return (
(spglib_point_group_symbol, not matching)
if return_periodicity_breaking
else (spglib_point_group_symbol)
)
# symm_ops approach failed, just use diagonal defect supercell approach:
if relaxed:
raise RuntimeError(
"Site symmetry could not be determined using the defect supercell, and so the relaxed site "
"symmetry cannot be automatically determined (set relaxed=False to obtain the (unrelaxed) "
"bulk site symmetry)."
)
defect_diagonal_supercell = defect_entry.defect.get_supercell_structure(
sc_mat=np.array([[2, 0, 0], [0, 2, 0], [0, 0, 2]]),
dummy_species="X",
) # create defect supercell, which is a diagonal expansion of the unit cell so that the defect
# periodic image retains the unit cell symmetry, in order not to affect the point group symmetry
sga = get_sga(defect_diagonal_supercell, symprec=symprec)
point_group = schoenflies_from_hermann(sga.get_point_group_symbol())
return (point_group, not matching) if return_periodicity_breaking else point_group
def _check_relaxed_defect_symmetry_determination(
defect_entry: DefectEntry,
unrelaxed_defect_structure: Structure = None,
symprec: float = 0.1,
verbose: bool = False,
):
defect_supercell_bulk_site_coords = _get_defect_supercell_bulk_site_coords(defect_entry, relaxed=False)
if defect_supercell_bulk_site_coords is None:
raise AttributeError(
"`defect_entry.defect_supercell_site` or `defect_entry.sc_defect_frac_coords` are not "
"defined! Needed to check defect supercell periodicity (for symmetry determination)"
)
if unrelaxed_defect_structure is None:
unrelaxed_defect_structure = _get_unrelaxed_defect_structure(defect_entry)
if unrelaxed_defect_structure is not None:
unrelaxed_spglib_point_group_symbol = schoenflies_from_hermann(
get_sga(unrelaxed_defect_structure, symprec=symprec).get_symmetry_dataset().pointgroup,
)
bulk_supercell = _get_bulk_supercell(defect_entry)
bulk_symm_ops = get_sga(bulk_supercell).get_symmetry_operations()
symm_dataset, _unique_sites = _get_symm_dataset_of_struc_with_all_equiv_sites(
defect_supercell_bulk_site_coords,
bulk_supercell,
symm_ops=bulk_symm_ops,
symprec=symprec,
dist_tol=symprec,
)
bulk_spglib_point_group_symbol = schoenflies_from_hermann(symm_dataset.site_symmetry_symbols[-1])
if bulk_spglib_point_group_symbol != unrelaxed_spglib_point_group_symbol:
if verbose:
warnings.warn(
"`relaxed` is set to True (i.e. get _relaxed_ defect symmetry), but doped has "
"detected that the defect supercell is likely a non-scalar matrix expansion which "
"could be breaking the cell periodicity and possibly preventing the correct _relaxed_ "
"point group symmetry from being automatically determined. You can set relaxed=False "
"to instead get the (unrelaxed) bulk site symmetry, and/or manually "
"check/set/edit the point symmetries and corresponding orientational degeneracy "
"factors by inspecting/editing the "
"calculation_metadata['relaxed point symmetry']/['bulk site symmetry'] and "
"degeneracy_factors['orientational degeneracy'] attributes."
)
return False
return True
return False # return False if symmetry couldn't be checked
[docs]
def point_symmetry_from_structure(
structure: Structure,
bulk_structure: Optional[Structure] = None,
symm_ops: Optional[list] = None,
symprec: Optional[float] = None,
relaxed: bool = True,
verbose: bool = True,
return_periodicity_breaking: bool = False,
):
r"""
Get the point symmetry of a given structure.
Note: For certain non-trivial supercell expansions, the broken cell
periodicity can break the site symmetry and lead to incorrect point
symmetry determination (particularly if using non-scalar supercell
matrices with high symmetry materials). If the unrelaxed bulk structure
(``bulk_structure``) is also supplied, then ``doped`` will determine
the defect site and then automatically check if this is the case,
and warn you if so.
This can also be checked by using this function on your doped `generated` defects:
.. code-block:: python
from doped.generation import get_defect_name_from_entry
for defect_name, defect_entry in defect_gen.items():
print(defect_name, get_defect_name_from_entry(defect_entry, relaxed=False),
get_defect_name_from_entry(defect_entry), "\n")
And if the point symmetries match in each case, then using this function
on your parsed `relaxed` ``DefectEntry`` objects should correctly determine
the final relaxed defect symmetry - otherwise periodicity-breaking prevents this.
If ``bulk_structure`` is supplied and ``relaxed`` is set to ``False``,
then returns the bulk site symmetry of the defect, which for
vacancies/substitutions is the symmetry of the corresponding bulk site,
while for interstitials it is the point symmetry of the `final relaxed`
interstitial site when placed in the (unrelaxed) bulk structure.
Args:
structure (Structure):
``Structure`` object for which to determine the point symmetry.
bulk_structure (Structure):
``Structure`` object of the bulk structure, if known. Default is None.
If provided and ``relaxed = True``, will be used to check if the supercell
is breaking the crystal periodicity (and thus preventing accurate
determination of the relaxed defect point symmetry) and warn you if so.
symm_ops (list):
List of symmetry operations of either the ``bulk_structure`` structure
(if ``relaxed=False``) or ``structure`` (if ``relaxed=True``), to avoid
re-calculating. Default is None (recalculates).
symprec (float):
Symmetry tolerance for ``spglib``. Default is 0.01 for unrelaxed structures,
0.1 for relaxed (to account for residual structural noise, matching that
used by the ``Materials Project``). You may want to adjust for your system
(e.g. if there are very slight octahedral distortions etc.).
relaxed (bool):
If ``False``, determines the site symmetry using the defect site `in the
unrelaxed bulk supercell` (i.e. the bulk site symmetry), otherwise
tries to determine the point symmetry of the relaxed defect in the
defect supercell. Default is True.
verbose (bool):
If True, prints a warning if the supercell is detected to break
the crystal periodicity (and hence not be able to return a reliable
`relaxed` point symmetry). Default is True.
return_periodicity_breaking (bool):
If True, also returns a boolean specifying if the supercell has been
detected to break the crystal periodicity (and hence not be able to
return a reliable `relaxed` point symmetry) or not. Default is False.
Returns:
str: Structure point symmetry (and if ``return_periodicity_breaking = True``,
a boolean specifying if the supercell has been detected to break the crystal
periodicity).
"""
if symprec is None:
symprec = 0.1 if relaxed else 0.01 # relaxed structures likely have structural noise
spglib_point_group_symbol = None
if relaxed and bulk_structure is None:
with contextlib.suppress(Exception):
spglib_point_group_symbol = schoenflies_from_hermann(
get_sga(structure, symprec=symprec).get_point_group_symbol()
)
if spglib_point_group_symbol is not None:
return (
(spglib_point_group_symbol, False)
if return_periodicity_breaking
else spglib_point_group_symbol
)
if bulk_structure is not None:
from doped.analysis import defect_from_structures
(
defect,
defect_site,
defect_site_in_bulk, # bulk site for vac/sub, relaxed defect site w/interstitials
defect_site_index, # in this initial_defect_structure
bulk_site_index,
guessed_initial_defect_structure,
unrelaxed_defect_structure,
bulk_voronoi_node_dict,
) = defect_from_structures(
bulk_structure,
structure,
oxi_state="Undetermined",
return_all_info=True,
)
defect_entry = DefectEntry(
# pmg attributes:
defect=defect, # this corresponds to _unrelaxed_ defect
charge_state=0,
sc_entry=ComputedStructureEntry(
structure=structure,
energy=0.0, # needs to be set, so set to 0.0
),
sc_defect_frac_coords=defect_site.frac_coords, # _relaxed_ defect site
bulk_entry=None,
# doped attributes:
defect_supercell_site=defect_site, # _relaxed_ defect site
defect_supercell=structure,
bulk_supercell=bulk_structure,
)
return point_symmetry_from_defect_entry(
defect_entry,
symm_ops=symm_ops,
symprec=symprec,
relaxed=relaxed,
verbose=verbose,
return_periodicity_breaking=return_periodicity_breaking,
)
# else bulk structure is None and normal relaxed structure symmetry determination failed
raise RuntimeError(
"Target site symmetry could not be determined using just the input structure. Please also supply "
"the unrelaxed bulk structure (`bulk_structure`)."
)
[docs]
def point_symmetry_from_site(
site: Union[PeriodicSite, np.ndarray, list],
structure: Structure,
coords_are_cartesian: bool = False,
symm_ops: Optional[list] = None,
symprec: float = 0.01,
):
r"""
Get the point symmetry of a site in a structure.
Args:
site (Union[PeriodicSite, np.ndarray, list]):
Site for which to determine the point symmetry. Can be a
``PeriodicSite`` object, or a list or numpy array of the
coordinates of the site (fractional coordinates by default,
or cartesian if ``coords_are_cartesian = True``).
structure (Structure):
``Structure`` object for which to determine the point symmetry
of the site.
coords_are_cartesian (bool):
If True, the site coordinates are assumed to be in cartesian
coordinates. Default is False.
symm_ops (list):
List of symmetry operations of the ``structure`` to avoid
re-calculating. Default is None (recalculates).
symprec (float):
Symmetry tolerance for ``spglib``. Default is 0.01. You may want
to adjust for your system (e.g. if there are very slight
octahedral distortions etc.).
Returns:
str: Site point symmetry.
"""
if isinstance(site, (np.ndarray, list)):
site = PeriodicSite(
species="X", coords=site, lattice=structure.lattice, coords_are_cartesian=coords_are_cartesian
)
symm_dataset, _unique_sites = _get_symm_dataset_of_struc_with_all_equiv_sites(
site.frac_coords, structure, symm_ops=symm_ops, symprec=symprec
)
return schoenflies_from_hermann(symm_dataset.site_symmetry_symbols[-1])
# Schoenflies, Hermann-Mauguin, spgid dict: (Taken from the excellent Abipy with GNU GPL License)
_PTG_IDS = [
("C1", "1", 1),
("Ci", "-1", 2),
("C2", "2", 3),
("Cs", "m", 6),
("C2h", "2/m", 10),
("D2", "222", 16),
("C2v", "mm2", 25),
("D2h", "mmm", 47),
("C4", "4", 75),
("S4", "-4", 81),
("C4h", "4/m", 83),
("D4", "422", 89),
("C4v", "4mm", 99),
("D2d", "-42m", 111),
("D4h", "4/mmm", 123),
("C3", "3", 143),
("C3i", "-3", 147),
("D3", "32", 149),
("C3v", "3m", 156),
("D3d", "-3m", 162),
("C6", "6", 168),
("C3h", "-6", 174),
("C6h", "6/m", 175),
("D6", "622", 177),
("C6v", "6mm", 183),
("D3h", "-6m2", 189),
("D6h", "6/mmm", 191),
("T", "23", 195),
("Th", "m-3", 200),
("O", "432", 207),
("Td", "-43m", 215),
("Oh", "m-3m", 221),
]
_SCH_to_HERM = {t[0]: t[1] for t in _PTG_IDS}
_HERM_to_SCH = {t[1]: t[0] for t in _PTG_IDS}
_SPGID_to_SCH = {t[2]: t[0] for t in _PTG_IDS}
_SCH_to_SPGID = {t[0]: t[2] for t in _PTG_IDS}
sch_symbols = list(_SCH_to_HERM.keys())
[docs]
def schoenflies_from_hermann(herm_symbol):
"""
Convert from Hermann-Mauguin to Schoenflies.
"""
herm_symbol = herm_symbol.replace(".", "")
schoenflies = _HERM_to_SCH.get(herm_symbol)
if schoenflies is None:
# try rearranging, symbols in spglib can be rearranged vs _HERM_to_SCH dict
# get _HERM_to_SCH key that has the same characters as herm_symbol
# (i.e. same characters, but possibly in a different order)
from collections import Counter
def find_matching_key(input_str, input_dict):
input_str_counter = Counter(input_str)
for key in input_dict:
if Counter(key) == input_str_counter:
return key
return None
herm_key = find_matching_key(herm_symbol, _HERM_to_SCH)
if herm_key is not None:
schoenflies = _HERM_to_SCH[herm_key]
return schoenflies
_point_group_order = {
"C1": 1,
"Ci": 2, # aka. S2, -1 in Hermann-Mauguin
"C2": 2,
"Cs": 2, # aka. C1h (m in Hermann-Mauguin)
"C3": 3,
"C4": 4,
"S4": 4, # C4 with improper rotation
"C2h": 4, # 2/m in Hermann-Mauguin
"D2": 4, # 222 in Hermann-Mauguin
"C2v": 4, # mm2 in Hermann-Mauguin
"C3i": 6, # aka. S6, -3 in Hermann-Mauguin
"C6": 6,
"C3h": 6,
"D3": 6, # 32 in Hermann-Mauguin
"C3v": 6, # 3m in Hermann-Mauguin
"D2h": 8, # mmm in Hermann-Mauguin
"C4h": 8, # 4/m in Hermann-Mauguin
"D4": 8, # 422 in Hermann-Mauguin
"C4v": 8, # 4mm in Hermann-Mauguin
"D2d": 8, # 42m in Hermann-Mauguin
"C6h": 12, # 6/m in Hermann-Mauguin
"T": 12, # 23 in Hermann-Mauguin
"D3d": 12, # 3m1 in Hermann-Mauguin
"D6": 12, # 622 in Hermann-Mauguin
"C6v": 12, # 6mm in Hermann-Mauguin
"D3h": 12, # 6m2 in Hermann-Mauguin
"D4h": 16, # 4/mmm in Hermann-Mauguin
"D6h": 24, # 6/mmm in Hermann-Mauguin
"Th": 24, # m3 in Hermann-Mauguin
"O": 24, # 432 in Hermann-Mauguin
"Td": 24, # 43m in Hermann-Mauguin
"Oh": 48, # m3m in Hermann-Mauguin
}
[docs]
def group_order_from_schoenflies(sch_symbol):
"""
Return the order of the point group from the Schoenflies symbol.
Useful for symmetry and orientational degeneracy analysis.
"""
return _point_group_order[sch_symbol]