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