Coverage for /builds/kinetik161/ase/ase/calculators/lammps/coordinatetransform.py: 95.83%
72 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +0000
1"""Prism"""
2import warnings
3from typing import Sequence
5import numpy as np
7from ase.geometry import wrap_positions
10def calc_box_parameters(cell: np.ndarray) -> np.ndarray:
11 """Calculate box parameters
13 https://docs.lammps.org/Howto_triclinic.html
14 """
15 ax = np.sqrt(cell[0] @ cell[0])
16 bx = cell[0] @ cell[1] / ax
17 by = np.sqrt(cell[1] @ cell[1] - bx ** 2)
18 cx = cell[0] @ cell[2] / ax
19 cy = (cell[1] @ cell[2] - bx * cx) / by
20 cz = np.sqrt(cell[2] @ cell[2] - cx ** 2 - cy ** 2)
21 return np.array((ax, by, cz, bx, cx, cy))
24def calc_rotated_cell(cell: np.ndarray) -> np.ndarray:
25 """Calculate rotated cell in LAMMPS coordinates
27 Parameters
28 ----------
29 cell : np.ndarray
30 Cell to be rotated.
32 Returns
33 -------
34 rotated_cell : np.ndarray
35 Rotated cell represented by a lower triangular matrix.
36 """
37 ax, by, cz, bx, cx, cy = calc_box_parameters(cell)
38 return np.array(((ax, 0.0, 0.0), (bx, by, 0.0), (cx, cy, cz)))
41def calc_reduced_cell(cell: np.ndarray, pbc: Sequence[bool]) -> np.ndarray:
42 """Calculate LAMMPS cell with short lattice basis vectors
44 The lengths of the second and the third lattice basis vectors, b and c, are
45 shortened with keeping the same periodicity of the system. b is modified by
46 adding multiple a vectors, and c is modified by adding first multiple b
47 vectors and then multiple a vectors.
49 Parameters
50 ----------
51 cell : np.ndarray
52 Cell to be reduced. This must be already a lower triangular matrix.
53 pbc : Sequence[bool]
54 True if the system is periodic along the corresponding direction.
56 Returns
57 -------
58 reduced_cell : np.ndarray
59 Reduced cell. `xx`, `yy`, `zz` are the same as the original cell,
60 and `abs(xy) <= xx`, `abs(xz) <= xx`, `abs(yz) <= yy`.
61 """
62 # cell 1 (reduced) <- cell 2 (original)
63 # o-----------------------------/==o-----------------------------/--o
64 # \ /--/ \ /--/
65 # \ /--/ \ /--/
66 # \ 1 /--/ \ 2 /--/
67 # \ /--/ \ /--/
68 # \ /--/ \ /--/
69 # o==/-----------------------------o--/
70 reduced_cell = cell.copy()
71 # Order in which off-diagonal elements are checked for strong tilt
72 # yz is updated before xz so that the latter does not affect the former
73 flip_order = ((1, 0), (2, 1), (2, 0))
74 for i, j in flip_order:
75 if not pbc[j]:
76 continue
77 ratio = reduced_cell[i, j] / reduced_cell[j, j]
78 if abs(ratio) > 0.5:
79 reduced_cell[i] -= reduced_cell[j] * np.round(ratio)
80 return reduced_cell
83class Prism:
84 """The representation of the unit cell in LAMMPS
86 The main purpose of the prism-object is to create suitable string
87 representations of prism limits and atom positions within the prism.
89 Parameters
90 ----------
91 cell : np.ndarray
92 Cell in ASE coordinate system.
93 pbc : one or three bool
94 Periodic boundary conditions flags.
95 reduce_cell : bool
96 If True, the LAMMPS cell is reduced for short lattice basis vectors.
97 The atomic positions are always wraped into the reduced cell,
98 regardress of `wrap` in `vector_to_lammps` and `vector_to_ase`.
99 tolerance : float
100 Precision for skewness test.
102 Methods
103 -------
104 vector_to_lammps
105 Rotate vectors from ASE to LAMMPS coordinates.
106 Positions can be further wrapped into the LAMMPS cell by `wrap=True`.
108 vector_to_ase
109 Rotate vectors from LAMMPS to ASE coordinates.
110 Positions can be further wrapped into the LAMMPS cell by `wrap=True`.
112 Notes
113 -----
114 LAMMPS prefers triangular matrixes without a strong tilt.
115 Therefore the 'Prism'-object contains three coordinate systems:
117 - ase_cell (the simulated system in the ASE coordination system)
118 - lammps_tilt (ase-cell rotated to be an lower triangular matrix)
119 - lammps_cell (same volume as tilted cell, but reduce edge length)
121 The translation between 'ase_cell' and 'lammps_tilt' is done with a
122 rotation matrix 'rot_mat'. (Mathematically this is a QR decomposition.)
124 The transformation between 'lammps_tilt' and 'lammps_cell' is done by
125 changing the off-diagonal elements.
127 Depending on the option `reduce`, vectors in ASE coordinates are
128 transformed either `lammps_tilt` or `lammps_cell`.
130 The vector conversion can fail as depending on the simulation run LAMMPS
131 might have changed the simulation box significantly. This is for example a
132 problem with hexagonal cells. LAMMPS might also wrap atoms across periodic
133 boundaries, which can lead to problems for example NEB calculations.
134 """
136 # !TODO: derive tolerance from cell-dimensions
137 def __init__(
138 self,
139 cell: np.ndarray,
140 pbc: bool = True,
141 reduce_cell: bool = False,
142 tolerance: float = 1.0e-8,
143 ):
144 # rot_mat * lammps_tilt^T = ase_cell^T
145 # => lammps_tilt * rot_mat^T = ase_cell
146 # => lammps_tilt = ase_cell * rot_mat
147 # LAMMPS requires positive diagonal elements of the triangular matrix.
148 # The diagonals of `lammps_tilt` are always positive by construction.
149 self.lammps_tilt = calc_rotated_cell(cell)
150 self.rot_mat = np.linalg.solve(self.lammps_tilt, cell).T
151 self.ase_cell = cell
152 self.tolerance = tolerance
153 self.pbc = np.zeros(3, bool) + pbc
154 self.lammps_cell = calc_reduced_cell(self.lammps_tilt, self.pbc)
155 self.is_reduced = reduce_cell
157 @property
158 def cell(self) -> np.ndarray:
159 return self.lammps_cell if self.is_reduced else self.lammps_tilt
161 def get_lammps_prism(self) -> np.ndarray:
162 """Return box parameters of the rotated cell in LAMMPS coordinates
164 Returns
165 -------
166 np.ndarray
167 xhi - xlo, yhi - ylo, zhi - zlo, xy, xz, yz
168 """
169 return self.cell[(0, 1, 2, 1, 2, 2), (0, 1, 2, 0, 0, 1)]
171 def update_cell(self, lammps_cell: np.ndarray) -> np.ndarray:
172 """Rotate new LAMMPS cell into ASE coordinate system
174 Parameters
175 ----------
176 lammps_cell : np.ndarray
177 New Cell in LAMMPS coordinates received after executing LAMMPS
179 Returns
180 -------
181 np.ndarray
182 New cell in ASE coordinates
183 """
184 # Transformation: integer matrix
185 # lammps_cell * transformation = lammps_tilt
186 transformation = np.linalg.solve(self.lammps_cell, self.lammps_tilt)
188 if self.is_reduced:
189 self.lammps_cell = lammps_cell
190 self.lammps_tilt = lammps_cell @ transformation
191 else:
192 self.lammps_tilt = lammps_cell
193 self.lammps_cell = calc_reduced_cell(self.lammps_tilt, self.pbc)
195 # try to detect potential flips in lammps
196 # (lammps minimizes the cell-vector lengths)
197 new_ase_cell = self.lammps_tilt @ self.rot_mat.T
198 # assuming the cell changes are mostly isotropic
199 new_vol = np.linalg.det(new_ase_cell)
200 old_vol = np.linalg.det(self.ase_cell)
201 test_residual = self.ase_cell.copy()
202 test_residual *= (new_vol / old_vol) ** (1.0 / 3.0)
203 test_residual -= new_ase_cell
204 if any(
205 np.linalg.norm(test_residual, axis=1)
206 > 0.5 * np.linalg.norm(self.ase_cell, axis=1)
207 ):
208 warnings.warn(
209 "Significant simulation cell changes from LAMMPS detected. "
210 "Backtransformation to ASE might fail!"
211 )
212 return new_ase_cell
214 def vector_to_lammps(
215 self,
216 vec: np.ndarray,
217 wrap: bool = False,
218 ) -> np.ndarray:
219 """Rotate vectors from ASE to LAMMPS coordinates
221 Parameters
222 ----------
223 vec : np.ndarray
224 Vectors in ASE coordinates to be rotated into LAMMPS coordinates
225 wrap : bool
226 If True, the vectors are wrapped into the cell
228 Returns
229 -------
230 np.array
231 Vectors in LAMMPS coordinates
232 """
233 # !TODO: right eps-limit
234 # lammps might not like atoms outside the cell
235 if wrap or self.is_reduced:
236 return wrap_positions(
237 vec @ self.rot_mat,
238 cell=self.cell,
239 pbc=self.pbc,
240 eps=1e-18,
241 )
242 return vec @ self.rot_mat
244 def vector_to_ase(
245 self,
246 vec: np.ndarray,
247 wrap: bool = False,
248 ) -> np.ndarray:
249 """Rotate vectors from LAMMPS to ASE coordinates
251 Parameters
252 ----------
253 vec : np.ndarray
254 Vectors in LAMMPS coordinates to be rotated into ASE coordinates
255 wrap : bool
256 If True, the vectors are wrapped into the cell
258 Returns
259 -------
260 np.ndarray
261 Vectors in ASE coordinates
262 """
263 if wrap or self.is_reduced:
264 # fractional in `lammps_tilt` (the same shape as ASE cell)
265 fractional = np.linalg.solve(self.lammps_tilt.T, vec.T).T
266 # wrap into 0 to 1 for periodic directions
267 fractional -= np.floor(fractional) * self.pbc
268 # Cartesian coordinates wrapped into `lammps_tilt`
269 vec = fractional @ self.lammps_tilt
270 # rotate back to the ASE cell
271 return vec @ self.rot_mat.T
273 def tensor2_to_ase(self, tensor: np.ndarray) -> np.ndarray:
274 """Rotate a second order tensor from LAMMPS to ASE coordinates
276 Parameters
277 ----------
278 tensor : np.ndarray
279 Tensor in LAMMPS coordinates to be rotated into ASE coordinates
281 Returns
282 -------
283 np.ndarray
284 Tensor in ASE coordinates
285 """
286 return self.rot_mat @ tensor @ self.rot_mat.T
288 def is_skewed(self) -> bool:
289 """Test if the lammps cell is skewed, i.e., monoclinic or triclinic.
291 Returns
292 -------
293 bool
294 True if the lammps cell is skewed.
295 """
296 cell_sq = self.cell ** 2
297 on_diag = np.sum(np.diag(cell_sq))
298 off_diag = np.sum(np.tril(cell_sq, -1))
299 return off_diag / on_diag > self.tolerance