1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
|
import numpy as np
from ase.geometry import wrap_positions
# Order in which off-diagonal elements are checked for strong tilt
FLIP_ORDER = [(1, 0, 0), (2, 0, 0), (2, 1, 1)]
class Prism:
"""The representation of the unit cell in LAMMPS
The main purpose of the prism-object is to create suitable
string representations of prism limits and atom positions
within the prism.
:param cell: cell in ase coordinate system
:param pbc: periodic boundaries
:param tolerance: precision for skewness test
**Implementation**
LAMMPS requires triangular matrixes without a strong tilt.
Therefore the 'Prism'-object contains three coordinate systems:
- ase_cell (the simulated system in the ASE coordination system)
- lammps_tilt (ase-cell rotated to be an lower triangular matrix)
- lammps_cell (same volume as tilted cell, but reduce edge length)
The translation between 'ase_cell' and 'lammps_cell' is down with a
rotation matrix 'rot_mat' obtained from a QR decomposition. The
transformation between 'lammps_tilt' and 'lammps_cell' is done by
changing the off-diagonal elements. 'flip' saves the modified
elements for later reversal. All vectors except positions are just
rotated with 'rot_mat'. Positions are rotated and wrapped into
'lammps_cell'. Translating results back from LAMMPS needs first the
unwrapping of positions. Then all vectors and the unit-cell are
rotated back into the ASE coordination system. This can fail as
depending on the simulation run LAMMPS might have changed the
simulation box significantly. This is for example a problem with
hexagonal cells. LAMMPS might also wrap atoms across periodic
boundaries, which can lead to problems for example NEB
calculations.
"""
# !TODO: derive tolerance from cell-dimensions
def __init__(self, cell, pbc=(True, True, True), tolerance=1.0e-8):
# Use QR decomposition to get the lammps cell
# rot_mat * lammps_tilt^T = ase_cell^T
# => lammps_tilt * rot_mat^T = ase_cell
# => lammps_tilt = ase_cell * rot_mat
qtrans, ltrans = np.linalg.qr(cell.T, mode="complete")
self.rot_mat = qtrans
self.lammps_tilt = ltrans.T
self.ase_cell = cell
self.tolerance = tolerance
self.pbc = pbc
# LAMMPS requires positive values on the diagonal of the
# triangluar matrix -> mirror if necessary
for i in range(3):
if self.lammps_tilt[i][i] < 0.0:
mirror_mat = np.eye(3)
mirror_mat[i][i] = -1.0
self.lammps_tilt = np.dot(mirror_mat, self.lammps_tilt.T).T
self.rot_mat = np.dot(self.rot_mat, mirror_mat)
self.lammps_cell = self.lammps_tilt.copy()
# LAMMPS minimizes the edge length of the parallelepiped
# What is ment with 'flip': cell 2 is transformed into cell 1
# cell 2 = 'lammps_tilt'; cell 1 = 'lammps_cell'
# o-----------------------------/==o-----------------------------/--o
# \ /--/ \ /--/
# \ /--/ \ /--/
# \ 1 /--/ \ 2 /--/
# \ /--/ \ /--/
# \ /--/ \ /--/
# o==/-----------------------------o--/
# !TODO: handle extreme tilt (= off-diagonal > 1.5)
self.flip = np.array(
[
abs(self.lammps_cell[i][j] / self.lammps_tilt[k][k]) > 0.5
and self.pbc[k]
for i, j, k in FLIP_ORDER
]
)
for iteri, (i, j, k) in enumerate(FLIP_ORDER):
if self.flip[iteri]:
change = self.lammps_cell[k][k]
change *= np.sign(self.lammps_cell[i][j])
self.lammps_cell[i][j] -= change
def get_lammps_prism(self):
"""Return into lammps coordination system rotated cell
:returns: lammps cell
:rtype: np.array
"""
return self.lammps_cell[(0, 1, 2, 1, 2, 2), (0, 1, 2, 0, 0, 1)]
# !TODO: detect flip in lammps
def update_cell(self, lammps_cell):
"""Rotate new lammps cell into ase coordinate system
:param lammps_cell: new lammps cell received after executing lammps
:returns: ase cell
:rtype: np.array
"""
self.lammps_cell = lammps_cell
self.lammps_tilt = self.lammps_cell.copy()
# reverse flip
for iteri, (i, j, k) in enumerate(FLIP_ORDER):
if self.flip[iteri]:
change = self.lammps_cell[k][k]
change *= np.sign(self.lammps_cell[i][j])
self.lammps_tilt[i][j] -= change
# try to detect potential flips in lammps
# (lammps minimizes the cell-vector lengths)
new_ase_cell = np.dot(self.lammps_tilt, self.rot_mat.T)
# assuming the cell changes are mostly isotropic
new_vol = np.linalg.det(new_ase_cell)
old_vol = np.linalg.det(self.ase_cell)
test_residual = self.ase_cell.copy()
test_residual *= (new_vol / old_vol) ** (1.0 / 3.0)
test_residual -= new_ase_cell
if any(
np.linalg.norm(test_residual, axis=1)
> 0.5 * np.linalg.norm(self.ase_cell, axis=1)
):
print(
"WARNING: Significant simulation cell changes from LAMMPS "
+ "detected.\n"
+ " " * 9
+ "Backtransformation to ASE might fail!"
)
return new_ase_cell
def vector_to_lammps(self, vec, wrap=False):
"""Rotate vector from ase coordinate system to lammps one
:param vec: to be rotated ase-vector
:returns: lammps-vector
:rtype: np.array
"""
# !TODO: right eps-limit
# lammps might not like atoms outside the cell
if wrap:
return wrap_positions(
np.dot(vec, self.rot_mat),
self.lammps_cell,
pbc=self.pbc,
eps=1e-18,
)
return np.dot(vec, self.rot_mat)
def vector_to_ase(self, vec, wrap=False):
"""Rotate vector from lammps coordinate system to ase one
:param vec: to be rotated lammps-vector
:param wrap: was vector wrapped into 'lammps_cell'
:returns: ase-vector
:rtype: np.array
"""
if wrap:
# trying to reverse wraping across periodic boundaries in the
# lammps coordination-system:
# translate: expresses lammps-coordinate system in the rotate, but
# without tilt removed system
# fractional: vector in tilted system
translate = np.linalg.solve(
self.lammps_tilt.T, self.lammps_cell.T
).T
fractional = np.linalg.solve(self.lammps_tilt.T, vec.T).T
# !TODO: make somehow nicer
# !TODO: handle extreme tilted cells
for ifrac in fractional:
for zyx in reversed(range(3)):
if ifrac[zyx] >= 1.0 and self.pbc[zyx]:
ifrac -= translate[zyx]
elif ifrac[zyx] < 0.0 and self.pbc[zyx]:
ifrac += translate[zyx]
vec = np.dot(fractional, self.lammps_tilt)
return np.dot(vec, self.rot_mat.T)
def is_skewed(self):
"""Test if a lammps cell is not tetragonal
:returns: bool
:rtype: bool
"""
cell_sq = self.lammps_cell ** 2
return (
np.sum(np.tril(cell_sq, -1)) / np.sum(np.diag(cell_sq))
> self.tolerance
)
|
