Coverage for /builds/ase/ase/ase/spacegroup/symmetrize.py: 92.06%

1# fmt: off 

2 

3""" 

4Provides utility functions for FixSymmetry class 

5""" 

6from collections.abc import MutableMapping 

7from typing import Optional 

8 

9import numpy as np 

10 

11from ase.utils import atoms_to_spglib_cell 

12 

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

14 

15 

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

17 """Temporary compatibility adapter around spglib dataset. 

18 

19 Return an object that allows attribute-based access 

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

21 older spglib versions. 

22 """ 

23 import spglib 

24 dataset = spglib.get_symmetry_dataset(*args, **kwargs) 

25 if dataset is None: 

26 return None 

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

28 return SpglibDatasetWrapper(dataset) 

29 return dataset # spglib >= 2.5.0 

30 

31 

32class SpglibDatasetWrapper(MutableMapping): 

33 # Spglib 2.5.0 returns SpglibDataset with deprecated __getitem__. 

34 # Spglib 2.4.0 and earlier return dict. 

35 # 

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

37 # work correctly. 

38 def __init__(self, spglib_dct): 

39 self._spglib_dct = spglib_dct 

40 

41 def __getattr__(self, attr): 

42 return self[attr] 

43 

44 def __getitem__(self, key): 

45 return self._spglib_dct[key] 

46 

47 def __len__(self): 

48 return len(self._spglib_dct) 

49 

50 def __iter__(self): 

51 return iter(self._spglib_dct) 

52 

53 def __setitem__(self, key, value): 

54 self._spglib_dct[key] = value 

55 

56 def __delitem__(self, item): 

57 del self._spglib_dct[item] 

58 

59 

60def print_symmetry(symprec, dataset): 

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

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

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

64 ", Hall ", dataset.hall) 

65 

66 

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

68 """ 

69 Refine symmetry of an Atoms object 

70 

71 Parameters 

72 ---------- 

73 atoms - input Atoms object 

74 symprec - symmetry precicion 

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

76 

77 Returns 

78 ------- 

79 

80 spglib dataset 

81 

82 """ 

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

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

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

86 

87 

88class IntermediateDatasetError(Exception): 

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

90 with the original symmetry dataset in `_check_and_symmetrize_cell`. 

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

92 

93 

94def get_symmetrized_atoms(atoms, 

95 symprec: float = 0.01, 

96 final_symprec: Optional[float] = None): 

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

98 

99 Checks internal consistency of the found symmetries. 

100 

101 Parameters 

102 ---------- 

103 atoms : Atoms 

104 Input atoms object. 

105 symprec : float 

106 Symmetry precision used to identify symmetries with spglib. 

107 final_symprec : float 

108 Symmetry precision used for testing the symmetrization. 

109 

110 Returns 

111 ------- 

112 symatoms : Atoms 

113 New atoms object symmetrized according to the input symprec. 

114 """ 

115 atoms = atoms.copy() 

116 original_dataset = _check_and_symmetrize_cell(atoms, symprec=symprec) 

117 intermediate_dataset = _check_and_symmetrize_positions( 

118 atoms, symprec=symprec) 

119 if intermediate_dataset.number != original_dataset.number: 

120 raise IntermediateDatasetError() 

121 final_symprec = final_symprec or symprec 

122 final_dataset = check_symmetry(atoms, symprec=final_symprec) 

123 assert final_dataset.number == original_dataset.number 

124 return atoms, final_dataset 

125 

126 

127def _check_and_symmetrize_cell(atoms, **kwargs): 

128 dataset = check_symmetry(atoms, **kwargs) 

129 _symmetrize_cell(atoms, dataset) 

130 return dataset 

131 

132 

133def _symmetrize_cell(atoms, dataset): 

134 # set actual cell to symmetrized cell vectors by copying 

135 # transformed and rotated standard cell 

136 std_cell = dataset.std_lattice 

137 trans_std_cell = dataset.transformation_matrix.T @ std_cell 

138 rot_trans_std_cell = trans_std_cell @ dataset.std_rotation_matrix 

139 atoms.set_cell(rot_trans_std_cell, True) 

140 

141 

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

143 import spglib 

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

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

146 # are compatible with std_lattice returned by get_symmetry_dataset 

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

148 _symmetrize_positions(atoms, dataset, res) 

149 return dataset 

150 

151 

152def _symmetrize_positions(atoms, dataset, primitive_spglib_cell): 

153 prim_cell, _prim_scaled_pos, _prim_types = primitive_spglib_cell 

154 

155 # calculate offset between standard cell and actual cell 

156 std_cell = dataset.std_lattice 

157 rot_std_cell = std_cell @ dataset.std_rotation_matrix 

158 rot_std_pos = dataset.std_positions @ rot_std_cell 

159 pos = atoms.get_positions() 

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

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

162 

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

164 rot_prim_cell = prim_cell @ dataset.std_rotation_matrix 

165 inv_rot_prim_cell = np.linalg.inv(rot_prim_cell) 

166 aligned_std_pos = rot_std_pos + dp0 

167 

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

169 # integer_vec . primitive cell vectors 

170 mapping_to_primitive = list(dataset.mapping_to_primitive) 

171 std_mapping_to_primitive = list(dataset.std_mapping_to_primitive) 

172 pos = atoms.get_positions() 

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

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

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

176 dp_s = dp @ inv_rot_prim_cell 

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

178 atoms.set_positions(pos) 

179 

180 

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

182 """ 

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

184 

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

186 """ 

187 dataset = spglib_get_symmetry_dataset(atoms_to_spglib_cell(atoms), 

188 symprec=symprec) 

189 if verbose: 

190 print_symmetry(symprec, dataset) 

191 return dataset 

192 

193 

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

195 """ 

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

197 """ 

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

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

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

201 break 

202 else: 

203 return False 

204 return True 

205 

206 

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

208 """ 

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

210 

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

212 """ 

213 dataset = spglib_get_symmetry_dataset(atoms_to_spglib_cell(atoms), 

214 symprec=symprec) 

215 if verbose: 

216 print_symmetry(symprec, dataset) 

217 rotations = dataset.rotations.copy() 

218 translations = dataset.translations.copy() 

219 symm_map = [] 

220 scaled_pos = atoms.get_scaled_positions() 

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

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

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

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

225 dp = scaled_pos - new_p 

226 dp -= np.round(dp) 

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

228 this_op_map[i_at] = i_at_map 

229 symm_map.append(this_op_map) 

230 return (rotations, translations, symm_map) 

231 

232 

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

234 """ 

235 Return symmetrized forces 

236 

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

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

239 """ 

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

241 

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

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

244 transformed_forces_T = np.dot(r, scaled_forces_T) 

245 scaled_symmetrized_forces_T[:, this_op_map] += transformed_forces_T 

246 scaled_symmetrized_forces_T /= len(rot) 

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

248 

249 return symmetrized_forces 

250 

251 

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

253 """ 

254 Return symmetrized stress 

255 

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

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

258 """ 

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

260 

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

262 for r in rot: 

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

264 symmetrized_scaled_stress /= len(rot) 

265 

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

267 return sym