r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

4"""

5Quasi-Newton algorithm

6"""

8__docformat__ = 'reStructuredText'

10import time

12import numpy as np

13from numpy.linalg import eigh

15from ase.optimize.optimize import Optimizer

19 b1 = b - lamda

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

21 return g

25 scale = 1.0

26# if(r<=0.01):

27# return scale

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

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

45 return scale

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

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

69 return scale

73 lowerlimit = upperlimit

74 step = 0.1

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

76 lowerlimit -= step

78 converged = False

80 while not converged:

82 midt = (upperlimit + lowerlimit) / 2.

83 lamda = midt

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

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

87 if fupper * fmidt < 0:

88 lowerlimit = midt

89 else:

90 upperlimit = midt

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

93 converged = True

95 return lamda

98class GoodOldQuasiNewton(Optimizer):

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

101 fmax=None, converged=None,

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

105 transitionstate=False, master=None):

106 """Parameters:

108 atoms: Atoms object

109 The Atoms object to relax.

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.

116 trajectory: string

117 File used to store trajectory of atomic movement.

119 maxstep: float

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

121 iteration (default value is 0.2 Angstroms).

124 logfile: file object or str

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

126 Use '-' for stdout.

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 """

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

135 self.eps = 1e-12

136 self.hessianupdate = hessianupdate

137 self.forcemin = forcemin

138 self.verbosity = verbosity

139 self.diagonal = diagonal

141 self.atoms = atoms

143 n = len(self.atoms) * 3

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

146 else:

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

151 else:

157 self.transitionstate = transitionstate

159 # check if this is a nudged elastic band calculation

160 if hasattr(atoms, 'springconstant'):

161 self.forcemin = False

163 self.t0 = time.time()

165 def initialize(self):

166 pass

168 def write_log(self, text):

169 if self.logfile is not None:

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

171 self.logfile.flush()

173 def set_hessian(self, hessian):

174 self.hessian = hessian

176 def get_hessian(self):

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

178 self.set_default_hessian()

179 return self.hessian

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)

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()

202 self.oldpos = copy.copy(pos)

203 self.oldG = copy.copy(G)

205 if self.verbosity:

206 print('hessian ', self.hessian)

208 def update_hessian_bfgs(self, pos, G):

209 n = len(self.hessian)

210 dgrad = G - self.oldG

211 dpos = pos - self.oldpos

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

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

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

231 tvecdpos = np.dot(tvec, dpos)

232 ddot = tvecdpos / absdpos

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

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

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

253 tvecdot = np.dot(tvec, tvec)

254 tvecdpos = np.dot(tvec, dpos)

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

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

258 coef3 = coef1 * tvecdpos / absdpos

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

269 def step(self, forces=None):

270 """ Do one QN step

271 """

273 if forces is None:

274 forces = self.atoms.get_forces()

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

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

278 energy = self.atoms.get_potential_energy()

280 if hasattr(self, 'oldenergy'):

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

283 ' ' + str(self.oldenergy))

285 if self.forcemin:

286 de = 1e-4

287 else:

288 de = 1e-2

290 if self.transitionstate:

291 de = 0.2

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

299 else:

300 self.update_hessian(pos, G)

301 de = energy - self.oldenergy

302 forces = 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 forces = abs((de / self.energy_estimate) - 1)

309 self.write_log('Energy prediction factor '

310 + str(forces))

311 # fg = self.get_force_prediction(G)

314 else:

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

318 fg = self.get_force_prediction(G)

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

324 else:

325 self.update_hessian(pos, G)

328 self.oldenergy = energy

330 b, V = eigh(self.hessian)

331 V = V.T.copy()

332 self.V = V

334 # calculate projection of G onto eigenvectors V

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

337 lamdas = self.get_lambdas(b, Gbar)

339 D = -Gbar / (b - lamdas)

340 n = len(D)

341 step = np.zeros((n))

342 for i in range(n):

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

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

346 pos += step

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

349 self.energy_estimate = energy_estimate

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

351 self.old_gbar = Gbar

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

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

357 de = 0.0

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

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

360 return de

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

363 gbar_est = (D * b) + Gbar

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

365 return gbar_est

367 def get_lambdas(self, b, Gbar):

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

370 D = -Gbar / b

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

373 eps = 1e-12

374 nminus = self.get_hessian_inertia(b)

377 if not self.transitionstate:

378 self.write_log('Newton step')

379 return lamdas

380 else:

381 if nminus == 1:

382 self.write_log('Newton step')

383 return lamdas

384 else:

385 self.write_log(

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

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

389 else:

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

392 if not self.transitionstate:

393 # upper limit

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

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

396 lamdas += lamda

397 else:

398 # upperlimit

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

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

401 lamdas += lamda

402 lamdas[0] -= 2 * lamda

404 return lamdas

406 def get_hessian_inertia(self, eigenvalues):

407 # return number of negative modes

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

409 eigenvalues[1],

410 eigenvalues[2]))

411 n = 0

412 while eigenvalues[n] < 0:

413 n += 1

414 return n

416 def get_force_prediction(self, G):

417 # return measure of how well the forces are predicted

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

419 dGbar_actual = Gbar - self.old_gbar

420 dGbar_predicted = Gbar - self.gbar_estimate

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

423 np.dot(dGbar_actual, dGbar_actual)

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

425 return f