1
|
|
|
"""Miscellaneous utility functions.""" |
2
|
|
|
|
3
|
|
|
from typing import Tuple |
4
|
|
|
|
5
|
|
|
import numpy as np |
6
|
|
|
import numba |
7
|
|
|
|
8
|
|
|
from ._types import FloatArray |
9
|
|
|
|
10
|
|
|
__all__ = [ |
11
|
|
|
"diameter", |
12
|
|
|
"cellpar_to_cell", |
13
|
|
|
"cellpar_to_cell_2D", |
14
|
|
|
"cell_to_cellpar", |
15
|
|
|
"cell_to_cellpar_2D", |
16
|
|
|
"random_cell", |
17
|
|
|
] |
18
|
|
|
|
19
|
|
|
|
20
|
|
|
@numba.njit(cache=True, fastmath=True) |
21
|
|
|
def diameter(cell: FloatArray) -> float: |
22
|
|
|
"""Return the diameter of a unit cell (given as a square matrix in |
23
|
|
|
orthogonal coordinates) in 3 or fewer dimensions. The diameter is |
24
|
|
|
the maxiumum distance between any two points in the unit cell. |
25
|
|
|
|
26
|
|
|
Parameters |
27
|
|
|
---------- |
28
|
|
|
cell : :class:`numpy.ndarray` |
29
|
|
|
Unit cell represented as a square matrix (in orthogonal |
30
|
|
|
coordinates). |
31
|
|
|
Returns |
32
|
|
|
------- |
33
|
|
|
diameter : float |
34
|
|
|
Diameter of the unit cell. |
35
|
|
|
""" |
36
|
|
|
|
37
|
|
|
dims = cell.shape[0] |
38
|
|
|
if dims == 1: |
39
|
|
|
return cell[0][0] |
40
|
|
|
if dims == 2: |
41
|
|
|
diagonals = np.empty((2, 2), dtype=np.float64) |
42
|
|
|
diagonals[0] = (cell[0] + cell[1]) ** 2 |
43
|
|
|
diagonals[1] = (cell[0] - cell[1]) ** 2 |
44
|
|
|
elif dims == 3: |
45
|
|
|
diagonals = np.empty((4, 3), dtype=np.float64) |
46
|
|
|
diagonals[0] = (cell[0] + cell[1] + cell[2]) ** 2 |
47
|
|
|
diagonals[1] = (cell[0] + cell[1] - cell[2]) ** 2 |
48
|
|
|
diagonals[2] = (cell[0] - cell[1] + cell[2]) ** 2 |
49
|
|
|
diagonals[3] = (-cell[0] + cell[1] + cell[2]) ** 2 |
50
|
|
|
else: |
51
|
|
|
raise ValueError("amd.diameter only implemented for dimensions <= 3") |
52
|
|
|
|
53
|
|
|
return np.sqrt(np.amax(np.sum(diagonals, axis=-1))) |
54
|
|
|
|
55
|
|
|
|
56
|
|
|
@numba.njit(cache=True, fastmath=True) |
57
|
|
|
def cellpar_to_cell(cellpar: FloatArray) -> FloatArray: |
58
|
|
|
"""Convert conventional unit cell parameters a,b,c,α,β,γ into a 3x3 |
59
|
|
|
:class:`numpy.ndarray` representing the unit cell in orthogonal |
60
|
|
|
coordinates. Numba-accelerated version of function from |
61
|
|
|
:mod:`ase.geometry` of the same name. |
62
|
|
|
|
63
|
|
|
Parameters |
64
|
|
|
---------- |
65
|
|
|
cellpar : :class:`numpy.ndarray` |
66
|
|
|
Vector of six unit cell parameters in the canonical order |
67
|
|
|
a,b,c,α,β,γ. |
68
|
|
|
|
69
|
|
|
Returns |
70
|
|
|
------- |
71
|
|
|
cell : :class:`numpy.ndarray` |
72
|
|
|
Unit cell represented as a 3x3 matrix in orthogonal coordinates. |
73
|
|
|
""" |
74
|
|
|
|
75
|
|
|
a, b, c, al, be, ga = cellpar |
76
|
|
|
eps = 2 * np.spacing(90) |
77
|
|
|
|
78
|
|
|
cos_alpha = 0.0 if abs(abs(al) - 90.0) < eps else np.cos(np.deg2rad(al)) |
79
|
|
|
cos_beta = 0.0 if abs(abs(be) - 90.0) < eps else np.cos(np.deg2rad(be)) |
80
|
|
|
cos_gamma = 0.0 if abs(abs(ga) - 90.0) < eps else np.cos(np.deg2rad(ga)) |
81
|
|
|
|
82
|
|
|
if abs(ga - 90.0) < eps: |
83
|
|
|
sin_gamma = 1.0 |
84
|
|
|
elif abs(ga + 90.0) < eps: |
85
|
|
|
sin_gamma = -1.0 |
86
|
|
|
else: |
87
|
|
|
sin_gamma = np.sin(np.deg2rad(ga)) |
88
|
|
|
|
89
|
|
|
cy = (cos_alpha - cos_beta * cos_gamma) / sin_gamma |
90
|
|
|
cz_sqr = 1.0 - cos_beta**2 - cy**2 |
91
|
|
|
if cz_sqr < 0: |
92
|
|
|
raise ValueError("Could not create a unit cell from given parameters.") |
93
|
|
|
|
94
|
|
|
cell = np.zeros((3, 3), dtype=np.float64) |
95
|
|
|
cell[0, 0] = a |
96
|
|
|
cell[1, 0] = b * cos_gamma |
97
|
|
|
cell[1, 1] = b * sin_gamma |
98
|
|
|
cell[2, 0] = c * cos_beta |
99
|
|
|
cell[2, 1] = c * cy |
100
|
|
|
cell[2, 2] = c * np.sqrt(cz_sqr) |
101
|
|
|
return cell |
102
|
|
|
|
103
|
|
|
|
104
|
|
|
@numba.njit(cache=True, fastmath=True) |
105
|
|
|
def cellpar_to_cell_2D(cellpar: FloatArray) -> FloatArray: |
106
|
|
|
"""Convert 3 parameters defining a 2D unit cell a,b,α into a 2x2 |
107
|
|
|
:class:`numpy.ndarray` representing the unit cell in orthogonal |
108
|
|
|
coordinates. |
109
|
|
|
|
110
|
|
|
Parameters |
111
|
|
|
---------- |
112
|
|
|
cellpar : :class:`numpy.ndarray` |
113
|
|
|
Vector of three unit cell parameters a,b,α, where a and b are |
114
|
|
|
lengths of sides of the unit cell and α is the angle between the |
115
|
|
|
sides. |
116
|
|
|
Returns |
117
|
|
|
------- |
118
|
|
|
cell : :class:`numpy.ndarray` |
119
|
|
|
Unit cell represented as a 2x2 matrix in orthogonal coordinates. |
120
|
|
|
""" |
121
|
|
|
|
122
|
|
|
a, b, alpha = cellpar |
123
|
|
|
cell = np.zeros((2, 2), dtype=np.float64) |
124
|
|
|
alpha_radians = np.deg2rad(alpha) |
125
|
|
|
cell[0, 0] = a |
126
|
|
|
cell[1, 0] = b * np.cos(alpha_radians) |
127
|
|
|
cell[1, 1] = b * np.sin(alpha_radians) |
128
|
|
|
return cell |
129
|
|
|
|
130
|
|
|
|
131
|
|
|
@numba.njit(cache=True, fastmath=True) |
132
|
|
|
def cell_to_cellpar(cell: FloatArray) -> FloatArray: |
133
|
|
|
"""Convert a 3x3 :class:`numpy.ndarray` representing a unit cell in |
134
|
|
|
orthogonal coordinates (as returned by |
135
|
|
|
:func:`cellpar_to_cell() <.utils.cellpar_to_cell>`) into a list of 3 |
136
|
|
|
lengths and 3 angles representing the unit cell. |
137
|
|
|
|
138
|
|
|
Parameters |
139
|
|
|
---------- |
140
|
|
|
cell : :class:`numpy.ndarray` |
141
|
|
|
Unit cell as a 3x3 matrix in orthogonal coordinates. |
142
|
|
|
|
143
|
|
|
Returns |
144
|
|
|
------- |
145
|
|
|
cellpar : :class:`numpy.ndarray` |
146
|
|
|
Vector of 6 unit cell parameters in the canonical order |
147
|
|
|
a,b,c,α,β,γ. |
148
|
|
|
""" |
149
|
|
|
|
150
|
|
|
cellpar = np.empty((6,), dtype=np.float64) |
151
|
|
|
cellpar[0] = np.linalg.norm(cell[0]) |
152
|
|
|
cellpar[1] = np.linalg.norm(cell[1]) |
153
|
|
|
cellpar[2] = np.linalg.norm(cell[2]) |
154
|
|
|
cellpar[3] = np.rad2deg( |
155
|
|
|
np.arccos(np.dot(cell[1], cell[2]) / (cellpar[1] * cellpar[2])) |
156
|
|
|
) |
157
|
|
|
cellpar[4] = np.rad2deg( |
158
|
|
|
np.arccos(np.dot(cell[0], cell[2]) / (cellpar[0] * cellpar[2])) |
159
|
|
|
) |
160
|
|
|
cellpar[5] = np.rad2deg( |
161
|
|
|
np.arccos(np.dot(cell[0], cell[1]) / (cellpar[0] * cellpar[1])) |
162
|
|
|
) |
163
|
|
|
return cellpar |
164
|
|
|
|
165
|
|
|
|
166
|
|
|
@numba.njit(cache=True, fastmath=True) |
167
|
|
|
def cell_to_cellpar_2D(cell: FloatArray) -> FloatArray: |
168
|
|
|
"""Convert a 2x2 :class:`numpy.ndarray` representing a unit cell in |
169
|
|
|
orthogonal coordinates (as returned by |
170
|
|
|
:func:`cellpar_to_cell_2D() <.utils.cellpar_to_cell_2D>`) into a |
171
|
|
|
list of 2 lengths and 1 angle representing the unit cell. |
172
|
|
|
|
173
|
|
|
Parameters |
174
|
|
|
---------- |
175
|
|
|
cell : :class:`numpy.ndarray` |
176
|
|
|
Unit cell as a 2x2 matrix in orthogonal coordinates. |
177
|
|
|
|
178
|
|
|
Returns |
179
|
|
|
------- |
180
|
|
|
cellpar : :class:`numpy.ndarray` |
181
|
|
|
Vector with 3 elements, containing 3 unit cell parameters. |
182
|
|
|
""" |
183
|
|
|
|
184
|
|
|
cellpar = np.empty((3,), dtype=np.float64) |
185
|
|
|
cellpar[0] = np.linalg.norm(cell[0]) |
186
|
|
|
cellpar[1] = np.linalg.norm(cell[1]) |
187
|
|
|
cellpar[2] = np.rad2deg( |
188
|
|
|
np.arccos(np.dot(cell[0], cell[1]) / (cellpar[0] * cellpar[1])) |
189
|
|
|
) |
190
|
|
|
return cellpar |
191
|
|
|
|
192
|
|
|
|
193
|
|
|
def random_cell( |
194
|
|
|
length_bounds: Tuple = (1.0, 2.0), |
195
|
|
|
angle_bounds: Tuple = (60.0, 120.0), |
196
|
|
|
dims: int = 3, |
197
|
|
|
) -> FloatArray: |
198
|
|
|
"""Return a random unit cell with uniformally chosen length and |
199
|
|
|
angle parameters between bounds. Dimensions 2 and 3 only. |
200
|
|
|
""" |
201
|
|
|
|
202
|
|
|
ll, lu = length_bounds |
203
|
|
|
al, au = angle_bounds |
204
|
|
|
|
205
|
|
|
if dims == 2: |
206
|
|
|
cellpar = [np.random.uniform(low=ll, high=lu) for _ in range(2)] |
207
|
|
|
cellpar.append(np.random.uniform(low=al, high=au)) |
208
|
|
|
cell = cellpar_to_cell_2D(np.array(cellpar)) |
209
|
|
|
elif dims == 3: |
210
|
|
|
while True: |
211
|
|
|
lengths = [np.random.uniform(low=ll, high=lu) for _ in range(dims)] |
212
|
|
|
angles = [np.random.uniform(low=al, high=au) for _ in range(dims)] |
213
|
|
|
cellpar = np.array(lengths + angles) |
214
|
|
|
try: |
215
|
|
|
cell = cellpar_to_cell(cellpar) |
216
|
|
|
break |
217
|
|
|
except RuntimeError: |
218
|
|
|
continue |
219
|
|
|
else: |
220
|
|
|
raise NotImplementedError( |
221
|
|
|
"amd.random_cell() only implemented for dimensions 2 and 3, " |
222
|
|
|
f"passed {dims}" |
223
|
|
|
) |
224
|
|
|
return cell |
225
|
|
|
|