Coverage for ase / constraints / fix_symmetry.py: 90.74%

1from warnings import warn 

2 

3import numpy as np 

4 

5from ase.constraints.constraint import FixConstraint 

6from ase.spacegroup.symmetrize import ( 

7 prep_symmetry, 

8 refine_symmetry, 

9 symmetrize_rank1, 

10 symmetrize_rank2, 

11) 

12from ase.stress import full_3x3_to_voigt_6_stress, voigt_6_to_full_3x3_stress 

13 

14 

15class FixSymmetry(FixConstraint): 

16 """ 

17 Constraint to preserve spacegroup symmetry during optimisation. 

18 

19 Requires spglib package to be available. 

20 """ 

21 

22 def __init__( 

23 self, 

24 atoms, 

25 symprec=0.01, 

26 adjust_positions=True, 

27 adjust_cell=True, 

28 verbose=False, 

29 ): 

30 self.atoms = atoms.copy() 

31 self.symprec = symprec 

32 self.verbose = verbose 

33 refine_symmetry(atoms, symprec, self.verbose) # refine initial symmetry 

34 sym = prep_symmetry(atoms, symprec, self.verbose) 

35 self.rotations, self.translations, self.symm_map = sym 

36 self.do_adjust_positions = adjust_positions 

37 self.do_adjust_cell = adjust_cell 

38 

39 def adjust_cell(self, atoms, cell): 

40 if not self.do_adjust_cell: 

41 return 

42 # stress should definitely be symmetrized as a rank 2 tensor 

43 # UnitCellFilter uses deformation gradient as cell DOF with steps 

44 # dF = stress.F^-T quantity that should be symmetrized is therefore dF . 

45 # F^T assume prev F = I, so just symmetrize dF 

46 cur_cell = atoms.get_cell() 

47 cur_cell_inv = atoms.cell.reciprocal().T 

48 

49 # F defined such that cell = cur_cell . F^T 

50 # assume prev F = I, so dF = F - I 

51 delta_deform_grad = np.dot(cur_cell_inv, cell).T - np.eye(3) 

52 

53 # symmetrization doesn't work properly with large steps, since 

54 # it depends on current cell, and cell is being changed by deformation 

55 # gradient 

56 max_delta_deform_grad = np.max(np.abs(delta_deform_grad)) 

57 if max_delta_deform_grad > 0.25: 

58 raise RuntimeError( 

59 'FixSymmetry adjust_cell does not work properly' 

60 ' with large deformation gradient step {} > 0.25'.format( 

61 max_delta_deform_grad 

62 ) 

63 ) 

64 elif max_delta_deform_grad > 0.15: 

65 warn( 

66 'FixSymmetry adjust_cell may be ill behaved with large ' 

67 'deformation gradient step {}'.format(max_delta_deform_grad) 

68 ) 

69 

70 symmetrized_delta_deform_grad = symmetrize_rank2( 

71 cur_cell, cur_cell_inv, delta_deform_grad, self.rotations 

72 ) 

73 cell[:] = np.dot( 

74 cur_cell, (symmetrized_delta_deform_grad + np.eye(3)).T 

75 ) 

76 

77 def adjust_positions(self, atoms, new): 

78 if not self.do_adjust_positions: 

79 return 

80 # symmetrize changes in position as rank 1 tensors 

81 step = new - atoms.positions 

82 symmetrized_step = symmetrize_rank1( 

83 atoms.get_cell(), 

84 atoms.cell.reciprocal().T, 

85 step, 

86 self.rotations, 

87 self.translations, 

88 self.symm_map, 

89 ) 

90 new[:] = atoms.positions + symmetrized_step 

91 

92 def adjust_forces(self, atoms, forces): 

93 # symmetrize forces as rank 1 tensors 

94 # print('adjusting forces') 

95 forces[:] = symmetrize_rank1( 

96 atoms.get_cell(), 

97 atoms.cell.reciprocal().T, 

98 forces, 

99 self.rotations, 

100 self.translations, 

101 self.symm_map, 

102 ) 

103 

104 def adjust_stress(self, atoms, stress): 

105 # symmetrize stress as rank 2 tensor 

106 raw_stress = voigt_6_to_full_3x3_stress(stress) 

107 symmetrized_stress = symmetrize_rank2( 

108 atoms.get_cell(), 

109 atoms.cell.reciprocal().T, 

110 raw_stress, 

111 self.rotations, 

112 ) 

113 stress[:] = full_3x3_to_voigt_6_stress(symmetrized_stress) 

114 

115 def index_shuffle(self, atoms, ind): 

116 if len(atoms) != len(ind) or len(set(ind)) != len(ind): 

117 raise RuntimeError( 

118 'FixSymmetry can only accomodate atom' 

119 ' permutions, and len(Atoms) == {} ' 

120 '!= len(ind) == {} or ind has duplicates'.format( 

121 len(atoms), len(ind) 

122 ) 

123 ) 

124 

125 ind_reversed = np.zeros((len(ind)), dtype=int) 

126 ind_reversed[ind] = range(len(ind)) 

127 new_symm_map = [] 

128 for sm in self.symm_map: 

129 new_sm = np.array([-1] * len(atoms)) 

130 for at_i in range(len(ind)): 

131 new_sm[ind_reversed[at_i]] = ind_reversed[sm[at_i]] 

132 new_symm_map.append(new_sm) 

133 

134 self.symm_map = new_symm_map 

135 

136 def todict(self): 

137 return { 

138 'name': 'FixSymmetry', 

139 'kwargs': { 

140 'atoms': self.atoms, 

141 'symprec': self.symprec, 

142 'adjust_positions': self.do_adjust_positions, 

143 'adjust_cell': self.do_adjust_cell, 

144 'verbose': self.verbose, 

145 }, 

146 }