1# Copyright (C) 2003 CAMP 

2# Please see the accompanying LICENSE file for further information. 

3 

4""" 

5Quasi-Newton algorithm 

6""" 

7 

8__docformat__ = 'reStructuredText' 

9 

10import time 

11 

12import numpy as np 

13from numpy.linalg import eigh 

14 

15from ase.optimize.optimize import Optimizer 

16 

17 

18def f(lamda, Gbar, b, radius): 

19 b1 = b - lamda 

20 g = radius**2 - np.dot(Gbar / b1, Gbar / b1) 

21 return g 

22 

23 

24def scale_radius_energy(f, r): 

25 scale = 1.0 

26# if(r<=0.01): 

27# return scale 

28 

29 if f < 0.01: 

30 scale *= 1.4 

31 if f < 0.05: 

32 scale *= 1.4 

33 if f < 0.10: 

34 scale *= 1.4 

35 if f < 0.40: 

36 scale *= 1.4 

37 

38 if f > 0.5: 

39 scale *= 1. / 1.4 

40 if f > 0.7: 

41 scale *= 1. / 1.4 

42 if f > 1.0: 

43 scale *= 1. / 1.4 

44 

45 return scale 

46 

47 

48def scale_radius_force(f, r): 

49 scale = 1.0 

50# if(r<=0.01): 

51# return scale 

52 g = abs(f - 1) 

53 if g < 0.01: 

54 scale *= 1.4 

55 if g < 0.05: 

56 scale *= 1.4 

57 if g < 0.10: 

58 scale *= 1.4 

59 if g < 0.40: 

60 scale *= 1.4 

61 

62 if g > 0.5: 

63 scale *= 1. / 1.4 

64 if g > 0.7: 

65 scale *= 1. / 1.4 

66 if g > 1.0: 

67 scale *= 1. / 1.4 

68 

69 return scale 

70 

71 

72def find_lamda(upperlimit, Gbar, b, radius): 

73 lowerlimit = upperlimit 

74 step = 0.1 

75 while f(lowerlimit, Gbar, b, radius) < 0: 

76 lowerlimit -= step 

77 

78 converged = False 

79 

80 while not converged: 

81 

82 midt = (upperlimit + lowerlimit) / 2. 

83 lamda = midt 

84 fmidt = f(midt, Gbar, b, radius) 

85 fupper = f(upperlimit, Gbar, b, radius) 

86 

87 if fupper * fmidt < 0: 

88 lowerlimit = midt 

89 else: 

90 upperlimit = midt 

91 

92 if abs(upperlimit - lowerlimit) < 1e-6: 

93 converged = True 

94 

95 return lamda 

96 

97 

98class GoodOldQuasiNewton(Optimizer): 

99 

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

101 fmax=None, converged=None, 

102 hessianupdate='BFGS', hessian=None, forcemin=True, 

103 verbosity=None, maxradius=None, 

104 diagonal=20., radius=None, 

105 transitionstate=False, master=None): 

106 """Parameters: 

107 

108 atoms: Atoms object 

109 The Atoms object to relax. 

110 

111 restart: string 

112 File used to store hessian matrix. If set, file with 

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

114 be used, if the file exists. 

115 

116 trajectory: string 

117 File used to store trajectory of atomic movement. 

118 

119 maxstep: float 

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

121 iteration (default value is 0.2 Angstroms). 

122 

123 

124 logfile: file object or str 

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

126 Use '-' for stdout. 

127 

128 master: boolean 

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

130 set to true, this rank will save files. 

131 """ 

132 

133 Optimizer.__init__(self, atoms, restart, logfile, trajectory, master) 

134 

135 self.eps = 1e-12 

136 self.hessianupdate = hessianupdate 

137 self.forcemin = forcemin 

138 self.verbosity = verbosity 

139 self.diagonal = diagonal 

140 

141 self.atoms = atoms 

142 

143 n = len(self.atoms) * 3 

144 if radius is None: 

145 self.radius = 0.05 * np.sqrt(n) / 10.0 

146 else: 

147 self.radius = radius 

148 

149 if maxradius is None: 

150 self.maxradius = 0.5 * np.sqrt(n) 

151 else: 

152 self.maxradius = maxradius 

153 

154 # 0.01 < radius < maxradius 

155 self.radius = max(min(self.radius, self.maxradius), 0.0001) 

156 

157 self.transitionstate = transitionstate 

158 

159 # check if this is a nudged elastic band calculation 

160 if hasattr(atoms, 'springconstant'): 

161 self.forcemin = False 

162 

163 self.t0 = time.time() 

164 

165 def initialize(self): 

166 pass 

167 

168 def write_log(self, text): 

169 if self.logfile is not None: 

170 self.logfile.write(text + '\n') 

171 self.logfile.flush() 

172 

173 def set_hessian(self, hessian): 

174 self.hessian = hessian 

175 

176 def get_hessian(self): 

177 if not hasattr(self, 'hessian'): 

178 self.set_default_hessian() 

179 return self.hessian 

180 

181 def set_default_hessian(self): 

182 # set unit matrix 

183 n = len(self.atoms) * 3 

184 hessian = np.zeros((n, n)) 

185 for i in range(n): 

186 hessian[i][i] = self.diagonal 

187 self.set_hessian(hessian) 

188 

189 def update_hessian(self, pos, G): 

190 import copy 

191 if hasattr(self, 'oldG'): 

192 if self.hessianupdate == 'BFGS': 

193 self.update_hessian_bfgs(pos, G) 

194 elif self.hessianupdate == 'Powell': 

195 self.update_hessian_powell(pos, G) 

196 else: 

197 self.update_hessian_bofill(pos, G) 

198 else: 

199 if not hasattr(self, 'hessian'): 

200 self.set_default_hessian() 

201 

202 self.oldpos = copy.copy(pos) 

203 self.oldG = copy.copy(G) 

204 

205 if self.verbosity: 

206 print('hessian ', self.hessian) 

207 

208 def update_hessian_bfgs(self, pos, G): 

209 n = len(self.hessian) 

210 dgrad = G - self.oldG 

211 dpos = pos - self.oldpos 

212 dotg = np.dot(dgrad, dpos) 

213 tvec = np.dot(dpos, self.hessian) 

214 dott = np.dot(dpos, tvec) 

215 if (abs(dott) > self.eps) and (abs(dotg) > self.eps): 

216 for i in range(n): 

217 for j in range(n): 

218 h = dgrad[i] * dgrad[j] / dotg - tvec[i] * tvec[j] / dott 

219 self.hessian[i][j] += h 

220 

221 def update_hessian_powell(self, pos, G): 

222 n = len(self.hessian) 

223 dgrad = G - self.oldG 

224 dpos = pos - self.oldpos 

225 absdpos = np.dot(dpos, dpos) 

226 if absdpos < self.eps: 

227 return 

228 

229 dotg = np.dot(dgrad, dpos) 

230 tvec = dgrad - np.dot(dpos, self.hessian) 

231 tvecdpos = np.dot(tvec, dpos) 

232 ddot = tvecdpos / absdpos 

233 

234 dott = np.dot(dpos, tvec) 

235 if (abs(dott) > self.eps) and (abs(dotg) > self.eps): 

236 for i in range(n): 

237 for j in range(n): 

238 h = tvec[i] * dpos[j] + dpos[i] * \ 

239 tvec[j] - ddot * dpos[i] * dpos[j] 

240 h *= 1. / absdpos 

241 self.hessian[i][j] += h 

242 

243 def update_hessian_bofill(self, pos, G): 

244 print('update Bofill') 

245 n = len(self.hessian) 

246 dgrad = G - self.oldG 

247 dpos = pos - self.oldpos 

248 absdpos = np.dot(dpos, dpos) 

249 if absdpos < self.eps: 

250 return 

251 dotg = np.dot(dgrad, dpos) 

252 tvec = dgrad - np.dot(dpos, self.hessian) 

253 tvecdot = np.dot(tvec, tvec) 

254 tvecdpos = np.dot(tvec, dpos) 

255 

256 coef1 = 1. - tvecdpos * tvecdpos / (absdpos * tvecdot) 

257 coef2 = (1. - coef1) * absdpos / tvecdpos 

258 coef3 = coef1 * tvecdpos / absdpos 

259 

260 dott = np.dot(dpos, tvec) 

261 if (abs(dott) > self.eps) and (abs(dotg) > self.eps): 

262 for i in range(n): 

263 for j in range(n): 

264 h = coef1 * (tvec[i] * dpos[j] + dpos[i] * tvec[j]) - \ 

265 dpos[i] * dpos[j] * coef3 + coef2 * tvec[i] * tvec[j] 

266 h *= 1. / absdpos 

267 self.hessian[i][j] += h 

268 

269 def step(self, f=None): 

270 """ Do one QN step 

271 """ 

272 

273 if f is None: 

274 f = self.atoms.get_forces() 

275 

276 pos = self.atoms.get_positions().ravel() 

277 G = -self.atoms.get_forces().ravel() 

278 energy = self.atoms.get_potential_energy() 

279 

280 if hasattr(self, 'oldenergy'): 

281 

282 self.write_log('energies ' + str(energy) + 

283 ' ' + str(self.oldenergy)) 

284 

285 if self.forcemin: 

286 de = 1e-4 

287 else: 

288 de = 1e-2 

289 

290 if self.transitionstate: 

291 de = 0.2 

292 

293 if (energy - self.oldenergy) > de: 

294 self.write_log('reject step') 

295 self.atoms.set_positions(self.oldpos.reshape((-1, 3))) 

296 G = self.oldG 

297 energy = self.oldenergy 

298 self.radius *= 0.5 

299 else: 

300 self.update_hessian(pos, G) 

301 de = energy - self.oldenergy 

302 f = 1.0 

303 if self.forcemin: 

304 self.write_log( 

305 "energy change; actual: %f estimated: %f " % 

306 (de, self.energy_estimate)) 

307 if abs(self.energy_estimate) > self.eps: 

308 f = abs((de / self.energy_estimate) - 1) 

309 self.write_log('Energy prediction factor ' + str(f)) 

310 # fg = self.get_force_prediction(G) 

311 self.radius *= scale_radius_energy(f, self.radius) 

312 

313 else: 

314 self.write_log("energy change; actual: %f " % (de)) 

315 self.radius *= 1.5 

316 

317 fg = self.get_force_prediction(G) 

318 self.write_log("Scale factors %f %f " % 

319 (scale_radius_energy(f, self.radius), 

320 scale_radius_force(fg, self.radius))) 

321 

322 self.radius = max(min(self.radius, self.maxradius), 0.0001) 

323 else: 

324 self.update_hessian(pos, G) 

325 

326 self.write_log("new radius %f " % (self.radius)) 

327 self.oldenergy = energy 

328 

329 b, V = eigh(self.hessian) 

330 V = V.T.copy() 

331 self.V = V 

332 

333 # calculate projection of G onto eigenvectors V 

334 Gbar = np.dot(G, np.transpose(V)) 

335 

336 lamdas = self.get_lambdas(b, Gbar) 

337 

338 D = -Gbar / (b - lamdas) 

339 n = len(D) 

340 step = np.zeros((n)) 

341 for i in range(n): 

342 step += D[i] * V[i] 

343 

344 pos = self.atoms.get_positions().ravel() 

345 pos += step 

346 

347 energy_estimate = self.get_energy_estimate(D, Gbar, b) 

348 self.energy_estimate = energy_estimate 

349 self.gbar_estimate = self.get_gbar_estimate(D, Gbar, b) 

350 self.old_gbar = Gbar 

351 

352 self.atoms.set_positions(pos.reshape((-1, 3))) 

353 

354 def get_energy_estimate(self, D, Gbar, b): 

355 

356 de = 0.0 

357 for n in range(len(D)): 

358 de += D[n] * Gbar[n] + 0.5 * D[n] * b[n] * D[n] 

359 return de 

360 

361 def get_gbar_estimate(self, D, Gbar, b): 

362 gbar_est = (D * b) + Gbar 

363 self.write_log('Abs Gbar estimate ' + str(np.dot(gbar_est, gbar_est))) 

364 return gbar_est 

365 

366 def get_lambdas(self, b, Gbar): 

367 lamdas = np.zeros((len(b))) 

368 

369 D = -Gbar / b 

370 absD = np.sqrt(np.dot(D, D)) 

371 

372 eps = 1e-12 

373 nminus = self.get_hessian_inertia(b) 

374 

375 if absD < self.radius: 

376 if not self.transitionstate: 

377 self.write_log('Newton step') 

378 return lamdas 

379 else: 

380 if nminus == 1: 

381 self.write_log('Newton step') 

382 return lamdas 

383 else: 

384 self.write_log( 

385 "Wrong inertia of Hessian matrix: %2.2f %2.2f " % 

386 (b[0], b[1])) 

387 

388 else: 

389 self.write_log("Corrected Newton step: abs(D) = %2.2f " % (absD)) 

390 

391 if not self.transitionstate: 

392 # upper limit 

393 upperlimit = min(0, b[0]) - eps 

394 lamda = find_lamda(upperlimit, Gbar, b, self.radius) 

395 lamdas += lamda 

396 else: 

397 # upperlimit 

398 upperlimit = min(-b[0], b[1], 0) - eps 

399 lamda = find_lamda(upperlimit, Gbar, b, self.radius) 

400 lamdas += lamda 

401 lamdas[0] -= 2 * lamda 

402 

403 return lamdas 

404 

405 def get_hessian_inertia(self, eigenvalues): 

406 # return number of negative modes 

407 self.write_log("eigenvalues %2.2f %2.2f %2.2f " % (eigenvalues[0], 

408 eigenvalues[1], 

409 eigenvalues[2])) 

410 n = 0 

411 while eigenvalues[n] < 0: 

412 n += 1 

413 return n 

414 

415 def get_force_prediction(self, G): 

416 # return measure of how well the forces are predicted 

417 Gbar = np.dot(G, np.transpose(self.V)) 

418 dGbar_actual = Gbar - self.old_gbar 

419 dGbar_predicted = Gbar - self.gbar_estimate 

420 

421 f = np.dot(dGbar_actual, dGbar_predicted) / \ 

422 np.dot(dGbar_actual, dGbar_actual) 

423 self.write_log('Force prediction factor ' + str(f)) 

424 return f