amd.periodicset.PeriodicSet._random() - Code Metrics - Inspection of "Update .scrutinizer.yml" - dwiddo/average-minimum-distance - Measure and Improve Code Quality continuously with Scrutinizer

Test Setup Failed

Push — master ( ad57f0...44572e )

by Daniel

created 2025-06-17 20:20 UTC

amd.periodicset.PeriodicSet._random() A

↳ Parent: amd.periodicset

Complexity

Conditions

Size

Total Lines	24
Code Lines	13

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
eloc	13
dl	0
loc	24
rs	9.75
c	0
b	0
f	0
cc	2
nop	4

"""Implements the :class:`PeriodicSet` class representing a periodic
set defined by a motif (collection of points) and unit cell (basis)
modelling a crystal with a point at the center of each atom.

:class:`PeriodicSet` s are yielded by the children of
:class:`amd.io._Reader <.io._Reader>` and are passed to
:func:`PDD() <.calculate.PDD>`, :func:`AMD() <.calculate.AMD>` etc to
calculate the invariants of crystals.
"""

from __future__ import annotations
from typing import Optional, List, NamedTuple
import math

import numpy as np
from scipy.special import factorial

from .utils import (
    cellpar_to_cell,
    cellpar_to_cell_2D,
    random_cell,
    cell_to_cellpar,
)
from ._types import FloatArray, UIntArray, AtomTypeArray
from .globals_ import ATOMIC_MASSES, ATOMIC_NUMBERS_TO_SYMS


__all__ = ["PeriodicSet", "DisorderGroup", "DisorderAssembly"]


class DisorderGroup(NamedTuple):
    indices: list[int]
    occupancy: float


class DisorderAssembly(NamedTuple):
    groups: list[DisorderGroup]
    is_substitutional: bool


class PeriodicSet:
    """A periodic set is a collection of points (motif) which
    periodically repeats according to a lattice (unit cell), which can
    model a crystal with a point in the centre of each atom.

    :class:`PeriodicSet` s are yielded by the children of
    :class:`amd.io._Reader <.io._Reader>` and are passed to
    :func:`PDD() <.calculate.PDD>`, :func:`AMD() <.calculate.AMD>` etc to
    calculate the invariants of crystals.

    Parameters
    ----------
    motif : :class:`numpy.ndarray`
        Collection of motif points in Cartesian (orthogonal)
        coordinates, shape ``(m, dims)``.
    cell : :class:`numpy.ndarray`
        Square array representing the unit cell in Cartesian
        (orthogonal) coordinates, shape ``(dims, dims)``. Use
        :func:`amd.cellpar_to_cell <.utils.cellpar_to_cell>` to convert
        a vector of 6 conventional cell parameters a,b,c,α,β,γ to a
        square matrix.
    name : str, optional
        Name of the periodic set.
    asym_unit : :class:`numpy.ndarray`, optional
        Indices of motif points constituting an asymmetric unit, so that
        ``motif[asym_unit]`` is an asymmetric unit. If given, the motif
        should be such that ``motif[asym_unit[i]:asym_unit[i+1]]`` are
        points corresponding to ``motif[asym_unit[i]``.
    multiplicities : :class:`numpy.ndarray`, optional
        Wyckoff multiplicities of points in the asymmetric unit (the
        number of unique sites generated under symmetries).
    types : list, optional
        Atomic numbers of atoms in the asymmetric unit.
    labels : list, optional
        Labels of atoms in the asymmetric unit.
    disorder : list, optional
        Contains disorder assemblies of the crystal.
    """

    def __init__(
        self,
        motif: FloatArray,
        cell: FloatArray,
        name: Optional[str] = None,
        asym_unit: Optional[UIntArray] = None,
        multiplicities: Optional[UIntArray] = None,
        types: Optional[AtomTypeArray] = None,
        labels: Optional[List[str]] = None,
        disorder: Optional[list[DisorderAssembly]] = None,
    ):
        self.motif = motif
        self.cell = cell
        self.name = name
        self._asym_unit = asym_unit
        self._multiplicities = multiplicities
        self.types = types
        self.labels = labels
        self.disorder = disorder

    @property
    def ndim(self) -> int:
        """Dimension of the ambient space."""
        return self.cell.shape[0]

    @property
    def m(self) -> int:
        """Number of motif points."""
        return self.motif.shape[0]

    @property
    def n_asym(self) -> int:
        """Number of points in the asymmetric unit."""
        if self.asym_unit is None:
            return self.m
        else:
            return self.asym_unit.shape[0]

    @property
    def asym_unit(self) -> np.ndarray:
        if self._asym_unit is None:
            return np.arange(self.m)
        else:
            return self._asym_unit

    @property
    def multiplicities(self) -> np.ndarray:
        if self._multiplicities is None:
            return np.ones((self.m,), dtype=np.uint64)
        else:
            return self._multiplicities

    @property
    def occupancies(self) -> np.ndarray:
        occs = np.ones((self.n_asym, ), dtype=np.float64)
        if self.disorder:
            for asm in self.disorder:
                for grp in asm.groups:
                    occs[grp.indices] = grp.occupancy
        return occs

    def __str__(self) -> str:
        m, n = self.motif.shape
        m_pl = "" if m == 1 else "s"
        n_pl = "" if n == 1 else "s"
        name_str = f"{self.name}: " if self.name is not None else ""
        return f"PeriodicSet({name_str}{m} point{m_pl} in {n} dimension{n_pl})"

    def __repr__(self) -> str:
        motif_str = str(self.motif).replace("\n ", "\n" + " " * 11)
        cell_str = str(self.cell).replace("\n ", "\n" + " " * 10)

        optional_attrs = []
        for attr in ("_asym_unit", "_multiplicities"):
            val = getattr(self, attr)
            if val is not None:
                st = str(val).replace("\n ", "\n" + " " * (len(attr) + 6))
                optional_attrs.append(f"{attr[1:]}={st}")

        for attr in ("types", "labels"):
            val = getattr(self, attr)
            if val is not None:
                st = str(val).replace("\n ", "\n" + " " * (len(attr) + 6))
                optional_attrs.append(f"{attr}={st}")

        optional_attrs_str = ",\n    " if optional_attrs else ""
        optional_attrs_str += ",\n    ".join(optional_attrs)

        return (
            f"PeriodicSet(name={self.name},\n"
            f"    motif={motif_str},\n"
            f"    cell={cell_str}{optional_attrs_str})"
        )

    def PPC(self):
        r"""Return the point packing coefficient (PPC) of the
        PeriodicSet.

        The PPC is a constant of any periodic set determining the
        asymptotic behaviour of its isometry invariants. As
        :math:`k \rightarrow \infty`, the ratio
        :math:`\text{AMD}_k / \sqrt[n]{k}` converges to the PPC, as does
        any row of its PDD.

        For a unit cell :math:`U` and :math:`m` motif points in
        :math:`n` dimensions,

        .. math::

            \text{PPC} = \sqrt[n]{\frac{\text{Vol}[U]}{m V_n}}

        where :math:`V_n` is the volume of a unit sphere in :math:`n`
        dimensions.

        Returns
        -------
        ppc : float
            The PPC of ``self``.
        """
        m, n = self.motif.shape
        t = int(n // 2)
        if n % 2 == 0:
            ball_vol = (np.pi ** t) / factorial(t)
        else:
            ball_vol = (2 * factorial(t) * (4 * np.pi) ** t) / factorial(n)
        return (np.abs(np.linalg.det(self.cell)) / (m * ball_vol)) ** (1.0 / n)

    def has_positional_disorder(self):
        if self.disorder is not None:
            return any(not asm.is_substitutional for asm in self.disorder)
        return False

    def has_substitutional_disorder(self):
        if self.disorder is not None:
            return any(asm.is_substitutional for asm in self.disorder)
        return False

    def density(self):
        """Return the physical density of the PeriodicSet representing a
        crystal in g/cm3.

        Returns
        -------
        density : float
            The physical density of ``self`` in g/cm3.
        """
        if self.types is None:
            raise ValueError(f"PeriodicSet(name={self.name}) has no types")
        cell_mass = sum(
            ATOMIC_MASSES[num] * tot_occ
            for num, tot_occ in zip(self.types, self.multiplicities * self.occupancies)
        )
        cell_vol = np.linalg.det(self.cell)
        return (cell_mass / cell_vol) * 1.66053907

    def formula(self):
        """Return the reduced formula of the PeriodicSet representing a
        crystal.

        Returns
        -------
        formula : str
            The reduced formula of ``self``.
        """
        counter = self._formula_dict()
        counter = {ATOMIC_NUMBERS_TO_SYMS[n]: count for n, count in counter.items()}
        return "".join(
            f"{lab}{num:.3g}" if num != 1 else lab for lab, num in counter.items()
        )

    def _formula_dict(self):
        """Return a dictionary of atomic types in the form
        {atomic_num: count}.
        """

        if self.types is None:
            raise ValueError(f"PeriodicSet(name={self.name}) has no types")

        counter = {}
        for num, mul, occ in zip(self.types, self.multiplicities, self.occupancies):
            if num not in counter:
                counter[num] = 0
            counter[num] += occ * mul

        if 0 in counter:
            del counter[0]

        gcd = math.gcd(*(int(v) for v in counter.values() if np.mod(v, 1.0) == 0.0))

        if gcd == 0:
            gcd = 1

        counter_ = {}
        if 6 in counter:
            counter_[6] = counter[6] / gcd
            del counter[6]
        if 1 in counter:
            counter_[1] = counter[1] / gcd
            del counter[1]

        for n, count in sorted(counter.items(), key=lambda x: x[0]):
            counter_[int(n)] = count / gcd

        return counter_

    def _equal_cell_and_motif(self, other: PeriodicSet, tol: float = 1e-8) -> bool:
        """Used for debugging. True if both 1) unit cells are (close to)
        identical, and 2) the motifs are the same shape, and
        every point in one motif has a (close to) identical point
        somewhere in the other, accounting for pbc.
        """

        if (
            self.cell.shape != other.cell.shape
            or self.motif.shape != other.motif.shape
            or not np.allclose(self.cell, other.cell)
        ):
            return False

        cell_inv = np.linalg.inv(self.cell)
        fm1 = np.mod(self.motif @ cell_inv, 1)
        fm2 = np.mod(other.motif @ cell_inv, 1)
        d1 = np.abs(fm2[:, None] - fm1)
        d2 = np.abs(d1 - 1)
        diffs = np.amax(np.minimum(d1, d2), axis=-1)

        if not np.all(
            (np.amin(diffs, axis=0) <= tol) | (np.amin(diffs, axis=-1) <= tol)
        ):
            return False
        return True

    @staticmethod
    def cubic(scale: float = 1.0, dims: int = 3) -> PeriodicSet:
        """Return a :class:`PeriodicSet` representing a cubic lattice."""
        return PeriodicSet(np.zeros((1, dims)), np.identity(dims) * scale)

    @staticmethod
    def hexagonal(scale: float = 1.0, dims: int = 3) -> PeriodicSet:
        """Return a :class:`PeriodicSet` representing a hexagonal
        lattice. Implemented for dimensions 2 and 3 only.
        """

        if dims == 3:
            cellpar = np.array([scale, scale, scale, 90.0, 90.0, 120.0])
            cell = cellpar_to_cell(cellpar)
        elif dims == 2:
            cell = cellpar_to_cell_2D(np.array([scale, scale, 60.0]))
        else:
            raise NotImplementedError(
                "amd.PeriodicSet.hexagonal() only implemented for dimensions "
                f"2 and 3, passed {dims}"
            )
        return PeriodicSet(np.zeros((1, dims)), cell)

    @staticmethod
    def random(
        n_points: int,
        length_bounds: tuple = (1.0, 2.0),
        angle_bounds: tuple = (60.0, 120.0),
        dims: int = 3,
    ) -> PeriodicSet:
        """Return a :class:`PeriodicSet` with a chosen number of
        randomly placed points in a random cell with edges between
        ``length_bounds`` and angles between ``angle_bounds``.
        Implemented for dimensions 2 and 3 only.
        """

        if dims not in (2, 3):
            raise NotImplementedError(
                "amd.PeriodicSet._random only implemented for dimensions 2 "
                f"and 3, passed {dims}"
            )

        cell = random_cell(
            length_bounds=length_bounds, angle_bounds=angle_bounds, dims=dims
        )
        frac_motif = np.random.uniform(size=(n_points, dims))
        return PeriodicSet(frac_motif @ cell, cell)

    def to_cif_str(self) -> str:
        """Return a CIF string representing the PeriodicSet. Does not
        write disorder data, atom labels or symmetries."""

        a, b, c, al, be, ga = cell_to_cellpar(self.cell)
        out = [
            f"data_{self.name}",
            "_symmetry_space_group_name_H-M   'P 1'",
            f"_cell_length_a   {a}",
            f"_cell_length_b   {b}",
            f"_cell_length_c   {c}",
            f"_cell_angle_alpha   {al}",
            f"_cell_angle_beta   {be}",
            f"_cell_angle_gamma   {ga}",
            "loop_",
            " _symmetry_equiv_pos_site_id",
            " _symmetry_equiv_pos_as_xyz",
            " 1  'x, y, z'",
            "loop_",
            " _atom_site_label",
            " _atom_site_type_symbol",
            " _atom_site_fract_x",
            " _atom_site_fract_y",
            " _atom_site_fract_z",
        ]

        frac_motif = self.motif @ np.linalg.inv(self.cell)
        ind = 0
        label_nums = {}
        for asym_ind, mul in enumerate(self.multiplicities):
            for _ in range(mul):

                sym = ATOMIC_NUMBERS_TO_SYMS[self.types[asym_ind]]
                if sym in label_nums:
                    label_nums[sym] += 1
                else:
                    label_nums[sym] = 1

                out.append(
                    "  ".join(
                        [
                            f"{sym}{label_nums[sym]}",
                            sym,
                            f"{frac_motif[ind, 0]:.4f}",
                            f"{frac_motif[ind, 1]:.4f}",
                            f"{frac_motif[ind, 2]:.4f}",
                        ]
                    )
                )
                ind += 1

        return "\n".join(out)

    def to_pymatgen_structure(self):

        if self.types is None:
            raise ValueError(
                "Cannot convert PeriodicSet to pymatgen Structure without "
                "atomic types."
            )

        from pymatgen.core import Structure

        species = [
            ATOMIC_NUMBERS_TO_SYMS[t]
            for m, t in zip(self.multiplicities, self.types)
            for _ in range(m)
        ]

        labels = None
        if self.labels is not None:
            labels = [
                f"{t}_{i+1}"
                for m, t in zip(self.multiplicities, self.labels)
                for i in range(m)
            ]

        return Structure(
            self.cell, species, self.motif, coords_are_cartesian=True, labels=labels
        )

    def majority_occupancy_config(self) -> PeriodicSet:

        if not self.disorder:
            return self

        asym_mask = np.full((self.n_asym, ), fill_value=True)
        motif_mask = np.full((self.m, ), fill_value=True)
        
        for asm in self.disorder:

            # find majority occupancy group
            majority_occ, majority_ind = 0, 0
            for i, grp in enumerate(asm.groups):
                if grp.occupancy > majority_occ:
                    majority_ind = i
                    majority_occ = grp.occupancy

            # mask all other groups
            for i, grp in enumerate(asm.groups):

                if i == majority_ind:
                    continue
                
                for j in grp.indices:
                    asym_mask[j] = False
                    m_i = self.asym_unit[j]
                    mul = self.multiplicities[j]
                    motif_mask[m_i : m_i + mul] = False

        muls = self.multiplicities[asym_mask]
        asym_unit = np.zeros_like(muls, dtype=np.int64)
        asym_unit[1:] = np.cumsum(muls)[:-1]
        masked_labels = [l for l, flag in zip(self.labels, asym_mask) if flag]

        return PeriodicSet(
            self.motif[motif_mask],
            self.cell,
            self.name,
            asym_unit,
            self.multiplicities[asym_mask],
            self.types[asym_mask],
            masked_labels
        )
    
    def compressed_dtypes(self) -> PeriodicSet:
        """Return a PeriodicSet using float32, uint32 and uint8 for
        atomic types."""
        
        if self._asym_unit is not None:
            asym_unit = np.array(self._asym_unit, dtype=np.uint32)
        else:
            asym_unit = None
        
        if self._multiplicities is not None:
            multiplicities = np.array(self._multiplicities, dtype=np.uint32)
        else:
            multiplicities = None
        
        if self.types is not None:
            types = np.array(self.types, dtype=np.uint8)
        else:
            types = None

        return PeriodicSet(
            np.array(self.motif, dtype=np.float32),
            np.array(self.cell, dtype=np.float32),
            name=self.name,
            asym_unit=asym_unit,
            multiplicities=multiplicities,
            types=types,
            labels=self.labels,
            disorder=self.disorder
        )


1			"""Implements the :class:`PeriodicSet` class representing a periodic
2			set defined by a motif (collection of points) and unit cell (basis)
3			modelling a crystal with a point at the center of each atom.
4
5			:class:`PeriodicSet` s are yielded by the children of
6			:class:`amd.io._Reader <.io._Reader>` and are passed to
7			:func:`PDD() <.calculate.PDD>`, :func:`AMD() <.calculate.AMD>` etc to
8			calculate the invariants of crystals.
9			"""
10
11			from __future__ import annotations
12			from typing import Optional, List, NamedTuple
13			import math
14
15			import numpy as np
16			from scipy.special import factorial
17
18			from .utils import (
19			cellpar_to_cell,
20			cellpar_to_cell_2D,
21			random_cell,
22			cell_to_cellpar,
23			)
24			from ._types import FloatArray, UIntArray, AtomTypeArray
25			from .globals_ import ATOMIC_MASSES, ATOMIC_NUMBERS_TO_SYMS
26
27
28			__all__ = ["PeriodicSet", "DisorderGroup", "DisorderAssembly"]
29
30
31			class DisorderGroup(NamedTuple):
32			indices: list[int]
33			occupancy: float
34
35
36			class DisorderAssembly(NamedTuple):
37			groups: list[DisorderGroup]
38			is_substitutional: bool
39
40
41			class PeriodicSet:
42			"""A periodic set is a collection of points (motif) which
43			periodically repeats according to a lattice (unit cell), which can
44			model a crystal with a point in the centre of each atom.
45
46			:class:`PeriodicSet` s are yielded by the children of
47			:class:`amd.io._Reader <.io._Reader>` and are passed to
48			:func:`PDD() <.calculate.PDD>`, :func:`AMD() <.calculate.AMD>` etc to
49			calculate the invariants of crystals.
50
51			Parameters
52			----------
53			motif : :class:`numpy.ndarray`
54			Collection of motif points in Cartesian (orthogonal)
55			coordinates, shape ``(m, dims)``.
56			cell : :class:`numpy.ndarray`
57			Square array representing the unit cell in Cartesian
58			(orthogonal) coordinates, shape ``(dims, dims)``. Use
59			:func:`amd.cellpar_to_cell <.utils.cellpar_to_cell>` to convert
60			a vector of 6 conventional cell parameters a,b,c,α,β,γ to a
61			square matrix.
62			name : str, optional
63			Name of the periodic set.
64			asym_unit : :class:`numpy.ndarray`, optional
65			Indices of motif points constituting an asymmetric unit, so that
66			``motif[asym_unit]`` is an asymmetric unit. If given, the motif
67			should be such that ``motif[asym_unit[i]:asym_unit[i+1]]`` are
68			points corresponding to ``motif[asym_unit[i]``.
69			multiplicities : :class:`numpy.ndarray`, optional
70			Wyckoff multiplicities of points in the asymmetric unit (the
71			number of unique sites generated under symmetries).
72			types : list, optional
73			Atomic numbers of atoms in the asymmetric unit.
74			labels : list, optional
75			Labels of atoms in the asymmetric unit.
76			disorder : list, optional
77			Contains disorder assemblies of the crystal.
78			"""
79
80			def __init__(
81			self,
82			motif: FloatArray,
83			cell: FloatArray,
84			name: Optional[str] = None,
85			asym_unit: Optional[UIntArray] = None,
86			multiplicities: Optional[UIntArray] = None,
87			types: Optional[AtomTypeArray] = None,
88			labels: Optional[List[str]] = None,
89			disorder: Optional[list[DisorderAssembly]] = None,
90			):
91			self.motif = motif
92			self.cell = cell
93			self.name = name
94			self._asym_unit = asym_unit
95			self._multiplicities = multiplicities
96			self.types = types
97			self.labels = labels
98			self.disorder = disorder
99
100			@property
101			def ndim(self) -> int:
102			"""Dimension of the ambient space."""
103			return self.cell.shape[0]
104
105			@property
106			def m(self) -> int:
107			"""Number of motif points."""
108			return self.motif.shape[0]
109
110			@property
111			def n_asym(self) -> int:
112			"""Number of points in the asymmetric unit."""
113			if self.asym_unit is None:
114			return self.m
115			else:
116			return self.asym_unit.shape[0]
117
118			@property
119			def asym_unit(self) -> np.ndarray:
120			if self._asym_unit is None:
121			return np.arange(self.m)
122			else:
123			return self._asym_unit
124
125			@property
126			def multiplicities(self) -> np.ndarray:
127			if self._multiplicities is None:
128			return np.ones((self.m,), dtype=np.uint64)
129			else:
130			return self._multiplicities
131
132			@property
133			def occupancies(self) -> np.ndarray:
134			occs = np.ones((self.n_asym, ), dtype=np.float64)
135			if self.disorder:
136			for asm in self.disorder:
137			for grp in asm.groups:
138			occs[grp.indices] = grp.occupancy
139			return occs
140
141			def __str__(self) -> str:
142			m, n = self.motif.shape
143			m_pl = "" if m == 1 else "s"
144			n_pl = "" if n == 1 else "s"
145			name_str = f"{self.name}: " if self.name is not None else ""
146			return f"PeriodicSet({name_str}{m} point{m_pl} in {n} dimension{n_pl})"
147
148			def __repr__(self) -> str:
149			motif_str = str(self.motif).replace("\n ", "\n" + " " * 11)
150			cell_str = str(self.cell).replace("\n ", "\n" + " " * 10)
151
152			optional_attrs = []
153			for attr in ("_asym_unit", "_multiplicities"):
154			val = getattr(self, attr)
155			if val is not None:
156			st = str(val).replace("\n ", "\n" + " " * (len(attr) + 6))
157			optional_attrs.append(f"{attr[1:]}={st}")
158
159			for attr in ("types", "labels"):
160			val = getattr(self, attr)
161			if val is not None:
162			st = str(val).replace("\n ", "\n" + " " * (len(attr) + 6))
163			optional_attrs.append(f"{attr}={st}")
164
165			optional_attrs_str = ",\n " if optional_attrs else ""
166			optional_attrs_str += ",\n ".join(optional_attrs)
167
168			return (
169			f"PeriodicSet(name={self.name},\n"
170			f" motif={motif_str},\n"
171			f" cell={cell_str}{optional_attrs_str})"
172			)
173
174			def PPC(self):
175			r"""Return the point packing coefficient (PPC) of the
176			PeriodicSet.
177
178			The PPC is a constant of any periodic set determining the
179			asymptotic behaviour of its isometry invariants. As
180			:math:`k \rightarrow \infty`, the ratio
181			:math:`\text{AMD}_k / \sqrt[n]{k}` converges to the PPC, as does
182			any row of its PDD.
183
184			For a unit cell :math:`U` and :math:`m` motif points in
185			:math:`n` dimensions,
186
187			.. math::
188
189			\text{PPC} = \sqrt[n]{\frac{\text{Vol}[U]}{m V_n}}
190
191			where :math:`V_n` is the volume of a unit sphere in :math:`n`
192			dimensions.
193
194			Returns
195			-------
196			ppc : float
197			The PPC of ``self``.
198			"""
199			m, n = self.motif.shape
200			t = int(n // 2)
201			if n % 2 == 0:
202			ball_vol = (np.pi ** t) / factorial(t)
203			else:
204			ball_vol = (2 * factorial(t) * (4 * np.pi) ** t) / factorial(n)
205			return (np.abs(np.linalg.det(self.cell)) / (m * ball_vol)) ** (1.0 / n)
206
207			def has_positional_disorder(self):
208			if self.disorder is not None:
209			return any(not asm.is_substitutional for asm in self.disorder)
210			return False
211
212			def has_substitutional_disorder(self):
213			if self.disorder is not None:
214			return any(asm.is_substitutional for asm in self.disorder)
215			return False
216
217			def density(self):
218			"""Return the physical density of the PeriodicSet representing a
219			crystal in g/cm3.
220
221			Returns
222			-------
223			density : float
224			The physical density of ``self`` in g/cm3.
225			"""
226			if self.types is None:
227			raise ValueError(f"PeriodicSet(name={self.name}) has no types")
228			cell_mass = sum(
229			ATOMIC_MASSES[num] * tot_occ
230			for num, tot_occ in zip(self.types, self.multiplicities * self.occupancies)
231			)
232			cell_vol = np.linalg.det(self.cell)
233			return (cell_mass / cell_vol) * 1.66053907
234
235			def formula(self):
236			"""Return the reduced formula of the PeriodicSet representing a
237			crystal.
238
239			Returns
240			-------
241			formula : str
242			The reduced formula of ``self``.
243			"""
244			counter = self._formula_dict()
245			counter = {ATOMIC_NUMBERS_TO_SYMS[n]: count for n, count in counter.items()}
246			return "".join(
247			f"{lab}{num:.3g}" if num != 1 else lab for lab, num in counter.items()
248			)
249
250			def _formula_dict(self):
251			"""Return a dictionary of atomic types in the form
252			{atomic_num: count}.
253			"""
254
255			if self.types is None:
256			raise ValueError(f"PeriodicSet(name={self.name}) has no types")
257
258			counter = {}
259			for num, mul, occ in zip(self.types, self.multiplicities, self.occupancies):
260			if num not in counter:
261			counter[num] = 0
262			counter[num] += occ * mul
263
264			if 0 in counter:
265			del counter[0]
266
267			gcd = math.gcd(*(int(v) for v in counter.values() if np.mod(v, 1.0) == 0.0))
			0 ignored issues – show Comprehensibility Best Practice introduced 2025-06-17 20:28 UTC by Report Bug Copy Issue Report The variable `int` does not seem to be defined. Loading history... Comprehensibility Best Practice introduced 2025-06-17 20:28 UTC by Report Bug Copy Issue Report The variable `v` does not seem to be defined. Loading history...
268			if gcd == 0:
269			gcd = 1
270
271			counter_ = {}
272			if 6 in counter:
273			counter_[6] = counter[6] / gcd
274			del counter[6]
275			if 1 in counter:
276			counter_[1] = counter[1] / gcd
277			del counter[1]
278
279			for n, count in sorted(counter.items(), key=lambda x: x[0]):
280			counter_[int(n)] = count / gcd
281
282			return counter_
283
284			def _equal_cell_and_motif(self, other: PeriodicSet, tol: float = 1e-8) -> bool:
285			"""Used for debugging. True if both 1) unit cells are (close to)
286			identical, and 2) the motifs are the same shape, and
287			every point in one motif has a (close to) identical point
288			somewhere in the other, accounting for pbc.
289			"""
290
291			if (
292			self.cell.shape != other.cell.shape
293			or self.motif.shape != other.motif.shape
294			or not np.allclose(self.cell, other.cell)
295			):
296			return False
297
298			cell_inv = np.linalg.inv(self.cell)
299			fm1 = np.mod(self.motif @ cell_inv, 1)
300			fm2 = np.mod(other.motif @ cell_inv, 1)
301			d1 = np.abs(fm2[:, None] - fm1)
302			d2 = np.abs(d1 - 1)
303			diffs = np.amax(np.minimum(d1, d2), axis=-1)
304
305			if not np.all(
306			(np.amin(diffs, axis=0) <= tol) \| (np.amin(diffs, axis=-1) <= tol)
307			):
308			return False
309			return True
310
311			@staticmethod
312			def cubic(scale: float = 1.0, dims: int = 3) -> PeriodicSet:
313			"""Return a :class:`PeriodicSet` representing a cubic lattice."""
314			return PeriodicSet(np.zeros((1, dims)), np.identity(dims) * scale)
315
316			@staticmethod
317			def hexagonal(scale: float = 1.0, dims: int = 3) -> PeriodicSet:
318			"""Return a :class:`PeriodicSet` representing a hexagonal
319			lattice. Implemented for dimensions 2 and 3 only.
320			"""
321
322			if dims == 3:
323			cellpar = np.array([scale, scale, scale, 90.0, 90.0, 120.0])
324			cell = cellpar_to_cell(cellpar)
325			elif dims == 2:
326			cell = cellpar_to_cell_2D(np.array([scale, scale, 60.0]))
327			else:
328			raise NotImplementedError(
329			"amd.PeriodicSet.hexagonal() only implemented for dimensions "
330			f"2 and 3, passed {dims}"
331			)
332			return PeriodicSet(np.zeros((1, dims)), cell)
333
334			@staticmethod
335			def random(
336			n_points: int,
337			length_bounds: tuple = (1.0, 2.0),
338			angle_bounds: tuple = (60.0, 120.0),
339			dims: int = 3,
340			) -> PeriodicSet:
341			"""Return a :class:`PeriodicSet` with a chosen number of
342			randomly placed points in a random cell with edges between
343			``length_bounds`` and angles between ``angle_bounds``.
344			Implemented for dimensions 2 and 3 only.
345			"""
346
347			if dims not in (2, 3):
348			raise NotImplementedError(
349			"amd.PeriodicSet._random only implemented for dimensions 2 "
350			f"and 3, passed {dims}"
351			)
352
353			cell = random_cell(
354			length_bounds=length_bounds, angle_bounds=angle_bounds, dims=dims
355			)
356			frac_motif = np.random.uniform(size=(n_points, dims))
357			return PeriodicSet(frac_motif @ cell, cell)
358
359			def to_cif_str(self) -> str:
360			"""Return a CIF string representing the PeriodicSet. Does not
361			write disorder data, atom labels or symmetries."""
362
363			a, b, c, al, be, ga = cell_to_cellpar(self.cell)
364			out = [
365			f"data_{self.name}",
366			"_symmetry_space_group_name_H-M 'P 1'",
367			f"_cell_length_a {a}",
368			f"_cell_length_b {b}",
369			f"_cell_length_c {c}",
370			f"_cell_angle_alpha {al}",
371			f"_cell_angle_beta {be}",
372			f"_cell_angle_gamma {ga}",
373			"loop_",
374			" _symmetry_equiv_pos_site_id",
375			" _symmetry_equiv_pos_as_xyz",
376			" 1 'x, y, z'",
377			"loop_",
378			" _atom_site_label",
379			" _atom_site_type_symbol",
380			" _atom_site_fract_x",
381			" _atom_site_fract_y",
382			" _atom_site_fract_z",
383			]
384
385			frac_motif = self.motif @ np.linalg.inv(self.cell)
386			ind = 0
387			label_nums = {}
388			for asym_ind, mul in enumerate(self.multiplicities):
389			for _ in range(mul):
390
391			sym = ATOMIC_NUMBERS_TO_SYMS[self.types[asym_ind]]
392			if sym in label_nums:
393			label_nums[sym] += 1
394			else:
395			label_nums[sym] = 1
396
397			out.append(
398			" ".join(
399			[
400			f"{sym}{label_nums[sym]}",
401			sym,
402			f"{frac_motif[ind, 0]:.4f}",
403			f"{frac_motif[ind, 1]:.4f}",
404			f"{frac_motif[ind, 2]:.4f}",
405			]
406			)
407			)
408			ind += 1
409
410			return "\n".join(out)
411
412			def to_pymatgen_structure(self):
413
414			if self.types is None:
415			raise ValueError(
416			"Cannot convert PeriodicSet to pymatgen Structure without "
417			"atomic types."
418			)
419
420			from pymatgen.core import Structure
421
422			species = [
423			ATOMIC_NUMBERS_TO_SYMS[t]
424			for m, t in zip(self.multiplicities, self.types)
425			for _ in range(m)
426			]
427
428			labels = None
429			if self.labels is not None:
430			labels = [
431			f"{t}_{i+1}"
432			for m, t in zip(self.multiplicities, self.labels)
433			for i in range(m)
434			]
435
436			return Structure(
437			self.cell, species, self.motif, coords_are_cartesian=True, labels=labels
438			)
439
440			def majority_occupancy_config(self) -> PeriodicSet:
441
442			if not self.disorder:
443			return self
444
445			asym_mask = np.full((self.n_asym, ), fill_value=True)
446			motif_mask = np.full((self.m, ), fill_value=True)
447
448			for asm in self.disorder:
449
450			# find majority occupancy group
451			majority_occ, majority_ind = 0, 0
452			for i, grp in enumerate(asm.groups):
453			if grp.occupancy > majority_occ:
454			majority_ind = i
455			majority_occ = grp.occupancy
456
457			# mask all other groups
458			for i, grp in enumerate(asm.groups):
459
460			if i == majority_ind:
461			continue
462
463			for j in grp.indices:
464			asym_mask[j] = False
465			m_i = self.asym_unit[j]
466			mul = self.multiplicities[j]
467			motif_mask[m_i : m_i + mul] = False
468
469			muls = self.multiplicities[asym_mask]
470			asym_unit = np.zeros_like(muls, dtype=np.int64)
471			asym_unit[1:] = np.cumsum(muls)[:-1]
472			masked_labels = [l for l, flag in zip(self.labels, asym_mask) if flag]
473
474			return PeriodicSet(
475			self.motif[motif_mask],
476			self.cell,
477			self.name,
478			asym_unit,
479			self.multiplicities[asym_mask],
480			self.types[asym_mask],
481			masked_labels
482			)
483
484			def compressed_dtypes(self) -> PeriodicSet:
485			"""Return a PeriodicSet using float32, uint32 and uint8 for
486			atomic types."""
487
488			if self._asym_unit is not None:
489			asym_unit = np.array(self._asym_unit, dtype=np.uint32)
490			else:
491			asym_unit = None
492
493			if self._multiplicities is not None:
494			multiplicities = np.array(self._multiplicities, dtype=np.uint32)
495			else:
496			multiplicities = None
497
498			if self.types is not None:
499			types = np.array(self.types, dtype=np.uint8)
500			else:
501			types = None
502
503			return PeriodicSet(
504			np.array(self.motif, dtype=np.float32),
505			np.array(self.cell, dtype=np.float32),
506			name=self.name,
507			asym_unit=asym_unit,
508			multiplicities=multiplicities,
509			types=types,
510			labels=self.labels,
511			disorder=self.disorder
512			)
513

dwiddo / average-minimum-distance

Push — master ( ad57f0...44572e )

amd.periodicset.PeriodicSet._random() A

Complexity

Size

Duplication

Importance

Duplication Side-by-Side

Filter issues like