Coverage for /builds/ase/ase/ase/optimize/oldqn.py: 72.96%

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 Optimizer.__init__(self, 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 self.logfile.flush() 

193 

194 def set_hessian(self, hessian): 

195 self.hessian = hessian 

196 

197 def get_hessian(self): 

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

199 self.set_default_hessian() 

200 return self.hessian 

201 

202 def set_default_hessian(self): 

203 # set unit matrix 

204 n = self.optimizable.ndofs() 

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

206 for i in range(n): 

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

208 self.set_hessian(hessian) 

209 

210 def update_hessian(self, pos, G): 

211 import copy 

212 if hasattr(self, 'oldG'): 

213 if self.hessianupdate == 'BFGS': 

214 self.update_hessian_bfgs(pos, G) 

215 elif self.hessianupdate == 'Powell': 

216 self.update_hessian_powell(pos, G) 

217 else: 

218 self.update_hessian_bofill(pos, G) 

219 else: 

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

221 self.set_default_hessian() 

222 

223 self.oldpos = copy.copy(pos) 

224 self.oldG = copy.copy(G) 

225 

226 if self.verbosity: 

227 print('hessian ', self.hessian) 

228 

229 def update_hessian_bfgs(self, pos, G): 

230 n = len(self.hessian) 

231 dgrad = G - self.oldG 

232 dpos = pos - self.oldpos 

233 dotg = np.dot(dgrad, dpos) 

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

235 dott = np.dot(dpos, tvec) 

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

237 for i in range(n): 

238 for j in range(n): 

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

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

241 

242 def update_hessian_powell(self, pos, G): 

243 n = len(self.hessian) 

244 dgrad = G - self.oldG 

245 dpos = pos - self.oldpos 

246 absdpos = np.dot(dpos, dpos) 

247 if absdpos < self.eps: 

248 return 

249 

250 dotg = np.dot(dgrad, dpos) 

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

252 tvecdpos = np.dot(tvec, dpos) 

253 ddot = tvecdpos / absdpos 

254 

255 dott = np.dot(dpos, tvec) 

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

257 for i in range(n): 

258 for j in range(n): 

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

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

261 h *= 1. / absdpos 

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

263 

264 def update_hessian_bofill(self, pos, G): 

265 print('update Bofill') 

266 n = len(self.hessian) 

267 dgrad = G - self.oldG 

268 dpos = pos - self.oldpos 

269 absdpos = np.dot(dpos, dpos) 

270 if absdpos < self.eps: 

271 return 

272 dotg = np.dot(dgrad, dpos) 

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

274 tvecdot = np.dot(tvec, tvec) 

275 tvecdpos = np.dot(tvec, dpos) 

276 

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

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

279 coef3 = coef1 * tvecdpos / absdpos 

280 

281 dott = np.dot(dpos, tvec) 

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

283 for i in range(n): 

284 for j in range(n): 

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

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

287 h *= 1. / absdpos 

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

289 

290 def step(self, forces=None): 

291 """ Do one QN step 

292 """ 

293 

294 if forces is None: 

295 forces = self.optimizable.get_gradient().reshape(-1, 3) 

296 

297 pos = self.optimizable.get_x() 

298 G = -self.optimizable.get_gradient() 

299 

300 energy = self.optimizable.get_value() 

301 

302 if hasattr(self, 'oldenergy'): 

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

304 ' ' + str(self.oldenergy)) 

305 

306 if self.forcemin: 

307 de = 1e-4 

308 else: 

309 de = 1e-2 

310 

311 if self.transitionstate: 

312 de = 0.2 

313 

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

315 self.write_log('reject step') 

316 self.optimizable.set_x(self.oldpos) 

317 G = self.oldG 

318 energy = self.oldenergy 

319 self.radius *= 0.5 

320 else: 

321 self.update_hessian(pos, G) 

322 de = energy - self.oldenergy 

323 forces = 1.0 

324 if self.forcemin: 

325 self.write_log( 

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

327 (de, self.energy_estimate)) 

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

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

330 self.write_log('Energy prediction factor ' 

331 + str(forces)) 

332 # fg = self.get_force_prediction(G) 

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

334 

335 else: 

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

337 self.radius *= 1.5 

338 

339 fg = self.get_force_prediction(G) 

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

341 (scale_radius_energy(forces, self.radius), 

342 scale_radius_force(fg, self.radius))) 

343 

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

345 else: 

346 self.update_hessian(pos, G) 

347 

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

349 self.oldenergy = energy 

350 

351 b, V = eigh(self.hessian) 

352 V = V.T.copy() 

353 self.V = V 

354 

355 # calculate projection of G onto eigenvectors V 

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

357 

358 lamdas = self.get_lambdas(b, Gbar) 

359 

360 D = -Gbar / (b - lamdas) 

361 n = len(D) 

362 step = np.zeros(n) 

363 for i in range(n): 

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

365 

366 pos = self.optimizable.get_x() 

367 pos += step 

368 

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

370 self.energy_estimate = energy_estimate 

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

372 self.old_gbar = Gbar 

373 

374 self.optimizable.set_x(pos) 

375 

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

377 

378 de = 0.0 

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

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

381 return de 

382 

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

384 gbar_est = (D * b) + Gbar 

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

386 return gbar_est 

387 

388 def get_lambdas(self, b, Gbar): 

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

390 

391 D = -Gbar / b 

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

393 

394 eps = 1e-12 

395 nminus = self.get_hessian_inertia(b) 

396 

397 if absD < self.radius: 

398 if not self.transitionstate: 

399 self.write_log('Newton step') 

400 return lamdas 

401 else: 

402 if nminus == 1: 

403 self.write_log('Newton step') 

404 return lamdas 

405 else: 

406 self.write_log( 

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

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

409 

410 else: 

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

412 

413 if not self.transitionstate: 

414 # upper limit 

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

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

417 lamdas += lamda 

418 else: 

419 # upperlimit 

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

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

422 lamdas += lamda 

423 lamdas[0] -= 2 * lamda 

424 

425 return lamdas 

426 

427 def get_hessian_inertia(self, eigenvalues): 

428 # return number of negative modes 

429 self.write_log("eigenvalues {:2.2f} {:2.2f} {:2.2f} ".format( 

430 eigenvalues[0], eigenvalues[1], eigenvalues[2])) 

431 n = 0 

432 while eigenvalues[n] < 0: 

433 n += 1 

434 return n 

435 

436 def get_force_prediction(self, G): 

437 # return measure of how well the forces are predicted 

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

439 dGbar_actual = Gbar - self.old_gbar 

440 dGbar_predicted = Gbar - self.gbar_estimate 

441 

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

443 np.dot(dGbar_actual, dGbar_actual) 

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

445 return f