1""" 

2Provides FixSymmetry class to preserve spacegroup symmetry during optimisation 

3""" 

4import warnings 

5import numpy as np 

6 

7from ase.constraints import FixConstraint 

8from ase.stress import voigt_6_to_full_3x3_stress, full_3x3_to_voigt_6_stress 

9from ase.utils import atoms_to_spglib_cell 

10 

11 

12__all__ = ['refine_symmetry', 'check_symmetry', 'FixSymmetry'] 

13 

14 

15def print_symmetry(symprec, dataset): 

16 print("ase.spacegroup.symmetrize: prec", symprec, 

17 "got symmetry group number", dataset["number"], 

18 ", international (Hermann-Mauguin)", dataset["international"], 

19 ", Hall ", dataset["hall"]) 

20 

21 

22def refine_symmetry(atoms, symprec=0.01, verbose=False): 

23 """ 

24 Refine symmetry of an Atoms object 

25 

26 Parameters 

27 ---------- 

28 atoms - input Atoms object 

29 symprec - symmetry precicion 

30 verbose - if True, print out symmetry information before and after 

31 

32 Returns 

33 ------- 

34 

35 spglib dataset 

36 

37 """ 

38 import spglib 

39 

40 # test orig config with desired tol 

41 dataset = check_symmetry(atoms, symprec, verbose=verbose) 

42 

43 # set actual cell to symmetrized cell vectors by copying 

44 # transformed and rotated standard cell 

45 std_cell = dataset['std_lattice'] 

46 trans_std_cell = dataset['transformation_matrix'].T @ std_cell 

47 rot_trans_std_cell = trans_std_cell @ dataset['std_rotation_matrix'] 

48 atoms.set_cell(rot_trans_std_cell, True) 

49 

50 # get new dataset and primitive cell 

51 dataset = check_symmetry(atoms, symprec=symprec, verbose=verbose) 

52 res = spglib.find_primitive(atoms_to_spglib_cell(atoms), symprec=symprec) 

53 prim_cell, prim_scaled_pos, prim_types = res 

54 

55 # calculate offset between standard cell and actual cell 

56 std_cell = dataset['std_lattice'] 

57 rot_std_cell = std_cell @ dataset['std_rotation_matrix'] 

58 rot_std_pos = dataset['std_positions'] @ rot_std_cell 

59 pos = atoms.get_positions() 

60 dp0 = (pos[list(dataset['mapping_to_primitive']).index(0)] - rot_std_pos[ 

61 list(dataset['std_mapping_to_primitive']).index(0)]) 

62 

63 # create aligned set of standard cell positions to figure out mapping 

64 rot_prim_cell = prim_cell @ dataset['std_rotation_matrix'] 

65 inv_rot_prim_cell = np.linalg.inv(rot_prim_cell) 

66 aligned_std_pos = rot_std_pos + dp0 

67 

68 # find ideal positions from position of corresponding std cell atom + 

69 # integer_vec . primitive cell vectors 

70 # here we are assuming that primitive vectors returned by find_primitive 

71 # are compatible with std_lattice returned by get_symmetry_dataset 

72 mapping_to_primitive = list(dataset['mapping_to_primitive']) 

73 std_mapping_to_primitive = list(dataset['std_mapping_to_primitive']) 

74 pos = atoms.get_positions() 

75 for i_at in range(len(atoms)): 

76 std_i_at = std_mapping_to_primitive.index(mapping_to_primitive[i_at]) 

77 dp = aligned_std_pos[std_i_at] - pos[i_at] 

78 dp_s = dp @ inv_rot_prim_cell 

79 pos[i_at] = (aligned_std_pos[std_i_at] - np.round(dp_s) @ rot_prim_cell) 

80 atoms.set_positions(pos) 

81 

82 # test final config with tight tol 

83 return check_symmetry(atoms, symprec=1e-4, verbose=verbose) 

84 

85 

86def check_symmetry(atoms, symprec=1.0e-6, verbose=False): 

87 """ 

88 Check symmetry of `atoms` with precision `symprec` using `spglib` 

89 

90 Prints a summary and returns result of `spglib.get_symmetry_dataset()` 

91 """ 

92 import spglib 

93 dataset = spglib.get_symmetry_dataset(atoms_to_spglib_cell(atoms), 

94 symprec=symprec) 

95 if verbose: 

96 print_symmetry(symprec, dataset) 

97 return dataset 

98 

99 

100def is_subgroup(sup_data, sub_data, tol=1e-10): 

101 """ 

102 Test if spglib dataset `sub_data` is a subgroup of dataset `sup_data` 

103 """ 

104 for rot1, trns1 in zip(sub_data['rotations'], sub_data['translations']): 

105 for rot2, trns2 in zip(sup_data['rotations'], sup_data['translations']): 

106 if np.all(rot1 == rot2) and np.linalg.norm(trns1 - trns2) < tol: 

107 break 

108 else: 

109 return False 

110 return True 

111 

112 

113def prep_symmetry(atoms, symprec=1.0e-6, verbose=False): 

114 """ 

115 Prepare `at` for symmetry-preserving minimisation at precision `symprec` 

116 

117 Returns a tuple `(rotations, translations, symm_map)` 

118 """ 

119 import spglib 

120 

121 dataset = spglib.get_symmetry_dataset(atoms_to_spglib_cell(atoms), 

122 symprec=symprec) 

123 if verbose: 

124 print_symmetry(symprec, dataset) 

125 rotations = dataset['rotations'].copy() 

126 translations = dataset['translations'].copy() 

127 symm_map = [] 

128 scaled_pos = atoms.get_scaled_positions() 

129 for (rot, trans) in zip(rotations, translations): 

130 this_op_map = [-1] * len(atoms) 

131 for i_at in range(len(atoms)): 

132 new_p = rot @ scaled_pos[i_at, :] + trans 

133 dp = scaled_pos - new_p 

134 dp -= np.round(dp) 

135 i_at_map = np.argmin(np.linalg.norm(dp, axis=1)) 

136 this_op_map[i_at] = i_at_map 

137 symm_map.append(this_op_map) 

138 return (rotations, translations, symm_map) 

139 

140 

141def symmetrize_rank1(lattice, inv_lattice, forces, rot, trans, symm_map): 

142 """ 

143 Return symmetrized forces 

144 

145 lattice vectors expected as row vectors (same as ASE get_cell() convention), 

146 inv_lattice is its matrix inverse (reciprocal().T) 

147 """ 

148 scaled_symmetrized_forces_T = np.zeros(forces.T.shape) 

149 

150 scaled_forces_T = np.dot(inv_lattice.T, forces.T) 

151 for (r, t, this_op_map) in zip(rot, trans, symm_map): 

152 transformed_forces_T = np.dot(r, scaled_forces_T) 

153 scaled_symmetrized_forces_T[:, this_op_map] += transformed_forces_T 

154 scaled_symmetrized_forces_T /= len(rot) 

155 symmetrized_forces = (lattice.T @ scaled_symmetrized_forces_T).T 

156 

157 return symmetrized_forces 

158 

159 

160def symmetrize_rank2(lattice, lattice_inv, stress_3_3, rot): 

161 """ 

162 Return symmetrized stress 

163 

164 lattice vectors expected as row vectors (same as ASE get_cell() convention), 

165 inv_lattice is its matrix inverse (reciprocal().T) 

166 """ 

167 scaled_stress = np.dot(np.dot(lattice, stress_3_3), lattice.T) 

168 

169 symmetrized_scaled_stress = np.zeros((3, 3)) 

170 for r in rot: 

171 symmetrized_scaled_stress += np.dot(np.dot(r.T, scaled_stress), r) 

172 symmetrized_scaled_stress /= len(rot) 

173 

174 sym = np.dot(np.dot(lattice_inv, symmetrized_scaled_stress), lattice_inv.T) 

175 return sym 

176 

177 

178class FixSymmetry(FixConstraint): 

179 """ 

180 Constraint to preserve spacegroup symmetry during optimisation. 

181 

182 Requires spglib package to be available. 

183 """ 

184 

185 def __init__(self, atoms, symprec=0.01, adjust_positions=True, 

186 adjust_cell=True, verbose=False): 

187 self.verbose = verbose 

188 refine_symmetry(atoms, symprec, self.verbose) # refine initial symmetry 

189 sym = prep_symmetry(atoms, symprec, self.verbose) 

190 self.rotations, self.translations, self.symm_map = sym 

191 self.do_adjust_positions = adjust_positions 

192 self.do_adjust_cell = adjust_cell 

193 

194 def adjust_cell(self, atoms, cell): 

195 if not self.do_adjust_cell: 

196 return 

197 # stress should definitely be symmetrized as a rank 2 tensor 

198 # UnitCellFilter uses deformation gradient as cell DOF with steps 

199 # dF = stress.F^-T quantity that should be symmetrized is therefore dF . 

200 # F^T assume prev F = I, so just symmetrize dF 

201 cur_cell = atoms.get_cell() 

202 cur_cell_inv = atoms.cell.reciprocal().T 

203 

204 # F defined such that cell = cur_cell . F^T 

205 # assume prev F = I, so dF = F - I 

206 delta_deform_grad = np.dot(cur_cell_inv, cell).T - np.eye(3) 

207 

208 # symmetrization doesn't work properly with large steps, since 

209 # it depends on current cell, and cell is being changed by deformation 

210 # gradient 

211 max_delta_deform_grad = np.max(np.abs(delta_deform_grad)) 

212 if max_delta_deform_grad > 0.25: 

213 raise RuntimeError('FixSymmetry adjust_cell does not work properly' 

214 ' with large deformation gradient step {} > 0.25' 

215 .format(max_delta_deform_grad)) 

216 elif max_delta_deform_grad > 0.15: 

217 warnings.warn('FixSymmetry adjust_cell may be ill behaved with' 

218 ' large deformation gradient step {}' 

219 .format(max_delta_deform_grad)) 

220 

221 symmetrized_delta_deform_grad = symmetrize_rank2(cur_cell, cur_cell_inv, 

222 delta_deform_grad, 

223 self.rotations) 

224 cell[:] = np.dot(cur_cell, 

225 (symmetrized_delta_deform_grad + np.eye(3)).T) 

226 

227 def adjust_positions(self, atoms, new): 

228 if not self.do_adjust_positions: 

229 return 

230 # symmetrize changes in position as rank 1 tensors 

231 step = new - atoms.positions 

232 symmetrized_step = symmetrize_rank1(atoms.get_cell(), 

233 atoms.cell.reciprocal().T, step, 

234 self.rotations, self.translations, 

235 self.symm_map) 

236 new[:] = atoms.positions + symmetrized_step 

237 

238 def adjust_forces(self, atoms, forces): 

239 # symmetrize forces as rank 1 tensors 

240 # print('adjusting forces') 

241 forces[:] = symmetrize_rank1(atoms.get_cell(), 

242 atoms.cell.reciprocal().T, forces, 

243 self.rotations, self.translations, 

244 self.symm_map) 

245 

246 def adjust_stress(self, atoms, stress): 

247 # symmetrize stress as rank 2 tensor 

248 raw_stress = voigt_6_to_full_3x3_stress(stress) 

249 symmetrized_stress = symmetrize_rank2(atoms.get_cell(), 

250 atoms.cell.reciprocal().T, 

251 raw_stress, self.rotations) 

252 stress[:] = full_3x3_to_voigt_6_stress(symmetrized_stress) 

253 

254 def index_shuffle(self, atoms, ind): 

255 if len(atoms) != len(ind) or len(set(ind)) != len(ind): 

256 raise RuntimeError("FixSymmetry can only accomodate atom" 

257 " permutions, and len(Atoms) == {} " 

258 "!= len(ind) == {} or ind has duplicates" 

259 .format(len(atoms), len(ind))) 

260 

261 ind_reversed = np.zeros((len(ind)), dtype=int) 

262 ind_reversed[ind] = range(len(ind)) 

263 new_symm_map = [] 

264 for sm in self.symm_map: 

265 new_sm = np.array([-1] * len(atoms)) 

266 for at_i in range(len(ind)): 

267 new_sm[ind_reversed[at_i]] = ind_reversed[sm[at_i]] 

268 new_symm_map.append(new_sm) 

269 

270 self.symm_map = new_symm_map