Coverage for ase / optimize / oldqn.py: 73.88%

1# fmt: off 

2 

3# Copyright (C) 2003 CAMP 

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

5 

6""" 

7Quasi-Newton algorithm 

8""" 

9 

10__docformat__ = 'reStructuredText' 

11 

12import time 

13from typing import IO, Optional, Union 

14 

15import numpy as np 

16from numpy.linalg import eigh 

17 

18from ase import Atoms 

19from ase.optimize.optimize import Optimizer 

20 

21 

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

23 b1 = b - lamda 

24 b1[abs(b1) < 1e-40] = 1e-40 # avoid divide-by-zero 

25 gbar_b_lamda = Gbar / b1 # only compute once 

26 g = radius**2 - np.dot(gbar_b_lamda, gbar_b_lamda) 

27 return g 

28 

29 

30def scale_radius_energy(f, r): 

31 scale = 1.0 

32# if(r<=0.01): 

33# return scale 

34 

35 if f < 0.01: 

36 scale *= 1.4 

37 if f < 0.05: 

38 scale *= 1.4 

39 if f < 0.10: 

40 scale *= 1.4 

41 if f < 0.40: 

42 scale *= 1.4 

43 

44 if f > 0.5: 

45 scale *= 1. / 1.4 

46 if f > 0.7: 

47 scale *= 1. / 1.4 

48 if f > 1.0: 

49 scale *= 1. / 1.4 

50 

51 return scale 

52 

53 

54def scale_radius_force(f, r): 

55 scale = 1.0 

56# if(r<=0.01): 

57# return scale 

58 g = abs(f - 1) 

59 if g < 0.01: 

60 scale *= 1.4 

61 if g < 0.05: 

62 scale *= 1.4 

63 if g < 0.10: 

64 scale *= 1.4 

65 if g < 0.40: 

66 scale *= 1.4 

67 

68 if g > 0.5: 

69 scale *= 1. / 1.4 

70 if g > 0.7: 

71 scale *= 1. / 1.4 

72 if g > 1.0: 

73 scale *= 1. / 1.4 

74 

75 return scale 

76 

77 

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

79 lowerlimit = upperlimit 

80 step = 0.1 

81 i = 0 

82 

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

84 lowerlimit -= step 

85 i += 1 

86 

87 # if many iterations are required, dynamically scale step for efficiency 

88 if i % 100 == 0: 

89 step *= 1.25 

90 

91 converged = False 

92 

93 while not converged: 

94 

95 midt = (upperlimit + lowerlimit) / 2. 

96 lamda = midt 

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

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

99 

100 if fupper * fmidt < 0: 

101 lowerlimit = midt 

102 else: 

103 upperlimit = midt 

104 

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

106 converged = True 

107 

108 return lamda 

109 

110 

111class GoodOldQuasiNewton(Optimizer): 

112 

113 def __init__( 

114 self, 

115 atoms: Atoms, 

116 restart: Optional[str] = None, 

117 logfile: Union[IO, str] = '-', 

118 trajectory: Optional[str] = None, 

119 fmax=None, 

120 converged=None, 

121 hessianupdate: str = 'BFGS', 

122 hessian=None, 

123 forcemin: bool = True, 

124 verbosity: bool = False, 

125 maxradius: Optional[float] = None, 

126 diagonal: float = 20.0, 

127 radius: Optional[float] = None, 

128 transitionstate: bool = False, 

129 **kwargs, 

130 ): 

131 """ 

132 

133 Parameters 

134 ---------- 

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

136 The Atoms object to relax. 

137 

138 restart: str 

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

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

141 be used, if the file exists. 

142 

143 trajectory: str 

144 File used to store trajectory of atomic movement. 

145 

146 maxstep: float 

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

148 iteration (default value is 0.2 Angstroms). 

149 

150 

151 logfile: file object or str 

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

153 Use '-' for stdout. 

154 

155 kwargs : dict, optional 

156 Extra arguments passed to 

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

158 

159 """ 

160 

161 super().__init__(atoms, restart, logfile, trajectory, **kwargs) 

162 

163 self.eps = 1e-12 

164 self.hessianupdate = hessianupdate 

165 self.forcemin = forcemin 

166 self.verbosity = verbosity 

167 self.diagonal = diagonal 

168 

169 n = self.optimizable.ndofs() 

170 if radius is None: 

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

172 else: 

173 self.radius = radius 

174 

175 if maxradius is None: 

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

177 else: 

178 self.maxradius = maxradius 

179 

180 # 0.01 < radius < maxradius 

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

182 

183 self.transitionstate = transitionstate 

184 self.t0 = time.time() 

185 

186 def initialize(self): 

187 pass 

188 

189 def write_log(self, text): 

190 if self.logfile is not None: 

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

192 

193 def set_hessian(self, hessian): 

194 self.hessian = hessian 

195 

196 def get_hessian(self): 

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

198 self.set_default_hessian() 

199 return self.hessian 

200 

201 def set_default_hessian(self): 

202 # set unit matrix 

203 n = self.optimizable.ndofs() 

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

205 for i in range(n): 

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

207 self.set_hessian(hessian) 

208 

209 def update_hessian(self, pos, G): 

210 import copy 

211 if hasattr(self, 'oldG'): 

212 if self.hessianupdate == 'BFGS': 

213 self.update_hessian_bfgs(pos, G) 

214 elif self.hessianupdate == 'Powell': 

215 self.update_hessian_powell(pos, G) 

216 else: 

217 self.update_hessian_bofill(pos, G) 

218 else: 

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

220 self.set_default_hessian() 

221 

222 self.oldpos = copy.copy(pos) 

223 self.oldG = copy.copy(G) 

224 

225 if self.verbosity: 

226 print('hessian ', self.hessian) 

227 

228 def update_hessian_bfgs(self, pos, G): 

229 n = len(self.hessian) 

230 dgrad = G - self.oldG 

231 dpos = pos - self.oldpos 

232 dotg = np.dot(dgrad, dpos) 

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

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 = dgrad[i] * dgrad[j] / dotg - tvec[i] * tvec[j] / dott 

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

240 

241 def update_hessian_powell(self, pos, G): 

242 n = len(self.hessian) 

243 dgrad = G - self.oldG 

244 dpos = pos - self.oldpos 

245 absdpos = np.dot(dpos, dpos) 

246 if absdpos < self.eps: 

247 return 

248 

249 dotg = np.dot(dgrad, dpos) 

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

251 tvecdpos = np.dot(tvec, dpos) 

252 ddot = tvecdpos / absdpos 

253 

254 dott = np.dot(dpos, tvec) 

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

256 for i in range(n): 

257 for j in range(n): 

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

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

260 h *= 1. / absdpos 

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

262 

263 def update_hessian_bofill(self, pos, G): 

264 print('update Bofill') 

265 n = len(self.hessian) 

266 dgrad = G - self.oldG 

267 dpos = pos - self.oldpos 

268 absdpos = np.dot(dpos, dpos) 

269 if absdpos < self.eps: 

270 return 

271 dotg = np.dot(dgrad, dpos) 

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

273 tvecdot = np.dot(tvec, tvec) 

274 tvecdpos = np.dot(tvec, dpos) 

275 

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

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

278 coef3 = coef1 * tvecdpos / absdpos 

279 

280 dott = np.dot(dpos, tvec) 

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

282 for i in range(n): 

283 for j in range(n): 

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

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

286 h *= 1. / absdpos 

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

288 

289 def step(self, forces=None): 

290 """ Do one QN step 

291 """ 

292 G = self._get_gradient(forces) 

293 pos = self.optimizable.get_x() 

294 energy = self.optimizable.get_value() 

295 

296 if hasattr(self, 'oldenergy'): 

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

298 ' ' + str(self.oldenergy)) 

299 

300 if self.forcemin: 

301 de = 1e-4 

302 else: 

303 de = 1e-2 

304 

305 if self.transitionstate: 

306 de = 0.2 

307 

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

309 self.write_log('reject step') 

310 self.optimizable.set_x(self.oldpos) 

311 G = self.oldG 

312 energy = self.oldenergy 

313 self.radius *= 0.5 

314 else: 

315 self.update_hessian(pos, G) 

316 de = energy - self.oldenergy 

317 forces = 1.0 

318 if self.forcemin: 

319 self.write_log( 

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

321 (de, self.energy_estimate)) 

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

323 forces = abs((de / self.energy_estimate) - 1) 

324 self.write_log('Energy prediction factor ' 

325 + str(forces)) 

326 # fg = self.get_force_prediction(G) 

327 self.radius *= scale_radius_energy(forces, self.radius) 

328 

329 else: 

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

331 self.radius *= 1.5 

332 

333 fg = self.get_force_prediction(G) 

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

335 (scale_radius_energy(forces, self.radius), 

336 scale_radius_force(fg, self.radius))) 

337 

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

339 else: 

340 self.update_hessian(pos, G) 

341 

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

343 self.oldenergy = energy 

344 

345 b, V = eigh(self.hessian) 

346 V = V.T.copy() 

347 self.V = V 

348 

349 # calculate projection of G onto eigenvectors V 

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

351 

352 lamdas = self.get_lambdas(b, Gbar) 

353 

354 D = -Gbar / (b - lamdas) 

355 n = len(D) 

356 step = np.zeros(n) 

357 for i in range(n): 

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

359 

360 pos = self.optimizable.get_x() 

361 pos += step 

362 

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

364 self.energy_estimate = energy_estimate 

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

366 self.old_gbar = Gbar 

367 

368 self.optimizable.set_x(pos) 

369 

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

371 

372 de = 0.0 

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

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

375 return de 

376 

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

378 gbar_est = (D * b) + Gbar 

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

380 return gbar_est 

381 

382 def get_lambdas(self, b, Gbar): 

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

384 

385 D = -Gbar / b 

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

387 

388 eps = 1e-12 

389 nminus = self.get_hessian_inertia(b) 

390 

391 if absD < self.radius: 

392 if not self.transitionstate: 

393 self.write_log('Newton step') 

394 return lamdas 

395 else: 

396 if nminus == 1: 

397 self.write_log('Newton step') 

398 return lamdas 

399 else: 

400 self.write_log( 

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

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

403 

404 else: 

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

406 

407 if not self.transitionstate: 

408 # upper limit 

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

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

411 lamdas += lamda 

412 else: 

413 # upperlimit 

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

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

416 lamdas += lamda 

417 lamdas[0] -= 2 * lamda 

418 

419 return lamdas 

420 

421 def get_hessian_inertia(self, eigenvalues): 

422 # return number of negative modes 

423 txt = ' '.join(f'{eig:2.2f}' for eig in eigenvalues[:3]) 

424 self.write_log(f'eigenvalues {txt}') 

425 n = 0 

426 while eigenvalues[n] < 0: 

427 n += 1 

428 return n 

429 

430 def get_force_prediction(self, G): 

431 # return measure of how well the forces are predicted 

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

433 dGbar_actual = Gbar - self.old_gbar 

434 dGbar_predicted = Gbar - self.gbar_estimate 

435 

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

437 np.dot(dGbar_actual, dGbar_actual) 

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

439 return f