Source code for ase.spacegroup.symmetrize

Provides FixSymmetry class to preserve spacegroup symmetry during optimisation
import warnings
import numpy as np

from ase.constraints import FixConstraint
from ase.stress import voigt_6_to_full_3x3_stress, full_3x3_to_voigt_6_stress
from ase.utils import atoms_to_spglib_cell

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

def print_symmetry(symprec, dataset):
    print("ase.spacegroup.symmetrize: prec", symprec,
          "got symmetry group number", dataset["number"],
          ", international (Hermann-Mauguin)", dataset["international"],
          ", Hall ", dataset["hall"])

[docs]def refine_symmetry(atoms, symprec=0.01, verbose=False): """ Refine symmetry of an Atoms object Parameters ---------- atoms - input Atoms object symprec - symmetry precicion verbose - if True, print out symmetry information before and after Returns ------- spglib dataset """ import spglib # test orig config with desired tol dataset = check_symmetry(atoms, symprec, verbose=verbose) # set actual cell to symmetrized cell vectors by copying # transformed and rotated standard cell std_cell = dataset['std_lattice'] trans_std_cell = dataset['transformation_matrix'].T @ std_cell rot_trans_std_cell = trans_std_cell @ dataset['std_rotation_matrix'] atoms.set_cell(rot_trans_std_cell, True) # get new dataset and primitive cell dataset = check_symmetry(atoms, symprec=symprec, verbose=verbose) res = spglib.find_primitive(atoms_to_spglib_cell(atoms), symprec=symprec) prim_cell, prim_scaled_pos, prim_types = res # calculate offset between standard cell and actual cell std_cell = dataset['std_lattice'] rot_std_cell = std_cell @ dataset['std_rotation_matrix'] rot_std_pos = dataset['std_positions'] @ rot_std_cell pos = atoms.get_positions() dp0 = (pos[list(dataset['mapping_to_primitive']).index(0)] - rot_std_pos[ list(dataset['std_mapping_to_primitive']).index(0)]) # create aligned set of standard cell positions to figure out mapping rot_prim_cell = prim_cell @ dataset['std_rotation_matrix'] inv_rot_prim_cell = np.linalg.inv(rot_prim_cell) aligned_std_pos = rot_std_pos + dp0 # find ideal positions from position of corresponding std cell atom + # integer_vec . primitive cell vectors # here we are assuming that primitive vectors returned by find_primitive # are compatible with std_lattice returned by get_symmetry_dataset mapping_to_primitive = list(dataset['mapping_to_primitive']) std_mapping_to_primitive = list(dataset['std_mapping_to_primitive']) pos = atoms.get_positions() for i_at in range(len(atoms)): std_i_at = std_mapping_to_primitive.index(mapping_to_primitive[i_at]) dp = aligned_std_pos[std_i_at] - pos[i_at] dp_s = dp @ inv_rot_prim_cell pos[i_at] = (aligned_std_pos[std_i_at] - np.round(dp_s) @ rot_prim_cell) atoms.set_positions(pos) # test final config with tight tol return check_symmetry(atoms, symprec=1e-4, verbose=verbose)
[docs]def check_symmetry(atoms, symprec=1.0e-6, verbose=False): """ Check symmetry of `atoms` with precision `symprec` using `spglib` Prints a summary and returns result of `spglib.get_symmetry_dataset()` """ import spglib dataset = spglib.get_symmetry_dataset(atoms_to_spglib_cell(atoms), symprec=symprec) if verbose: print_symmetry(symprec, dataset) return dataset
def is_subgroup(sup_data, sub_data, tol=1e-10): """ Test if spglib dataset `sub_data` is a subgroup of dataset `sup_data` """ for rot1, trns1 in zip(sub_data['rotations'], sub_data['translations']): for rot2, trns2 in zip(sup_data['rotations'], sup_data['translations']): if np.all(rot1 == rot2) and np.linalg.norm(trns1 - trns2) < tol: break else: return False return True def prep_symmetry(atoms, symprec=1.0e-6, verbose=False): """ Prepare `at` for symmetry-preserving minimisation at precision `symprec` Returns a tuple `(rotations, translations, symm_map)` """ import spglib dataset = spglib.get_symmetry_dataset(atoms_to_spglib_cell(atoms), symprec=symprec) if verbose: print_symmetry(symprec, dataset) rotations = dataset['rotations'].copy() translations = dataset['translations'].copy() symm_map = [] scaled_pos = atoms.get_scaled_positions() for (rot, trans) in zip(rotations, translations): this_op_map = [-1] * len(atoms) for i_at in range(len(atoms)): new_p = rot @ scaled_pos[i_at, :] + trans dp = scaled_pos - new_p dp -= np.round(dp) i_at_map = np.argmin(np.linalg.norm(dp, axis=1)) this_op_map[i_at] = i_at_map symm_map.append(this_op_map) return (rotations, translations, symm_map) def symmetrize_rank1(lattice, inv_lattice, forces, rot, trans, symm_map): """ Return symmetrized forces lattice vectors expected as row vectors (same as ASE get_cell() convention), inv_lattice is its matrix inverse (reciprocal().T) """ scaled_symmetrized_forces_T = np.zeros(forces.T.shape) scaled_forces_T =, forces.T) for (r, t, this_op_map) in zip(rot, trans, symm_map): transformed_forces_T =, scaled_forces_T) scaled_symmetrized_forces_T[:, this_op_map] += transformed_forces_T scaled_symmetrized_forces_T /= len(rot) symmetrized_forces = (lattice.T @ scaled_symmetrized_forces_T).T return symmetrized_forces def symmetrize_rank2(lattice, lattice_inv, stress_3_3, rot): """ Return symmetrized stress lattice vectors expected as row vectors (same as ASE get_cell() convention), inv_lattice is its matrix inverse (reciprocal().T) """ scaled_stress =, stress_3_3), lattice.T) symmetrized_scaled_stress = np.zeros((3, 3)) for r in rot: symmetrized_scaled_stress +=, scaled_stress), r) symmetrized_scaled_stress /= len(rot) sym =, symmetrized_scaled_stress), lattice_inv.T) return sym
[docs]class FixSymmetry(FixConstraint): """ Constraint to preserve spacegroup symmetry during optimisation. Requires spglib package to be available. """ def __init__(self, atoms, symprec=0.01, adjust_positions=True, adjust_cell=True, verbose=False): self.verbose = verbose refine_symmetry(atoms, symprec, self.verbose) # refine initial symmetry sym = prep_symmetry(atoms, symprec, self.verbose) self.rotations, self.translations, self.symm_map = sym self.do_adjust_positions = adjust_positions self.do_adjust_cell = adjust_cell def adjust_cell(self, atoms, cell): if not self.do_adjust_cell: return # stress should definitely be symmetrized as a rank 2 tensor # UnitCellFilter uses deformation gradient as cell DOF with steps # dF = stress.F^-T quantity that should be symmetrized is therefore dF . # F^T assume prev F = I, so just symmetrize dF cur_cell = atoms.get_cell() cur_cell_inv = atoms.cell.reciprocal().T # F defined such that cell = cur_cell . F^T # assume prev F = I, so dF = F - I delta_deform_grad =, cell).T - np.eye(3) # symmetrization doesn't work properly with large steps, since # it depends on current cell, and cell is being changed by deformation # gradient max_delta_deform_grad = np.max(np.abs(delta_deform_grad)) if max_delta_deform_grad > 0.25: raise RuntimeError('FixSymmetry adjust_cell does not work properly' ' with large deformation gradient step {} > 0.25' .format(max_delta_deform_grad)) elif max_delta_deform_grad > 0.15: warnings.warn('FixSymmetry adjust_cell may be ill behaved with' ' large deformation gradient step {}' .format(max_delta_deform_grad)) symmetrized_delta_deform_grad = symmetrize_rank2(cur_cell, cur_cell_inv, delta_deform_grad, self.rotations) cell[:] =, (symmetrized_delta_deform_grad + np.eye(3)).T) def adjust_positions(self, atoms, new): if not self.do_adjust_positions: return # symmetrize changes in position as rank 1 tensors step = new - atoms.positions symmetrized_step = symmetrize_rank1(atoms.get_cell(), atoms.cell.reciprocal().T, step, self.rotations, self.translations, self.symm_map) new[:] = atoms.positions + symmetrized_step def adjust_forces(self, atoms, forces): # symmetrize forces as rank 1 tensors # print('adjusting forces') forces[:] = symmetrize_rank1(atoms.get_cell(), atoms.cell.reciprocal().T, forces, self.rotations, self.translations, self.symm_map) def adjust_stress(self, atoms, stress): # symmetrize stress as rank 2 tensor raw_stress = voigt_6_to_full_3x3_stress(stress) symmetrized_stress = symmetrize_rank2(atoms.get_cell(), atoms.cell.reciprocal().T, raw_stress, self.rotations) stress[:] = full_3x3_to_voigt_6_stress(symmetrized_stress) def index_shuffle(self, atoms, ind): if len(atoms) != len(ind) or len(set(ind)) != len(ind): raise RuntimeError("FixSymmetry can only accomodate atom" " permutions, and len(Atoms) == {} " "!= len(ind) == {} or ind has duplicates" .format(len(atoms), len(ind))) ind_reversed = np.zeros((len(ind)), dtype=int) ind_reversed[ind] = range(len(ind)) new_symm_map = [] for sm in self.symm_map: new_sm = np.array([-1] * len(atoms)) for at_i in range(len(ind)): new_sm[ind_reversed[at_i]] = ind_reversed[sm[at_i]] new_symm_map.append(new_sm) self.symm_map = new_symm_map