Coverage for ase / constraints / mirror_torque.py: 80.00%

1import numpy as np 

2 

3from ase.constraints.constraint import FixConstraint, slice2enlist 

4 

5 

6class MirrorTorque(FixConstraint): 

7 """Constraint object for mirroring the torque acting on a dihedral 

8 angle defined by four atoms. 

9 

10 This class is designed to find a transition state with the help of a 

11 single optimization. It can be used if the transition state belongs to a 

12 cis-trans-isomerization with a change of dihedral angle. First the given 

13 dihedral angle will be fixed until all other degrees of freedom are 

14 optimized, then the torque acting on the dihedral angle will be mirrored 

15 to find the transition state. Transition states in 

16 dependence of the force can be obtained by stretching the molecule and 

17 fixing its total length with *FixBondLength* or by using *ExternalForce* 

18 during the optimization with *MirrorTorque*. 

19 

20 This constraint can be used to find 

21 transition states of cis-trans-isomerization. 

22 

23 a1 a4 

24 | | 

25 a2 __ a3 

26 

27 Parameters 

28 ---------- 

29 a1: int 

30 First atom index. 

31 a2: int 

32 Second atom index. 

33 a3: int 

34 Third atom index. 

35 a4: int 

36 Fourth atom index. 

37 max_angle: float 

38 Upper limit of the dihedral angle interval where the transition state 

39 can be found. 

40 min_angle: float 

41 Lower limit of the dihedral angle interval where the transition state 

42 can be found. 

43 fmax: float 

44 Maximum force used for the optimization. 

45 

46 Notes 

47 ----- 

48 You can combine this constraint for example with FixBondLength but make 

49 sure that *MirrorTorque* comes first in the list if there are overlaps 

50 between atom1-4 and atom5-6: 

51 

52 >>> from ase.build import bulk 

53 

54 >>> atoms = bulk('Cu', 'fcc', a=3.6) 

55 >>> atom1, atom2, atom3, atom4, atom5, atom6 = atoms[:6] 

56 >>> con1 = MirrorTorque(atom1, atom2, atom3, atom4) 

57 >>> con2 = FixBondLength(atom5, atom6) 

58 >>> atoms.set_constraint([con1, con2]) 

59 

60 """ 

61 

62 def __init__( 

63 self, a1, a2, a3, a4, max_angle=2 * np.pi, min_angle=0.0, fmax=0.1 

64 ): 

65 self.indices = [a1, a2, a3, a4] 

66 self.min_angle = min_angle 

67 self.max_angle = max_angle 

68 self.fmax = fmax 

69 

70 def adjust_positions(self, atoms, new): 

71 pass 

72 

73 def adjust_forces(self, atoms, forces): 

74 angle = atoms.get_dihedral( 

75 self.indices[0], self.indices[1], self.indices[2], self.indices[3] 

76 ) 

77 angle *= np.pi / 180.0 

78 if (angle < self.min_angle) or (angle > self.max_angle): 

79 # Stop structure optimization 

80 forces[:] *= 0 

81 return 

82 p = atoms.positions[self.indices] 

83 f = forces[self.indices] 

84 

85 f0 = (f[1] + f[2]) / 2.0 

86 ff = f - f0 

87 p0 = (p[2] + p[1]) / 2.0 

88 m0 = np.cross(p[1] - p0, ff[1]) / (p[1] - p0).dot(p[1] - p0) 

89 fff = ff - np.cross(m0, p - p0) 

90 d1 = np.cross(np.cross(p[1] - p0, p[0] - p[1]), p[1] - p0) / ( 

91 p[1] - p0 

92 ).dot(p[1] - p0) 

93 d2 = np.cross(np.cross(p[2] - p0, p[3] - p[2]), p[2] - p0) / ( 

94 p[2] - p0 

95 ).dot(p[2] - p0) 

96 omegap1 = (np.cross(d1, fff[0]) / d1.dot(d1)).dot( 

97 p[1] - p0 

98 ) / np.linalg.norm(p[1] - p0) 

99 omegap2 = (np.cross(d2, fff[3]) / d2.dot(d2)).dot( 

100 p[2] - p0 

101 ) / np.linalg.norm(p[2] - p0) 

102 omegap = omegap1 + omegap2 

103 con_saved = atoms.constraints 

104 try: 

105 con = [ 

106 con for con in con_saved if not isinstance(con, MirrorTorque) 

107 ] 

108 atoms.set_constraint(con) 

109 forces_copy = atoms.get_forces() 

110 finally: 

111 atoms.set_constraint(con_saved) 

112 df1 = ( 

113 -1 

114 / 2.0 

115 * omegap 

116 * np.cross(p[1] - p0, d1) 

117 / np.linalg.norm(p[1] - p0) 

118 ) 

119 df2 = ( 

120 -1 

121 / 2.0 

122 * omegap 

123 * np.cross(p[2] - p0, d2) 

124 / np.linalg.norm(p[2] - p0) 

125 ) 

126 forces_copy[self.indices] += ( 

127 df1, 

128 [0.0, 0.0, 0.0], 

129 [0.0, 0.0, 0.0], 

130 df2, 

131 ) 

132 # Check if forces would be converged if the dihedral angle with 

133 # mirrored torque would also be fixed 

134 if (forces_copy**2).sum(axis=1).max() < self.fmax**2: 

135 factor = 1.0 

136 else: 

137 factor = 0.0 

138 df1 = ( 

139 -(1 + factor) 

140 / 2.0 

141 * omegap 

142 * np.cross(p[1] - p0, d1) 

143 / np.linalg.norm(p[1] - p0) 

144 ) 

145 df2 = ( 

146 -(1 + factor) 

147 / 2.0 

148 * omegap 

149 * np.cross(p[2] - p0, d2) 

150 / np.linalg.norm(p[2] - p0) 

151 ) 

152 forces[self.indices] += (df1, [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], df2) 

153 

154 def index_shuffle(self, atoms, ind): 

155 # See docstring of superclass 

156 indices = [] 

157 for new, old in slice2enlist(ind, len(atoms)): 

158 if old in self.indices: 

159 indices.append(new) 

160 if len(indices) == 0: 

161 raise IndexError('All indices in MirrorTorque not part of slice') 

162 self.indices = np.asarray(indices, int) 

163 

164 def __repr__(self): 

165 return 'MirrorTorque(%d, %d, %d, %d, %f, %f, %f)' % ( 

166 self.indices[0], 

167 self.indices[1], 

168 self.indices[2], 

169 self.indices[3], 

170 self.max_angle, 

171 self.min_angle, 

172 self.fmax, 

173 ) 

174 

175 def todict(self): 

176 return { 

177 'name': 'MirrorTorque', 

178 'kwargs': { 

179 'a1': self.indices[0], 

180 'a2': self.indices[1], 

181 'a3': self.indices[2], 

182 'a4': self.indices[3], 

183 'max_angle': self.max_angle, 

184 'min_angle': self.min_angle, 

185 'fmax': self.fmax, 

186 }, 

187 }