Coverage for /builds/debichem-team/python-ase/ase/optimize/oldqn.py: 72.18%

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 **kwargs, 

119 ): 

120 """ 

121 

122 Parameters 

123 ---------- 

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

125 The Atoms object to relax. 

126 

127 restart: str 

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

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

130 be used, if the file exists. 

131 

132 trajectory: str 

133 File used to store trajectory of atomic movement. 

134 

135 maxstep: float 

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

137 iteration (default value is 0.2 Angstroms). 

138 

139 

140 logfile: file object or str 

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

142 Use '-' for stdout. 

143 

144 kwargs : dict, optional 

145 Extra arguments passed to 

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

147 

148 """ 

149 

150 Optimizer.__init__(self, atoms, restart, logfile, trajectory, **kwargs) 

151 

152 self.eps = 1e-12 

153 self.hessianupdate = hessianupdate 

154 self.forcemin = forcemin 

155 self.verbosity = verbosity 

156 self.diagonal = diagonal 

157 

158 n = len(self.optimizable) * 3 

159 if radius is None: 

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

161 else: 

162 self.radius = radius 

163 

164 if maxradius is None: 

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

166 else: 

167 self.maxradius = maxradius 

168 

169 # 0.01 < radius < maxradius 

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

171 

172 self.transitionstate = transitionstate 

173 

174 if self.optimizable.is_neb(): 

175 self.forcemin = False 

176 

177 self.t0 = time.time() 

178 

179 def initialize(self): 

180 pass 

181 

182 def write_log(self, text): 

183 if self.logfile is not None: 

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

185 self.logfile.flush() 

186 

187 def set_hessian(self, hessian): 

188 self.hessian = hessian 

189 

190 def get_hessian(self): 

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

192 self.set_default_hessian() 

193 return self.hessian 

194 

195 def set_default_hessian(self): 

196 # set unit matrix 

197 n = len(self.optimizable) * 3 

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

199 for i in range(n): 

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

201 self.set_hessian(hessian) 

202 

203 def update_hessian(self, pos, G): 

204 import copy 

205 if hasattr(self, 'oldG'): 

206 if self.hessianupdate == 'BFGS': 

207 self.update_hessian_bfgs(pos, G) 

208 elif self.hessianupdate == 'Powell': 

209 self.update_hessian_powell(pos, G) 

210 else: 

211 self.update_hessian_bofill(pos, G) 

212 else: 

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

214 self.set_default_hessian() 

215 

216 self.oldpos = copy.copy(pos) 

217 self.oldG = copy.copy(G) 

218 

219 if self.verbosity: 

220 print('hessian ', self.hessian) 

221 

222 def update_hessian_bfgs(self, pos, G): 

223 n = len(self.hessian) 

224 dgrad = G - self.oldG 

225 dpos = pos - self.oldpos 

226 dotg = np.dot(dgrad, dpos) 

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

228 dott = np.dot(dpos, tvec) 

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

230 for i in range(n): 

231 for j in range(n): 

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

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

234 

235 def update_hessian_powell(self, pos, G): 

236 n = len(self.hessian) 

237 dgrad = G - self.oldG 

238 dpos = pos - self.oldpos 

239 absdpos = np.dot(dpos, dpos) 

240 if absdpos < self.eps: 

241 return 

242 

243 dotg = np.dot(dgrad, dpos) 

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

245 tvecdpos = np.dot(tvec, dpos) 

246 ddot = tvecdpos / absdpos 

247 

248 dott = np.dot(dpos, tvec) 

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

250 for i in range(n): 

251 for j in range(n): 

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

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

254 h *= 1. / absdpos 

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

256 

257 def update_hessian_bofill(self, pos, G): 

258 print('update Bofill') 

259 n = len(self.hessian) 

260 dgrad = G - self.oldG 

261 dpos = pos - self.oldpos 

262 absdpos = np.dot(dpos, dpos) 

263 if absdpos < self.eps: 

264 return 

265 dotg = np.dot(dgrad, dpos) 

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

267 tvecdot = np.dot(tvec, tvec) 

268 tvecdpos = np.dot(tvec, dpos) 

269 

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

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

272 coef3 = coef1 * tvecdpos / absdpos 

273 

274 dott = np.dot(dpos, tvec) 

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

276 for i in range(n): 

277 for j in range(n): 

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

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

280 h *= 1. / absdpos 

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

282 

283 def step(self, forces=None): 

284 """ Do one QN step 

285 """ 

286 

287 if forces is None: 

288 forces = self.optimizable.get_forces() 

289 

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

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

292 

293 energy = self.optimizable.get_potential_energy() 

294 

295 if hasattr(self, 'oldenergy'): 

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

297 ' ' + str(self.oldenergy)) 

298 

299 if self.forcemin: 

300 de = 1e-4 

301 else: 

302 de = 1e-2 

303 

304 if self.transitionstate: 

305 de = 0.2 

306 

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

308 self.write_log('reject step') 

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

310 G = self.oldG 

311 energy = self.oldenergy 

312 self.radius *= 0.5 

313 else: 

314 self.update_hessian(pos, G) 

315 de = energy - self.oldenergy 

316 forces = 1.0 

317 if self.forcemin: 

318 self.write_log( 

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

320 (de, self.energy_estimate)) 

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

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

323 self.write_log('Energy prediction factor ' 

324 + str(forces)) 

325 # fg = self.get_force_prediction(G) 

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

327 

328 else: 

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

330 self.radius *= 1.5 

331 

332 fg = self.get_force_prediction(G) 

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

334 (scale_radius_energy(forces, self.radius), 

335 scale_radius_force(fg, self.radius))) 

336 

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

338 else: 

339 self.update_hessian(pos, G) 

340 

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

342 self.oldenergy = energy 

343 

344 b, V = eigh(self.hessian) 

345 V = V.T.copy() 

346 self.V = V 

347 

348 # calculate projection of G onto eigenvectors V 

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

350 

351 lamdas = self.get_lambdas(b, Gbar) 

352 

353 D = -Gbar / (b - lamdas) 

354 n = len(D) 

355 step = np.zeros(n) 

356 for i in range(n): 

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

358 

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

360 pos += step 

361 

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

363 self.energy_estimate = energy_estimate 

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

365 self.old_gbar = Gbar 

366 

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

368 

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

370 

371 de = 0.0 

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

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

374 return de 

375 

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

377 gbar_est = (D * b) + Gbar 

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

379 return gbar_est 

380 

381 def get_lambdas(self, b, Gbar): 

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

383 

384 D = -Gbar / b 

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

386 

387 eps = 1e-12 

388 nminus = self.get_hessian_inertia(b) 

389 

390 if absD < self.radius: 

391 if not self.transitionstate: 

392 self.write_log('Newton step') 

393 return lamdas 

394 else: 

395 if nminus == 1: 

396 self.write_log('Newton step') 

397 return lamdas 

398 else: 

399 self.write_log( 

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

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

402 

403 else: 

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

405 

406 if not self.transitionstate: 

407 # upper limit 

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

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

410 lamdas += lamda 

411 else: 

412 # upperlimit 

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

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

415 lamdas += lamda 

416 lamdas[0] -= 2 * lamda 

417 

418 return lamdas 

419 

420 def get_hessian_inertia(self, eigenvalues): 

421 # return number of negative modes 

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

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

424 n = 0 

425 while eigenvalues[n] < 0: 

426 n += 1 

427 return n 

428 

429 def get_force_prediction(self, G): 

430 # return measure of how well the forces are predicted 

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

432 dGbar_actual = Gbar - self.old_gbar 

433 dGbar_predicted = Gbar - self.gbar_estimate 

434 

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

436 np.dot(dGbar_actual, dGbar_actual) 

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

438 return f