Passed
Push — master ( ae6bce...0853a7 )
by Daniel
03:53
created

amd.periodicset.PeriodicSet.density()   A

Complexity

Conditions 2

Size

Total Lines 7
Code Lines 7

Duplication

Lines 0
Ratio 0 %

Importance

Changes 0
Metric Value
eloc 7
dl 0
loc 7
rs 10
c 0
b 0
f 0
cc 2
nop 1
1
"""Implements the :class:`PeriodicSet` class representing a periodic
2
set, defined by a motif and unit cell. This models a crystal with a
3
point at the center of each atom.
4
5
This is the type yielded by the readers 
6
:class:`amd.CifReader <.io.CifReader>` and 
7
:class:`amd.CSDReader <.io.CSDReader>`. A :class:`PeriodicSet` can be
8
passed as the first argument to :func:`amd.AMD() <.calculate.AMD>` or
9
:func:`amd.PDD() <.calculate.PDD>` to calculate its invariants.
10
"""
11
12
from typing import Optional
13
import numpy as np
14
15
from .utils import (
16
    cellpar_to_cell,
17
    cellpar_to_cell_2D,
18
    cell_to_cellpar,
19
    cell_to_cellpar_2D,
20
    random_cell,
21
)
22
23
from ase.data import atomic_masses
24
25
26
class PeriodicSet:
27
    """A periodic set representats a crystal by putting a point in the
28
    center of each atom. Mathematically it's defined by a basis (unit 
29
    cell) and collection of points (motif) which repeats according to
30
    the basis. Periodic sets are defined for any dimension.
31
32
    :class:`PeriodicSet` s are returned by the readers in the
33
    :mod:`.io` module. They can be passed to
34
    :func:`amd.AMD() <.calculate.AMD>` or 
35
    :func:`amd.PDD() <.calculate.PDD>` to calculate the invariants.
36
37
    Parameters
38
    ----------
39
    motif : :class:`numpy.ndarray`
40
        Cartesian (orthogonal) coordinates of the motif, shape (no
41
        points, dims).
42
    cell : :class:`numpy.ndarray`
43
        Cartesian (orthogonal) square array representing the unit cell,
44
        shape (dims, dims). Use
45
        :func:`amd.cellpar_to_cell <.utils.cellpar_to_cell>` to convert
46
        6 cell parameters to an orthogonal square matrix.
47
    name : str, optional
48
        Name of the periodic set.
49
    asymmetric_unit : :class:`numpy.ndarray`, optional
50
        Indices for the asymmetric unit, pointing to the motif. Useful
51
        in invariant calculations.
52
    wyckoff_multiplicities : :class:`numpy.ndarray`, optional
53
        Wyckoff multiplicities of each atom in the asymmetric unit
54
        (number of unique sites generated under all symmetries). Useful
55
        in invariant calculations.
56
    types : :class:`numpy.ndarray`, optional
57
        Array of atomic numbers of motif points.
58
    """
59
60
    def __init__(
61
            self,
62
            motif: np.ndarray,
63
            cell: np.ndarray,
64
            name: Optional[str] = None,
65
            asymmetric_unit: Optional[np.ndarray] = None,
66
            wyckoff_multiplicities: Optional[np.ndarray] = None,
67
            types: Optional[np.ndarray] = None
68
    ):
69
70
        self.motif = motif
71
        self.cell = cell
72
        self.name = name
73
        self.asymmetric_unit = asymmetric_unit
74
        self.wyckoff_multiplicities = wyckoff_multiplicities
75
        self.types = types
76
77
    @property
78
    def ndim(self):
79
        return self.cell.shape[0]
80
81
    @property
82
    def density(self):
83
        if self.types is None:
84
            raise ValueError('')
85
        cell_mass = sum(atomic_masses[n] for n in self.types)
86
        cell_vol = np.linalg.det(self.cell)
87
        return cell_mass / cell_vol
88
89
    def __str__(self):
90
        """Returns a string representation of the format
91
        PeriodicSet(name, motif (m, d), t asym sites, 
92
        abcαβγ=cellpar).
93
        """
94
95
        if self.asymmetric_unit is None:
96
            n_asym_sites = len(self.motif)
97
        else:
98
            n_asym_sites = len(self.asymmetric_unit)
99
100
        cellpar_str = None
101
        if self.ndim == 1:
102
            cellpar_str = f', cell={self.cell[0][0]}'
103
        if self.ndim == 2:
104
            cellpar = np.round(cell_to_cellpar_2D(self.cell), 2)
105
            cellpar_str = ','.join(str(round(p, 2)) for p in cellpar)
106
            cellpar_str = f', abα={cellpar_str}'
107
        elif self.ndim == 3:
108
            cellpar = np.round(cell_to_cellpar(self.cell), 2)
109
            cellpar_str = ','.join(str(round(p, 2)) for p in cellpar)
110
            cellpar_str = f', abcαβγ={cellpar_str}'
111
        else:
112
            cellpar_str = f''
113
114
        s = f'PeriodicSet(name={self.name}, ' \
115
            f'motif {self.motif.shape} ({n_asym_sites} asym sites)' \
116
            f'{cellpar_str})'
117
118
        return s
119
120
    def __repr__(self):
121
122
        if self.asymmetric_unit is None:
123
            n_asym_sites = len(self.motif)
124
        else:
125
            n_asym_sites = len(self.asymmetric_unit)
126
127
        s = f'PeriodicSet(name={self.name}, ' \
128
            f'motif={self.motif} ({n_asym_sites} asym sites), ' \
129
            f'cell={self.cell})'
130
131
        return s
132
133
    
134
    def __eq__(self, other):
135
        """Used for debugging/tests. Returns True if:
136
        1. The unit cells are identical (element for element within tol)
137
        2. The motifs are the same shape, and when transformed into
138
        fractional coordinates every point in one motif has an identical
139
        (within tol) point somewhere in the other, accounting for pbc.
140
        """
141
142
        if self.cell.shape != other.cell.shape or \
143
           self.motif.shape != other.motif.shape or \
144
           not np.allclose(self.cell, other.cell):
145
            return False
146
147
        fm1 = np.mod(self.motif @ np.linalg.inv(self.cell), 1)
148
        fm2 = np.mod(other.motif @ np.linalg.inv(other.cell), 1)
149
        d1 = np.abs(fm2[:, None] - fm1)
150
        d2 = np.abs(d1 - 1)
151
        diffs = np.amax(np.minimum(d1, d2), axis=-1)
152
153
        if not np.all((np.amin(diffs, axis=0) <= 1e-8) | (np.amin(diffs, axis=-1) <= 1e-8)):
154
            return False
155
156
        return True
157
158
    def __ne__(self, other):
159
        return not self.__eq__(other)
160
161
    @staticmethod
162
    def cubic(scale=1, dims=3):
163
        """Return a :class:`PeriodicSet` representing a cubic lattice.
164
        """
165
166
        cell = np.identity(dims) * scale
167
        return PeriodicSet(np.zeros((1, dims)), cell)
168
169
    @staticmethod
170
    def hexagonal(scale=1, dims=3):
171
        """Dimensions 2 and 3 only. Return a :class:`PeriodicSet`
172
        representing a hexagonal lattice.
173
        """
174
175
        if dims == 3:
176
            cell = cellpar_to_cell(scale, scale, scale, 90, 90, 120)
177
        elif dims == 2:
178
            cell = cellpar_to_cell_2D(scale, scale, 60)
179
        else:
180
            msg = 'PeriodicSet.hexagonal() only implemented for dimensions ' \
181
                  f'2 and 3 (passed {dims})'
182
            raise NotImplementedError(msg)
183
184
        return PeriodicSet(np.zeros((1, dims)), cell)
185
186
    @staticmethod
187
    def _random(
188
        n_points,
189
        length_bounds=(1, 2),
190
        angle_bounds=(60, 120),
191
        dims=3
192
    ):
193
        """Dimensions 2 and 3 only. Return a :class:`PeriodicSet` with a
194
        chosen number of randomly placed points, in random cell with
195
        edges between length_bounds and angles between angle_bounds.
196
        """
197
198
        cell = random_cell(
199
            length_bounds=length_bounds, 
200
            angle_bounds=angle_bounds, 
201
            dims=dims
202
        )
203
        frac_motif = np.random.uniform(size=(n_points, dims))
204
        return PeriodicSet(frac_motif @ cell, cell)
205