Coverage for /builds/ase/ase/ase/optimize/bfgs.py: 78.41%

1# fmt: off 

2 

3import warnings 

4from pathlib import Path 

5from typing import IO, Optional, Union 

6 

7import numpy as np 

8from numpy.linalg import eigh 

9 

10from ase import Atoms 

11from ase.optimize.optimize import Optimizer, UnitCellFilter 

12 

13 

14class BFGS(Optimizer): 

15 # default parameters 

16 defaults = {**Optimizer.defaults, 'alpha': 70.0} 

17 

18 def __init__( 

19 self, 

20 atoms: Atoms, 

21 restart: Optional[str] = None, 

22 logfile: Optional[Union[IO, str, Path]] = '-', 

23 trajectory: Optional[Union[str, Path]] = None, 

24 append_trajectory: bool = False, 

25 maxstep: Optional[float] = None, 

26 alpha: Optional[float] = None, 

27 **kwargs, 

28 ): 

29 """BFGS optimizer. 

30 

31 Parameters 

32 ---------- 

33 atoms: :class:`~ase.Atoms` 

34 The Atoms object to relax. 

35 

36 restart: str 

37 JSON file used to store hessian matrix. If set, file with 

38 such a name will be searched and hessian matrix stored will 

39 be used, if the file exists. 

40 

41 trajectory: str or Path 

42 Trajectory file used to store optimisation path. 

43 

44 logfile: file object, Path, or str 

45 If *logfile* is a string, a file with that name will be opened. 

46 Use '-' for stdout. 

47 

48 maxstep: float 

49 Used to set the maximum distance an atom can move per 

50 iteration (default value is 0.2 Å). 

51 

52 alpha: float 

53 Initial guess for the Hessian (curvature of energy surface). A 

54 conservative value of 70.0 is the default, but number of needed 

55 steps to converge might be less if a lower value is used. However, 

56 a lower value also means risk of instability. 

57 

58 kwargs : dict, optional 

59 Extra arguments passed to 

60 :class:`~ase.optimize.optimize.Optimizer`. 

61 

62 """ 

63 if maxstep is None: 

64 self.maxstep = self.defaults['maxstep'] 

65 else: 

66 self.maxstep = maxstep 

67 

68 if self.maxstep > 1.0: 

69 warnings.warn('You are using a *very* large value for ' 

70 'the maximum step size: %.1f Å' % self.maxstep) 

71 

72 self.alpha = alpha 

73 if self.alpha is None: 

74 self.alpha = self.defaults['alpha'] 

75 Optimizer.__init__(self, atoms=atoms, restart=restart, 

76 logfile=logfile, trajectory=trajectory, 

77 append_trajectory=append_trajectory, 

78 **kwargs) 

79 

80 def initialize(self): 

81 # initial hessian 

82 self.H0 = np.eye(self.optimizable.ndofs()) * self.alpha 

83 

84 self.H = None 

85 self.pos0 = None 

86 self.forces0 = None 

87 

88 def read(self): 

89 file = self.load() 

90 if len(file) == 5: 

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

92 self.atoms.orig_cell) = file 

93 else: 

94 self.H, self.pos0, self.forces0, self.maxstep = file 

95 

96 def step(self, gradient=None): 

97 optimizable = self.optimizable 

98 

99 if gradient is None: 

100 gradient = optimizable.get_gradient() 

101 

102 pos = optimizable.get_x() 

103 dpos, steplengths = self.prepare_step(pos, gradient) 

104 dpos = self.determine_step(dpos, steplengths) 

105 optimizable.set_x(pos + dpos) 

106 if isinstance(self.atoms, UnitCellFilter): 

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

108 self.atoms.orig_cell)) 

109 else: 

110 self.dump((self.H, self.pos0, self.forces0, self.maxstep)) 

111 

112 def prepare_step(self, pos, gradient): 

113 pos = pos.ravel() 

114 gradient = gradient.ravel() 

115 self.update(pos, gradient, self.pos0, self.forces0) 

116 omega, V = eigh(self.H) 

117 

118 # FUTURE: Log this properly 

119 # # check for negative eigenvalues of the hessian 

120 # if any(omega < 0): 

121 # n_negative = len(omega[omega < 0]) 

122 # msg = '\n** BFGS Hessian has {} negative eigenvalues.'.format( 

123 # n_negative 

124 # ) 

125 # print(msg, flush=True) 

126 # if self.logfile is not None: 

127 # self.logfile.write(msg) 

128 # self.logfile.flush() 

129 

130 dpos = np.dot(V, np.dot(gradient, V) / np.fabs(omega)) 

131 # XXX Here we are calling gradient_norm() on some positions. 

132 # Should there be a general norm concept 

133 steplengths = self.optimizable.gradient_norm(dpos) 

134 self.pos0 = pos 

135 self.forces0 = gradient.copy() 

136 return dpos, steplengths 

137 

138 def determine_step(self, dpos, steplengths): 

139 """Determine step to take according to maxstep 

140 

141 Normalize all steps as the largest step. This way 

142 we still move along the direction. 

143 """ 

144 maxsteplength = np.max(steplengths) 

145 if maxsteplength >= self.maxstep: 

146 scale = self.maxstep / maxsteplength 

147 # FUTURE: Log this properly 

148 # msg = '\n** scale step by {:.3f} to be shorter than {}'.format( 

149 # scale, self.maxstep 

150 # ) 

151 # print(msg, flush=True) 

152 

153 dpos *= scale 

154 return dpos 

155 

156 def update(self, pos, forces, pos0, forces0): 

157 if self.H is None: 

158 self.H = self.H0 

159 return 

160 dpos = pos - pos0 

161 

162 if np.abs(dpos).max() < 1e-7: 

163 # Same configuration again (maybe a restart): 

164 return 

165 

166 dforces = forces - forces0 

167 a = np.dot(dpos, dforces) 

168 dg = np.dot(self.H, dpos) 

169 b = np.dot(dpos, dg) 

170 self.H -= np.outer(dforces, dforces) / a + np.outer(dg, dg) / b 

171 

172 def replay_trajectory(self, traj): 

173 """Initialize hessian from old trajectory.""" 

174 if isinstance(traj, str): 

175 from ase.io.trajectory import Trajectory 

176 traj = Trajectory(traj, 'r') 

177 self.H = None 

178 atoms = traj[0] 

179 pos0 = atoms.get_positions().ravel() 

180 forces0 = atoms.get_forces().ravel() 

181 for atoms in traj: 

182 pos = atoms.get_positions().ravel() 

183 forces = atoms.get_forces().ravel() 

184 self.update(pos, forces, pos0, forces0) 

185 pos0 = pos 

186 forces0 = forces 

187 

188 self.pos0 = pos0 

189 self.forces0 = forces0 

190 

191 

192class oldBFGS(BFGS): 

193 def determine_step(self, dpos, steplengths): 

194 """Old BFGS behaviour for scaling step lengths 

195 

196 This keeps the behaviour of truncating individual steps. Some might 

197 depend of this as some absurd kind of stimulated annealing to find the 

198 global minimum. 

199 """ 

200 dpos /= np.maximum(steplengths / self.maxstep, 1.0).reshape(-1, 1) 

201 return dpos