Coverage for /builds/alexhroom/ase/ase/spacegroup/symmetrize.py: 99.04%
104 statements
« prev ^ index » next coverage.py v7.5.3, created at 2024-08-05 14:37 +0000
1"""
2Provides utility functions for FixSymmetry class
3"""
4from typing import Optional
6import numpy as np
8from ase.utils import atoms_to_spglib_cell
10__all__ = ['refine_symmetry', 'check_symmetry']
13def print_symmetry(symprec, dataset):
14 print("ase.spacegroup.symmetrize: prec", symprec,
15 "got symmetry group number", dataset["number"],
16 ", international (Hermann-Mauguin)", dataset["international"],
17 ", Hall ", dataset["hall"])
20def refine_symmetry(atoms, symprec=0.01, verbose=False):
21 """
22 Refine symmetry of an Atoms object
24 Parameters
25 ----------
26 atoms - input Atoms object
27 symprec - symmetry precicion
28 verbose - if True, print out symmetry information before and after
30 Returns
31 -------
33 spglib dataset
35 """
36 _check_and_symmetrize_cell(atoms, symprec=symprec, verbose=verbose)
37 _check_and_symmetrize_positions(atoms, symprec=symprec, verbose=verbose)
38 return check_symmetry(atoms, symprec=1e-4, verbose=verbose)
41class IntermediateDatasetError(Exception):
42 """The symmetry dataset in `_check_and_symmetrize_positions` can be at odds
43 with the original symmetry dataset in `_check_and_symmetrize_cell`.
44 This implies a faulty partial symmetrization if not handled by exception."""
47def get_symmetrized_atoms(atoms,
48 symprec: float = 0.01,
49 final_symprec: Optional[float] = None):
50 """Get new Atoms object with refined symmetries.
52 Checks internal consistency of the found symmetries.
54 Parameters
55 ----------
56 atoms : Atoms
57 Input atoms object.
58 symprec : float
59 Symmetry precision used to identify symmetries with spglib.
60 final_symprec : float
61 Symmetry precision used for testing the symmetrization.
63 Returns
64 -------
65 symatoms : Atoms
66 New atoms object symmetrized according to the input symprec.
67 """
68 atoms = atoms.copy()
69 original_dataset = _check_and_symmetrize_cell(atoms, symprec=symprec)
70 intermediate_dataset = _check_and_symmetrize_positions(
71 atoms, symprec=symprec)
72 if intermediate_dataset['number'] != original_dataset['number']:
73 raise IntermediateDatasetError()
74 final_symprec = final_symprec or symprec
75 final_dataset = check_symmetry(atoms, symprec=final_symprec)
76 assert final_dataset['number'] == original_dataset['number']
77 return atoms, final_dataset
80def _check_and_symmetrize_cell(atoms, **kwargs):
81 dataset = check_symmetry(atoms, **kwargs)
82 _symmetrize_cell(atoms, dataset)
83 return dataset
86def _symmetrize_cell(atoms, dataset):
87 # set actual cell to symmetrized cell vectors by copying
88 # transformed and rotated standard cell
89 std_cell = dataset['std_lattice']
90 trans_std_cell = dataset['transformation_matrix'].T @ std_cell
91 rot_trans_std_cell = trans_std_cell @ dataset['std_rotation_matrix']
92 atoms.set_cell(rot_trans_std_cell, True)
95def _check_and_symmetrize_positions(atoms, *, symprec, **kwargs):
96 import spglib
97 dataset = check_symmetry(atoms, symprec=symprec, **kwargs)
98 # here we are assuming that primitive vectors returned by find_primitive
99 # are compatible with std_lattice returned by get_symmetry_dataset
100 res = spglib.find_primitive(atoms_to_spglib_cell(atoms), symprec=symprec)
101 _symmetrize_positions(atoms, dataset, res)
102 return dataset
105def _symmetrize_positions(atoms, dataset, primitive_spglib_cell):
106 prim_cell, prim_scaled_pos, prim_types = primitive_spglib_cell
108 # calculate offset between standard cell and actual cell
109 std_cell = dataset['std_lattice']
110 rot_std_cell = std_cell @ dataset['std_rotation_matrix']
111 rot_std_pos = dataset['std_positions'] @ rot_std_cell
112 pos = atoms.get_positions()
113 dp0 = (pos[list(dataset['mapping_to_primitive']).index(0)] - rot_std_pos[
114 list(dataset['std_mapping_to_primitive']).index(0)])
116 # create aligned set of standard cell positions to figure out mapping
117 rot_prim_cell = prim_cell @ dataset['std_rotation_matrix']
118 inv_rot_prim_cell = np.linalg.inv(rot_prim_cell)
119 aligned_std_pos = rot_std_pos + dp0
121 # find ideal positions from position of corresponding std cell atom +
122 # integer_vec . primitive cell vectors
123 mapping_to_primitive = list(dataset['mapping_to_primitive'])
124 std_mapping_to_primitive = list(dataset['std_mapping_to_primitive'])
125 pos = atoms.get_positions()
126 for i_at in range(len(atoms)):
127 std_i_at = std_mapping_to_primitive.index(mapping_to_primitive[i_at])
128 dp = aligned_std_pos[std_i_at] - pos[i_at]
129 dp_s = dp @ inv_rot_prim_cell
130 pos[i_at] = (aligned_std_pos[std_i_at] - np.round(dp_s) @ rot_prim_cell)
131 atoms.set_positions(pos)
134def check_symmetry(atoms, symprec=1.0e-6, verbose=False):
135 """
136 Check symmetry of `atoms` with precision `symprec` using `spglib`
138 Prints a summary and returns result of `spglib.get_symmetry_dataset()`
139 """
140 import spglib
141 dataset = spglib.get_symmetry_dataset(atoms_to_spglib_cell(atoms),
142 symprec=symprec)
143 if verbose:
144 print_symmetry(symprec, dataset)
145 return dataset
148def is_subgroup(sup_data, sub_data, tol=1e-10):
149 """
150 Test if spglib dataset `sub_data` is a subgroup of dataset `sup_data`
151 """
152 for rot1, trns1 in zip(sub_data['rotations'], sub_data['translations']):
153 for rot2, trns2 in zip(sup_data['rotations'], sup_data['translations']):
154 if np.all(rot1 == rot2) and np.linalg.norm(trns1 - trns2) < tol:
155 break
156 else:
157 return False
158 return True
161def prep_symmetry(atoms, symprec=1.0e-6, verbose=False):
162 """
163 Prepare `at` for symmetry-preserving minimisation at precision `symprec`
165 Returns a tuple `(rotations, translations, symm_map)`
166 """
167 import spglib
169 dataset = spglib.get_symmetry_dataset(atoms_to_spglib_cell(atoms),
170 symprec=symprec)
171 if verbose:
172 print_symmetry(symprec, dataset)
173 rotations = dataset['rotations'].copy()
174 translations = dataset['translations'].copy()
175 symm_map = []
176 scaled_pos = atoms.get_scaled_positions()
177 for (rot, trans) in zip(rotations, translations):
178 this_op_map = [-1] * len(atoms)
179 for i_at in range(len(atoms)):
180 new_p = rot @ scaled_pos[i_at, :] + trans
181 dp = scaled_pos - new_p
182 dp -= np.round(dp)
183 i_at_map = np.argmin(np.linalg.norm(dp, axis=1))
184 this_op_map[i_at] = i_at_map
185 symm_map.append(this_op_map)
186 return (rotations, translations, symm_map)
189def symmetrize_rank1(lattice, inv_lattice, forces, rot, trans, symm_map):
190 """
191 Return symmetrized forces
193 lattice vectors expected as row vectors (same as ASE get_cell() convention),
194 inv_lattice is its matrix inverse (reciprocal().T)
195 """
196 scaled_symmetrized_forces_T = np.zeros(forces.T.shape)
198 scaled_forces_T = np.dot(inv_lattice.T, forces.T)
199 for (r, t, this_op_map) in zip(rot, trans, symm_map):
200 transformed_forces_T = np.dot(r, scaled_forces_T)
201 scaled_symmetrized_forces_T[:, this_op_map] += transformed_forces_T
202 scaled_symmetrized_forces_T /= len(rot)
203 symmetrized_forces = (lattice.T @ scaled_symmetrized_forces_T).T
205 return symmetrized_forces
208def symmetrize_rank2(lattice, lattice_inv, stress_3_3, rot):
209 """
210 Return symmetrized stress
212 lattice vectors expected as row vectors (same as ASE get_cell() convention),
213 inv_lattice is its matrix inverse (reciprocal().T)
214 """
215 scaled_stress = np.dot(np.dot(lattice, stress_3_3), lattice.T)
217 symmetrized_scaled_stress = np.zeros((3, 3))
218 for r in rot:
219 symmetrized_scaled_stress += np.dot(np.dot(r.T, scaled_stress), r)
220 symmetrized_scaled_stress /= len(rot)
222 sym = np.dot(np.dot(lattice_inv, symmetrized_scaled_stress), lattice_inv.T)
223 return sym