1from numpy.linalg import norm 

2from ase.optimize.bfgs import BFGS 

3from ase.constraints import FixInternals 

4 

5 

6class BFGSClimbFixInternals(BFGS): 

7 """Class for transition state search and optimization 

8 

9 Climbs the 1D reaction coordinate defined as constrained internal coordinate 

10 via the :class:`~ase.constraints.FixInternals` class while minimizing all 

11 remaining degrees of freedom. 

12 

13 Details: Two optimizers, 'A' and 'B', are applied orthogonal to each other. 

14 Optimizer 'A' climbs the constrained coordinate while optimizer 'B' 

15 optimizes the remaining degrees of freedom after each climbing step. 

16 Optimizer 'A' uses the BFGS algorithm to climb along the projected force of 

17 the selected constraint. Optimizer 'B' can be any ASE optimizer 

18 (default: BFGS). 

19 

20 In combination with other constraints, the order of constraints matters. 

21 Generally, the FixInternals constraint should come first in the list of 

22 constraints. 

23 This has been tested with the :class:`~ase.constraints.FixAtoms` constraint. 

24 

25 Inspired by concepts described by P. N. Plessow. [1]_ 

26 

27 .. [1] Plessow, P. N., Efficient Transition State Optimization of Periodic 

28 Structures through Automated Relaxed Potential Energy Surface Scans. 

29 J. Chem. Theory Comput. 2018, 14 (2), 981–990. 

30 https://doi.org/10.1021/acs.jctc.7b01070. 

31 

32 .. note:: 

33 Convergence is based on 'fmax' of the total forces which is the sum of 

34 the projected forces and the forces of the remaining degrees of freedom. 

35 This value is logged in the 'logfile'. Optimizer 'B' logs 'fmax' of the 

36 remaining degrees of freedom without the projected forces. The projected 

37 forces can be inspected using the :meth:`get_projected_forces` method: 

38 

39 >>> for _ in dyn.irun(): 

40 ... projected_forces = dyn.get_projected_forces() 

41 

42 Example 

43 ------- 

44 .. literalinclude:: ../../ase/test/optimize/test_climb_fix_internals.py 

45 :end-before: # end example for documentation 

46 """ 

47 def __init__(self, atoms, restart=None, logfile='-', trajectory=None, 

48 maxstep=None, master=None, alpha=None, 

49 climb_coordinate=None, 

50 optB=BFGS, optB_kwargs=None, optB_fmax=0.05, 

51 optB_fmax_scaling=0.0): 

52 """Allowed parameters are similar to the parent class 

53 :class:`~ase.optimize.bfgs.BFGS` with the following additions: 

54 

55 Parameters 

56 ---------- 

57 climb_coordinate: list 

58 Specifies which subconstraint of the 

59 :class:`~ase.constraints.FixInternals` constraint is to be climbed. 

60 Provide the corresponding nested list of indices 

61 (including coefficients in the case of Combo constraints). 

62 See the example above. 

63 

64 optB: any ASE optimizer, optional 

65 Optimizer 'B' for optimization of the remaining degrees of freedom 

66 after each climbing step. 

67 

68 optB_kwargs: dict, optional 

69 Specifies keyword arguments to be passed to optimizer 'B' at its 

70 initialization. By default, optimizer 'B' writes a logfile and 

71 trajectory (optB_{...}.log, optB_{...}.traj) where {...} is the 

72 current value of the ``climb_coordinate``. Set ``logfile`` to '-' 

73 for console output. Set ``trajectory`` to 'None' to suppress 

74 writing of the trajectory file. 

75 

76 optB_fmax: float, optional 

77 Specifies the convergence criterion 'fmax' of optimizer 'B'. 

78 

79 optB_fmax_scaling: float, optional 

80 Scaling factor to dynamically tighten 'fmax' of optimizer 'B' to 

81 the value of ``optB_fmax`` when close to convergence. 

82 Can speed up the climbing process. The scaling formula is 

83 

84 'fmax' = ``optB_fmax`` + ``optB_fmax_scaling`` 

85 :math:`\\cdot` norm_of_projected_forces 

86 

87 The final optimization with optimizer 'B' is 

88 performed with ``optB_fmax`` independent of ``optB_fmax_scaling``. 

89 """ 

90 self.targetvalue = None # may be assigned during restart in self.read() 

91 super().__init__(atoms, restart=restart, logfile=logfile, 

92 trajectory=trajectory, maxstep=maxstep, master=master, 

93 alpha=alpha) 

94 

95 self.constr2climb = self.get_constr2climb(self.atoms, climb_coordinate) 

96 self.targetvalue = self.targetvalue or self.constr2climb.targetvalue 

97 

98 self.optB = optB 

99 self.optB_kwargs = optB_kwargs or {} 

100 self.optB_fmax = optB_fmax 

101 self.scaling = optB_fmax_scaling 

102 # log optimizer 'B' in logfiles named after current value of constraint 

103 self.autolog = 'logfile' not in self.optB_kwargs 

104 self.autotraj = 'trajectory' not in self.optB_kwargs 

105 

106 def get_constr2climb(self, atoms, climb_coordinate): 

107 """Get pointer to the subconstraint that is to be climbed. 

108 Identification by its definition via indices (and coefficients).""" 

109 constr = self.get_fixinternals(atoms) 

110 return constr.get_subconstraint(atoms, climb_coordinate) 

111 

112 def get_fixinternals(self, atoms): 

113 """Get pointer to the FixInternals constraint on the atoms object.""" 

114 all_constr_types = list(map(type, atoms.constraints)) 

115 index = all_constr_types.index(FixInternals) # locate constraint 

116 return atoms.constraints[index] 

117 

118 def read(self): 

119 (self.H, self.pos0, self.forces0, self.maxstep, 

120 self.targetvalue) = self.load() 

121 

122 def step(self): 

123 self.relax_remaining_dof() # optimization with optimizer 'B' 

124 

125 pos, dpos = self.pretend2climb() # with optimizer 'A' 

126 self.update_positions_and_targetvalue(pos, dpos) # obey other constr. 

127 

128 self.dump((self.H, self.pos0, self.forces0, self.maxstep, 

129 self.targetvalue)) 

130 

131 def pretend2climb(self): 

132 """Get directions for climbing and climb with optimizer 'A'.""" 

133 proj_forces = self.get_projected_forces() 

134 pos = self.atoms.get_positions() 

135 dpos, steplengths = self.prepare_step(pos, proj_forces) 

136 dpos = self.determine_step(dpos, steplengths) 

137 return pos, dpos 

138 

139 def update_positions_and_targetvalue(self, pos, dpos): 

140 """Adjust constrained targetvalue of constraint and update positions.""" 

141 self.constr2climb.adjust_positions(pos, pos + dpos) # update sigma 

142 self.targetvalue += self.constr2climb.sigma # climb constraint 

143 self.constr2climb.targetvalue = self.targetvalue # adjust positions 

144 self.atoms.set_positions(self.atoms.get_positions()) # to targetvalue 

145 

146 def relax_remaining_dof(self): 

147 """Optimize remaining degrees of freedom with optimizer 'B'.""" 

148 if self.autolog: 

149 self.optB_kwargs['logfile'] = f'optB_{self.targetvalue}.log' 

150 if self.autotraj: 

151 self.optB_kwargs['trajectory'] = f'optB_{self.targetvalue}.traj' 

152 fmax = self.get_scaled_fmax() 

153 with self.optB(self.atoms, **self.optB_kwargs) as opt: 

154 opt.run(fmax) # optimize with scaled fmax 

155 if self.converged() and fmax > self.optB_fmax: 

156 # (final) optimization with desired fmax 

157 opt.run(self.optB_fmax) 

158 

159 def get_scaled_fmax(self): 

160 """Return the adaptive 'fmax' based on the estimated distance to the 

161 transition state.""" 

162 return (self.optB_fmax + 

163 self.scaling * norm(self.constr2climb.projected_forces)) 

164 

165 def get_projected_forces(self): 

166 """Return the projected forces along the constrained coordinate in 

167 uphill direction (negative sign).""" 

168 forces = self.constr2climb.projected_forces 

169 forces = -forces.reshape(self.atoms.positions.shape) 

170 return forces 

171 

172 def get_total_forces(self): 

173 """Return forces obeying all constraints plus projected forces.""" 

174 return self.atoms.get_forces() + self.get_projected_forces() 

175 

176 def converged(self, forces=None): 

177 """Did the optimization converge based on the total forces?""" 

178 forces = forces or self.get_total_forces() 

179 return super().converged(forces=forces) 

180 

181 def log(self, forces=None): 

182 forces = forces or self.get_total_forces() 

183 super().log(forces=forces)