Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1import sys 

2import threading 

3import warnings 

4from abc import ABC, abstractmethod 

5import time 

6 

7import numpy as np 

8 

9from scipy.interpolate import CubicSpline 

10from scipy.integrate import cumtrapz 

11 

12import ase.parallel 

13from ase.build import minimize_rotation_and_translation 

14from ase.calculators.calculator import Calculator 

15from ase.calculators.singlepoint import SinglePointCalculator 

16from ase.optimize import MDMin 

17from ase.optimize.optimize import Optimizer 

18from ase.optimize.sciopt import OptimizerConvergenceError 

19from ase.geometry import find_mic 

20from ase.utils import lazyproperty, deprecated 

21from ase.utils.forcecurve import fit_images 

22from ase.optimize.precon import Precon, PreconImages 

23from ase.optimize.ode import ode12r 

24 

25 

26class Spring: 

27 def __init__(self, atoms1, atoms2, energy1, energy2, k): 

28 self.atoms1 = atoms1 

29 self.atoms2 = atoms2 

30 self.energy1 = energy1 

31 self.energy2 = energy2 

32 self.k = k 

33 

34 def _find_mic(self): 

35 pos1 = self.atoms1.get_positions() 

36 pos2 = self.atoms2.get_positions() 

37 # XXX If we want variable cells we will need to edit this. 

38 mic, _ = find_mic(pos2 - pos1, self.atoms1.cell, self.atoms1.pbc) 

39 return mic 

40 

41 @lazyproperty 

42 def t(self): 

43 return self._find_mic() 

44 

45 @lazyproperty 

46 def nt(self): 

47 return np.linalg.norm(self.t) 

48 

49 

50class NEBState: 

51 def __init__(self, neb, images, energies): 

52 self.neb = neb 

53 self.images = images 

54 self.energies = energies 

55 

56 def spring(self, i): 

57 return Spring(self.images[i], self.images[i + 1], 

58 self.energies[i], self.energies[i + 1], 

59 self.neb.k[i]) 

60 

61 @lazyproperty 

62 def imax(self): 

63 return 1 + np.argsort(self.energies[1:-1])[-1] 

64 

65 @property 

66 def emax(self): 

67 return self.energies[self.imax] 

68 

69 @lazyproperty 

70 def eqlength(self): 

71 images = self.images 

72 beeline = (images[self.neb.nimages - 1].get_positions() - 

73 images[0].get_positions()) 

74 beelinelength = np.linalg.norm(beeline) 

75 return beelinelength / (self.neb.nimages - 1) 

76 

77 @lazyproperty 

78 def nimages(self): 

79 return len(self.images) 

80 

81 @property 

82 def precon(self): 

83 return self.neb.precon 

84 

85 

86class NEBMethod(ABC): 

87 def __init__(self, neb): 

88 self.neb = neb 

89 

90 @abstractmethod 

91 def get_tangent(self, state, spring1, spring2, i): 

92 ... 

93 

94 @abstractmethod 

95 def add_image_force(self, state, tangential_force, tangent, imgforce, 

96 spring1, spring2, i): 

97 ... 

98 

99 def adjust_positions(self, positions): 

100 return positions 

101 

102 

103class ImprovedTangentMethod(NEBMethod): 

104 """ 

105 Tangent estimates are improved according to Eqs. 8-11 in paper I. 

106 Tangents are weighted at extrema to ensure smooth transitions between 

107 the positive and negative tangents. 

108 """ 

109 

110 def get_tangent(self, state, spring1, spring2, i): 

111 energies = state.energies 

112 if energies[i + 1] > energies[i] > energies[i - 1]: 

113 tangent = spring2.t.copy() 

114 elif energies[i + 1] < energies[i] < energies[i - 1]: 

115 tangent = spring1.t.copy() 

116 else: 

117 deltavmax = max(abs(energies[i + 1] - energies[i]), 

118 abs(energies[i - 1] - energies[i])) 

119 deltavmin = min(abs(energies[i + 1] - energies[i]), 

120 abs(energies[i - 1] - energies[i])) 

121 if energies[i + 1] > energies[i - 1]: 

122 tangent = spring2.t * deltavmax + spring1.t * deltavmin 

123 else: 

124 tangent = spring2.t * deltavmin + spring1.t * deltavmax 

125 # Normalize the tangent vector 

126 tangent /= np.linalg.norm(tangent) 

127 return tangent 

128 

129 def add_image_force(self, state, tangential_force, tangent, imgforce, 

130 spring1, spring2, i): 

131 imgforce -= tangential_force * tangent 

132 # Improved parallel spring force (formula 12 of paper I) 

133 imgforce += (spring2.nt * spring2.k - spring1.nt * spring1.k) * tangent 

134 

135 

136class ASENEBMethod(NEBMethod): 

137 """ 

138 Standard NEB implementation in ASE. The tangent of each image is 

139 estimated from the spring closest to the saddle point in each 

140 spring pair. 

141 """ 

142 

143 def get_tangent(self, state, spring1, spring2, i): 

144 imax = self.neb.imax 

145 if i < imax: 

146 tangent = spring2.t 

147 elif i > imax: 

148 tangent = spring1.t 

149 else: 

150 tangent = spring1.t + spring2.t 

151 return tangent 

152 

153 def add_image_force(self, state, tangential_force, tangent, imgforce, 

154 spring1, spring2, i): 

155 # Magnitude for normalizing. Ensure it is not 0 

156 tangent_mag = np.vdot(tangent, tangent) or 1 

157 factor = tangent / tangent_mag 

158 imgforce -= tangential_force * factor 

159 imgforce -= np.vdot( 

160 spring1.t * spring1.k - 

161 spring2.t * spring2.k, tangent) * factor 

162 

163 

164class FullSpringMethod(NEBMethod): 

165 """ 

166 Elastic band method. The full spring force is included. 

167 """ 

168 

169 def get_tangent(self, state, spring1, spring2, i): 

170 # Tangents are bisections of spring-directions 

171 # (formula C8 of paper III) 

172 tangent = spring1.t / spring1.nt + spring2.t / spring2.nt 

173 tangent /= np.linalg.norm(tangent) 

174 return tangent 

175 

176 def add_image_force(self, state, tangential_force, tangent, imgforce, 

177 spring1, spring2, i): 

178 imgforce -= tangential_force * tangent 

179 energies = state.energies 

180 # Spring forces 

181 # Eqs. C1, C5, C6 and C7 in paper III) 

182 f1 = -(spring1.nt - 

183 state.eqlength) * spring1.t / spring1.nt * spring1.k 

184 f2 = (spring2.nt - state.eqlength) * spring2.t / spring2.nt * spring2.k 

185 if self.neb.climb and abs(i - self.neb.imax) == 1: 

186 deltavmax = max(abs(energies[i + 1] - energies[i]), 

187 abs(energies[i - 1] - energies[i])) 

188 deltavmin = min(abs(energies[i + 1] - energies[i]), 

189 abs(energies[i - 1] - energies[i])) 

190 imgforce += (f1 + f2) * deltavmin / deltavmax 

191 else: 

192 imgforce += f1 + f2 

193 

194 

195class BaseSplineMethod(NEBMethod): 

196 """ 

197 Base class for SplineNEB and String methods 

198 

199 Can optionally be preconditioned, as described in the following article: 

200 

201 S. Makri, C. Ortner and J. R. Kermode, J. Chem. Phys. 

202 150, 094109 (2019) 

203 https://dx.doi.org/10.1063/1.5064465 

204 """ 

205 

206 def __init__(self, neb): 

207 NEBMethod.__init__(self, neb) 

208 

209 def get_tangent(self, state, spring1, spring2, i): 

210 return state.precon.get_tangent(i) 

211 

212 def add_image_force(self, state, tangential_force, tangent, imgforce, 

213 spring1, spring2, i): 

214 # project out tangential component (Eqs 6 and 7 in Paper IV) 

215 imgforce -= tangential_force * tangent 

216 

217 

218class SplineMethod(BaseSplineMethod): 

219 """ 

220 NEB using spline interpolation, plus optional preconditioning 

221 """ 

222 

223 def add_image_force(self, state, tangential_force, tangent, imgforce, 

224 spring1, spring2, i): 

225 super().add_image_force(state, tangential_force, 

226 tangent, imgforce, spring1, spring2, i) 

227 eta = state.precon.get_spring_force(i, spring1.k, spring2.k, tangent) 

228 imgforce += eta 

229 

230 

231class StringMethod(BaseSplineMethod): 

232 """ 

233 String method using spline interpolation, plus optional preconditioning 

234 """ 

235 

236 def adjust_positions(self, positions): 

237 # fit cubic spline to positions, reinterpolate to equispace images 

238 # note this uses the preconditioned distance metric. 

239 fit = self.neb.spline_fit(positions) 

240 new_s = np.linspace(0.0, 1.0, self.neb.nimages) 

241 new_positions = fit.x(new_s[1:-1]).reshape(-1, 3) 

242 return new_positions 

243 

244 

245def get_neb_method(neb, method): 

246 if method == 'eb': 

247 return FullSpringMethod(neb) 

248 elif method == 'aseneb': 

249 return ASENEBMethod(neb) 

250 elif method == 'improvedtangent': 

251 return ImprovedTangentMethod(neb) 

252 elif method == 'spline': 

253 return SplineMethod(neb) 

254 elif method == 'string': 

255 return StringMethod(neb) 

256 else: 

257 raise ValueError(f'Bad method: {method}') 

258 

259 

260class BaseNEB: 

261 def __init__(self, images, k=0.1, climb=False, parallel=False, 

262 remove_rotation_and_translation=False, world=None, 

263 method='aseneb', allow_shared_calculator=False, precon=None): 

264 

265 self.images = images 

266 self.climb = climb 

267 self.parallel = parallel 

268 self.allow_shared_calculator = allow_shared_calculator 

269 

270 for img in images: 

271 if len(img) != self.natoms: 

272 raise ValueError('Images have different numbers of atoms') 

273 if np.any(img.pbc != images[0].pbc): 

274 raise ValueError('Images have different boundary conditions') 

275 if np.any(img.get_atomic_numbers() != 

276 images[0].get_atomic_numbers()): 

277 raise ValueError('Images have atoms in different orders') 

278 # check periodic cell directions 

279 cell_ok = True 

280 for pbc, vc, vc0 in zip(img.pbc, img.cell, images[0].cell): 

281 if pbc and np.any(np.abs(vc - vc0) > 1e-8): 

282 cell_ok = False 

283 if not cell_ok: 

284 raise NotImplementedError( 

285 "Variable cell in periodic directions " 

286 "is not implemented yet for NEB") 

287 

288 self.emax = np.nan 

289 

290 self.remove_rotation_and_translation = remove_rotation_and_translation 

291 

292 if method in ['aseneb', 'eb', 'improvedtangent', 'spline', 'string']: 

293 self.method = method 

294 else: 

295 raise NotImplementedError(method) 

296 

297 if precon is not None and method not in ['spline', 'string']: 

298 raise NotImplementedError(f'no precon implemented: {method}') 

299 self.precon = precon 

300 

301 self.neb_method = get_neb_method(self, method) 

302 if isinstance(k, (float, int)): 

303 k = [k] * (self.nimages - 1) 

304 self.k = list(k) 

305 

306 if world is None: 

307 world = ase.parallel.world 

308 self.world = world 

309 

310 if parallel: 

311 if self.allow_shared_calculator: 

312 raise RuntimeError( 

313 "Cannot use shared calculators in parallel in NEB.") 

314 self.real_forces = None # ndarray of shape (nimages, natom, 3) 

315 self.energies = None # ndarray of shape (nimages,) 

316 self.residuals = None # ndarray of shape (nimages,) 

317 

318 @property 

319 def natoms(self): 

320 return len(self.images[0]) 

321 

322 @property 

323 def nimages(self): 

324 return len(self.images) 

325 

326 @staticmethod 

327 def freeze_results_on_image(atoms: ase.Atoms, 

328 **results_to_include): 

329 atoms.calc = SinglePointCalculator(atoms=atoms, **results_to_include) 

330 

331 def interpolate(self, method='linear', mic=False, apply_constraint=None): 

332 """Interpolate the positions of the interior images between the 

333 initial state (image 0) and final state (image -1). 

334 

335 method: str 

336 Method by which to interpolate: 'linear' or 'idpp'. 

337 linear provides a standard straight-line interpolation, while 

338 idpp uses an image-dependent pair potential. 

339 mic: bool 

340 Use the minimum-image convention when interpolating. 

341 apply_constraint: bool 

342 Controls if the constraints attached to the images 

343 are ignored or applied when setting the interpolated positions. 

344 Default value is None, in this case the resulting constrained 

345 positions (apply_constraint=True) are compared with unconstrained 

346 positions (apply_constraint=False), 

347 if the positions are not the same 

348 the user is required to specify the desired behaviour 

349 by setting up apply_constraint keyword argument to False or True. 

350 """ 

351 if self.remove_rotation_and_translation: 

352 minimize_rotation_and_translation(self.images[0], self.images[-1]) 

353 

354 interpolate(self.images, mic, apply_constraint=apply_constraint) 

355 

356 if method == 'idpp': 

357 idpp_interpolate(images=self, traj=None, log=None, mic=mic) 

358 

359 @deprecated("Please use NEB's interpolate(method='idpp') method or " 

360 "directly call the idpp_interpolate function from ase.neb") 

361 def idpp_interpolate(self, traj='idpp.traj', log='idpp.log', fmax=0.1, 

362 optimizer=MDMin, mic=False, steps=100): 

363 idpp_interpolate(self, traj=traj, log=log, fmax=fmax, 

364 optimizer=optimizer, mic=mic, steps=steps) 

365 

366 def get_positions(self): 

367 positions = np.empty(((self.nimages - 2) * self.natoms, 3)) 

368 n1 = 0 

369 for image in self.images[1:-1]: 

370 n2 = n1 + self.natoms 

371 positions[n1:n2] = image.get_positions() 

372 n1 = n2 

373 return positions 

374 

375 def set_positions(self, positions, adjust_positions=True): 

376 if adjust_positions: 

377 # optional reparameterisation step: some NEB methods need to adjust 

378 # positions e.g. string method does this to equispace the images) 

379 positions = self.neb_method.adjust_positions(positions) 

380 n1 = 0 

381 for image in self.images[1:-1]: 

382 n2 = n1 + self.natoms 

383 image.set_positions(positions[n1:n2]) 

384 n1 = n2 

385 

386 def get_forces(self): 

387 """Evaluate and return the forces.""" 

388 images = self.images 

389 

390 if not self.allow_shared_calculator: 

391 calculators = [image.calc for image in images 

392 if image.calc is not None] 

393 if len(set(calculators)) != len(calculators): 

394 msg = ('One or more NEB images share the same calculator. ' 

395 'Each image must have its own calculator. ' 

396 'You may wish to use the ase.neb.SingleCalculatorNEB ' 

397 'class instead, although using separate calculators ' 

398 'is recommended.') 

399 raise ValueError(msg) 

400 

401 forces = np.empty(((self.nimages - 2), self.natoms, 3)) 

402 energies = np.empty(self.nimages) 

403 

404 if self.remove_rotation_and_translation: 

405 for i in range(1, self.nimages): 

406 minimize_rotation_and_translation(images[i - 1], images[i]) 

407 

408 if self.method != 'aseneb': 

409 energies[0] = images[0].get_potential_energy() 

410 energies[-1] = images[-1].get_potential_energy() 

411 

412 if not self.parallel: 

413 # Do all images - one at a time: 

414 for i in range(1, self.nimages - 1): 

415 forces[i - 1] = images[i].get_forces() 

416 energies[i] = images[i].get_potential_energy() 

417 

418 elif self.world.size == 1: 

419 def run(image, energies, forces): 

420 forces[:] = image.get_forces() 

421 energies[:] = image.get_potential_energy() 

422 

423 threads = [threading.Thread(target=run, 

424 args=(images[i], 

425 energies[i:i + 1], 

426 forces[i - 1:i])) 

427 for i in range(1, self.nimages - 1)] 

428 for thread in threads: 

429 thread.start() 

430 for thread in threads: 

431 thread.join() 

432 else: 

433 # Parallelize over images: 

434 i = self.world.rank * (self.nimages - 2) // self.world.size + 1 

435 try: 

436 forces[i - 1] = images[i].get_forces() 

437 energies[i] = images[i].get_potential_energy() 

438 except Exception: 

439 # Make sure other images also fail: 

440 error = self.world.sum(1.0) 

441 raise 

442 else: 

443 error = self.world.sum(0.0) 

444 if error: 

445 raise RuntimeError('Parallel NEB failed!') 

446 

447 for i in range(1, self.nimages - 1): 

448 root = (i - 1) * self.world.size // (self.nimages - 2) 

449 self.world.broadcast(energies[i:i + 1], root) 

450 self.world.broadcast(forces[i - 1], root) 

451 

452 # if this is the first force call, we need to build the preconditioners 

453 if (self.precon is None or isinstance(self.precon, str) or 

454 isinstance(self.precon, Precon)): 

455 self.precon = PreconImages(self.precon, images) 

456 

457 # apply preconditioners to transform forces 

458 # for the default IdentityPrecon this does not change their values 

459 precon_forces = self.precon.apply(forces, index=slice(1, -1)) 

460 

461 # Save for later use in iterimages: 

462 self.energies = energies 

463 self.real_forces = np.zeros((self.nimages, self.natoms, 3)) 

464 self.real_forces[1:-1] = forces 

465 

466 state = NEBState(self, images, energies) 

467 

468 # Can we get rid of self.energies, self.imax, self.emax etc.? 

469 self.imax = state.imax 

470 self.emax = state.emax 

471 

472 spring1 = state.spring(0) 

473 

474 self.residuals = [] 

475 for i in range(1, self.nimages - 1): 

476 spring2 = state.spring(i) 

477 tangent = self.neb_method.get_tangent(state, spring1, spring2, i) 

478 

479 # Get overlap between full PES-derived force and tangent 

480 tangential_force = np.vdot(forces[i - 1], tangent) 

481 

482 # from now on we use the preconditioned forces (equal for precon=ID) 

483 imgforce = precon_forces[i - 1] 

484 

485 if i == self.imax and self.climb: 

486 """The climbing image, imax, is not affected by the spring 

487 forces. This image feels the full PES-derived force, 

488 but the tangential component is inverted: 

489 see Eq. 5 in paper II.""" 

490 if self.method == 'aseneb': 

491 tangent_mag = np.vdot(tangent, tangent) # For normalizing 

492 imgforce -= 2 * tangential_force / tangent_mag * tangent 

493 else: 

494 imgforce -= 2 * tangential_force * tangent 

495 else: 

496 self.neb_method.add_image_force(state, tangential_force, 

497 tangent, imgforce, spring1, 

498 spring2, i) 

499 # compute the residual - with ID precon, this is just max force 

500 residual = self.precon.get_residual(i, imgforce) 

501 self.residuals.append(residual) 

502 

503 spring1 = spring2 

504 

505 return precon_forces.reshape((-1, 3)) 

506 

507 def get_residual(self): 

508 """Return residual force along the band. 

509 

510 Typically this the maximum force component on any image. For 

511 non-trivial preconditioners, the appropriate preconditioned norm 

512 is used to compute the residual. 

513 """ 

514 if self.residuals is None: 

515 raise RuntimeError("get_residual() called before get_forces()") 

516 return np.max(self.residuals) 

517 

518 def get_potential_energy(self, force_consistent=False): 

519 """Return the maximum potential energy along the band. 

520 Note that the force_consistent keyword is ignored and is only 

521 present for compatibility with ase.Atoms.get_potential_energy.""" 

522 return self.emax 

523 

524 def set_calculators(self, calculators): 

525 """Set new calculators to the images. 

526 

527 Parameters 

528 ---------- 

529 calculators : Calculator / list(Calculator) 

530 calculator(s) to attach to images 

531 - single calculator, only if allow_shared_calculator=True 

532 list of calculators if length: 

533 - length nimages, set to all images 

534 - length nimages-2, set to non-end images only 

535 """ 

536 

537 if not isinstance(calculators, list): 

538 if self.allow_shared_calculator: 

539 calculators = [calculators] * self.nimages 

540 else: 

541 raise RuntimeError("Cannot set shared calculator to NEB " 

542 "with allow_shared_calculator=False") 

543 

544 n = len(calculators) 

545 if n == self.nimages: 

546 for i in range(self.nimages): 

547 self.images[i].calc = calculators[i] 

548 elif n == self.nimages - 2: 

549 for i in range(1, self.nimages - 1): 

550 self.images[i].calc = calculators[i - 1] 

551 else: 

552 raise RuntimeError( 

553 'len(calculators)=%d does not fit to len(images)=%d' 

554 % (n, self.nimages)) 

555 

556 def __len__(self): 

557 # Corresponds to number of optimizable degrees of freedom, i.e. 

558 # virtual atom count for the optimization algorithm. 

559 return (self.nimages - 2) * self.natoms 

560 

561 def iterimages(self): 

562 # Allows trajectory to convert NEB into several images 

563 for i, atoms in enumerate(self.images): 

564 if i == 0 or i == self.nimages - 1: 

565 yield atoms 

566 else: 

567 atoms = atoms.copy() 

568 self.freeze_results_on_image( 

569 atoms, energy=self.energies[i], 

570 forces=self.real_forces[i]) 

571 

572 yield atoms 

573 

574 def spline_fit(self, positions=None, norm='precon'): 

575 """ 

576 Fit a cubic spline to this NEB 

577 

578 Args: 

579 norm (str, optional): Norm to use: 'precon' (default) or 'euclidean' 

580 

581 Returns: 

582 fit: ase.precon.precon.SplineFit instance 

583 """ 

584 if norm == 'precon': 

585 if self.precon is None or isinstance(self.precon, str): 

586 self.precon = PreconImages(self.precon, self.images) 

587 precon = self.precon 

588 # if this is the first call, we need to build the preconditioners 

589 elif norm == 'euclidean': 

590 precon = PreconImages('ID', self.images) 

591 else: 

592 raise ValueError(f'unsupported norm {norm}') 

593 return precon.spline_fit(positions) 

594 

595 def integrate_forces(self, spline_points=1000, bc_type='not-a-knot'): 

596 """Use spline fit to integrate forces along MEP to approximate 

597 energy differences using the virtual work approach. 

598 

599 Args: 

600 spline_points (int, optional): Number of points. Defaults to 1000. 

601 bc_type (str, optional): Boundary conditions, default 'not-a-knot'. 

602 

603 Returns: 

604 s: reaction coordinate in range [0, 1], with `spline_points` entries 

605 E: result of integrating forces, on the same grid as `s`. 

606 F: projected forces along MEP 

607 """ 

608 # note we use standard Euclidean rather than preconditioned norm 

609 # to compute the virtual work 

610 fit = self.spline_fit(norm='euclidean') 

611 forces = np.array([image.get_forces().reshape(-1) 

612 for image in self.images]) 

613 f = CubicSpline(fit.s, forces, bc_type=bc_type) 

614 

615 s = np.linspace(0.0, 1.0, spline_points, endpoint=True) 

616 dE = f(s) * fit.dx_ds(s) 

617 F = dE.sum(axis=1) 

618 E = -cumtrapz(F, s, initial=0.0) 

619 return s, E, F 

620 

621 

622class DyNEB(BaseNEB): 

623 def __init__(self, images, k=0.1, fmax=0.05, climb=False, parallel=False, 

624 remove_rotation_and_translation=False, world=None, 

625 dynamic_relaxation=True, scale_fmax=0., method='aseneb', 

626 allow_shared_calculator=False, precon=None): 

627 """ 

628 Subclass of NEB that allows for scaled and dynamic optimizations of 

629 images. This method, which only works in series, does not perform 

630 force calls on images that are below the convergence criterion. 

631 The convergence criteria can be scaled with a displacement metric 

632 to focus the optimization on the saddle point region. 

633 

634 'Scaled and Dynamic Optimizations of Nudged Elastic Bands', 

635 P. Lindgren, G. Kastlunger and A. A. Peterson, 

636 J. Chem. Theory Comput. 15, 11, 5787-5793 (2019). 

637 

638 dynamic_relaxation: bool 

639 True skips images with forces below the convergence criterion. 

640 This is updated after each force call; if a previously converged 

641 image goes out of tolerance (due to spring adjustments between 

642 the image and its neighbors), it will be optimized again. 

643 False reverts to the default NEB implementation. 

644 

645 fmax: float 

646 Must be identical to the fmax of the optimizer. 

647 

648 scale_fmax: float 

649 Scale convergence criteria along band based on the distance between 

650 an image and the image with the highest potential energy. This 

651 keyword determines how rapidly the convergence criteria are scaled. 

652 """ 

653 super().__init__( 

654 images, k=k, climb=climb, parallel=parallel, 

655 remove_rotation_and_translation=remove_rotation_and_translation, 

656 world=world, method=method, 

657 allow_shared_calculator=allow_shared_calculator, precon=precon) 

658 self.fmax = fmax 

659 self.dynamic_relaxation = dynamic_relaxation 

660 self.scale_fmax = scale_fmax 

661 

662 if not self.dynamic_relaxation and self.scale_fmax: 

663 msg = ('Scaled convergence criteria only implemented in series ' 

664 'with dynamic relaxation.') 

665 raise ValueError(msg) 

666 

667 def set_positions(self, positions): 

668 if not self.dynamic_relaxation: 

669 return super().set_positions(positions) 

670 

671 n1 = 0 

672 for i, image in enumerate(self.images[1:-1]): 

673 if self.parallel: 

674 msg = ('Dynamic relaxation does not work efficiently ' 

675 'when parallelizing over images. Try AutoNEB ' 

676 'routine for freezing images in parallel.') 

677 raise ValueError(msg) 

678 else: 

679 forces_dyn = self._fmax_all(self.images) 

680 if forces_dyn[i] < self.fmax: 

681 n1 += self.natoms 

682 else: 

683 n2 = n1 + self.natoms 

684 image.set_positions(positions[n1:n2]) 

685 n1 = n2 

686 

687 def _fmax_all(self, images): 

688 """Store maximum force acting on each image in list. This is used in 

689 the dynamic optimization routine in the set_positions() function.""" 

690 n = self.natoms 

691 forces = self.get_forces() 

692 fmax_images = [ 

693 np.sqrt((forces[n * i:n + n * i] ** 2).sum(axis=1)).max() 

694 for i in range(self.nimages - 2)] 

695 return fmax_images 

696 

697 def get_forces(self): 

698 forces = super().get_forces() 

699 if not self.dynamic_relaxation: 

700 return forces 

701 

702 """Get NEB forces and scale the convergence criteria to focus 

703 optimization on saddle point region. The keyword scale_fmax 

704 determines the rate of convergence scaling.""" 

705 n = self.natoms 

706 for i in range(self.nimages - 2): 

707 n1 = n * i 

708 n2 = n1 + n 

709 force = np.sqrt((forces[n1:n2] ** 2.).sum(axis=1)).max() 

710 n_imax = (self.imax - 1) * n # Image with highest energy. 

711 

712 positions = self.get_positions() 

713 pos_imax = positions[n_imax:n_imax + n] 

714 

715 """Scale convergence criteria based on distance between an 

716 image and the image with the highest potential energy.""" 

717 rel_pos = np.sqrt(((positions[n1:n2] - pos_imax) ** 2).sum()) 

718 if force < self.fmax * (1 + rel_pos * self.scale_fmax): 

719 if i == self.imax - 1: 

720 # Keep forces at saddle point for the log file. 

721 pass 

722 else: 

723 # Set forces to zero before they are sent to optimizer. 

724 forces[n1:n2, :] = 0 

725 return forces 

726 

727 

728def _check_deprecation(keyword, kwargs): 

729 if keyword in kwargs: 

730 warnings.warn(f'Keyword {keyword} of NEB is deprecated. ' 

731 'Please use the DyNEB class instead for dynamic ' 

732 'relaxation', FutureWarning) 

733 

734 

735class NEB(DyNEB): 

736 def __init__(self, images, k=0.1, climb=False, parallel=False, 

737 remove_rotation_and_translation=False, world=None, 

738 method='aseneb', allow_shared_calculator=False, 

739 precon=None, **kwargs): 

740 """Nudged elastic band. 

741 

742 Paper I: 

743 

744 G. Henkelman and H. Jonsson, Chem. Phys, 113, 9978 (2000). 

745 :doi:`10.1063/1.1323224` 

746 

747 Paper II: 

748 

749 G. Henkelman, B. P. Uberuaga, and H. Jonsson, Chem. Phys, 

750 113, 9901 (2000). 

751 :doi:`10.1063/1.1329672` 

752 

753 Paper III: 

754 

755 E. L. Kolsbjerg, M. N. Groves, and B. Hammer, J. Chem. Phys, 

756 145, 094107 (2016) 

757 :doi:`10.1063/1.4961868` 

758 

759 Paper IV: 

760 

761 S. Makri, C. Ortner and J. R. Kermode, J. Chem. Phys. 

762 150, 094109 (2019) 

763 https://dx.doi.org/10.1063/1.5064465 

764 

765 images: list of Atoms objects 

766 Images defining path from initial to final state. 

767 k: float or list of floats 

768 Spring constant(s) in eV/Ang. One number or one for each spring. 

769 climb: bool 

770 Use a climbing image (default is no climbing image). 

771 parallel: bool 

772 Distribute images over processors. 

773 remove_rotation_and_translation: bool 

774 TRUE actives NEB-TR for removing translation and 

775 rotation during NEB. By default applied non-periodic 

776 systems 

777 method: string of method 

778 Choice betweeen five methods: 

779 

780 * aseneb: standard ase NEB implementation 

781 * improvedtangent: Paper I NEB implementation 

782 * eb: Paper III full spring force implementation 

783 * spline: Paper IV spline interpolation (supports precon) 

784 * string: Paper IV string method (supports precon) 

785 allow_shared_calculator: bool 

786 Allow images to share the same calculator between them. 

787 Incompatible with parallelisation over images. 

788 precon: string, :class:`ase.optimize.precon.Precon` instance or list of 

789 instances. If present, enable preconditioing as in Paper IV. This is 

790 possible using the 'spline' or 'string' methods only. 

791 Default is no preconditioning (precon=None), which is converted to 

792 a list of :class:`ase.precon.precon.IdentityPrecon` instances. 

793 """ 

794 for keyword in 'dynamic_relaxation', 'fmax', 'scale_fmax': 

795 _check_deprecation(keyword, kwargs) 

796 defaults = dict(dynamic_relaxation=False, 

797 fmax=0.05, 

798 scale_fmax=0.0) 

799 defaults.update(kwargs) 

800 # Only reason for separating BaseNEB/NEB is that we are 

801 # deprecating dynamic_relaxation. 

802 # 

803 # We can turn BaseNEB into NEB once we get rid of the 

804 # deprecated variables. 

805 # 

806 # Then we can also move DyNEB into ase.dyneb without cyclic imports. 

807 # We can do that in ase-3.22 or 3.23. 

808 super().__init__( 

809 images, k=k, climb=climb, parallel=parallel, 

810 remove_rotation_and_translation=remove_rotation_and_translation, 

811 world=world, method=method, 

812 allow_shared_calculator=allow_shared_calculator, 

813 precon=precon, 

814 **defaults) 

815 

816 

817class NEBOptimizer(Optimizer): 

818 """ 

819 This optimizer applies an adaptive ODE solver to a NEB 

820 

821 Details of the adaptive ODE solver are described in paper IV 

822 """ 

823 

824 def __init__(self, 

825 neb, 

826 restart=None, logfile='-', trajectory=None, 

827 master=None, 

828 append_trajectory=False, 

829 method='ODE', 

830 alpha=0.01, 

831 verbose=0, 

832 rtol=0.1, 

833 C1=1e-2, 

834 C2=2.0): 

835 

836 super().__init__(atoms=neb, restart=restart, 

837 logfile=logfile, trajectory=trajectory, 

838 master=master, 

839 append_trajectory=append_trajectory, 

840 force_consistent=False) 

841 self.neb = neb 

842 

843 method = method.lower() 

844 methods = ['ode', 'static', 'krylov'] 

845 if method not in methods: 

846 raise ValueError(f'method must be one of {methods}') 

847 self.method = method 

848 

849 self.alpha = alpha 

850 self.verbose = verbose 

851 self.rtol = rtol 

852 self.C1 = C1 

853 self.C2 = C2 

854 

855 def force_function(self, X): 

856 positions = X.reshape((self.neb.nimages - 2) * 

857 self.neb.natoms, 3) 

858 self.neb.set_positions(positions) 

859 forces = self.neb.get_forces().reshape(-1) 

860 return forces 

861 

862 def get_residual(self, F=None, X=None): 

863 return self.neb.get_residual() 

864 

865 def log(self): 

866 fmax = self.get_residual() 

867 T = time.localtime() 

868 if self.logfile is not None: 

869 name = f'{self.__class__.__name__}[{self.method}]' 

870 if self.nsteps == 0: 

871 args = (" " * len(name), "Step", "Time", "fmax") 

872 msg = "%s %4s %8s %12s\n" % args 

873 self.logfile.write(msg) 

874 

875 args = (name, self.nsteps, T[3], T[4], T[5], fmax) 

876 msg = "%s: %3d %02d:%02d:%02d %12.4f\n" % args 

877 self.logfile.write(msg) 

878 self.logfile.flush() 

879 

880 def callback(self, X, F=None): 

881 self.log() 

882 self.call_observers() 

883 self.nsteps += 1 

884 

885 def run_ode(self, fmax): 

886 try: 

887 ode12r(self.force_function, 

888 self.neb.get_positions().reshape(-1), 

889 fmax=fmax, 

890 rtol=self.rtol, 

891 C1=self.C1, 

892 C2=self.C2, 

893 steps=self.max_steps, 

894 verbose=self.verbose, 

895 callback=self.callback, 

896 residual=self.get_residual) 

897 return True 

898 except OptimizerConvergenceError: 

899 return False 

900 

901 def run_static(self, fmax): 

902 X = self.neb.get_positions().reshape(-1) 

903 for step in range(self.max_steps): 

904 F = self.force_function(X) 

905 if self.neb.get_residual() <= fmax: 

906 return True 

907 X += self.alpha * F 

908 self.callback(X) 

909 return False 

910 

911 def run(self, fmax=0.05, steps=None, method=None): 

912 """ 

913 Optimize images to obtain the minimum energy path 

914 

915 Parameters 

916 ---------- 

917 fmax - desired force tolerance 

918 steps - maximum number of steps 

919 """ 

920 if steps: 

921 self.max_steps = steps 

922 if method is None: 

923 method = self.method 

924 if method == 'ode': 

925 return self.run_ode(fmax) 

926 elif method == 'static': 

927 return self.run_static(fmax) 

928 else: 

929 raise ValueError(f'unknown method: {self.method}') 

930 

931 

932class IDPP(Calculator): 

933 """Image dependent pair potential. 

934 

935 See: 

936 Improved initial guess for minimum energy path calculations. 

937 Søren Smidstrup, Andreas Pedersen, Kurt Stokbro and Hannes Jónsson 

938 Chem. Phys. 140, 214106 (2014) 

939 """ 

940 

941 implemented_properties = ['energy', 'forces'] 

942 

943 def __init__(self, target, mic): 

944 Calculator.__init__(self) 

945 self.target = target 

946 self.mic = mic 

947 

948 def calculate(self, atoms, properties, system_changes): 

949 Calculator.calculate(self, atoms, properties, system_changes) 

950 

951 P = atoms.get_positions() 

952 d = [] 

953 D = [] 

954 for p in P: 

955 Di = P - p 

956 if self.mic: 

957 Di, di = find_mic(Di, atoms.get_cell(), atoms.get_pbc()) 

958 else: 

959 di = np.sqrt((Di ** 2).sum(1)) 

960 d.append(di) 

961 D.append(Di) 

962 d = np.array(d) 

963 D = np.array(D) 

964 

965 dd = d - self.target 

966 d.ravel()[::len(d) + 1] = 1 # avoid dividing by zero 

967 d4 = d ** 4 

968 e = 0.5 * (dd ** 2 / d4).sum() 

969 f = -2 * ((dd * (1 - 2 * dd / d) / d ** 5)[..., np.newaxis] * D).sum( 

970 0) 

971 self.results = {'energy': e, 'forces': f} 

972 

973 

974@deprecated("SingleCalculatorNEB is deprecated. " 

975 "Please use NEB(allow_shared_calculator=True) instead.") 

976class SingleCalculatorNEB(NEB): 

977 def __init__(self, images, *args, **kwargs): 

978 kwargs["allow_shared_calculator"] = True 

979 super().__init__(images, *args, **kwargs) 

980 

981 

982def interpolate(images, mic=False, interpolate_cell=False, 

983 use_scaled_coord=False, apply_constraint=None): 

984 """Given a list of images, linearly interpolate the positions of the 

985 interior images. 

986 

987 mic: bool 

988 Map movement into the unit cell by using the minimum image convention. 

989 interpolate_cell: bool 

990 Interpolate the three cell vectors linearly just like the atomic 

991 positions. Not implemented for NEB calculations! 

992 use_scaled_coord: bool 

993 Use scaled/internal/fractional coordinates instead of real ones for the 

994 interpolation. Not implemented for NEB calculations! 

995 apply_constraint: bool 

996 Controls if the constraints attached to the images 

997 are ignored or applied when setting the interpolated positions. 

998 Default value is None, in this case the resulting constrained positions 

999 (apply_constraint=True) are compared with unconstrained positions 

1000 (apply_constraint=False), if the positions are not the same 

1001 the user is required to specify the desired behaviour 

1002 by setting up apply_constraint keyword argument to False or True. 

1003 """ 

1004 if use_scaled_coord: 

1005 pos1 = images[0].get_scaled_positions(wrap=mic) 

1006 pos2 = images[-1].get_scaled_positions(wrap=mic) 

1007 else: 

1008 pos1 = images[0].get_positions() 

1009 pos2 = images[-1].get_positions() 

1010 d = pos2 - pos1 

1011 if not use_scaled_coord and mic: 

1012 d = find_mic(d, images[0].get_cell(), images[0].pbc)[0] 

1013 d /= (len(images) - 1.0) 

1014 if interpolate_cell: 

1015 cell1 = images[0].get_cell() 

1016 cell2 = images[-1].get_cell() 

1017 cell_diff = cell2 - cell1 

1018 cell_diff /= (len(images) - 1.0) 

1019 for i in range(1, len(images) - 1): 

1020 # first the new cell, otherwise scaled positions are wrong 

1021 if interpolate_cell: 

1022 images[i].set_cell(cell1 + i * cell_diff) 

1023 new_pos = pos1 + i * d 

1024 if use_scaled_coord: 

1025 images[i].set_scaled_positions(new_pos) 

1026 else: 

1027 if apply_constraint is None: 

1028 unconstrained_image = images[i].copy() 

1029 unconstrained_image.set_positions(new_pos, 

1030 apply_constraint=False) 

1031 images[i].set_positions(new_pos, apply_constraint=True) 

1032 try: 

1033 np.testing.assert_allclose(unconstrained_image.positions, 

1034 images[i].positions) 

1035 except AssertionError: 

1036 raise RuntimeError(f"Constraint(s) in image number {i} \n" 

1037 f"affect the interpolation results.\n" 

1038 "Please specify if you want to \n" 

1039 "apply or ignore the constraints \n" 

1040 "during the interpolation \n" 

1041 "with apply_constraint argument.") 

1042 else: 

1043 images[i].set_positions(new_pos, 

1044 apply_constraint=apply_constraint) 

1045 

1046 

1047def idpp_interpolate(images, traj='idpp.traj', log='idpp.log', fmax=0.1, 

1048 optimizer=MDMin, mic=False, steps=100): 

1049 """Interpolate using the IDPP method. 'images' can either be a plain 

1050 list of images or an NEB object (containing a list of images).""" 

1051 if hasattr(images, 'interpolate'): 

1052 neb = images 

1053 else: 

1054 neb = NEB(images) 

1055 

1056 d1 = neb.images[0].get_all_distances(mic=mic) 

1057 d2 = neb.images[-1].get_all_distances(mic=mic) 

1058 d = (d2 - d1) / (neb.nimages - 1) 

1059 real_calcs = [] 

1060 for i, image in enumerate(neb.images): 

1061 real_calcs.append(image.calc) 

1062 image.calc = IDPP(d1 + i * d, mic=mic) 

1063 

1064 with optimizer(neb, trajectory=traj, logfile=log) as opt: 

1065 opt.run(fmax=fmax, steps=steps) 

1066 

1067 for image, calc in zip(neb.images, real_calcs): 

1068 image.calc = calc 

1069 

1070 

1071class NEBTools: 

1072 """Class to make many of the common tools for NEB analysis available to 

1073 the user. Useful for scripting the output of many jobs. Initialize with 

1074 list of images which make up one or more band of the NEB relaxation.""" 

1075 

1076 def __init__(self, images): 

1077 self.images = images 

1078 

1079 @deprecated('NEBTools.get_fit() is deprecated. ' 

1080 'Please use ase.utils.forcecurve.fit_images(images).') 

1081 def get_fit(self): 

1082 return fit_images(self.images) 

1083 

1084 def get_barrier(self, fit=True, raw=False): 

1085 """Returns the barrier estimate from the NEB, along with the 

1086 Delta E of the elementary reaction. If fit=True, the barrier is 

1087 estimated based on the interpolated fit to the images; if 

1088 fit=False, the barrier is taken as the maximum-energy image 

1089 without interpolation. Set raw=True to get the raw energy of the 

1090 transition state instead of the forward barrier.""" 

1091 forcefit = fit_images(self.images) 

1092 energies = forcefit.energies 

1093 fit_energies = forcefit.fit_energies 

1094 dE = energies[-1] - energies[0] 

1095 if fit: 

1096 barrier = max(fit_energies) 

1097 else: 

1098 barrier = max(energies) 

1099 if raw: 

1100 barrier += self.images[0].get_potential_energy() 

1101 return barrier, dE 

1102 

1103 def get_fmax(self, **kwargs): 

1104 """Returns fmax, as used by optimizers with NEB.""" 

1105 neb = NEB(self.images, **kwargs) 

1106 forces = neb.get_forces() 

1107 return np.sqrt((forces ** 2).sum(axis=1).max()) 

1108 

1109 def plot_band(self, ax=None): 

1110 """Plots the NEB band on matplotlib axes object 'ax'. If ax=None 

1111 returns a new figure object.""" 

1112 forcefit = fit_images(self.images) 

1113 ax = forcefit.plot(ax=ax) 

1114 return ax.figure 

1115 

1116 def plot_bands(self, constant_x=False, constant_y=False, 

1117 nimages=None, label='nebplots'): 

1118 """Given a trajectory containing many steps of a NEB, makes 

1119 plots of each band in the series in a single PDF. 

1120 

1121 constant_x: bool 

1122 Use the same x limits on all plots. 

1123 constant_y: bool 

1124 Use the same y limits on all plots. 

1125 nimages: int 

1126 Number of images per band. Guessed if not supplied. 

1127 label: str 

1128 Name for the output file. .pdf will be appended. 

1129 """ 

1130 from matplotlib import pyplot 

1131 from matplotlib.backends.backend_pdf import PdfPages 

1132 if nimages is None: 

1133 nimages = self._guess_nimages() 

1134 nebsteps = len(self.images) // nimages 

1135 if constant_x or constant_y: 

1136 sys.stdout.write('Scaling axes.\n') 

1137 sys.stdout.flush() 

1138 # Plot all to one plot, then pull its x and y range. 

1139 fig, ax = pyplot.subplots() 

1140 for index in range(nebsteps): 

1141 images = self.images[index * nimages:(index + 1) * nimages] 

1142 NEBTools(images).plot_band(ax=ax) 

1143 xlim = ax.get_xlim() 

1144 ylim = ax.get_ylim() 

1145 pyplot.close(fig) # Reference counting "bug" in pyplot. 

1146 with PdfPages(label + '.pdf') as pdf: 

1147 for index in range(nebsteps): 

1148 sys.stdout.write('\rProcessing band {:10d} / {:10d}' 

1149 .format(index, nebsteps)) 

1150 sys.stdout.flush() 

1151 fig, ax = pyplot.subplots() 

1152 images = self.images[index * nimages:(index + 1) * nimages] 

1153 NEBTools(images).plot_band(ax=ax) 

1154 if constant_x: 

1155 ax.set_xlim(xlim) 

1156 if constant_y: 

1157 ax.set_ylim(ylim) 

1158 pdf.savefig(fig) 

1159 pyplot.close(fig) # Reference counting "bug" in pyplot. 

1160 sys.stdout.write('\n') 

1161 

1162 def _guess_nimages(self): 

1163 """Attempts to guess the number of images per band from 

1164 a trajectory, based solely on the repetition of the 

1165 potential energy of images. This should also work for symmetric 

1166 cases.""" 

1167 e_first = self.images[0].get_potential_energy() 

1168 nimages = None 

1169 for index, image in enumerate(self.images[1:], start=1): 

1170 e = image.get_potential_energy() 

1171 if e == e_first: 

1172 # Need to check for symmetric case when e_first = e_last. 

1173 try: 

1174 e_next = self.images[index + 1].get_potential_energy() 

1175 except IndexError: 

1176 pass 

1177 else: 

1178 if e_next == e_first: 

1179 nimages = index + 1 # Symmetric 

1180 break 

1181 nimages = index # Normal 

1182 break 

1183 if nimages is None: 

1184 sys.stdout.write('Appears to be only one band in the images.\n') 

1185 return len(self.images) 

1186 # Sanity check that the energies of the last images line up too. 

1187 e_last = self.images[nimages - 1].get_potential_energy() 

1188 e_nextlast = self.images[2 * nimages - 1].get_potential_energy() 

1189 if not (e_last == e_nextlast): 

1190 raise RuntimeError('Could not guess number of images per band.') 

1191 sys.stdout.write('Number of images per band guessed to be {:d}.\n' 

1192 .format(nimages)) 

1193 return nimages 

1194 

1195 

1196class NEBtools(NEBTools): 

1197 @deprecated('NEBtools has been renamed; please use NEBTools.') 

1198 def __init__(self, images): 

1199 NEBTools.__init__(self, images) 

1200 

1201 

1202@deprecated('Please use NEBTools.plot_band_from_fit.') 

1203def plot_band_from_fit(s, E, Sfit, Efit, lines, ax=None): 

1204 NEBTools.plot_band_from_fit(s, E, Sfit, Efit, lines, ax=None) 

1205 

1206 

1207def fit0(*args, **kwargs): 

1208 raise DeprecationWarning('fit0 is deprecated. Use `fit_raw` from ' 

1209 '`ase.utils.forcecurve` instead.')