Coverage for /builds/kinetik161/ase/ase/optimize/oldqn.py: 72.28%

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 

11from typing import IO, Optional, Union 

12 

13import numpy as np 

14from numpy.linalg import eigh 

15 

16from ase import Atoms 

17from ase.optimize.optimize import Optimizer 

18 

19 

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

21 b1 = b - lamda 

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

23 return g 

24 

25 

26def scale_radius_energy(f, r): 

27 scale = 1.0 

28# if(r<=0.01): 

29# return scale 

30 

31 if f < 0.01: 

32 scale *= 1.4 

33 if f < 0.05: 

34 scale *= 1.4 

35 if f < 0.10: 

36 scale *= 1.4 

37 if f < 0.40: 

38 scale *= 1.4 

39 

40 if f > 0.5: 

41 scale *= 1. / 1.4 

42 if f > 0.7: 

43 scale *= 1. / 1.4 

44 if f > 1.0: 

45 scale *= 1. / 1.4 

46 

47 return scale 

48 

49 

50def scale_radius_force(f, r): 

51 scale = 1.0 

52# if(r<=0.01): 

53# return scale 

54 g = abs(f - 1) 

55 if g < 0.01: 

56 scale *= 1.4 

57 if g < 0.05: 

58 scale *= 1.4 

59 if g < 0.10: 

60 scale *= 1.4 

61 if g < 0.40: 

62 scale *= 1.4 

63 

64 if g > 0.5: 

65 scale *= 1. / 1.4 

66 if g > 0.7: 

67 scale *= 1. / 1.4 

68 if g > 1.0: 

69 scale *= 1. / 1.4 

70 

71 return scale 

72 

73 

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

75 lowerlimit = upperlimit 

76 step = 0.1 

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

78 lowerlimit -= step 

79 

80 converged = False 

81 

82 while not converged: 

83 

84 midt = (upperlimit + lowerlimit) / 2. 

85 lamda = midt 

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

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

88 

89 if fupper * fmidt < 0: 

90 lowerlimit = midt 

91 else: 

92 upperlimit = midt 

93 

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

95 converged = True 

96 

97 return lamda 

98 

99 

100class GoodOldQuasiNewton(Optimizer): 

101 

102 def __init__( 

103 self, 

104 atoms: Atoms, 

105 restart: Optional[str] = None, 

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

107 trajectory: Optional[str] = None, 

108 fmax=None, 

109 converged=None, 

110 hessianupdate: str = 'BFGS', 

111 hessian=None, 

112 forcemin: bool = True, 

113 verbosity: bool = False, 

114 maxradius: Optional[float] = None, 

115 diagonal: float = 20.0, 

116 radius: Optional[float] = None, 

117 transitionstate: bool = False, 

118 master: Optional[bool] = None, 

119 ): 

120 """Parameters: 

121 

122 atoms: Atoms object 

123 The Atoms object to relax. 

124 

125 restart: string 

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

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

128 be used, if the file exists. 

129 

130 trajectory: string 

131 File used to store trajectory of atomic movement. 

132 

133 maxstep: float 

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

135 iteration (default value is 0.2 Angstroms). 

136 

137 

138 logfile: file object or str 

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

140 Use '-' for stdout. 

141 

142 master: boolean 

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

144 set to true, this rank will save files. 

145 """ 

146 

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

148 

149 self.eps = 1e-12 

150 self.hessianupdate = hessianupdate 

151 self.forcemin = forcemin 

152 self.verbosity = verbosity 

153 self.diagonal = diagonal 

154 

155 n = len(self.optimizable) * 3 

156 if radius is None: 

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

158 else: 

159 self.radius = radius 

160 

161 if maxradius is None: 

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

163 else: 

164 self.maxradius = maxradius 

165 

166 # 0.01 < radius < maxradius 

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

168 

169 self.transitionstate = transitionstate 

170 

171 if self.optimizable.is_neb(): 

172 self.forcemin = False 

173 

174 self.t0 = time.time() 

175 

176 def initialize(self): 

177 pass 

178 

179 def write_log(self, text): 

180 if self.logfile is not None: 

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

182 self.logfile.flush() 

183 

184 def set_hessian(self, hessian): 

185 self.hessian = hessian 

186 

187 def get_hessian(self): 

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

189 self.set_default_hessian() 

190 return self.hessian 

191 

192 def set_default_hessian(self): 

193 # set unit matrix 

194 n = len(self.optimizable) * 3 

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

196 for i in range(n): 

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

198 self.set_hessian(hessian) 

199 

200 def update_hessian(self, pos, G): 

201 import copy 

202 if hasattr(self, 'oldG'): 

203 if self.hessianupdate == 'BFGS': 

204 self.update_hessian_bfgs(pos, G) 

205 elif self.hessianupdate == 'Powell': 

206 self.update_hessian_powell(pos, G) 

207 else: 

208 self.update_hessian_bofill(pos, G) 

209 else: 

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

211 self.set_default_hessian() 

212 

213 self.oldpos = copy.copy(pos) 

214 self.oldG = copy.copy(G) 

215 

216 if self.verbosity: 

217 print('hessian ', self.hessian) 

218 

219 def update_hessian_bfgs(self, pos, G): 

220 n = len(self.hessian) 

221 dgrad = G - self.oldG 

222 dpos = pos - self.oldpos 

223 dotg = np.dot(dgrad, dpos) 

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

225 dott = np.dot(dpos, tvec) 

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

227 for i in range(n): 

228 for j in range(n): 

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

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

231 

232 def update_hessian_powell(self, pos, G): 

233 n = len(self.hessian) 

234 dgrad = G - self.oldG 

235 dpos = pos - self.oldpos 

236 absdpos = np.dot(dpos, dpos) 

237 if absdpos < self.eps: 

238 return 

239 

240 dotg = np.dot(dgrad, dpos) 

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

242 tvecdpos = np.dot(tvec, dpos) 

243 ddot = tvecdpos / absdpos 

244 

245 dott = np.dot(dpos, tvec) 

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

247 for i in range(n): 

248 for j in range(n): 

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

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

251 h *= 1. / absdpos 

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

253 

254 def update_hessian_bofill(self, pos, G): 

255 print('update Bofill') 

256 n = len(self.hessian) 

257 dgrad = G - self.oldG 

258 dpos = pos - self.oldpos 

259 absdpos = np.dot(dpos, dpos) 

260 if absdpos < self.eps: 

261 return 

262 dotg = np.dot(dgrad, dpos) 

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

264 tvecdot = np.dot(tvec, tvec) 

265 tvecdpos = np.dot(tvec, dpos) 

266 

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

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

269 coef3 = coef1 * tvecdpos / absdpos 

270 

271 dott = np.dot(dpos, tvec) 

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

273 for i in range(n): 

274 for j in range(n): 

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

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

277 h *= 1. / absdpos 

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

279 

280 def step(self, forces=None): 

281 """ Do one QN step 

282 """ 

283 

284 if forces is None: 

285 forces = self.optimizable.get_forces() 

286 

287 pos = self.optimizable.get_positions().ravel() 

288 G = -self.optimizable.get_forces().ravel() 

289 

290 energy = self.optimizable.get_potential_energy() 

291 

292 if hasattr(self, 'oldenergy'): 

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

294 ' ' + str(self.oldenergy)) 

295 

296 if self.forcemin: 

297 de = 1e-4 

298 else: 

299 de = 1e-2 

300 

301 if self.transitionstate: 

302 de = 0.2 

303 

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

305 self.write_log('reject step') 

306 self.optimizable.set_positions(self.oldpos.reshape((-1, 3))) 

307 G = self.oldG 

308 energy = self.oldenergy 

309 self.radius *= 0.5 

310 else: 

311 self.update_hessian(pos, G) 

312 de = energy - self.oldenergy 

313 forces = 1.0 

314 if self.forcemin: 

315 self.write_log( 

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

317 (de, self.energy_estimate)) 

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

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

320 self.write_log('Energy prediction factor ' 

321 + str(forces)) 

322 # fg = self.get_force_prediction(G) 

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

324 

325 else: 

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

327 self.radius *= 1.5 

328 

329 fg = self.get_force_prediction(G) 

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

331 (scale_radius_energy(forces, self.radius), 

332 scale_radius_force(fg, self.radius))) 

333 

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

335 else: 

336 self.update_hessian(pos, G) 

337 

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

339 self.oldenergy = energy 

340 

341 b, V = eigh(self.hessian) 

342 V = V.T.copy() 

343 self.V = V 

344 

345 # calculate projection of G onto eigenvectors V 

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

347 

348 lamdas = self.get_lambdas(b, Gbar) 

349 

350 D = -Gbar / (b - lamdas) 

351 n = len(D) 

352 step = np.zeros(n) 

353 for i in range(n): 

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

355 

356 pos = self.optimizable.get_positions().ravel() 

357 pos += step 

358 

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

360 self.energy_estimate = energy_estimate 

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

362 self.old_gbar = Gbar 

363 

364 self.optimizable.set_positions(pos.reshape((-1, 3))) 

365 

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

367 

368 de = 0.0 

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

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

371 return de 

372 

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

374 gbar_est = (D * b) + Gbar 

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

376 return gbar_est 

377 

378 def get_lambdas(self, b, Gbar): 

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

380 

381 D = -Gbar / b 

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

383 

384 eps = 1e-12 

385 nminus = self.get_hessian_inertia(b) 

386 

387 if absD < self.radius: 

388 if not self.transitionstate: 

389 self.write_log('Newton step') 

390 return lamdas 

391 else: 

392 if nminus == 1: 

393 self.write_log('Newton step') 

394 return lamdas 

395 else: 

396 self.write_log( 

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

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

399 

400 else: 

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

402 

403 if not self.transitionstate: 

404 # upper limit 

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

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

407 lamdas += lamda 

408 else: 

409 # upperlimit 

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

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

412 lamdas += lamda 

413 lamdas[0] -= 2 * lamda 

414 

415 return lamdas 

416 

417 def get_hessian_inertia(self, eigenvalues): 

418 # return number of negative modes 

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

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

421 n = 0 

422 while eigenvalues[n] < 0: 

423 n += 1 

424 return n 

425 

426 def get_force_prediction(self, G): 

427 # return measure of how well the forces are predicted 

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

429 dGbar_actual = Gbar - self.old_gbar 

430 dGbar_predicted = Gbar - self.gbar_estimate 

431 

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

433 np.dot(dGbar_actual, dGbar_actual) 

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

435 return f