Coverage for /builds/kinetik161/ase/ase/optimize/minimahopping.py: 57.57%

1import os

3import numpy as np

5from ase import io, units

6from ase.md import MDLogger, VelocityVerlet

7from ase.md.velocitydistribution import MaxwellBoltzmannDistribution

8from ase.optimize import QuasiNewton

9from ase.parallel import paropen, world

12class MinimaHopping:

13 """Implements the minima hopping method of global optimization outlined

14 by S. Goedecker, J. Chem. Phys. 120: 9911 (2004). Initialize with an

15 ASE atoms object. Optional parameters are fed through keywords.

16 To run multiple searches in parallel, specify the minima_traj keyword,

17 and have each run point to the same path.

18 """

20 _default_settings = {

21 'T0': 1000., # K, initial MD 'temperature'

22 'beta1': 1.1, # temperature adjustment parameter

23 'beta2': 1.1, # temperature adjustment parameter

24 'beta3': 1. / 1.1, # temperature adjustment parameter

25 'Ediff0': 0.5, # eV, initial energy acceptance threshold

26 'alpha1': 0.98, # energy threshold adjustment parameter

27 'alpha2': 1. / 0.98, # energy threshold adjustment parameter

28 'mdmin': 2, # criteria to stop MD simulation (no. of minima)

29 'logfile': 'hop.log', # text log

30 'minima_threshold': 0.5, # A, threshold for identical configs

31 'timestep': 1.0, # fs, timestep for MD simulations

32 'optimizer': QuasiNewton, # local optimizer to use

33 'minima_traj': 'minima.traj', # storage file for minima list

34 'fmax': 0.05} # eV/A, max force for optimizations

36 def __init__(self, atoms, **kwargs):

37 """Initialize with an ASE atoms object and keyword arguments."""

38 self._atoms = atoms

39 for key in kwargs:

40 if key not in self._default_settings:

41 raise RuntimeError(f'Unknown keyword: {key}')

42 for k, v in self._default_settings.items():

43 setattr(self, f'_{k}', kwargs.pop(k, v))

45 # when a MD sim. has passed a local minimum:

46 self._passedminimum = PassedMinimum()

48 # Misc storage.

49 self._previous_optimum = None

50 self._previous_energy = None

51 self._temperature = self._T0

52 self._Ediff = self._Ediff0

54 def __call__(self, totalsteps=None, maxtemp=None):

55 """Run the minima hopping algorithm. Can specify stopping criteria

56 with total steps allowed or maximum searching temperature allowed.

57 If neither is specified, runs indefinitely (or until stopped by

58 batching software)."""

59 self._startup()

60 while True:

61 if (totalsteps and self._counter >= totalsteps):

62 self._log('msg', 'Run terminated. Step #%i reached of '

63 '%i allowed. Increase totalsteps if resuming.'

64 % (self._counter, totalsteps))

65 return

66 if (maxtemp and self._temperature >= maxtemp):

67 self._log('msg', 'Run terminated. Temperature is %.2f K;'

68 ' max temperature allowed %.2f K.'

69 % (self._temperature, maxtemp))

70 return

72 self._previous_optimum = self._atoms.copy()

73 self._previous_energy = self._atoms.get_potential_energy()

74 self._molecular_dynamics()

75 self._optimize()

76 self._counter += 1

77 self._check_results()

79 def _startup(self):

80 """Initiates a run, and determines if running from previous data or

81 a fresh run."""

83 status = np.array(-1.)

84 exists = self._read_minima()

85 if world.rank == 0:

86 if not exists:

87 # Fresh run with new minima file.

88 status = np.array(0.)

89 elif not os.path.exists(self._logfile):

90 # Fresh run with existing or shared minima file.

91 status = np.array(1.)

92 else:

93 # Must be resuming from within a working directory.

94 status = np.array(2.)

95 world.barrier()

96 world.broadcast(status, 0)

98 if status == 2.:

99 self._resume()

100 else:

101 self._counter = 0

102 self._log('init')

103 self._log('msg', 'Performing initial optimization.')

104 if status == 1.:

105 self._log('msg', 'Using existing minima file with %i prior '

106 'minima: %s' % (len(self._minima),

107 self._minima_traj))

108 self._optimize()

109 self._check_results()

110 self._counter += 1

111

112 def _resume(self):

113 """Attempt to resume a run, based on information in the log

114 file. Note it will almost always be interrupted in the middle of

115 either a qn or md run or when exceeding totalsteps, so it only has

116 been tested in those cases currently."""

117 f = paropen(self._logfile, 'r')

118 lines = f.read().splitlines()

119 f.close()

120 self._log('msg', 'Attempting to resume stopped run.')

121 self._log('msg', 'Using existing minima file with %i prior '

122 'minima: %s' % (len(self._minima), self._minima_traj))

123 mdcount, qncount = 0, 0

124 for line in lines:

125 if (line[:4] == 'par:') and ('Ediff' not in line):

126 self._temperature = float(line.split()[1])

127 self._Ediff = float(line.split()[2])

128 elif line[:18] == 'msg: Optimization:':

129 qncount = int(line[19:].split('qn')[1])

130 elif line[:24] == 'msg: Molecular dynamics:':

131 mdcount = int(line[25:].split('md')[1])

132 self._counter = max((mdcount, qncount))

133 if qncount == mdcount:

134 # Either stopped during local optimization or terminated due to

135 # max steps.

136 self._log('msg', 'Attempting to resume at qn%05i' % qncount)

137 if qncount > 0:

138 atoms = io.read('qn%05i.traj' % (qncount - 1), index=-1)

139 self._previous_optimum = atoms.copy()

140 self._previous_energy = atoms.get_potential_energy()

141 if os.path.getsize('qn%05i.traj' % qncount) > 0:

142 atoms = io.read('qn%05i.traj' % qncount, index=-1)

143 else:

144 atoms = io.read('md%05i.traj' % qncount, index=-3)

145 self._atoms.positions = atoms.get_positions()

146 fmax = np.sqrt((atoms.get_forces() ** 2).sum(axis=1).max())

147 if fmax < self._fmax:

148 # Stopped after a qn finished.

149 self._log('msg', 'qn%05i fmax already less than fmax=%.3f'

150 % (qncount, self._fmax))

151 self._counter += 1

152 return

153 self._optimize()

154 self._counter += 1

155 if qncount > 0:

156 self._check_results()

157 else:

158 self._record_minimum()

159 self._log('msg', 'Found a new minimum.')

160 self._log('msg', 'Accepted new minimum.')

161 self._log('par')

162 elif qncount < mdcount:

163 # Probably stopped during molecular dynamics.

164 self._log('msg', 'Attempting to resume at md%05i.' % mdcount)

165 atoms = io.read('qn%05i.traj' % qncount, index=-1)

166 self._previous_optimum = atoms.copy()

167 self._previous_energy = atoms.get_potential_energy()

168 self._molecular_dynamics(resume=mdcount)

169 self._optimize()

170 self._counter += 1

171 self._check_results()

172

173 def _check_results(self):

174 """Adjusts parameters and positions based on outputs."""

175

176 # No prior minima found?

177 self._read_minima()

178 if len(self._minima) == 0:

179 self._log('msg', 'Found a new minimum.')

180 self._log('msg', 'Accepted new minimum.')

181 self._record_minimum()

182 self._log('par')

183 return

184 # Returned to starting position?

185 if self._previous_optimum:

186 compare = ComparePositions(translate=False)

187 dmax = compare(self._atoms, self._previous_optimum)

188 self._log('msg', 'Max distance to last minimum: %.3f A' % dmax)

189 if dmax < self._minima_threshold:

190 self._log('msg', 'Re-found last minimum.')

191 self._temperature *= self._beta1

192 self._log('par')

193 return

194 # In a previously found position?

195 unique, dmax_closest = self._unique_minimum_position()

196 self._log('msg', 'Max distance to closest minimum: %.3f A' %

197 dmax_closest)

198 if not unique:

199 self._temperature *= self._beta2

200 self._log('msg', 'Found previously found minimum.')

201 self._log('par')

202 if self._previous_optimum:

203 self._log('msg', 'Restoring last minimum.')

204 self._atoms.positions = self._previous_optimum.positions

205 return

206 # Must have found a unique minimum.

207 self._temperature *= self._beta3

208 self._log('msg', 'Found a new minimum.')

209 self._log('par')

210 if (self._previous_energy is None or

211 (self._atoms.get_potential_energy() <

212 self._previous_energy + self._Ediff)):

213 self._log('msg', 'Accepted new minimum.')

214 self._Ediff *= self._alpha1

215 self._log('par')

216 self._record_minimum()

217 else:

218 self._log('msg', 'Rejected new minimum due to energy. '

219 'Restoring last minimum.')

220 self._atoms.positions = self._previous_optimum.positions

221 self._Ediff *= self._alpha2

222 self._log('par')

223

224 def _log(self, cat='msg', message=None):

225 """Records the message as a line in the log file."""

226 if cat == 'init':

227 if world.rank == 0:

228 if os.path.exists(self._logfile):

229 raise RuntimeError(f'File exists: {self._logfile}')

230 fd = paropen(self._logfile, 'w')

231 fd.write('par: %12s %12s %12s\n' % ('T (K)', 'Ediff (eV)',

232 'mdmin'))

233 fd.write('ene: %12s %12s %12s\n' % ('E_current', 'E_previous',

234 'Difference'))

235 fd.close()

236 return

237 fd = paropen(self._logfile, 'a')

238 if cat == 'msg':

239 line = f'msg: {message}'

240 elif cat == 'par':

241 line = ('par: %12.4f %12.4f %12i' %

242 (self._temperature, self._Ediff, self._mdmin))

243 elif cat == 'ene':

244 current = self._atoms.get_potential_energy()

245 if self._previous_optimum:

246 previous = self._previous_energy

247 line = ('ene: %12.5f %12.5f %12.5f' %

248 (current, previous, current - previous))

249 else:

250 line = ('ene: %12.5f' % current)

251 fd.write(line + '\n')

252 fd.close()

253

254 def _optimize(self):

255 """Perform an optimization."""

256 self._atoms.set_momenta(np.zeros(self._atoms.get_momenta().shape))

257 with self._optimizer(self._atoms,

258 trajectory='qn%05i.traj' % self._counter,

259 logfile='qn%05i.log' % self._counter) as opt:

260 self._log('msg', 'Optimization: qn%05i' % self._counter)

261 opt.run(fmax=self._fmax)

262 self._log('ene')

263

264 def _record_minimum(self):

265 """Adds the current atoms configuration to the minima list."""

266 with io.Trajectory(self._minima_traj, 'a') as traj:

267 traj.write(self._atoms)

268 self._read_minima()

269 self._log('msg', 'Recorded minima #%i.' % (len(self._minima) - 1))

270

271 def _read_minima(self):

272 """Reads in the list of minima from the minima file."""

273 exists = os.path.exists(self._minima_traj)

274 if exists:

275 empty = os.path.getsize(self._minima_traj) == 0

276 if not empty:

277 with io.Trajectory(self._minima_traj, 'r') as traj:

278 self._minima = [atoms for atoms in traj]

279 else:

280 self._minima = []

281 return True

282 else:

283 self._minima = []

284 return False

285

286 def _molecular_dynamics(self, resume=None):

287 """Performs a molecular dynamics simulation, until mdmin is

288 exceeded. If resuming, the file number (md%05i) is expected."""

289 self._log('msg', 'Molecular dynamics: md%05i' % self._counter)

290 mincount = 0

291 energies, oldpositions = [], []

292 thermalized = False

293 if resume:

294 self._log('msg', 'Resuming MD from md%05i.traj' % resume)

295 if os.path.getsize('md%05i.traj' % resume) == 0:

296 self._log('msg', 'md%05i.traj is empty. Resuming from '

297 'qn%05i.traj.' % (resume, resume - 1))

298 atoms = io.read('qn%05i.traj' % (resume - 1), index=-1)

299 else:

300 with io.Trajectory('md%05i.traj' % resume, 'r') as images:

301 for atoms in images:

302 energies.append(atoms.get_potential_energy())

303 oldpositions.append(atoms.positions.copy())

304 passedmin = self._passedminimum(energies)

305 if passedmin:

306 mincount += 1

307 self._atoms.set_momenta(atoms.get_momenta())

308 thermalized = True

309 self._atoms.positions = atoms.get_positions()

310 self._log('msg', 'Starting MD with %i existing energies.' %

311 len(energies))

312 if not thermalized:

313 MaxwellBoltzmannDistribution(self._atoms,

314 temperature_K=self._temperature,

315 force_temp=True)

316 traj = io.Trajectory('md%05i.traj' % self._counter, 'a',

317 self._atoms)

318 dyn = VelocityVerlet(self._atoms, timestep=self._timestep * units.fs)

319 log = MDLogger(dyn, self._atoms, 'md%05i.log' % self._counter,

320 header=True, stress=False, peratom=False)

321

322 with traj, dyn, log:

323 dyn.attach(log, interval=1)

324 dyn.attach(traj, interval=1)

325 while mincount < self._mdmin:

326 dyn.run(1)

327 energies.append(self._atoms.get_potential_energy())

328 passedmin = self._passedminimum(energies)

329 if passedmin:

330 mincount += 1

331 oldpositions.append(self._atoms.positions.copy())

332 # Reset atoms to minimum point.

333 self._atoms.positions = oldpositions[passedmin[0]]

334

335 def _unique_minimum_position(self):

336 """Identifies if the current position of the atoms, which should be

337 a local minima, has been found before."""

338 unique = True

339 dmax_closest = 99999.

340 compare = ComparePositions(translate=True)

341 self._read_minima()

342 for minimum in self._minima:

343 dmax = compare(minimum, self._atoms)

344 if dmax < self._minima_threshold:

345 unique = False

346 if dmax < dmax_closest:

347 dmax_closest = dmax

348 return unique, dmax_closest

349

350

351class ComparePositions:

352 """Class that compares the atomic positions between two ASE atoms

353 objects. Returns the maximum distance that any atom has moved, assuming

354 all atoms of the same element are indistinguishable. If translate is

355 set to True, allows for arbitrary translations within the unit cell,

356 as well as translations across any periodic boundary conditions. When

357 called, returns the maximum displacement of any one atom."""

358

359 def __init__(self, translate=True):

360 self._translate = translate

361

362 def __call__(self, atoms1, atoms2):

363 atoms1 = atoms1.copy()

364 atoms2 = atoms2.copy()

365 if not self._translate:

366 dmax = self. _indistinguishable_compare(atoms1, atoms2)

367 else:

368 dmax = self._translated_compare(atoms1, atoms2)

369 return dmax

370

371 def _translated_compare(self, atoms1, atoms2):

372 """Moves the atoms around and tries to pair up atoms, assuming any

373 atoms with the same symbol are indistinguishable, and honors

374 periodic boundary conditions (for example, so that an atom at

375 (0.1, 0., 0.) correctly is found to be close to an atom at

376 (7.9, 0., 0.) if the atoms are in an orthorhombic cell with

377 x-dimension of 8. Returns dmax, the maximum distance between any

378 two atoms in the optimal configuration."""

379 atoms1.set_constraint()

380 atoms2.set_constraint()

381 for index in range(3):

382 assert atoms1.pbc[index] == atoms2.pbc[index]

383 least = self._get_least_common(atoms1)

384 indices1 = [atom.index for atom in atoms1 if atom.symbol == least[0]]

385 indices2 = [atom.index for atom in atoms2 if atom.symbol == least[0]]

386 # Make comparison sets from atoms2, which contain repeated atoms in

387 # all pbc's and bring the atom listed in indices2 to (0,0,0)

388 comparisons = []

389 repeat = []

390 for bc in atoms2.pbc:

391 if bc:

392 repeat.append(3)

393 else:

394 repeat.append(1)

395 repeated = atoms2.repeat(repeat)

396 moved_cell = atoms2.cell * atoms2.pbc

397 for moved in moved_cell:

398 repeated.translate(-moved)

399 repeated.set_cell(atoms2.cell)

400 for index in indices2:

401 comparison = repeated.copy()

402 comparison.translate(-atoms2[index].position)

403 comparisons.append(comparison)

404 # Bring the atom listed in indices1 to (0,0,0) [not whole list]

405 standard = atoms1.copy()

406 standard.translate(-atoms1[indices1[0]].position)

407 # Compare the standard to the comparison sets.

408 dmaxes = []

409 for comparison in comparisons:

410 dmax = self._indistinguishable_compare(standard, comparison)

411 dmaxes.append(dmax)

412 return min(dmaxes)

413

414 def _get_least_common(self, atoms):

415 """Returns the least common element in atoms. If more than one,

416 returns the first encountered."""

417 symbols = [atom.symbol for atom in atoms]

418 least = ['', np.inf]

419 for element in set(symbols):

420 count = symbols.count(element)

421 if count < least[1]:

422 least = [element, count]

423 return least

424

425 def _indistinguishable_compare(self, atoms1, atoms2):

426 """Finds each atom in atoms1's nearest neighbor with the same

427 chemical symbol in atoms2. Return dmax, the farthest distance an

428 individual atom differs by."""

429 atoms2 = atoms2.copy() # allow deletion

430 atoms2.set_constraint()

431 dmax = 0.

432 for atom1 in atoms1:

433 closest = [np.nan, np.inf]

434 for index, atom2 in enumerate(atoms2):

435 if atom2.symbol == atom1.symbol:

436 d = np.linalg.norm(atom1.position - atom2.position)

437 if d < closest[1]:

438 closest = [index, d]

439 if closest[1] > dmax:

440 dmax = closest[1]

441 del atoms2[closest[0]]

442 return dmax

443

444

445class PassedMinimum:

446 """Simple routine to find if a minimum in the potential energy surface

447 has been passed. In its default settings, a minimum is found if the

448 sequence ends with two downward points followed by two upward points.

449 Initialize with n_down and n_up, integer values of the number of up and

450 down points. If it has successfully determined it passed a minimum, it

451 returns the value (energy) of that minimum and the number of positions

452 back it occurred, otherwise returns None."""

453

454 def __init__(self, n_down=2, n_up=2):

455 self._ndown = n_down

456 self._nup = n_up

457

458 def __call__(self, energies):

459 if len(energies) < (self._nup + self._ndown + 1):

460 return None

461 status = True

462 index = -1

463 for i_up in range(self._nup):

464 if energies[index] < energies[index - 1]:

465 status = False

466 index -= 1

467 for i_down in range(self._ndown):

468 if energies[index] > energies[index - 1]:

469 status = False

470 index -= 1

471 if status:

472 return (-self._nup - 1), energies[-self._nup - 1]

473

474

475class MHPlot:

476 """Makes a plot summarizing the output of the MH algorithm from the

477 specified rundirectory. If no rundirectory is supplied, uses the

478 current directory."""

479

480 def __init__(self, rundirectory=None, logname='hop.log'):

481 if not rundirectory:

482 rundirectory = os.getcwd()

483 self._rundirectory = rundirectory

484 self._logname = logname

485 self._read_log()

486 self._fig, self._ax = self._makecanvas()

487 self._plot_data()

488

489 def get_figure(self):

490 """Returns the matplotlib figure object."""

491 return self._fig

492

493 def save_figure(self, filename):

494 """Saves the file to the specified path, with any allowed

495 matplotlib extension (e.g., .pdf, .png, etc.)."""

496 self._fig.savefig(filename)

497

498 def _read_log(self):

499 """Reads relevant parts of the log file."""

500 data = [] # format: [energy, status, temperature, ediff]

501

502 with open(os.path.join(self._rundirectory, self._logname)) as fd:

503 lines = fd.read().splitlines()

504

505 step_almost_over = False

506 step_over = False

507 for line in lines:

508 if line.startswith('msg: Molecular dynamics:'):

509 status = 'performing MD'

510 elif line.startswith('msg: Optimization:'):

511 status = 'performing QN'

512 elif line.startswith('ene:'):

513 status = 'local optimum reached'

514 energy = floatornan(line.split()[1])

515 elif line.startswith('msg: Accepted new minimum.'):

516 status = 'accepted'

517 step_almost_over = True

518 elif line.startswith('msg: Found previously found minimum.'):

519 status = 'previously found minimum'

520 step_almost_over = True

521 elif line.startswith('msg: Re-found last minimum.'):

522 status = 'previous minimum'

523 step_almost_over = True

524 elif line.startswith('msg: Rejected new minimum'):

525 status = 'rejected'

526 step_almost_over = True

527 elif line.startswith('par: '):

528 temperature = floatornan(line.split()[1])

529 ediff = floatornan(line.split()[2])

530 if step_almost_over:

531 step_over = True

532 step_almost_over = False

533 if step_over:

534 data.append([energy, status, temperature, ediff])

535 step_over = False

536 if data[-1][1] != status:

537 data.append([np.nan, status, temperature, ediff])

538 self._data = data

539

540 def _makecanvas(self):

541 from matplotlib import pyplot

542 from matplotlib.ticker import ScalarFormatter

543 fig = pyplot.figure(figsize=(6., 8.))

544 lm, rm, bm, tm = 0.22, 0.02, 0.05, 0.04

545 vg1 = 0.01 # between adjacent energy plots

546 vg2 = 0.03 # between different types of plots

547 ratio = 2. # size of an energy plot to a parameter plot

548 figwidth = 1. - lm - rm

549 totalfigheight = 1. - bm - tm - vg1 - 2. * vg2

550 parfigheight = totalfigheight / (2. * ratio + 2)

551 epotheight = ratio * parfigheight

552 ax1 = fig.add_axes((lm, bm, figwidth, epotheight))

553 ax2 = fig.add_axes((lm, bm + epotheight + vg1,

554 figwidth, epotheight))

555 for ax in [ax1, ax2]:

556 ax.yaxis.set_major_formatter(ScalarFormatter(useOffset=False))

557 ediffax = fig.add_axes((lm, bm + 2. * epotheight + vg1 + vg2,

558 figwidth, parfigheight))

559 tempax = fig.add_axes((lm, (bm + 2 * epotheight + vg1 + 2 * vg2 +

560 parfigheight), figwidth, parfigheight))

561 for ax in [ax2, tempax, ediffax]:

562 ax.set_xticklabels([])

563 ax1.set_xlabel('step')

564 tempax.set_ylabel('$T$, K')

565 ediffax.set_ylabel(r'$E_\mathrm{diff}$, eV')

566 for ax in [ax1, ax2]:

567 ax.set_ylabel(r'$E_\mathrm{pot}$, eV')

568 ax = CombinedAxis(ax1, ax2, tempax, ediffax)

569 self._set_zoomed_range(ax)

570 ax1.spines['top'].set_visible(False)

571 ax2.spines['bottom'].set_visible(False)

572 return fig, ax

573

574 def _set_zoomed_range(self, ax):

575 """Try to intelligently set the range for the zoomed-in part of the

576 graph."""

577 energies = [line[0] for line in self._data

578 if not np.isnan(line[0])]

579 dr = max(energies) - min(energies)

580 if dr == 0.:

581 dr = 1.

582 ax.set_ax1_range((min(energies) - 0.2 * dr,

583 max(energies) + 0.2 * dr))

584

585 def _plot_data(self):

586 for step, line in enumerate(self._data):

587 self._plot_energy(step, line)

588 self._plot_qn(step, line)

589 self._plot_md(step, line)

590 self._plot_parameters()

591 self._ax.set_xlim(self._ax.ax1.get_xlim())

592

593 def _plot_energy(self, step, line):

594 """Plots energy and annotation for acceptance."""

595 energy, status = line[0], line[1]

596 if np.isnan(energy):

597 return

598 self._ax.plot([step, step + 0.5], [energy] * 2, '-',

599 color='k', linewidth=2.)

600 if status == 'accepted':

601 self._ax.text(step + 0.51, energy, r'$\checkmark$')

602 elif status == 'rejected':

603 self._ax.text(step + 0.51, energy, r'$\Uparrow$', color='red')

604 elif status == 'previously found minimum':

605 self._ax.text(step + 0.51, energy, r'$\hookleftarrow$',

606 color='red', va='center')

607 elif status == 'previous minimum':

608 self._ax.text(step + 0.51, energy, r'$\leftarrow$',

609 color='red', va='center')

610

611 def _plot_md(self, step, line):

612 """Adds a curved plot of molecular dynamics trajectory."""

613 if step == 0:

614 return

615 energies = [self._data[step - 1][0]]

616 file = os.path.join(self._rundirectory, 'md%05i.traj' % step)

617 with io.Trajectory(file, 'r') as traj:

618 for atoms in traj:

619 energies.append(atoms.get_potential_energy())

620 xi = step - 1 + .5

621 if len(energies) > 2:

622 xf = xi + (step + 0.25 - xi) * len(energies) / (len(energies) - 2.)

623 else:

624 xf = step

625 if xf > (step + .75):

626 xf = step

627 self._ax.plot(np.linspace(xi, xf, num=len(energies)), energies,

628 '-k')

629

630 def _plot_qn(self, index, line):

631 """Plots a dashed vertical line for the optimization."""

632 if line[1] == 'performing MD':

633 return

634 file = os.path.join(self._rundirectory, 'qn%05i.traj' % index)

635 if os.path.getsize(file) == 0:

636 return

637 with io.Trajectory(file, 'r') as traj:

638 energies = [traj[0].get_potential_energy(),

639 traj[-1].get_potential_energy()]

640 if index > 0:

641 file = os.path.join(self._rundirectory, 'md%05i.traj' % index)

642 atoms = io.read(file, index=-3)

643 energies[0] = atoms.get_potential_energy()

644 self._ax.plot([index + 0.25] * 2, energies, ':k')

645

646 def _plot_parameters(self):

647 """Adds a plot of temperature and Ediff to the plot."""

648 steps, Ts, ediffs = [], [], []

649 for step, line in enumerate(self._data):

650 steps.extend([step + 0.5, step + 1.5])

651 Ts.extend([line[2]] * 2)

652 ediffs.extend([line[3]] * 2)

653 self._ax.tempax.plot(steps, Ts)

654 self._ax.ediffax.plot(steps, ediffs)

655

656 for ax in [self._ax.tempax, self._ax.ediffax]:

657 ylim = ax.get_ylim()

658 yrange = ylim[1] - ylim[0]

659 ax.set_ylim((ylim[0] - 0.1 * yrange, ylim[1] + 0.1 * yrange))

660

661

662def floatornan(value):

663 """Converts the argument into a float if possible, np.nan if not."""

664 try:

665 output = float(value)

666 except ValueError:

667 output = np.nan

668 return output

669

670

671class CombinedAxis:

672 """Helper class for MHPlot to plot on split y axis and adjust limits

673 simultaneously."""

674

675 def __init__(self, ax1, ax2, tempax, ediffax):

676 self.ax1 = ax1

677 self.ax2 = ax2

678 self.tempax = tempax

679 self.ediffax = ediffax

680 self._ymax = -np.inf

681

682 def set_ax1_range(self, ylim):

683 self._ax1_ylim = ylim

684 self.ax1.set_ylim(ylim)

685

686 def plot(self, *args, **kwargs):

687 self.ax1.plot(*args, **kwargs)

688 self.ax2.plot(*args, **kwargs)

689 # Re-adjust yrange

690 for yvalue in args[1]:

691 if yvalue > self._ymax:

692 self._ymax = yvalue

693 self.ax1.set_ylim(self._ax1_ylim)

694 self.ax2.set_ylim((self._ax1_ylim[1], self._ymax))

695

696 def set_xlim(self, *args):

697 self.ax1.set_xlim(*args)

698 self.ax2.set_xlim(*args)

699 self.tempax.set_xlim(*args)

700 self.ediffax.set_xlim(*args)

701

702 def text(self, *args, **kwargs):

703 y = args[1]

704 if y < self._ax1_ylim[1]:

705 ax = self.ax1

706 else:

707 ax = self.ax2

708 ax.text(*args, **kwargs)