Coverage for ase / constraints / mirror_force.py: 73.33%

1import numpy as np 

2 

3from ase.constraints.constraint import FixConstraint, slice2enlist 

4 

5 

6class MirrorForce(FixConstraint): 

7 """Constraint object for mirroring the force between two atoms. 

8 

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

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

11 bond breaking reaction. First the given bond length will be fixed until 

12 all other degrees of freedom are optimized, then the forces of the two 

13 atoms will be mirrored to find the transition state. The mirror plane is 

14 perpendicular to the connecting line of the atoms. Transition states in 

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

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

17 during the optimization with *MirrorForce*. 

18 

19 Parameters 

20 ---------- 

21 a1: int 

22 First atom index. 

23 a2: int 

24 Second atom index. 

25 max_dist: float 

26 Upper limit of the bond length interval where the transition state 

27 can be found. 

28 min_dist: float 

29 Lower limit of the bond length interval where the transition state 

30 can be found. 

31 fmax: float 

32 Maximum force used for the optimization. 

33 

34 Notes 

35 ----- 

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

37 sure that *MirrorForce* comes first in the list if there are overlaps 

38 between atom1-2 and atom3-4: 

39 

40 >>> from ase.build import bulk 

41 

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

43 >>> atom1, atom2, atom3, atom4 = atoms[:4] 

44 >>> con1 = MirrorForce(atom1, atom2) 

45 >>> con2 = FixBondLength(atom3, atom4) 

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

47 

48 """ 

49 

50 def __init__(self, a1, a2, max_dist=2.5, min_dist=1.0, fmax=0.1): 

51 self.indices = [a1, a2] 

52 self.min_dist = min_dist 

53 self.max_dist = max_dist 

54 self.fmax = fmax 

55 

56 def adjust_positions(self, atoms, new): 

57 pass 

58 

59 def adjust_forces(self, atoms, forces): 

60 dist = np.subtract.reduce(atoms.positions[self.indices]) 

61 d = np.linalg.norm(dist) 

62 if (d < self.min_dist) or (d > self.max_dist): 

63 # Stop structure optimization 

64 forces[:] *= 0 

65 return 

66 dist /= d 

67 df = np.subtract.reduce(forces[self.indices]) 

68 f = df.dot(dist) 

69 con_saved = atoms.constraints 

70 try: 

71 con = [con for con in con_saved if not isinstance(con, MirrorForce)] 

72 atoms.set_constraint(con) 

73 forces_copy = atoms.get_forces() 

74 finally: 

75 atoms.set_constraint(con_saved) 

76 df1 = -1 / 2.0 * f * dist 

77 forces_copy[self.indices] += (df1, -df1) 

78 # Check if forces would be converged if the bond with mirrored forces 

79 # would also be fixed 

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

81 factor = 1.0 

82 else: 

83 factor = 0.0 

84 df1 = -(1 + factor) / 2.0 * f * dist 

85 forces[self.indices] += (df1, -df1) 

86 

87 def index_shuffle(self, atoms, ind): 

88 """Shuffle the indices of the two atoms in this constraint""" 

89 newa = [-1, -1] # Signal error 

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

91 for i, a in enumerate(self.indices): 

92 if old == a: 

93 newa[i] = new 

94 if newa[0] == -1 or newa[1] == -1: 

95 raise IndexError('Constraint not part of slice') 

96 self.indices = newa 

97 

98 def __repr__(self): 

99 return 'MirrorForce(%d, %d, %f, %f, %f)' % ( 

100 self.indices[0], 

101 self.indices[1], 

102 self.max_dist, 

103 self.min_dist, 

104 self.fmax, 

105 ) 

106 

107 def todict(self): 

108 return { 

109 'name': 'MirrorForce', 

110 'kwargs': { 

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

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

113 'max_dist': self.max_dist, 

114 'min_dist': self.min_dist, 

115 'fmax': self.fmax, 

116 }, 

117 }