 r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1"""Determine symmetry equivalence of two structures.

2Based on the recipe from Comput. Phys. Commun. 183, 690-697 (2012)."""

3from collections import Counter

4from itertools import combinations, product, filterfalse

6import numpy as np

7from scipy.spatial import cKDTree as KDTree

9from ase import Atom, Atoms

10from ase.build.tools import niggli_reduce

13def normalize(cell):

14 for i in range(3):

15 cell[i] /= np.linalg.norm(cell[i])

18class SpgLibNotFoundError(Exception):

21 def __init__(self, msg):

22 super(SpgLibNotFoundError, self).__init__(msg)

25class SymmetryEquivalenceCheck:

26 """Compare two structures to determine if they are symmetry equivalent.

28 Based on the recipe from Comput. Phys. Commun. 183, 690-697 (2012).

30 Parameters:

32 angle_tol: float

33 angle tolerance for the lattice vectors in degrees

35 ltol: float

36 relative tolerance for the length of the lattice vectors (per atom)

38 stol: float

39 position tolerance for the site comparison in units of

40 (V/N)^(1/3) (average length between atoms)

42 vol_tol: float

43 volume tolerance in angstrom cubed to compare the volumes of

44 the two structures

46 scale_volume: bool

47 if True the volumes of the two structures are scaled to be equal

49 to_primitive: bool

50 if True the structures are reduced to their primitive cells

51 note that this feature requires spglib to installed

53 Examples:

55 >>> from ase.build import bulk

56 >>> from ase.utils.structure_comparator import SymmetryEquivalenceCheck

57 >>> comp = SymmetryEquivalenceCheck()

59 Compare a cell with a rotated version

61 >>> a = bulk('Al', orthorhombic=True)

62 >>> b = a.copy()

63 >>> b.rotate(60, 'x', rotate_cell=True)

64 >>> comp.compare(a, b)

65 True

67 Transform to the primitive cell and then compare

69 >>> pa = bulk('Al')

70 >>> comp.compare(a, pa)

71 False

72 >>> comp = SymmetryEquivalenceCheck(to_primitive=True)

73 >>> comp.compare(a, pa)

74 True

76 Compare one structure with a list of other structures

78 >>> import numpy as np

79 >>> from ase import Atoms

80 >>> s1 = Atoms('H3', positions=[[0.5, 0.5, 0],

81 ... [0.5, 1.5, 0],

82 ... [1.5, 1.5, 0]],

83 ... cell=[2, 2, 2], pbc=True)

84 >>> comp = SymmetryEquivalenceCheck(stol=0.068)

85 >>> s2_list = []

86 >>> for d in np.linspace(0.1, 1.0, 5):

87 ... s2 = s1.copy()

88 ... s2.positions += [d, 0, 0]

89 ... s2_list.append(s2)

90 >>> comp.compare(s1, s2_list[:-1])

91 False

92 >>> comp.compare(s1, s2_list)

93 True

95 """

97 def __init__(self, angle_tol=1.0, ltol=0.05, stol=0.05, vol_tol=0.1,

98 scale_volume=False, to_primitive=False):

99 self.angle_tol = angle_tol * np.pi / 180.0 # convert to radians

100 self.scale_volume = scale_volume

101 self.stol = stol

102 self.ltol = ltol

103 self.vol_tol = vol_tol

104 self.position_tolerance = 0.0

105 self.to_primitive = to_primitive

107 # Variables to be used in the compare function

108 self.s1 = None

109 self.s2 = None

110 self.expanded_s1 = None

111 self.expanded_s2 = None

112 self.least_freq_element = None

114 def _niggli_reduce(self, atoms):

115 """Reduce to niggli cells.

117 Reduce the atoms to niggli cells, then rotates the niggli cells to

118 the so called "standard" orientation with one lattice vector along the

119 x-axis and a second vector in the xy plane.

120 """

121 niggli_reduce(atoms)

122 self._standarize_cell(atoms)

124 def _standarize_cell(self, atoms):

125 """Rotate the first vector such that it points along the x-axis.

126 Then rotate around the first vector so the second vector is in the

127 xy plane.

128 """

129 # Rotate first vector to x axis

130 cell = atoms.get_cell().T

131 total_rot_mat = np.eye(3)

132 v1 = cell[:, 0]

133 l1 = np.sqrt(v1**2 + v1**2)

134 angle = np.abs(np.arcsin(v1 / l1))

135 if (v1 < 0.0 and v1 > 0.0):

136 angle = np.pi - angle

137 elif (v1 < 0.0 and v1 < 0.0):

138 angle = np.pi + angle

139 elif (v1 > 0.0 and v1 < 0.0):

140 angle = -angle

141 ca = np.cos(angle)

142 sa = np.sin(angle)

143 rotmat = np.array([[ca, 0.0, sa], [0.0, 1.0, 0.0], [-sa, 0.0, ca]])

144 total_rot_mat = rotmat.dot(total_rot_mat)

145 cell = rotmat.dot(cell)

147 v1 = cell[:, 0]

148 l1 = np.sqrt(v1**2 + v1**2)

149 angle = np.abs(np.arcsin(v1 / l1))

150 if (v1 < 0.0 and v1 > 0.0):

151 angle = np.pi - angle

152 elif (v1 < 0.0 and v1 < 0.0):

153 angle = np.pi + angle

154 elif (v1 > 0.0 and v1 < 0.0):

155 angle = -angle

156 ca = np.cos(angle)

157 sa = np.sin(angle)

158 rotmat = np.array([[ca, sa, 0.0], [-sa, ca, 0.0], [0.0, 0.0, 1.0]])

159 total_rot_mat = rotmat.dot(total_rot_mat)

160 cell = rotmat.dot(cell)

162 # Rotate around x axis such that the second vector is in the xy plane

163 v2 = cell[:, 1]

164 l2 = np.sqrt(v2**2 + v2**2)

165 angle = np.abs(np.arcsin(v2 / l2))

166 if (v2 < 0.0 and v2 > 0.0):

167 angle = np.pi - angle

168 elif (v2 < 0.0 and v2 < 0.0):

169 angle = np.pi + angle

170 elif (v2 > 0.0 and v2 < 0.0):

171 angle = -angle

172 ca = np.cos(angle)

173 sa = np.sin(angle)

174 rotmat = np.array([[1.0, 0.0, 0.0], [0.0, ca, sa], [0.0, -sa, ca]])

175 total_rot_mat = rotmat.dot(total_rot_mat)

176 cell = rotmat.dot(cell)

178 atoms.set_cell(cell.T)

179 atoms.set_positions(total_rot_mat.dot(atoms.get_positions().T).T)

180 atoms.wrap(pbc=[1, 1, 1])

181 return atoms

183 def _get_element_count(self, struct):

184 """Count the number of elements in each of the structures."""

185 return Counter(struct.numbers)

187 def _get_angles(self, cell):

188 """Get the internal angles of the unit cell."""

189 cell = cell.copy()

191 normalize(cell)

193 dot = cell.dot(cell.T)

195 # Extract only the relevant dot products

196 dot = [dot[0, 1], dot[0, 2], dot[1, 2]]

198 # Return angles

199 return np.arccos(dot)

201 def _has_same_elements(self):

202 """Check if two structures have same elements."""

203 elem1 = self._get_element_count(self.s1)

204 return elem1 == self._get_element_count(self.s2)

206 def _has_same_angles(self):

207 """Check that the Niggli unit vectors has the same internal angles."""

208 ang1 = np.sort(self._get_angles(self.s1.get_cell()))

209 ang2 = np.sort(self._get_angles(self.s2.get_cell()))

211 return np.allclose(ang1, ang2, rtol=0, atol=self.angle_tol)

213 def _has_same_volume(self):

214 vol1 = self.s1.get_volume()

215 vol2 = self.s2.get_volume()

216 return np.abs(vol1 - vol2) < self.vol_tol

218 def _scale_volumes(self):

219 """Scale the cell of s2 to have the same volume as s1."""

220 cell2 = self.s2.get_cell()

221 # Get the volumes

222 v2 = np.linalg.det(cell2)

223 v1 = np.linalg.det(self.s1.get_cell())

225 # Scale the cells

226 coordinate_scaling = (v1 / v2)**(1.0 / 3.0)

227 cell2 *= coordinate_scaling

228 self.s2.set_cell(cell2, scale_atoms=True)

230 def compare(self, s1, s2):

231 """Compare the two structures.

233 Return *True* if the two structures are equivalent, *False* otherwise.

235 Parameters:

237 s1: Atoms object.

238 Transformation matrices are calculated based on this structure.

240 s2: Atoms or list

241 s1 can be compared to one structure or many structures supplied in

242 a list. If s2 is a list it returns True if any structure in s2

243 matches s1, False otherwise.

244 """

245 if self.to_primitive:

246 s1 = self._reduce_to_primitive(s1)

247 self._set_least_frequent_element(s1)

248 self._least_frequent_element_to_origin(s1)

249 self.s1 = s1.copy()

250 vol = self.s1.get_volume()

251 self.expanded_s1 = None

252 s1_niggli_reduced = False

254 if isinstance(s2, Atoms):

255 # Just make it a list of length 1

256 s2 = [s2]

258 matrices = None

259 translations = None

260 transposed_matrices = None

261 for struct in s2:

262 self.s2 = struct.copy()

263 self.expanded_s2 = None

265 if self.to_primitive:

266 self.s2 = self._reduce_to_primitive(self.s2)

268 # Compare number of elements in structures

269 if len(self.s1) != len(self.s2):

270 continue

272 # Compare chemical formulae

273 if not self._has_same_elements():

274 continue

276 # Compare angles

277 if not s1_niggli_reduced:

278 self._niggli_reduce(self.s1)

279 self._niggli_reduce(self.s2)

280 if not self._has_same_angles():

281 continue

283 # Compare volumes

284 if self.scale_volume:

285 self._scale_volumes()

286 if not self._has_same_volume():

287 continue

289 if matrices is None:

290 matrices = self._get_rotation_reflection_matrices()

291 if matrices is None:

292 continue

294 if translations is None:

295 translations = self._get_least_frequent_positions(self.s1)

297 # After the candidate translation based on s1 has been computed

298 # we need potentially to swap s1 and s2 for robust comparison

299 self._least_frequent_element_to_origin(self.s2)

300 switch = self._switch_reference_struct()

301 if switch:

302 # Remember the matrices and translations used before

303 old_matrices = matrices

304 old_translations = translations

306 # If a s1 and s2 has been switched we need to use the

307 # transposed version of the matrices to map atoms the

308 # other way

309 if transposed_matrices is None:

310 transposed_matrices = np.transpose(matrices,

311 axes=[0, 2, 1])

312 matrices = transposed_matrices

313 translations = self._get_least_frequent_positions(self.s1)

315 # Calculate tolerance on positions

316 self.position_tolerance = \

317 self.stol * (vol / len(self.s2))**(1.0 / 3.0)

319 if self._positions_match(matrices, translations):

320 return True

322 # Set the reference structure back to its original

323 self.s1 = s1.copy()

324 if switch:

325 self.expanded_s1 = self.expanded_s2

326 matrices = old_matrices

327 translations = old_translations

328 return False

330 def _set_least_frequent_element(self, atoms):

331 """Save the atomic number of the least frequent element."""

332 elem1 = self._get_element_count(atoms)

333 self.least_freq_element = elem1.most_common()[-1]

335 def _get_least_frequent_positions(self, atoms):

336 """Get the positions of the least frequent element in atoms."""

337 pos = atoms.get_positions(wrap=True)

338 return pos[atoms.numbers == self.least_freq_element]

340 def _get_only_least_frequent_of(self, struct):

341 """Get the atoms object with all other elements than the least frequent

342 one removed. Wrap the positions to get everything in the cell."""

343 pos = struct.get_positions(wrap=True)

345 indices = struct.numbers == self.least_freq_element

346 least_freq_struct = struct[indices]

347 least_freq_struct.set_positions(pos[indices])

349 return least_freq_struct

351 def _switch_reference_struct(self):

352 """There is an intrinsic assymetry in the system because

353 one of the atoms are being expanded, while the other is not.

354 This can cause the algorithm to return different result

355 depending on which structure is passed first.

356 We adopt the convention of using the atoms object

357 having the fewest atoms in its expanded cell as the

358 reference object.

359 We return True if a switch of structures has been performed."""

361 # First expand the cells

362 if self.expanded_s1 is None:

363 self.expanded_s1 = self._expand(self.s1)

364 if self.expanded_s2 is None:

365 self.expanded_s2 = self._expand(self.s2)

367 exp1 = self.expanded_s1

368 exp2 = self.expanded_s2

369 if len(exp1) < len(exp2):

370 # s1 should be the reference structure

371 # We have to swap s1 and s2

372 s1_temp = self.s1.copy()

373 self.s1 = self.s2

374 self.s2 = s1_temp

375 exp1_temp = self.expanded_s1.copy()

376 self.expanded_s1 = self.expanded_s2

377 self.expanded_s2 = exp1_temp

378 return True

379 return False

381 def _positions_match(self, rotation_reflection_matrices, translations):

382 """Check if the position and elements match.

384 Note that this function changes self.s1 and self.s2 to the rotation and

385 translation that matches best. Hence, it is crucial that this function

386 calls the element comparison, not the other way around.

387 """

388 pos1_ref = self.s1.get_positions(wrap=True)

390 # Get the expanded reference object

391 exp2 = self.expanded_s2

392 # Build a KD tree to enable fast look-up of nearest neighbours

393 tree = KDTree(exp2.get_positions())

394 for i in range(translations.shape):

395 # Translate

396 pos1_trans = pos1_ref - translations[i]

397 for matrix in rotation_reflection_matrices:

398 # Rotate

399 pos1 = matrix.dot(pos1_trans.T).T

401 # Update the atoms positions

402 self.s1.set_positions(pos1)

403 self.s1.wrap(pbc=[1, 1, 1])

404 if self._elements_match(self.s1, exp2, tree):

405 return True

406 return False

408 def _expand(self, ref_atoms, tol=0.0001):

409 """If an atom is closer to a boundary than tol it is repeated at the

410 opposite boundaries.

412 This ensures that atoms having crossed the cell boundaries due to

413 numerical noise are properly detected.

415 The distance between a position and cell boundary is calculated as:

416 dot(position, (b_vec x c_vec) / (|b_vec| |c_vec|) ), where x is the

417 cross product.

418 """

419 syms = ref_atoms.get_chemical_symbols()

420 cell = ref_atoms.get_cell()

421 positions = ref_atoms.get_positions(wrap=True)

422 expanded_atoms = ref_atoms.copy()

424 # Calculate normal vectors to the unit cell faces

425 normal_vectors = np.array([np.cross(cell[1, :], cell[2, :]),

426 np.cross(cell[0, :], cell[2, :]),

427 np.cross(cell[0, :], cell[1, :])])

428 normalize(normal_vectors)

430 # Get the distance to the unit cell faces from each atomic position

431 pos2faces = np.abs(positions.dot(normal_vectors.T))

433 # And the opposite faces

434 pos2oppofaces = np.abs(np.dot(positions - np.sum(cell, axis=0),

435 normal_vectors.T))

437 for i, i2face in enumerate(pos2faces):

438 # Append indices for positions close to the other faces

439 # and convert to boolean array signifying if the position at

440 # index i is close to the faces bordering origo (0, 1, 2) or

441 # the opposite faces (3, 4, 5)

442 i_close2face = np.append(i2face, pos2oppofaces[i]) < tol

443 # For each position i.e. row it holds that

444 # 1 x True -> close to face -> 1 extra atom at opposite face

445 # 2 x True -> close to edge -> 3 extra atoms at opposite edges

446 # 3 x True -> close to corner -> 7 extra atoms opposite corners

447 # E.g. to add atoms at all corners we need to use the cell

448 # vectors: (a, b, c, a + b, a + c, b + c, a + b + c), we use

449 # itertools.combinations to get them all

450 for j in range(sum(i_close2face)):

451 for c in combinations(np.nonzero(i_close2face), j + 1):

452 # Get the displacement vectors by adding the corresponding

453 # cell vectors, if the atom is close to an opposite face

454 # i.e. k > 2 subtract the cell vector

455 disp_vec = np.zeros(3)

456 for k in c:

457 disp_vec += cell[k % 3] * (int(k < 3) * 2 - 1)

458 pos = positions[i] + disp_vec

459 expanded_atoms.append(Atom(syms[i], position=pos))

460 return expanded_atoms

462 def _equal_elements_in_array(self, arr):

463 s = np.sort(arr)

464 return np.any(s[1:] == s[:-1])

466 def _elements_match(self, s1, s2, kdtree):

467 """Check if all the elements in s1 match corresponding position in s2

469 NOTE: The unit cells may be in different octants

470 Hence, try all cyclic permutations of x,y and z

471 """

472 pos1 = s1.get_positions()

473 for order in range(1): # Is the order still needed?

474 pos_order = [order, (order + 1) % 3, (order + 2) % 3]

475 pos = pos1[:, np.argsort(pos_order)]

476 dists, closest_in_s2 = kdtree.query(pos)

478 # Check if the elements are the same

479 if not np.all(s2.numbers[closest_in_s2] == s1.numbers):

480 return False

482 # Check if any distance is too large

483 if np.any(dists > self.position_tolerance):

484 return False

486 # Check for duplicates in what atom is closest

487 if self._equal_elements_in_array(closest_in_s2):

488 return False

490 return True

492 def _least_frequent_element_to_origin(self, atoms):

493 """Put one of the least frequent elements at the origin."""

494 least_freq_pos = self._get_least_frequent_positions(atoms)

495 cell_diag = np.sum(atoms.get_cell(), axis=0)

496 d = least_freq_pos - 1e-6 * cell_diag

497 atoms.positions -= d

498 atoms.wrap(pbc=[1, 1, 1])

500 def _get_rotation_reflection_matrices(self):

501 """Compute candidates for the transformation matrix."""

502 atoms1_ref = self._get_only_least_frequent_of(self.s1)

503 cell = self.s1.get_cell().T

504 cell_diag = np.sum(cell, axis=1)

505 angle_tol = self.angle_tol

508 # there always is an atom at the origin

509 delta_vec = 1E-6 * cell_diag

511 # Store three reference vectors and their lengths

512 ref_vec = self.s2.get_cell()

513 ref_vec_lengths = np.linalg.norm(ref_vec, axis=1)

515 # Compute ref vec angles

516 # ref_angles are arranged as [angle12, angle13, angle23]

517 ref_angles = np.array(self._get_angles(ref_vec))

518 large_angles = ref_angles > np.pi / 2.0

519 ref_angles[large_angles] = np.pi - ref_angles[large_angles]

521 # Translate by one cell diagonal so that a central cell is

522 # surrounded by cells in all directions

523 sc_atom_search = atoms1_ref * (3, 3, 3)

524 new_sc_pos = sc_atom_search.get_positions()

525 new_sc_pos -= new_sc_pos + cell_diag - delta_vec

527 lengths = np.linalg.norm(new_sc_pos, axis=1)

529 candidate_indices = []

530 rtol = self.ltol / len(self.s1)

531 for k in range(3):

533 ref_vec_lengths[k],

534 rtol=rtol, atol=0)

535 # The first vector is not interesting

538 # If no trial vectors can be found (for any direction)

539 # then the candidates are different and we return None

541 return None

545 # Now we calculate all relevant angles in one step. The relevant angles

546 # are the ones made by the current candidates. We will have to keep

547 # track of the indices in the angles matrix and the indices in the

548 # position and length arrays.

550 # Get all candidate indices (aci), only unique values

551 aci = np.sort(list(set().union(*candidate_indices)))

553 # Make a dictionary from original positions and lengths index to

554 # index in angle matrix

555 i2ang = dict(zip(aci, range(len(aci))))

557 # Calculate the dot product divided by the lengths:

558 # cos(angle) = dot(vec1, vec2) / |vec1| |vec2|

559 cosa = np.inner(new_sc_pos[aci],

560 new_sc_pos[aci]) / np.outer(lengths[aci],

561 lengths[aci])

562 # Make sure the inverse cosine will work

563 cosa[cosa > 1] = 1

564 cosa[cosa < -1] = -1

565 angles = np.arccos(cosa)

566 # Do trick for enantiomorphic structures

567 angles[angles > np.pi / 2] = np.pi - angles[angles > np.pi / 2]

569 # Check which angles match the reference angles

570 # Test for all combinations on candidates. filterfalse makes sure

571 # that there are no duplicate candidates. product is the same as

572 # nested for loops.

573 refined_candidate_list = []

574 for p in filterfalse(self._equal_elements_in_array,

575 product(*candidate_indices)):

576 a = np.array([angles[i2ang[p], i2ang[p]],

577 angles[i2ang[p], i2ang[p]],

578 angles[i2ang[p], i2ang[p]]])

580 if np.allclose(a, ref_angles, atol=angle_tol, rtol=0):

581 refined_candidate_list.append(new_sc_pos[np.array(p)].T)

583 # Get the rotation/reflection matrix [R] by:

584 # [R] = [V][T]^-1, where [V] is the reference vectors and

585 # [T] is the trial vectors

586 # XXX What do we know about the length/shape of refined_candidate_list?

587 if len(refined_candidate_list) == 0:

588 return None

589 else:

590 inverted_trial = np.linalg.inv(refined_candidate_list)

592 # Equivalent to np.matmul(ref_vec.T, inverted_trial)

593 candidate_trans_mat = np.dot(ref_vec.T, inverted_trial.T).T

594 return candidate_trans_mat

596 def _reduce_to_primitive(self, structure):

597 """Reduce the two structure to their primitive type"""

598 try:

599 import spglib

600 except ImportError:

601 raise SpgLibNotFoundError(

602 "SpgLib is required if to_primitive=True")

603 cell = (structure.get_cell()).tolist()

604 pos = structure.get_scaled_positions().tolist()

605 numbers = structure.get_atomic_numbers()

607 cell, scaled_pos, numbers = spglib.standardize_cell(

608 (cell, pos, numbers), to_primitive=True)

610 atoms = Atoms(

611 scaled_positions=scaled_pos,

612 numbers=numbers,

613 cell=cell,

614 pbc=True)

615 return atoms