1import numpy as np 

2 

3from ase.optimize.optimize import Optimizer 

4from ase.utils.linesearch import LineSearch 

5 

6 

7class LBFGS(Optimizer): 

8 """Limited memory BFGS optimizer. 

9 

10 A limited memory version of the bfgs algorithm. Unlike the bfgs algorithm 

11 used in bfgs.py, the inverse of Hessian matrix is updated. The inverse 

12 Hessian is represented only as a diagonal matrix to save memory 

13 

14 """ 

15 

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

17 maxstep=None, memory=100, damping=1.0, alpha=70.0, 

18 use_line_search=False, master=None, 

19 force_consistent=None): 

20 """Parameters: 

21 

22 atoms: Atoms object 

23 The Atoms object to relax. 

24 

25 restart: string 

26 Pickle file used to store vectors for updating the inverse of 

27 Hessian matrix. If set, file with such a name will be searched 

28 and information stored will be used, if the file exists. 

29 

30 logfile: file object or str 

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

32 Use '-' for stdout. 

33 

34 trajectory: string 

35 Pickle file used to store trajectory of atomic movement. 

36 

37 maxstep: float 

38 How far is a single atom allowed to move. This is useful for DFT 

39 calculations where wavefunctions can be reused if steps are small. 

40 Default is 0.2 Angstrom. 

41 

42 memory: int 

43 Number of steps to be stored. Default value is 100. Three numpy 

44 arrays of this length containing floats are stored. 

45 

46 damping: float 

47 The calculated step is multiplied with this number before added to 

48 the positions. 

49 

50 alpha: float 

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

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

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

54 a lower value also means risk of instability. 

55 

56 master: boolean 

57 Defaults to None, which causes only rank 0 to save files. If 

58 set to true, this rank will save files. 

59 

60 force_consistent: boolean or None 

61 Use force-consistent energy calls (as opposed to the energy 

62 extrapolated to 0 K). By default (force_consistent=None) uses 

63 force-consistent energies if available in the calculator, but 

64 falls back to force_consistent=False if not. 

65 """ 

66 Optimizer.__init__(self, atoms, restart, logfile, trajectory, master, 

67 force_consistent=force_consistent) 

68 

69 if maxstep is not None: 

70 self.maxstep = maxstep 

71 else: 

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

73 

74 if self.maxstep > 1.0: 

75 raise ValueError('You are using a much too large value for ' + 

76 'the maximum step size: %.1f Angstrom' % 

77 maxstep) 

78 

79 self.memory = memory 

80 # Initial approximation of inverse Hessian 1./70. is to emulate the 

81 # behaviour of BFGS. Note that this is never changed! 

82 self.H0 = 1. / alpha 

83 self.damping = damping 

84 self.use_line_search = use_line_search 

85 self.p = None 

86 self.function_calls = 0 

87 self.force_calls = 0 

88 

89 def initialize(self): 

90 """Initialize everything so no checks have to be done in step""" 

91 self.iteration = 0 

92 self.s = [] 

93 self.y = [] 

94 # Store also rho, to avoid calculating the dot product again and 

95 # again. 

96 self.rho = [] 

97 

98 self.r0 = None 

99 self.f0 = None 

100 self.e0 = None 

101 self.task = 'START' 

102 self.load_restart = False 

103 

104 def read(self): 

105 """Load saved arrays to reconstruct the Hessian""" 

106 self.iteration, self.s, self.y, self.rho, \ 

107 self.r0, self.f0, self.e0, self.task = self.load() 

108 self.load_restart = True 

109 

110 def step(self, forces=None): 

111 """Take a single step 

112 

113 Use the given forces, update the history and calculate the next step -- 

114 then take it""" 

115 

116 if forces is None: 

117 forces = self.atoms.get_forces() 

118 

119 pos = self.atoms.get_positions() 

120 

121 self.update(pos, forces, self.r0, self.f0) 

122 

123 s = self.s 

124 y = self.y 

125 rho = self.rho 

126 H0 = self.H0 

127 

128 loopmax = np.min([self.memory, self.iteration]) 

129 a = np.empty((loopmax,), dtype=np.float64) 

130 

131 # ## The algorithm itself: 

132 q = -forces.reshape(-1) 

133 for i in range(loopmax - 1, -1, -1): 

134 a[i] = rho[i] * np.dot(s[i], q) 

135 q -= a[i] * y[i] 

136 z = H0 * q 

137 

138 for i in range(loopmax): 

139 b = rho[i] * np.dot(y[i], z) 

140 z += s[i] * (a[i] - b) 

141 

142 self.p = - z.reshape((-1, 3)) 

143 # ## 

144 

145 g = -forces 

146 if self.use_line_search is True: 

147 e = self.func(pos) 

148 self.line_search(pos, g, e) 

149 dr = (self.alpha_k * self.p).reshape(len(self.atoms), -1) 

150 else: 

151 self.force_calls += 1 

152 self.function_calls += 1 

153 dr = self.determine_step(self.p) * self.damping 

154 self.atoms.set_positions(pos + dr) 

155 

156 self.iteration += 1 

157 self.r0 = pos 

158 self.f0 = -g 

159 self.dump((self.iteration, self.s, self.y, 

160 self.rho, self.r0, self.f0, self.e0, self.task)) 

161 

162 def determine_step(self, dr): 

163 """Determine step to take according to maxstep 

164 

165 Normalize all steps as the largest step. This way 

166 we still move along the eigendirection. 

167 """ 

168 steplengths = (dr**2).sum(1)**0.5 

169 longest_step = np.max(steplengths) 

170 if longest_step >= self.maxstep: 

171 dr *= self.maxstep / longest_step 

172 

173 return dr 

174 

175 def update(self, pos, forces, r0, f0): 

176 """Update everything that is kept in memory 

177 

178 This function is mostly here to allow for replay_trajectory. 

179 """ 

180 if self.iteration > 0: 

181 s0 = pos.reshape(-1) - r0.reshape(-1) 

182 self.s.append(s0) 

183 

184 # We use the gradient which is minus the force! 

185 y0 = f0.reshape(-1) - forces.reshape(-1) 

186 self.y.append(y0) 

187 

188 rho0 = 1.0 / np.dot(y0, s0) 

189 self.rho.append(rho0) 

190 

191 if self.iteration > self.memory: 

192 self.s.pop(0) 

193 self.y.pop(0) 

194 self.rho.pop(0) 

195 

196 def replay_trajectory(self, traj): 

197 """Initialize history from old trajectory.""" 

198 if isinstance(traj, str): 

199 from ase.io.trajectory import Trajectory 

200 traj = Trajectory(traj, 'r') 

201 r0 = None 

202 f0 = None 

203 # The last element is not added, as we get that for free when taking 

204 # the first qn-step after the replay 

205 for i in range(0, len(traj) - 1): 

206 pos = traj[i].get_positions() 

207 forces = traj[i].get_forces() 

208 self.update(pos, forces, r0, f0) 

209 r0 = pos.copy() 

210 f0 = forces.copy() 

211 self.iteration += 1 

212 self.r0 = r0 

213 self.f0 = f0 

214 

215 def func(self, x): 

216 """Objective function for use of the optimizers""" 

217 self.atoms.set_positions(x.reshape(-1, 3)) 

218 self.function_calls += 1 

219 return self.atoms.get_potential_energy( 

220 force_consistent=self.force_consistent) 

221 

222 def fprime(self, x): 

223 """Gradient of the objective function for use of the optimizers""" 

224 self.atoms.set_positions(x.reshape(-1, 3)) 

225 self.force_calls += 1 

226 # Remember that forces are minus the gradient! 

227 return - self.atoms.get_forces().reshape(-1) 

228 

229 def line_search(self, r, g, e): 

230 self.p = self.p.ravel() 

231 p_size = np.sqrt((self.p**2).sum()) 

232 if p_size <= np.sqrt(len(self.atoms) * 1e-10): 

233 self.p /= (p_size / np.sqrt(len(self.atoms) * 1e-10)) 

234 g = g.ravel() 

235 r = r.ravel() 

236 ls = LineSearch() 

237 self.alpha_k, e, self.e0, self.no_update = \ 

238 ls._line_search(self.func, self.fprime, r, self.p, g, e, self.e0, 

239 maxstep=self.maxstep, c1=.23, 

240 c2=.46, stpmax=50.) 

241 if self.alpha_k is None: 

242 raise RuntimeError('LineSearch failed!') 

243 

244 

245class LBFGSLineSearch(LBFGS): 

246 """This optimizer uses the LBFGS algorithm, but does a line search that 

247 fulfills the Wolff conditions. 

248 """ 

249 

250 def __init__(self, *args, **kwargs): 

251 kwargs['use_line_search'] = True 

252 LBFGS.__init__(self, *args, **kwargs) 

253 

254# """Modified version of LBFGS. 

255# 

256# This optimizer uses the LBFGS algorithm, but does a line search for the 

257# minimum along the search direction. This is done by issuing an additional 

258# force call for each step, thus doubling the number of calculations. 

259# 

260# Additionally the Hessian is reset if the new guess is not sufficiently 

261# better than the old one. 

262# """ 

263# def __init__(self, *args, **kwargs): 

264# self.dR = kwargs.pop('dR', 0.1) 

265# LBFGS.__init__(self, *args, **kwargs) 

266# 

267# def update(self, r, f, r0, f0): 

268# """Update everything that is kept in memory 

269# 

270# This function is mostly here to allow for replay_trajectory. 

271# """ 

272# if self.iteration > 0: 

273# a1 = abs(np.dot(f.reshape(-1), f0.reshape(-1))) 

274# a2 = np.dot(f0.reshape(-1), f0.reshape(-1)) 

275# if not (a1 <= 0.5 * a2 and a2 != 0): 

276# # Reset optimization 

277# self.initialize() 

278# 

279# # Note that the reset above will set self.iteration to 0 again 

280# # which is why we should check again 

281# if self.iteration > 0: 

282# s0 = r.reshape(-1) - r0.reshape(-1) 

283# self.s.append(s0) 

284# 

285# # We use the gradient which is minus the force! 

286# y0 = f0.reshape(-1) - f.reshape(-1) 

287# self.y.append(y0) 

288# 

289# rho0 = 1.0 / np.dot(y0, s0) 

290# self.rho.append(rho0) 

291# 

292# if self.iteration > self.memory: 

293# self.s.pop(0) 

294# self.y.pop(0) 

295# self.rho.pop(0) 

296# 

297# def determine_step(self, dr): 

298# f = self.atoms.get_forces() 

299# 

300# # Unit-vector along the search direction 

301# du = dr / np.sqrt(np.dot(dr.reshape(-1), dr.reshape(-1))) 

302# 

303# # We keep the old step determination before we figure 

304# # out what is the best to do. 

305# maxstep = self.maxstep * np.sqrt(3 * len(self.atoms)) 

306# 

307# # Finite difference step using temporary point 

308# self.atoms.positions += (du * self.dR) 

309# # Decide how much to move along the line du 

310# Fp1 = np.dot(f.reshape(-1), du.reshape(-1)) 

311# Fp2 = np.dot(self.atoms.get_forces().reshape(-1), du.reshape(-1)) 

312# CR = (Fp1 - Fp2) / self.dR 

313# #RdR = Fp1*0.1 

314# if CR < 0.0: 

315# #print "negcurve" 

316# RdR = maxstep 

317# #if(abs(RdR) > maxstep): 

318# # RdR = self.sign(RdR) * maxstep 

319# else: 

320# Fp = (Fp1 + Fp2) * 0.5 

321# RdR = Fp / CR 

322# if abs(RdR) > maxstep: 

323# RdR = np.sign(RdR) * maxstep 

324# else: 

325# RdR += self.dR * 0.5 

326# return du * RdR