1import numpy as np 

2from ase.geometry import wrap_positions 

3 

4# Order in which off-diagonal elements are checked for strong tilt 

5FLIP_ORDER = [(1, 0, 0), (2, 0, 0), (2, 1, 1)] 

6 

7 

8class Prism: 

9 """The representation of the unit cell in LAMMPS 

10 

11 The main purpose of the prism-object is to create suitable 

12 string representations of prism limits and atom positions 

13 within the prism. 

14 

15 :param cell: cell in ase coordinate system 

16 :param pbc: periodic boundaries 

17 :param tolerance: precision for skewness test 

18 

19 **Implementation** 

20 

21 LAMMPS requires triangular matrixes without a strong tilt. 

22 Therefore the 'Prism'-object contains three coordinate systems: 

23 

24 - ase_cell (the simulated system in the ASE coordination system) 

25 - lammps_tilt (ase-cell rotated to be an lower triangular matrix) 

26 - lammps_cell (same volume as tilted cell, but reduce edge length) 

27 

28 The translation between 'ase_cell' and 'lammps_cell' is down with a 

29 rotation matrix 'rot_mat' obtained from a QR decomposition. The 

30 transformation between 'lammps_tilt' and 'lammps_cell' is done by 

31 changing the off-diagonal elements. 'flip' saves the modified 

32 elements for later reversal. All vectors except positions are just 

33 rotated with 'rot_mat'. Positions are rotated and wrapped into 

34 'lammps_cell'. Translating results back from LAMMPS needs first the 

35 unwrapping of positions. Then all vectors and the unit-cell are 

36 rotated back into the ASE coordination system. This can fail as 

37 depending on the simulation run LAMMPS might have changed the 

38 simulation box significantly. This is for example a problem with 

39 hexagonal cells. LAMMPS might also wrap atoms across periodic 

40 boundaries, which can lead to problems for example NEB 

41 calculations. 

42 """ 

43 

44 # !TODO: derive tolerance from cell-dimensions 

45 def __init__(self, cell, pbc=(True, True, True), tolerance=1.0e-8): 

46 # Use QR decomposition to get the lammps cell 

47 # rot_mat * lammps_tilt^T = ase_cell^T 

48 # => lammps_tilt * rot_mat^T = ase_cell 

49 # => lammps_tilt = ase_cell * rot_mat 

50 qtrans, ltrans = np.linalg.qr(cell.T, mode="complete") 

51 self.rot_mat = qtrans 

52 self.lammps_tilt = ltrans.T 

53 self.ase_cell = cell 

54 self.tolerance = tolerance 

55 self.pbc = pbc 

56 

57 # LAMMPS requires positive values on the diagonal of the 

58 # triangluar matrix -> mirror if necessary 

59 for i in range(3): 

60 if self.lammps_tilt[i][i] < 0.0: 

61 mirror_mat = np.eye(3) 

62 mirror_mat[i][i] = -1.0 

63 self.lammps_tilt = np.dot(mirror_mat, self.lammps_tilt.T).T 

64 self.rot_mat = np.dot(self.rot_mat, mirror_mat) 

65 

66 self.lammps_cell = self.lammps_tilt.copy() 

67 

68 # LAMMPS minimizes the edge length of the parallelepiped 

69 # What is ment with 'flip': cell 2 is transformed into cell 1 

70 # cell 2 = 'lammps_tilt'; cell 1 = 'lammps_cell' 

71 # o-----------------------------/==o-----------------------------/--o 

72 # \ /--/ \ /--/ 

73 # \ /--/ \ /--/ 

74 # \ 1 /--/ \ 2 /--/ 

75 # \ /--/ \ /--/ 

76 # \ /--/ \ /--/ 

77 # o==/-----------------------------o--/ 

78 # !TODO: handle extreme tilt (= off-diagonal > 1.5) 

79 self.flip = np.array( 

80 [ 

81 abs(self.lammps_cell[i][j] / self.lammps_tilt[k][k]) > 0.5 

82 and self.pbc[k] 

83 for i, j, k in FLIP_ORDER 

84 ] 

85 ) 

86 for iteri, (i, j, k) in enumerate(FLIP_ORDER): 

87 if self.flip[iteri]: 

88 change = self.lammps_cell[k][k] 

89 change *= np.sign(self.lammps_cell[i][j]) 

90 self.lammps_cell[i][j] -= change 

91 

92 def get_lammps_prism(self): 

93 """Return into lammps coordination system rotated cell 

94 

95 :returns: lammps cell 

96 :rtype: np.array 

97 

98 """ 

99 return self.lammps_cell[(0, 1, 2, 1, 2, 2), (0, 1, 2, 0, 0, 1)] 

100 

101 # !TODO: detect flip in lammps 

102 def update_cell(self, lammps_cell): 

103 """Rotate new lammps cell into ase coordinate system 

104 

105 :param lammps_cell: new lammps cell received after executing lammps 

106 :returns: ase cell 

107 :rtype: np.array 

108 

109 """ 

110 self.lammps_cell = lammps_cell 

111 self.lammps_tilt = self.lammps_cell.copy() 

112 

113 # reverse flip 

114 for iteri, (i, j, k) in enumerate(FLIP_ORDER): 

115 if self.flip[iteri]: 

116 change = self.lammps_cell[k][k] 

117 change *= np.sign(self.lammps_cell[i][j]) 

118 self.lammps_tilt[i][j] -= change 

119 

120 # try to detect potential flips in lammps 

121 # (lammps minimizes the cell-vector lengths) 

122 new_ase_cell = np.dot(self.lammps_tilt, self.rot_mat.T) 

123 # assuming the cell changes are mostly isotropic 

124 new_vol = np.linalg.det(new_ase_cell) 

125 old_vol = np.linalg.det(self.ase_cell) 

126 test_residual = self.ase_cell.copy() 

127 test_residual *= (new_vol / old_vol) ** (1.0 / 3.0) 

128 test_residual -= new_ase_cell 

129 if any( 

130 np.linalg.norm(test_residual, axis=1) 

131 > 0.5 * np.linalg.norm(self.ase_cell, axis=1) 

132 ): 

133 print( 

134 "WARNING: Significant simulation cell changes from LAMMPS " 

135 + "detected.\n" 

136 + " " * 9 

137 + "Backtransformation to ASE might fail!" 

138 ) 

139 return new_ase_cell 

140 

141 def vector_to_lammps(self, vec, wrap=False): 

142 """Rotate vector from ase coordinate system to lammps one 

143 

144 :param vec: to be rotated ase-vector 

145 :returns: lammps-vector 

146 :rtype: np.array 

147 

148 """ 

149 # !TODO: right eps-limit 

150 # lammps might not like atoms outside the cell 

151 if wrap: 

152 return wrap_positions( 

153 np.dot(vec, self.rot_mat), 

154 self.lammps_cell, 

155 pbc=self.pbc, 

156 eps=1e-18, 

157 ) 

158 

159 return np.dot(vec, self.rot_mat) 

160 

161 def vector_to_ase(self, vec, wrap=False): 

162 """Rotate vector from lammps coordinate system to ase one 

163 

164 :param vec: to be rotated lammps-vector 

165 :param wrap: was vector wrapped into 'lammps_cell' 

166 :returns: ase-vector 

167 :rtype: np.array 

168 

169 """ 

170 if wrap: 

171 # trying to reverse wraping across periodic boundaries in the 

172 # lammps coordination-system: 

173 # translate: expresses lammps-coordinate system in the rotate, but 

174 # without tilt removed system 

175 # fractional: vector in tilted system 

176 translate = np.linalg.solve( 

177 self.lammps_tilt.T, self.lammps_cell.T 

178 ).T 

179 fractional = np.linalg.solve(self.lammps_tilt.T, vec.T).T 

180 

181 # !TODO: make somehow nicer 

182 # !TODO: handle extreme tilted cells 

183 for ifrac in fractional: 

184 for zyx in reversed(range(3)): 

185 if ifrac[zyx] >= 1.0 and self.pbc[zyx]: 

186 ifrac -= translate[zyx] 

187 elif ifrac[zyx] < 0.0 and self.pbc[zyx]: 

188 ifrac += translate[zyx] 

189 vec = np.dot(fractional, self.lammps_tilt) 

190 

191 return np.dot(vec, self.rot_mat.T) 

192 

193 def is_skewed(self): 

194 """Test if a lammps cell is not tetragonal 

195 

196 :returns: bool 

197 :rtype: bool 

198 

199 """ 

200 cell_sq = self.lammps_cell ** 2 

201 return ( 

202 np.sum(np.tril(cell_sq, -1)) / np.sum(np.diag(cell_sq)) 

203 > self.tolerance 

204 )