Coverage for /builds/kinetik161/ase/ase/spacegroup/symmetrize.py: 95.38%
130 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +0000
1"""
2Provides FixSymmetry class to preserve spacegroup symmetry during optimisation
3"""
4import warnings
6import numpy as np
8from ase.constraints import FixConstraint
9from ase.stress import full_3x3_to_voigt_6_stress, voigt_6_to_full_3x3_stress
10from ase.utils import atoms_to_spglib_cell
12__all__ = ['refine_symmetry', 'check_symmetry', 'FixSymmetry']
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"])
22def refine_symmetry(atoms, symprec=0.01, verbose=False):
23 """
24 Refine symmetry of an Atoms object
26 Parameters
27 ----------
28 atoms - input Atoms object
29 symprec - symmetry precicion
30 verbose - if True, print out symmetry information before and after
32 Returns
33 -------
35 spglib dataset
37 """
38 import spglib
40 # test orig config with desired tol
41 dataset = check_symmetry(atoms, symprec, verbose=verbose)
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)
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
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)])
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
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)
82 # test final config with tight tol
83 return check_symmetry(atoms, symprec=1e-4, verbose=verbose)
86def check_symmetry(atoms, symprec=1.0e-6, verbose=False):
87 """
88 Check symmetry of `atoms` with precision `symprec` using `spglib`
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
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
113def prep_symmetry(atoms, symprec=1.0e-6, verbose=False):
114 """
115 Prepare `at` for symmetry-preserving minimisation at precision `symprec`
117 Returns a tuple `(rotations, translations, symm_map)`
118 """
119 import spglib
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)
141def symmetrize_rank1(lattice, inv_lattice, forces, rot, trans, symm_map):
142 """
143 Return symmetrized forces
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)
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
157 return symmetrized_forces
160def symmetrize_rank2(lattice, lattice_inv, stress_3_3, rot):
161 """
162 Return symmetrized stress
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)
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)
174 sym = np.dot(np.dot(lattice_inv, symmetrized_scaled_stress), lattice_inv.T)
175 return sym
178class FixSymmetry(FixConstraint):
179 """
180 Constraint to preserve spacegroup symmetry during optimisation.
182 Requires spglib package to be available.
183 """
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
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
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)
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))
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)
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
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)
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)
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)))
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)
270 self.symm_map = new_symm_map