Coverage for /builds/ase/ase/ase/neb.py : 83.25%
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
7import numpy as np
9from scipy.interpolate import CubicSpline
10from scipy.integrate import cumtrapz
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
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
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
41 @lazyproperty
42 def t(self):
43 return self._find_mic()
45 @lazyproperty
46 def nt(self):
47 return np.linalg.norm(self.t)
50class NEBState:
51 def __init__(self, neb, images, energies):
52 self.neb = neb
53 self.images = images
54 self.energies = energies
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])
61 @lazyproperty
62 def imax(self):
63 return 1 + np.argsort(self.energies[1:-1])[-1]
65 @property
66 def emax(self):
67 return self.energies[self.imax]
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)
77 @lazyproperty
78 def nimages(self):
79 return len(self.images)
81 @property
82 def precon(self):
83 return self.neb.precon
86class NEBMethod(ABC):
87 def __init__(self, neb):
88 self.neb = neb
90 @abstractmethod
91 def get_tangent(self, state, spring1, spring2, i):
92 ...
94 @abstractmethod
95 def add_image_force(self, state, tangential_force, tangent, imgforce,
96 spring1, spring2, i):
97 ...
99 def adjust_positions(self, positions):
100 return positions
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 """
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
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
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 """
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
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
164class FullSpringMethod(NEBMethod):
165 """
166 Elastic band method. The full spring force is included.
167 """
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
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
195class BaseSplineMethod(NEBMethod):
196 """
197 Base class for SplineNEB and String methods
199 Can optionally be preconditioned, as described in the following article:
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 """
206 def __init__(self, neb):
207 NEBMethod.__init__(self, neb)
209 def get_tangent(self, state, spring1, spring2, i):
210 return state.precon.get_tangent(i)
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
218class SplineMethod(BaseSplineMethod):
219 """
220 NEB using spline interpolation, plus optional preconditioning
221 """
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
231class StringMethod(BaseSplineMethod):
232 """
233 String method using spline interpolation, plus optional preconditioning
234 """
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
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}')
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):
265 self.images = images
266 self.climb = climb
267 self.parallel = parallel
268 self.allow_shared_calculator = allow_shared_calculator
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")
288 self.emax = np.nan
290 self.remove_rotation_and_translation = remove_rotation_and_translation
292 if method in ['aseneb', 'eb', 'improvedtangent', 'spline', 'string']:
293 self.method = method
294 else:
295 raise NotImplementedError(method)
297 if precon is not None and method not in ['spline', 'string']:
298 raise NotImplementedError(f'no precon implemented: {method}')
299 self.precon = precon
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)
306 if world is None:
307 world = ase.parallel.world
308 self.world = world
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,)
318 @property
319 def natoms(self):
320 return len(self.images[0])
322 @property
323 def nimages(self):
324 return len(self.images)
326 @staticmethod
327 def freeze_results_on_image(atoms: ase.Atoms,
328 **results_to_include):
329 atoms.calc = SinglePointCalculator(atoms=atoms, **results_to_include)
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).
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])
354 interpolate(self.images, mic, apply_constraint=apply_constraint)
356 if method == 'idpp':
357 idpp_interpolate(images=self, traj=None, log=None, mic=mic)
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)
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
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
386 def get_forces(self):
387 """Evaluate and return the forces."""
388 images = self.images
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)
401 forces = np.empty(((self.nimages - 2), self.natoms, 3))
402 energies = np.empty(self.nimages)
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])
408 if self.method != 'aseneb':
409 energies[0] = images[0].get_potential_energy()
410 energies[-1] = images[-1].get_potential_energy()
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()
418 elif self.world.size == 1:
419 def run(image, energies, forces):
420 forces[:] = image.get_forces()
421 energies[:] = image.get_potential_energy()
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!')
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)
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) or
455 isinstance(self.precon, list)):
456 self.precon = PreconImages(self.precon, images)
458 # apply preconditioners to transform forces
459 # for the default IdentityPrecon this does not change their values
460 precon_forces = self.precon.apply(forces, index=slice(1, -1))
462 # Save for later use in iterimages:
463 self.energies = energies
464 self.real_forces = np.zeros((self.nimages, self.natoms, 3))
465 self.real_forces[1:-1] = forces
467 state = NEBState(self, images, energies)
469 # Can we get rid of self.energies, self.imax, self.emax etc.?
470 self.imax = state.imax
471 self.emax = state.emax
473 spring1 = state.spring(0)
475 self.residuals = []
476 for i in range(1, self.nimages - 1):
477 spring2 = state.spring(i)
478 tangent = self.neb_method.get_tangent(state, spring1, spring2, i)
480 # Get overlap between full PES-derived force and tangent
481 tangential_force = np.vdot(forces[i - 1], tangent)
483 # from now on we use the preconditioned forces (equal for precon=ID)
484 imgforce = precon_forces[i - 1]
486 if i == self.imax and self.climb:
487 """The climbing image, imax, is not affected by the spring
488 forces. This image feels the full PES-derived force,
489 but the tangential component is inverted:
490 see Eq. 5 in paper II."""
491 if self.method == 'aseneb':
492 tangent_mag = np.vdot(tangent, tangent) # For normalizing
493 imgforce -= 2 * tangential_force / tangent_mag * tangent
494 else:
495 imgforce -= 2 * tangential_force * tangent
496 else:
497 self.neb_method.add_image_force(state, tangential_force,
498 tangent, imgforce, spring1,
499 spring2, i)
500 # compute the residual - with ID precon, this is just max force
501 residual = self.precon.get_residual(i, imgforce)
502 self.residuals.append(residual)
504 spring1 = spring2
506 return precon_forces.reshape((-1, 3))
508 def get_residual(self):
509 """Return residual force along the band.
511 Typically this the maximum force component on any image. For
512 non-trivial preconditioners, the appropriate preconditioned norm
513 is used to compute the residual.
514 """
515 if self.residuals is None:
516 raise RuntimeError("get_residual() called before get_forces()")
517 return np.max(self.residuals)
519 def get_potential_energy(self, force_consistent=False):
520 """Return the maximum potential energy along the band.
521 Note that the force_consistent keyword is ignored and is only
522 present for compatibility with ase.Atoms.get_potential_energy."""
523 return self.emax
525 def set_calculators(self, calculators):
526 """Set new calculators to the images.
528 Parameters
529 ----------
530 calculators : Calculator / list(Calculator)
531 calculator(s) to attach to images
532 - single calculator, only if allow_shared_calculator=True
533 list of calculators if length:
534 - length nimages, set to all images
535 - length nimages-2, set to non-end images only
536 """
538 if not isinstance(calculators, list):
539 if self.allow_shared_calculator:
540 calculators = [calculators] * self.nimages
541 else:
542 raise RuntimeError("Cannot set shared calculator to NEB "
543 "with allow_shared_calculator=False")
545 n = len(calculators)
546 if n == self.nimages:
547 for i in range(self.nimages):
548 self.images[i].calc = calculators[i]
549 elif n == self.nimages - 2:
550 for i in range(1, self.nimages - 1):
551 self.images[i].calc = calculators[i - 1]
552 else:
553 raise RuntimeError(
554 'len(calculators)=%d does not fit to len(images)=%d'
555 % (n, self.nimages))
557 def __len__(self):
558 # Corresponds to number of optimizable degrees of freedom, i.e.
559 # virtual atom count for the optimization algorithm.
560 return (self.nimages - 2) * self.natoms
562 def iterimages(self):
563 # Allows trajectory to convert NEB into several images
564 for i, atoms in enumerate(self.images):
565 if i == 0 or i == self.nimages - 1:
566 yield atoms
567 else:
568 atoms = atoms.copy()
569 self.freeze_results_on_image(
570 atoms, energy=self.energies[i],
571 forces=self.real_forces[i])
573 yield atoms
575 def spline_fit(self, positions=None, norm='precon'):
576 """
577 Fit a cubic spline to this NEB
579 Args:
580 norm (str, optional): Norm to use: 'precon' (default) or 'euclidean'
582 Returns:
583 fit: ase.precon.precon.SplineFit instance
584 """
585 if norm == 'precon':
586 if self.precon is None or isinstance(self.precon, str):
587 self.precon = PreconImages(self.precon, self.images)
588 precon = self.precon
589 # if this is the first call, we need to build the preconditioners
590 elif norm == 'euclidean':
591 precon = PreconImages('ID', self.images)
592 else:
593 raise ValueError(f'unsupported norm {norm}')
594 return precon.spline_fit(positions)
596 def integrate_forces(self, spline_points=1000, bc_type='not-a-knot'):
597 """Use spline fit to integrate forces along MEP to approximate
598 energy differences using the virtual work approach.
600 Args:
601 spline_points (int, optional): Number of points. Defaults to 1000.
602 bc_type (str, optional): Boundary conditions, default 'not-a-knot'.
604 Returns:
605 s: reaction coordinate in range [0, 1], with `spline_points` entries
606 E: result of integrating forces, on the same grid as `s`.
607 F: projected forces along MEP
608 """
609 # note we use standard Euclidean rather than preconditioned norm
610 # to compute the virtual work
611 fit = self.spline_fit(norm='euclidean')
612 forces = np.array([image.get_forces().reshape(-1)
613 for image in self.images])
614 f = CubicSpline(fit.s, forces, bc_type=bc_type)
616 s = np.linspace(0.0, 1.0, spline_points, endpoint=True)
617 dE = f(s) * fit.dx_ds(s)
618 F = dE.sum(axis=1)
619 E = -cumtrapz(F, s, initial=0.0)
620 return s, E, F
623class DyNEB(BaseNEB):
624 def __init__(self, images, k=0.1, fmax=0.05, climb=False, parallel=False,
625 remove_rotation_and_translation=False, world=None,
626 dynamic_relaxation=True, scale_fmax=0., method='aseneb',
627 allow_shared_calculator=False, precon=None):
628 """
629 Subclass of NEB that allows for scaled and dynamic optimizations of
630 images. This method, which only works in series, does not perform
631 force calls on images that are below the convergence criterion.
632 The convergence criteria can be scaled with a displacement metric
633 to focus the optimization on the saddle point region.
635 'Scaled and Dynamic Optimizations of Nudged Elastic Bands',
636 P. Lindgren, G. Kastlunger and A. A. Peterson,
637 J. Chem. Theory Comput. 15, 11, 5787-5793 (2019).
639 dynamic_relaxation: bool
640 True skips images with forces below the convergence criterion.
641 This is updated after each force call; if a previously converged
642 image goes out of tolerance (due to spring adjustments between
643 the image and its neighbors), it will be optimized again.
644 False reverts to the default NEB implementation.
646 fmax: float
647 Must be identical to the fmax of the optimizer.
649 scale_fmax: float
650 Scale convergence criteria along band based on the distance between
651 an image and the image with the highest potential energy. This
652 keyword determines how rapidly the convergence criteria are scaled.
653 """
654 super().__init__(
655 images, k=k, climb=climb, parallel=parallel,
656 remove_rotation_and_translation=remove_rotation_and_translation,
657 world=world, method=method,
658 allow_shared_calculator=allow_shared_calculator, precon=precon)
659 self.fmax = fmax
660 self.dynamic_relaxation = dynamic_relaxation
661 self.scale_fmax = scale_fmax
663 if not self.dynamic_relaxation and self.scale_fmax:
664 msg = ('Scaled convergence criteria only implemented in series '
665 'with dynamic relaxation.')
666 raise ValueError(msg)
668 def set_positions(self, positions):
669 if not self.dynamic_relaxation:
670 return super().set_positions(positions)
672 n1 = 0
673 for i, image in enumerate(self.images[1:-1]):
674 if self.parallel:
675 msg = ('Dynamic relaxation does not work efficiently '
676 'when parallelizing over images. Try AutoNEB '
677 'routine for freezing images in parallel.')
678 raise ValueError(msg)
679 else:
680 forces_dyn = self._fmax_all(self.images)
681 if forces_dyn[i] < self.fmax:
682 n1 += self.natoms
683 else:
684 n2 = n1 + self.natoms
685 image.set_positions(positions[n1:n2])
686 n1 = n2
688 def _fmax_all(self, images):
689 """Store maximum force acting on each image in list. This is used in
690 the dynamic optimization routine in the set_positions() function."""
691 n = self.natoms
692 forces = self.get_forces()
693 fmax_images = [
694 np.sqrt((forces[n * i:n + n * i] ** 2).sum(axis=1)).max()
695 for i in range(self.nimages - 2)]
696 return fmax_images
698 def get_forces(self):
699 forces = super().get_forces()
700 if not self.dynamic_relaxation:
701 return forces
703 """Get NEB forces and scale the convergence criteria to focus
704 optimization on saddle point region. The keyword scale_fmax
705 determines the rate of convergence scaling."""
706 n = self.natoms
707 for i in range(self.nimages - 2):
708 n1 = n * i
709 n2 = n1 + n
710 force = np.sqrt((forces[n1:n2] ** 2.).sum(axis=1)).max()
711 n_imax = (self.imax - 1) * n # Image with highest energy.
713 positions = self.get_positions()
714 pos_imax = positions[n_imax:n_imax + n]
716 """Scale convergence criteria based on distance between an
717 image and the image with the highest potential energy."""
718 rel_pos = np.sqrt(((positions[n1:n2] - pos_imax) ** 2).sum())
719 if force < self.fmax * (1 + rel_pos * self.scale_fmax):
720 if i == self.imax - 1:
721 # Keep forces at saddle point for the log file.
722 pass
723 else:
724 # Set forces to zero before they are sent to optimizer.
725 forces[n1:n2, :] = 0
726 return forces
729def _check_deprecation(keyword, kwargs):
730 if keyword in kwargs:
731 warnings.warn(f'Keyword {keyword} of NEB is deprecated. '
732 'Please use the DyNEB class instead for dynamic '
733 'relaxation', FutureWarning)
736class NEB(DyNEB):
737 def __init__(self, images, k=0.1, climb=False, parallel=False,
738 remove_rotation_and_translation=False, world=None,
739 method='aseneb', allow_shared_calculator=False,
740 precon=None, **kwargs):
741 """Nudged elastic band.
743 Paper I:
745 G. Henkelman and H. Jonsson, Chem. Phys, 113, 9978 (2000).
746 :doi:`10.1063/1.1323224`
748 Paper II:
750 G. Henkelman, B. P. Uberuaga, and H. Jonsson, Chem. Phys,
751 113, 9901 (2000).
752 :doi:`10.1063/1.1329672`
754 Paper III:
756 E. L. Kolsbjerg, M. N. Groves, and B. Hammer, J. Chem. Phys,
757 145, 094107 (2016)
758 :doi:`10.1063/1.4961868`
760 Paper IV:
762 S. Makri, C. Ortner and J. R. Kermode, J. Chem. Phys.
763 150, 094109 (2019)
764 https://dx.doi.org/10.1063/1.5064465
766 images: list of Atoms objects
767 Images defining path from initial to final state.
768 k: float or list of floats
769 Spring constant(s) in eV/Ang. One number or one for each spring.
770 climb: bool
771 Use a climbing image (default is no climbing image).
772 parallel: bool
773 Distribute images over processors.
774 remove_rotation_and_translation: bool
775 TRUE actives NEB-TR for removing translation and
776 rotation during NEB. By default applied non-periodic
777 systems
778 method: string of method
779 Choice betweeen five methods:
781 * aseneb: standard ase NEB implementation
782 * improvedtangent: Paper I NEB implementation
783 * eb: Paper III full spring force implementation
784 * spline: Paper IV spline interpolation (supports precon)
785 * string: Paper IV string method (supports precon)
786 allow_shared_calculator: bool
787 Allow images to share the same calculator between them.
788 Incompatible with parallelisation over images.
789 precon: string, :class:`ase.optimize.precon.Precon` instance or list of
790 instances. If present, enable preconditioing as in Paper IV. This is
791 possible using the 'spline' or 'string' methods only.
792 Default is no preconditioning (precon=None), which is converted to
793 a list of :class:`ase.precon.precon.IdentityPrecon` instances.
794 """
795 for keyword in 'dynamic_relaxation', 'fmax', 'scale_fmax':
796 _check_deprecation(keyword, kwargs)
797 defaults = dict(dynamic_relaxation=False,
798 fmax=0.05,
799 scale_fmax=0.0)
800 defaults.update(kwargs)
801 # Only reason for separating BaseNEB/NEB is that we are
802 # deprecating dynamic_relaxation.
803 #
804 # We can turn BaseNEB into NEB once we get rid of the
805 # deprecated variables.
806 #
807 # Then we can also move DyNEB into ase.dyneb without cyclic imports.
808 # We can do that in ase-3.22 or 3.23.
809 super().__init__(
810 images, k=k, climb=climb, parallel=parallel,
811 remove_rotation_and_translation=remove_rotation_and_translation,
812 world=world, method=method,
813 allow_shared_calculator=allow_shared_calculator,
814 precon=precon,
815 **defaults)
818class NEBOptimizer(Optimizer):
819 """
820 This optimizer applies an adaptive ODE solver to a NEB
822 Details of the adaptive ODE solver are described in paper IV
823 """
825 def __init__(self,
826 neb,
827 restart=None, logfile='-', trajectory=None,
828 master=None,
829 append_trajectory=False,
830 method='ODE',
831 alpha=0.01,
832 verbose=0,
833 rtol=0.1,
834 C1=1e-2,
835 C2=2.0):
837 super().__init__(atoms=neb, restart=restart,
838 logfile=logfile, trajectory=trajectory,
839 master=master,
840 append_trajectory=append_trajectory,
841 force_consistent=False)
842 self.neb = neb
844 method = method.lower()
845 methods = ['ode', 'static', 'krylov']
846 if method not in methods:
847 raise ValueError(f'method must be one of {methods}')
848 self.method = method
850 self.alpha = alpha
851 self.verbose = verbose
852 self.rtol = rtol
853 self.C1 = C1
854 self.C2 = C2
856 def force_function(self, X):
857 positions = X.reshape((self.neb.nimages - 2) *
858 self.neb.natoms, 3)
859 self.neb.set_positions(positions)
860 forces = self.neb.get_forces().reshape(-1)
861 return forces
863 def get_residual(self, F=None, X=None):
864 return self.neb.get_residual()
866 def log(self):
867 fmax = self.get_residual()
868 T = time.localtime()
869 if self.logfile is not None:
870 name = f'{self.__class__.__name__}[{self.method}]'
871 if self.nsteps == 0:
872 args = (" " * len(name), "Step", "Time", "fmax")
873 msg = "%s %4s %8s %12s\n" % args
874 self.logfile.write(msg)
876 args = (name, self.nsteps, T[3], T[4], T[5], fmax)
877 msg = "%s: %3d %02d:%02d:%02d %12.4f\n" % args
878 self.logfile.write(msg)
879 self.logfile.flush()
881 def callback(self, X, F=None):
882 self.log()
883 self.call_observers()
884 self.nsteps += 1
886 def run_ode(self, fmax):
887 try:
888 ode12r(self.force_function,
889 self.neb.get_positions().reshape(-1),
890 fmax=fmax,
891 rtol=self.rtol,
892 C1=self.C1,
893 C2=self.C2,
894 steps=self.max_steps,
895 verbose=self.verbose,
896 callback=self.callback,
897 residual=self.get_residual)
898 return True
899 except OptimizerConvergenceError:
900 return False
902 def run_static(self, fmax):
903 X = self.neb.get_positions().reshape(-1)
904 for step in range(self.max_steps):
905 F = self.force_function(X)
906 if self.neb.get_residual() <= fmax:
907 return True
908 X += self.alpha * F
909 self.callback(X)
910 return False
912 def run(self, fmax=0.05, steps=None, method=None):
913 """
914 Optimize images to obtain the minimum energy path
916 Parameters
917 ----------
918 fmax - desired force tolerance
919 steps - maximum number of steps
920 """
921 if steps:
922 self.max_steps = steps
923 if method is None:
924 method = self.method
925 if method == 'ode':
926 return self.run_ode(fmax)
927 elif method == 'static':
928 return self.run_static(fmax)
929 else:
930 raise ValueError(f'unknown method: {self.method}')
933class IDPP(Calculator):
934 """Image dependent pair potential.
936 See:
937 Improved initial guess for minimum energy path calculations.
938 Søren Smidstrup, Andreas Pedersen, Kurt Stokbro and Hannes Jónsson
939 Chem. Phys. 140, 214106 (2014)
940 """
942 implemented_properties = ['energy', 'forces']
944 def __init__(self, target, mic):
945 Calculator.__init__(self)
946 self.target = target
947 self.mic = mic
949 def calculate(self, atoms, properties, system_changes):
950 Calculator.calculate(self, atoms, properties, system_changes)
952 P = atoms.get_positions()
953 d = []
954 D = []
955 for p in P:
956 Di = P - p
957 if self.mic:
958 Di, di = find_mic(Di, atoms.get_cell(), atoms.get_pbc())
959 else:
960 di = np.sqrt((Di ** 2).sum(1))
961 d.append(di)
962 D.append(Di)
963 d = np.array(d)
964 D = np.array(D)
966 dd = d - self.target
967 d.ravel()[::len(d) + 1] = 1 # avoid dividing by zero
968 d4 = d ** 4
969 e = 0.5 * (dd ** 2 / d4).sum()
970 f = -2 * ((dd * (1 - 2 * dd / d) / d ** 5)[..., np.newaxis] * D).sum(
971 0)
972 self.results = {'energy': e, 'forces': f}
975@deprecated("SingleCalculatorNEB is deprecated. "
976 "Please use NEB(allow_shared_calculator=True) instead.")
977class SingleCalculatorNEB(NEB):
978 def __init__(self, images, *args, **kwargs):
979 kwargs["allow_shared_calculator"] = True
980 super().__init__(images, *args, **kwargs)
983def interpolate(images, mic=False, interpolate_cell=False,
984 use_scaled_coord=False, apply_constraint=None):
985 """Given a list of images, linearly interpolate the positions of the
986 interior images.
988 mic: bool
989 Map movement into the unit cell by using the minimum image convention.
990 interpolate_cell: bool
991 Interpolate the three cell vectors linearly just like the atomic
992 positions. Not implemented for NEB calculations!
993 use_scaled_coord: bool
994 Use scaled/internal/fractional coordinates instead of real ones for the
995 interpolation. Not implemented for NEB calculations!
996 apply_constraint: bool
997 Controls if the constraints attached to the images
998 are ignored or applied when setting the interpolated positions.
999 Default value is None, in this case the resulting constrained positions
1000 (apply_constraint=True) are compared with unconstrained positions
1001 (apply_constraint=False), if the positions are not the same
1002 the user is required to specify the desired behaviour
1003 by setting up apply_constraint keyword argument to False or True.
1004 """
1005 if use_scaled_coord:
1006 pos1 = images[0].get_scaled_positions(wrap=mic)
1007 pos2 = images[-1].get_scaled_positions(wrap=mic)
1008 else:
1009 pos1 = images[0].get_positions()
1010 pos2 = images[-1].get_positions()
1011 d = pos2 - pos1
1012 if not use_scaled_coord and mic:
1013 d = find_mic(d, images[0].get_cell(), images[0].pbc)[0]
1014 d /= (len(images) - 1.0)
1015 if interpolate_cell:
1016 cell1 = images[0].get_cell()
1017 cell2 = images[-1].get_cell()
1018 cell_diff = cell2 - cell1
1019 cell_diff /= (len(images) - 1.0)
1020 for i in range(1, len(images) - 1):
1021 # first the new cell, otherwise scaled positions are wrong
1022 if interpolate_cell:
1023 images[i].set_cell(cell1 + i * cell_diff)
1024 new_pos = pos1 + i * d
1025 if use_scaled_coord:
1026 images[i].set_scaled_positions(new_pos)
1027 else:
1028 if apply_constraint is None:
1029 unconstrained_image = images[i].copy()
1030 unconstrained_image.set_positions(new_pos,
1031 apply_constraint=False)
1032 images[i].set_positions(new_pos, apply_constraint=True)
1033 try:
1034 np.testing.assert_allclose(unconstrained_image.positions,
1035 images[i].positions)
1036 except AssertionError:
1037 raise RuntimeError(f"Constraint(s) in image number {i} \n"
1038 f"affect the interpolation results.\n"
1039 "Please specify if you want to \n"
1040 "apply or ignore the constraints \n"
1041 "during the interpolation \n"
1042 "with apply_constraint argument.")
1043 else:
1044 images[i].set_positions(new_pos,
1045 apply_constraint=apply_constraint)
1048def idpp_interpolate(images, traj='idpp.traj', log='idpp.log', fmax=0.1,
1049 optimizer=MDMin, mic=False, steps=100):
1050 """Interpolate using the IDPP method. 'images' can either be a plain
1051 list of images or an NEB object (containing a list of images)."""
1052 if hasattr(images, 'interpolate'):
1053 neb = images
1054 else:
1055 neb = NEB(images)
1057 d1 = neb.images[0].get_all_distances(mic=mic)
1058 d2 = neb.images[-1].get_all_distances(mic=mic)
1059 d = (d2 - d1) / (neb.nimages - 1)
1060 real_calcs = []
1061 for i, image in enumerate(neb.images):
1062 real_calcs.append(image.calc)
1063 image.calc = IDPP(d1 + i * d, mic=mic)
1065 with optimizer(neb, trajectory=traj, logfile=log) as opt:
1066 opt.run(fmax=fmax, steps=steps)
1068 for image, calc in zip(neb.images, real_calcs):
1069 image.calc = calc
1072class NEBTools:
1073 """Class to make many of the common tools for NEB analysis available to
1074 the user. Useful for scripting the output of many jobs. Initialize with
1075 list of images which make up one or more band of the NEB relaxation."""
1077 def __init__(self, images):
1078 self.images = images
1080 @deprecated('NEBTools.get_fit() is deprecated. '
1081 'Please use ase.utils.forcecurve.fit_images(images).')
1082 def get_fit(self):
1083 return fit_images(self.images)
1085 def get_barrier(self, fit=True, raw=False):
1086 """Returns the barrier estimate from the NEB, along with the
1087 Delta E of the elementary reaction. If fit=True, the barrier is
1088 estimated based on the interpolated fit to the images; if
1089 fit=False, the barrier is taken as the maximum-energy image
1090 without interpolation. Set raw=True to get the raw energy of the
1091 transition state instead of the forward barrier."""
1092 forcefit = fit_images(self.images)
1093 energies = forcefit.energies
1094 fit_energies = forcefit.fit_energies
1095 dE = energies[-1] - energies[0]
1096 if fit:
1097 barrier = max(fit_energies)
1098 else:
1099 barrier = max(energies)
1100 if raw:
1101 barrier += self.images[0].get_potential_energy()
1102 return barrier, dE
1104 def get_fmax(self, **kwargs):
1105 """Returns fmax, as used by optimizers with NEB."""
1106 neb = NEB(self.images, **kwargs)
1107 forces = neb.get_forces()
1108 return np.sqrt((forces ** 2).sum(axis=1).max())
1110 def plot_band(self, ax=None):
1111 """Plots the NEB band on matplotlib axes object 'ax'. If ax=None
1112 returns a new figure object."""
1113 forcefit = fit_images(self.images)
1114 ax = forcefit.plot(ax=ax)
1115 return ax.figure
1117 def plot_bands(self, constant_x=False, constant_y=False,
1118 nimages=None, label='nebplots'):
1119 """Given a trajectory containing many steps of a NEB, makes
1120 plots of each band in the series in a single PDF.
1122 constant_x: bool
1123 Use the same x limits on all plots.
1124 constant_y: bool
1125 Use the same y limits on all plots.
1126 nimages: int
1127 Number of images per band. Guessed if not supplied.
1128 label: str
1129 Name for the output file. .pdf will be appended.
1130 """
1131 from matplotlib import pyplot
1132 from matplotlib.backends.backend_pdf import PdfPages
1133 if nimages is None:
1134 nimages = self._guess_nimages()
1135 nebsteps = len(self.images) // nimages
1136 if constant_x or constant_y:
1137 sys.stdout.write('Scaling axes.\n')
1138 sys.stdout.flush()
1139 # Plot all to one plot, then pull its x and y range.
1140 fig, ax = pyplot.subplots()
1141 for index in range(nebsteps):
1142 images = self.images[index * nimages:(index + 1) * nimages]
1143 NEBTools(images).plot_band(ax=ax)
1144 xlim = ax.get_xlim()
1145 ylim = ax.get_ylim()
1146 pyplot.close(fig) # Reference counting "bug" in pyplot.
1147 with PdfPages(label + '.pdf') as pdf:
1148 for index in range(nebsteps):
1149 sys.stdout.write('\rProcessing band {:10d} / {:10d}'
1150 .format(index, nebsteps))
1151 sys.stdout.flush()
1152 fig, ax = pyplot.subplots()
1153 images = self.images[index * nimages:(index + 1) * nimages]
1154 NEBTools(images).plot_band(ax=ax)
1155 if constant_x:
1156 ax.set_xlim(xlim)
1157 if constant_y:
1158 ax.set_ylim(ylim)
1159 pdf.savefig(fig)
1160 pyplot.close(fig) # Reference counting "bug" in pyplot.
1161 sys.stdout.write('\n')
1163 def _guess_nimages(self):
1164 """Attempts to guess the number of images per band from
1165 a trajectory, based solely on the repetition of the
1166 potential energy of images. This should also work for symmetric
1167 cases."""
1168 e_first = self.images[0].get_potential_energy()
1169 nimages = None
1170 for index, image in enumerate(self.images[1:], start=1):
1171 e = image.get_potential_energy()
1172 if e == e_first:
1173 # Need to check for symmetric case when e_first = e_last.
1174 try:
1175 e_next = self.images[index + 1].get_potential_energy()
1176 except IndexError:
1177 pass
1178 else:
1179 if e_next == e_first:
1180 nimages = index + 1 # Symmetric
1181 break
1182 nimages = index # Normal
1183 break
1184 if nimages is None:
1185 sys.stdout.write('Appears to be only one band in the images.\n')
1186 return len(self.images)
1187 # Sanity check that the energies of the last images line up too.
1188 e_last = self.images[nimages - 1].get_potential_energy()
1189 e_nextlast = self.images[2 * nimages - 1].get_potential_energy()
1190 if not (e_last == e_nextlast):
1191 raise RuntimeError('Could not guess number of images per band.')
1192 sys.stdout.write('Number of images per band guessed to be {:d}.\n'
1193 .format(nimages))
1194 return nimages
1197class NEBtools(NEBTools):
1198 @deprecated('NEBtools has been renamed; please use NEBTools.')
1199 def __init__(self, images):
1200 NEBTools.__init__(self, images)
1203@deprecated('Please use NEBTools.plot_band_from_fit.')
1204def plot_band_from_fit(s, E, Sfit, Efit, lines, ax=None):
1205 NEBTools.plot_band_from_fit(s, E, Sfit, Efit, lines, ax=None)
1208def fit0(*args, **kwargs):
1209 raise DeprecationWarning('fit0 is deprecated. Use `fit_raw` from '
1210 '`ase.utils.forcecurve` instead.')