Coverage for ase / spacegroup / symmetrize.py: 92.00%

1# fmt: off 

2 

3""" 

4Provides utility functions for FixSymmetry class 

5""" 

6from collections.abc import MutableMapping 

7 

8import numpy as np 

9 

10from ase.utils import atoms_to_spglib_cell, spglib_new_errorhandling 

11 

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

13 

14 

15def spglib_get_symmetry_dataset(*args, **kwargs): 

16 """Temporary compatibility adapter around spglib dataset. 

17 

18 Return an object that allows attribute-based access 

19 in line with recent spglib. This allows ASE code to not care about 

20 older spglib versions. 

21 """ 

22 import spglib 

23 

24 dataset = spglib_new_errorhandling(spglib.get_symmetry_dataset)( 

25 *args, **kwargs 

26 ) 

27 if dataset is None: 

28 return None 

29 if isinstance(dataset, dict): # spglib < 2.5.0 

30 return SpglibDatasetWrapper(dataset) 

31 return dataset # spglib >= 2.5.0 

32 

33 

34class SpglibDatasetWrapper(MutableMapping): 

35 # Spglib 2.5.0 returns SpglibDataset with deprecated __getitem__. 

36 # Spglib 2.4.0 and earlier return dict. 

37 # 

38 # We use this object to wrap dictionaries such that both types of access 

39 # work correctly. 

40 def __init__(self, spglib_dct): 

41 self._spglib_dct = spglib_dct 

42 

43 def __getattr__(self, attr): 

44 return self[attr] 

45 

46 def __getitem__(self, key): 

47 return self._spglib_dct[key] 

48 

49 def __len__(self): 

50 return len(self._spglib_dct) 

51 

52 def __iter__(self): 

53 return iter(self._spglib_dct) 

54 

55 def __setitem__(self, key, value): 

56 self._spglib_dct[key] = value 

57 

58 def __delitem__(self, item): 

59 del self._spglib_dct[item] 

60 

61 

62def print_symmetry(symprec, dataset): 

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

64 "got symmetry group number", dataset.number, 

65 ", international (Hermann-Mauguin)", dataset.international, 

66 ", Hall ", dataset.hall) 

67 

68 

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

70 """ 

71 Refine symmetry of an Atoms object 

72 

73 Parameters 

74 ---------- 

75 atoms - input Atoms object 

76 symprec - symmetry precicion 

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

78 

79 Returns 

80 ------- 

81 

82 spglib dataset 

83 

84 """ 

85 _check_and_symmetrize_cell(atoms, symprec=symprec, verbose=verbose) 

86 _check_and_symmetrize_positions(atoms, symprec=symprec, verbose=verbose) 

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

88 

89 

90class IntermediateDatasetError(Exception): 

91 """The symmetry dataset in `_check_and_symmetrize_positions` can be at odds 

92 with the original symmetry dataset in `_check_and_symmetrize_cell`. 

93 This implies a faulty partial symmetrization if not handled by exception.""" 

94 

95 

96def get_symmetrized_atoms(atoms, 

97 symprec: float = 0.01, 

98 final_symprec: float | None = None): 

99 """Get new Atoms object with refined symmetries. 

100 

101 Checks internal consistency of the found symmetries. 

102 

103 Parameters 

104 ---------- 

105 atoms : Atoms 

106 Input atoms object. 

107 symprec : float 

108 Symmetry precision used to identify symmetries with spglib. 

109 final_symprec : float 

110 Symmetry precision used for testing the symmetrization. 

111 

112 Returns 

113 ------- 

114 symatoms : Atoms 

115 New atoms object symmetrized according to the input symprec. 

116 """ 

117 atoms = atoms.copy() 

118 original_dataset = _check_and_symmetrize_cell(atoms, symprec=symprec) 

119 intermediate_dataset = _check_and_symmetrize_positions( 

120 atoms, symprec=symprec) 

121 if intermediate_dataset.number != original_dataset.number: 

122 raise IntermediateDatasetError() 

123 final_symprec = final_symprec or symprec 

124 final_dataset = check_symmetry(atoms, symprec=final_symprec) 

125 assert final_dataset.number == original_dataset.number 

126 return atoms, final_dataset 

127 

128 

129def _check_and_symmetrize_cell(atoms, **kwargs): 

130 dataset = check_symmetry(atoms, **kwargs) 

131 _symmetrize_cell(atoms, dataset) 

132 return dataset 

133 

134 

135def _symmetrize_cell(atoms, dataset): 

136 # set actual cell to symmetrized cell vectors by copying 

137 # transformed and rotated standard cell 

138 std_cell = dataset.std_lattice 

139 trans_std_cell = dataset.transformation_matrix.T @ std_cell 

140 rot_trans_std_cell = trans_std_cell @ dataset.std_rotation_matrix 

141 atoms.set_cell(rot_trans_std_cell, True) 

142 

143 

144def _check_and_symmetrize_positions(atoms, *, symprec, **kwargs): 

145 import spglib 

146 dataset = check_symmetry(atoms, symprec=symprec, **kwargs) 

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

148 # are compatible with std_lattice returned by get_symmetry_dataset 

149 

150 res = spglib_new_errorhandling(spglib.find_primitive)( 

151 atoms_to_spglib_cell(atoms), symprec=symprec) 

152 _symmetrize_positions(atoms, dataset, res) 

153 return dataset 

154 

155 

156def _symmetrize_positions(atoms, dataset, primitive_spglib_cell): 

157 prim_cell, _prim_scaled_pos, _prim_types = primitive_spglib_cell 

158 

159 # calculate offset between standard cell and actual cell 

160 std_cell = dataset.std_lattice 

161 rot_std_cell = std_cell @ dataset.std_rotation_matrix 

162 rot_std_pos = dataset.std_positions @ rot_std_cell 

163 pos = atoms.get_positions() 

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

165 list(dataset.std_mapping_to_primitive).index(0)]) 

166 

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

168 rot_prim_cell = prim_cell @ dataset.std_rotation_matrix 

169 inv_rot_prim_cell = np.linalg.inv(rot_prim_cell) 

170 aligned_std_pos = rot_std_pos + dp0 

171 

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

173 # integer_vec . primitive cell vectors 

174 mapping_to_primitive = list(dataset.mapping_to_primitive) 

175 std_mapping_to_primitive = list(dataset.std_mapping_to_primitive) 

176 pos = atoms.get_positions() 

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

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

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

180 dp_s = dp @ inv_rot_prim_cell 

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

182 atoms.set_positions(pos) 

183 

184 

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

186 """ 

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

188 

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

190 """ 

191 dataset = spglib_get_symmetry_dataset(atoms_to_spglib_cell(atoms), 

192 symprec=symprec) 

193 if verbose: 

194 print_symmetry(symprec, dataset) 

195 return dataset 

196 

197 

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

199 """ 

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

201 """ 

202 for rot1, trns1 in zip(sub_data.rotations, sub_data.translations): 

203 for rot2, trns2 in zip(sup_data.rotations, sup_data.translations): 

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

205 break 

206 else: 

207 return False 

208 return True 

209 

210 

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

212 """ 

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

214 

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

216 """ 

217 dataset = spglib_get_symmetry_dataset(atoms_to_spglib_cell(atoms), 

218 symprec=symprec) 

219 if verbose: 

220 print_symmetry(symprec, dataset) 

221 rotations = dataset.rotations.copy() 

222 translations = dataset.translations.copy() 

223 symm_map = [] 

224 scaled_pos = atoms.get_scaled_positions() 

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

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

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

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

229 dp = scaled_pos - new_p 

230 dp -= np.round(dp) 

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

232 this_op_map[i_at] = i_at_map 

233 symm_map.append(this_op_map) 

234 return (rotations, translations, symm_map) 

235 

236 

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

238 """ 

239 Return symmetrized forces 

240 

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

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

243 """ 

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

245 

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

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

248 transformed_forces_T = np.dot(r, scaled_forces_T) 

249 scaled_symmetrized_forces_T[:, this_op_map] += transformed_forces_T 

250 scaled_symmetrized_forces_T /= len(rot) 

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

252 

253 return symmetrized_forces 

254 

255 

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

257 """ 

258 Return symmetrized stress 

259 

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

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

262 """ 

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

264 

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

266 for r in rot: 

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

268 symmetrized_scaled_stress /= len(rot) 

269 

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

271 return sym