Coverage for /builds/ase/ase/ase/mep/neb.py: 83.54%
638 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +0000
1# fmt: off
3import sys
4import threading
5import time
6import warnings
7from abc import ABC, abstractmethod
8from functools import cached_property
10import numpy as np
11from scipy.integrate import cumulative_trapezoid
12from scipy.interpolate import CubicSpline
14import ase.parallel
15from ase.build import minimize_rotation_and_translation
16from ase.calculators.calculator import Calculator
17from ase.calculators.singlepoint import SinglePointCalculator
18from ase.geometry import find_mic
19from ase.optimize import MDMin
20from ase.optimize.ode import ode12r
21from ase.optimize.optimize import DEFAULT_MAX_STEPS, Optimizer
22from ase.optimize.precon import Precon, PreconImages
23from ase.optimize.sciopt import OptimizerConvergenceError
24from ase.utils import deprecated
25from ase.utils.abc import Optimizable
26from ase.utils.forcecurve import fit_images
29class Spring:
30 def __init__(self, atoms1, atoms2, energy1, energy2, k):
31 self.atoms1 = atoms1
32 self.atoms2 = atoms2
33 self.energy1 = energy1
34 self.energy2 = energy2
35 self.k = k
37 def _find_mic(self):
38 pos1 = self.atoms1.get_positions()
39 pos2 = self.atoms2.get_positions()
40 # XXX If we want variable cells we will need to edit this.
41 mic, _ = find_mic(pos2 - pos1, self.atoms1.cell, self.atoms1.pbc)
42 return mic
44 @cached_property
45 def t(self):
46 return self._find_mic()
48 @cached_property
49 def nt(self):
50 return np.linalg.norm(self.t)
53class NEBState:
54 def __init__(self, neb, images, energies):
55 self.neb = neb
56 self.images = images
57 self.energies = energies
59 def spring(self, i):
60 return Spring(self.images[i], self.images[i + 1],
61 self.energies[i], self.energies[i + 1],
62 self.neb.k[i])
64 @cached_property
65 def imax(self):
66 return 1 + np.argsort(self.energies[1:-1])[-1]
68 @property
69 def emax(self):
70 return self.energies[self.imax]
72 @cached_property
73 def eqlength(self):
74 images = self.images
75 beeline = (images[self.neb.nimages - 1].get_positions() -
76 images[0].get_positions())
77 beelinelength = np.linalg.norm(beeline)
78 return beelinelength / (self.neb.nimages - 1)
80 @cached_property
81 def nimages(self):
82 return len(self.images)
84 @property
85 def precon(self):
86 return self.neb.precon
89class NEBMethod(ABC):
90 def __init__(self, neb):
91 self.neb = neb
93 @abstractmethod
94 def get_tangent(self, state, spring1, spring2, i):
95 ...
97 @abstractmethod
98 def add_image_force(self, state, tangential_force, tangent, imgforce,
99 spring1, spring2, i):
100 ...
102 def adjust_positions(self, positions):
103 return positions
106class ImprovedTangentMethod(NEBMethod):
107 """
108 Tangent estimates are improved according to Eqs. 8-11 in paper I.
109 Tangents are weighted at extrema to ensure smooth transitions between
110 the positive and negative tangents.
111 """
113 def get_tangent(self, state, spring1, spring2, i):
114 energies = state.energies
115 if energies[i + 1] > energies[i] > energies[i - 1]:
116 tangent = spring2.t.copy()
117 elif energies[i + 1] < energies[i] < energies[i - 1]:
118 tangent = spring1.t.copy()
119 else:
120 deltavmax = max(abs(energies[i + 1] - energies[i]),
121 abs(energies[i - 1] - energies[i]))
122 deltavmin = min(abs(energies[i + 1] - energies[i]),
123 abs(energies[i - 1] - energies[i]))
124 if energies[i + 1] > energies[i - 1]:
125 tangent = spring2.t * deltavmax + spring1.t * deltavmin
126 else:
127 tangent = spring2.t * deltavmin + spring1.t * deltavmax
128 # Normalize the tangent vector
129 tangent /= np.linalg.norm(tangent)
130 return tangent
132 def add_image_force(self, state, tangential_force, tangent, imgforce,
133 spring1, spring2, i):
134 imgforce -= tangential_force * tangent
135 # Improved parallel spring force (formula 12 of paper I)
136 imgforce += (spring2.nt * spring2.k - spring1.nt * spring1.k) * tangent
139class ASENEBMethod(NEBMethod):
140 """
141 Standard NEB implementation in ASE. The tangent of each image is
142 estimated from the spring closest to the saddle point in each
143 spring pair.
144 """
146 def get_tangent(self, state, spring1, spring2, i):
147 imax = self.neb.imax
148 if i < imax:
149 tangent = spring2.t
150 elif i > imax:
151 tangent = spring1.t
152 else:
153 tangent = spring1.t + spring2.t
154 return tangent
156 def add_image_force(self, state, tangential_force, tangent, imgforce,
157 spring1, spring2, i):
158 # Magnitude for normalizing. Ensure it is not 0
159 tangent_mag = np.vdot(tangent, tangent) or 1
160 factor = tangent / tangent_mag
161 imgforce -= tangential_force * factor
162 imgforce -= np.vdot(
163 spring1.t * spring1.k -
164 spring2.t * spring2.k, tangent) * factor
167class FullSpringMethod(NEBMethod):
168 """
169 Elastic band method. The full spring force is included.
170 """
172 def get_tangent(self, state, spring1, spring2, i):
173 # Tangents are bisections of spring-directions
174 # (formula C8 of paper III)
175 tangent = spring1.t / spring1.nt + spring2.t / spring2.nt
176 tangent /= np.linalg.norm(tangent)
177 return tangent
179 def add_image_force(self, state, tangential_force, tangent, imgforce,
180 spring1, spring2, i):
181 imgforce -= tangential_force * tangent
182 energies = state.energies
183 # Spring forces
184 # Eqs. C1, C5, C6 and C7 in paper III)
185 f1 = -(spring1.nt -
186 state.eqlength) * spring1.t / spring1.nt * spring1.k
187 f2 = (spring2.nt - state.eqlength) * spring2.t / spring2.nt * spring2.k
188 if self.neb.climb and abs(i - self.neb.imax) == 1:
189 deltavmax = max(abs(energies[i + 1] - energies[i]),
190 abs(energies[i - 1] - energies[i]))
191 deltavmin = min(abs(energies[i + 1] - energies[i]),
192 abs(energies[i - 1] - energies[i]))
193 imgforce += (f1 + f2) * deltavmin / deltavmax
194 else:
195 imgforce += f1 + f2
198class BaseSplineMethod(NEBMethod):
199 """
200 Base class for SplineNEB and String methods
202 Can optionally be preconditioned, as described in the following article:
204 S. Makri, C. Ortner and J. R. Kermode, J. Chem. Phys.
205 150, 094109 (2019)
206 https://dx.doi.org/10.1063/1.5064465
207 """
209 def __init__(self, neb):
210 NEBMethod.__init__(self, neb)
212 def get_tangent(self, state, spring1, spring2, i):
213 return state.precon.get_tangent(i)
215 def add_image_force(self, state, tangential_force, tangent, imgforce,
216 spring1, spring2, i):
217 # project out tangential component (Eqs 6 and 7 in Paper IV)
218 imgforce -= tangential_force * tangent
221class SplineMethod(BaseSplineMethod):
222 """
223 NEB using spline interpolation, plus optional preconditioning
224 """
226 def add_image_force(self, state, tangential_force, tangent, imgforce,
227 spring1, spring2, i):
228 super().add_image_force(state, tangential_force,
229 tangent, imgforce, spring1, spring2, i)
230 eta = state.precon.get_spring_force(i, spring1.k, spring2.k, tangent)
231 imgforce += eta
234class StringMethod(BaseSplineMethod):
235 """
236 String method using spline interpolation, plus optional preconditioning
237 """
239 def adjust_positions(self, positions):
240 # fit cubic spline to positions, reinterpolate to equispace images
241 # note this uses the preconditioned distance metric.
242 fit = self.neb.spline_fit(positions)
243 new_s = np.linspace(0.0, 1.0, self.neb.nimages)
244 new_positions = fit.x(new_s[1:-1]).reshape(-1, 3)
245 return new_positions
248def get_neb_method(neb, method):
249 if method == 'eb':
250 return FullSpringMethod(neb)
251 elif method == 'aseneb':
252 return ASENEBMethod(neb)
253 elif method == 'improvedtangent':
254 return ImprovedTangentMethod(neb)
255 elif method == 'spline':
256 return SplineMethod(neb)
257 elif method == 'string':
258 return StringMethod(neb)
259 else:
260 raise ValueError(f'Bad method: {method}')
263class NEBOptimizable(Optimizable):
264 def __init__(self, neb):
265 self.neb = neb
267 def get_gradient(self):
268 return self.neb.get_forces().ravel()
270 def get_value(self):
271 return self.neb.get_potential_energy()
273 def get_x(self):
274 return self.neb.get_positions().ravel()
276 def set_x(self, x):
277 self.neb.set_positions(x.reshape(-1, 3))
279 def ndofs(self):
280 return 3 * len(self.neb)
282 def iterimages(self):
283 return self.neb.iterimages()
286class BaseNEB:
287 def __init__(self, images, k=0.1, climb=False, parallel=False,
288 remove_rotation_and_translation=False, world=None,
289 method='aseneb', allow_shared_calculator=False, precon=None):
291 self.images = images
292 self.climb = climb
293 self.parallel = parallel
294 self.allow_shared_calculator = allow_shared_calculator
296 for img in images:
297 if len(img) != self.natoms:
298 raise ValueError('Images have different numbers of atoms')
299 if np.any(img.pbc != images[0].pbc):
300 raise ValueError('Images have different boundary conditions')
301 if np.any(img.get_atomic_numbers() !=
302 images[0].get_atomic_numbers()):
303 raise ValueError('Images have atoms in different orders')
304 # check periodic cell directions
305 cell_ok = True
306 for pbc, vc, vc0 in zip(img.pbc, img.cell, images[0].cell):
307 if pbc and np.any(np.abs(vc - vc0) > 1e-8):
308 cell_ok = False
309 if not cell_ok:
310 raise NotImplementedError(
311 "Variable cell in periodic directions "
312 "is not implemented yet for NEB")
314 self.emax = np.nan
316 self.remove_rotation_and_translation = remove_rotation_and_translation
318 if method in ['aseneb', 'eb', 'improvedtangent', 'spline', 'string']:
319 self.method = method
320 else:
321 raise NotImplementedError(method)
323 if precon is not None and method not in ['spline', 'string']:
324 raise NotImplementedError(f'no precon implemented: {method}')
325 self.precon = precon
327 self.neb_method = get_neb_method(self, method)
328 if isinstance(k, (float, int)):
329 k = [k] * (self.nimages - 1)
330 self.k = list(k)
332 if world is None:
333 world = ase.parallel.world
334 self.world = world
336 if parallel:
337 if self.allow_shared_calculator:
338 raise RuntimeError(
339 "Cannot use shared calculators in parallel in NEB.")
340 self.real_forces = None # ndarray of shape (nimages, natom, 3)
341 self.energies = None # ndarray of shape (nimages,)
342 self.residuals = None # ndarray of shape (nimages,)
344 def __ase_optimizable__(self):
345 return NEBOptimizable(self)
347 @property
348 def natoms(self):
349 return len(self.images[0])
351 @property
352 def nimages(self):
353 return len(self.images)
355 @staticmethod
356 def freeze_results_on_image(atoms: ase.Atoms,
357 **results_to_include):
358 atoms.calc = SinglePointCalculator(atoms=atoms, **results_to_include)
360 def interpolate(self, method='linear', mic=False, apply_constraint=None):
361 """Interpolate the positions of the interior images between the
362 initial state (image 0) and final state (image -1).
364 method: str
365 Method by which to interpolate: 'linear' or 'idpp'.
366 linear provides a standard straight-line interpolation, while
367 idpp uses an image-dependent pair potential.
368 mic: bool
369 Use the minimum-image convention when interpolating.
370 apply_constraint: bool
371 Controls if the constraints attached to the images
372 are ignored or applied when setting the interpolated positions.
373 Default value is None, in this case the resulting constrained
374 positions (apply_constraint=True) are compared with unconstrained
375 positions (apply_constraint=False),
376 if the positions are not the same
377 the user is required to specify the desired behaviour
378 by setting up apply_constraint keyword argument to False or True.
379 """
380 if self.remove_rotation_and_translation:
381 minimize_rotation_and_translation(self.images[0], self.images[-1])
383 interpolate(self.images, mic, apply_constraint=apply_constraint)
385 if method == 'idpp':
386 idpp_interpolate(images=self, traj=None, log=None, mic=mic)
388 @deprecated("Please use NEB's interpolate(method='idpp') method or "
389 "directly call the idpp_interpolate function from ase.mep")
390 def idpp_interpolate(self, traj='idpp.traj', log='idpp.log', fmax=0.1,
391 optimizer=MDMin, mic=False, steps=100):
392 """
393 .. deprecated:: 3.23.0
394 Please use :class:`~ase.mep.NEB`'s ``interpolate(method='idpp')``
395 method
396 """
397 idpp_interpolate(self, traj=traj, log=log, fmax=fmax,
398 optimizer=optimizer, mic=mic, steps=steps)
400 def get_positions(self):
401 positions = np.empty(((self.nimages - 2) * self.natoms, 3))
402 n1 = 0
403 for image in self.images[1:-1]:
404 n2 = n1 + self.natoms
405 positions[n1:n2] = image.get_positions()
406 n1 = n2
407 return positions
409 def set_positions(self, positions, adjust_positions=True):
410 if adjust_positions:
411 # optional reparameterisation step: some NEB methods need to adjust
412 # positions e.g. string method does this to equispace the images)
413 positions = self.neb_method.adjust_positions(positions)
414 n1 = 0
415 for image in self.images[1:-1]:
416 n2 = n1 + self.natoms
417 image.set_positions(positions[n1:n2])
418 n1 = n2
420 def get_forces(self):
421 """Evaluate and return the forces."""
422 images = self.images
424 if not self.allow_shared_calculator:
425 calculators = [image.calc for image in images
426 if image.calc is not None]
427 if len(set(calculators)) != len(calculators):
428 msg = ('One or more NEB images share the same calculator. '
429 'Each image must have its own calculator. '
430 'You may wish to use the ase.mep.SingleCalculatorNEB '
431 'class instead, although using separate calculators '
432 'is recommended.')
433 raise ValueError(msg)
435 forces = np.empty(((self.nimages - 2), self.natoms, 3))
436 energies = np.empty(self.nimages)
438 if self.remove_rotation_and_translation:
439 for i in range(1, self.nimages):
440 minimize_rotation_and_translation(images[i - 1], images[i])
442 if self.method != 'aseneb':
443 energies[0] = images[0].get_potential_energy()
444 energies[-1] = images[-1].get_potential_energy()
446 if not self.parallel:
447 # Do all images - one at a time:
448 for i in range(1, self.nimages - 1):
449 forces[i - 1] = images[i].get_forces()
450 energies[i] = images[i].get_potential_energy()
452 elif self.world.size == 1:
453 def run(image, energies, forces):
454 forces[:] = image.get_forces()
455 energies[:] = image.get_potential_energy()
457 threads = [threading.Thread(target=run,
458 args=(images[i],
459 energies[i:i + 1],
460 forces[i - 1:i]))
461 for i in range(1, self.nimages - 1)]
462 for thread in threads:
463 thread.start()
464 for thread in threads:
465 thread.join()
466 else:
467 # Parallelize over images:
468 i = self.world.rank * (self.nimages - 2) // self.world.size + 1
469 try:
470 forces[i - 1] = images[i].get_forces()
471 energies[i] = images[i].get_potential_energy()
472 except Exception:
473 # Make sure other images also fail:
474 error = self.world.sum(1.0)
475 raise
476 else:
477 error = self.world.sum(0.0)
478 if error:
479 raise RuntimeError('Parallel NEB failed!')
481 for i in range(1, self.nimages - 1):
482 root = (i - 1) * self.world.size // (self.nimages - 2)
483 self.world.broadcast(energies[i:i + 1], root)
484 self.world.broadcast(forces[i - 1], root)
486 # if this is the first force call, we need to build the preconditioners
487 if self.precon is None or isinstance(self.precon, (str, Precon, list)):
488 self.precon = PreconImages(self.precon, images)
490 # apply preconditioners to transform forces
491 # for the default IdentityPrecon this does not change their values
492 precon_forces = self.precon.apply(forces, index=slice(1, -1))
494 # Save for later use in iterimages:
495 self.energies = energies
496 self.real_forces = np.zeros((self.nimages, self.natoms, 3))
497 self.real_forces[1:-1] = forces
499 state = NEBState(self, images, energies)
501 # Can we get rid of self.energies, self.imax, self.emax etc.?
502 self.imax = state.imax
503 self.emax = state.emax
505 spring1 = state.spring(0)
507 self.residuals = []
508 for i in range(1, self.nimages - 1):
509 spring2 = state.spring(i)
510 tangent = self.neb_method.get_tangent(state, spring1, spring2, i)
512 # Get overlap between full PES-derived force and tangent
513 tangential_force = np.vdot(forces[i - 1], tangent)
515 # from now on we use the preconditioned forces (equal for precon=ID)
516 imgforce = precon_forces[i - 1]
518 if i == self.imax and self.climb:
519 """The climbing image, imax, is not affected by the spring
520 forces. This image feels the full PES-derived force,
521 but the tangential component is inverted:
522 see Eq. 5 in paper II."""
523 if self.method == 'aseneb':
524 tangent_mag = np.vdot(tangent, tangent) # For normalizing
525 imgforce -= 2 * tangential_force / tangent_mag * tangent
526 else:
527 imgforce -= 2 * tangential_force * tangent
528 else:
529 self.neb_method.add_image_force(state, tangential_force,
530 tangent, imgforce, spring1,
531 spring2, i)
532 # compute the residual - with ID precon, this is just max force
533 residual = self.precon.get_residual(i, imgforce)
534 self.residuals.append(residual)
536 spring1 = spring2
538 return precon_forces.reshape((-1, 3))
540 def get_residual(self):
541 """Return residual force along the band.
543 Typically this the maximum force component on any image. For
544 non-trivial preconditioners, the appropriate preconditioned norm
545 is used to compute the residual.
546 """
547 if self.residuals is None:
548 raise RuntimeError("get_residual() called before get_forces()")
549 return np.max(self.residuals)
551 def get_potential_energy(self):
552 """Return the maximum potential energy along the band."""
553 return self.emax
555 def set_calculators(self, calculators):
556 """Set new calculators to the images.
558 Parameters
559 ----------
560 calculators : Calculator / list(Calculator)
561 calculator(s) to attach to images
562 - single calculator, only if allow_shared_calculator=True
563 list of calculators if length:
564 - length nimages, set to all images
565 - length nimages-2, set to non-end images only
566 """
568 if not isinstance(calculators, list):
569 if self.allow_shared_calculator:
570 calculators = [calculators] * self.nimages
571 else:
572 raise RuntimeError("Cannot set shared calculator to NEB "
573 "with allow_shared_calculator=False")
575 n = len(calculators)
576 if n == self.nimages:
577 for i in range(self.nimages):
578 self.images[i].calc = calculators[i]
579 elif n == self.nimages - 2:
580 for i in range(1, self.nimages - 1):
581 self.images[i].calc = calculators[i - 1]
582 else:
583 raise RuntimeError(
584 'len(calculators)=%d does not fit to len(images)=%d'
585 % (n, self.nimages))
587 def __len__(self):
588 # Corresponds to number of optimizable degrees of freedom, i.e.
589 # virtual atom count for the optimization algorithm.
590 return (self.nimages - 2) * self.natoms
592 def iterimages(self):
593 # Allows trajectory to convert NEB into several images
594 for i, atoms in enumerate(self.images):
595 if i == 0 or i == self.nimages - 1:
596 yield atoms
597 else:
598 atoms = atoms.copy()
599 self.freeze_results_on_image(
600 atoms, energy=self.energies[i],
601 forces=self.real_forces[i])
603 yield atoms
605 def spline_fit(self, positions=None, norm='precon'):
606 """
607 Fit a cubic spline to this NEB
609 Args:
610 norm (str, optional): Norm to use: 'precon' (default) or 'euclidean'
612 Returns:
613 fit: ase.precon.precon.SplineFit instance
614 """
615 if norm == 'precon':
616 if self.precon is None or isinstance(self.precon, str):
617 self.precon = PreconImages(self.precon, self.images)
618 precon = self.precon
619 # if this is the first call, we need to build the preconditioners
620 elif norm == 'euclidean':
621 precon = PreconImages('ID', self.images)
622 else:
623 raise ValueError(f'unsupported norm {norm}')
624 return precon.spline_fit(positions)
626 def integrate_forces(self, spline_points=1000, bc_type='not-a-knot'):
627 """Use spline fit to integrate forces along MEP to approximate
628 energy differences using the virtual work approach.
630 Args:
631 spline_points (int, optional): Number of points. Defaults to 1000.
632 bc_type (str, optional): Boundary conditions, default 'not-a-knot'.
634 Returns:
635 s: reaction coordinate in range [0, 1], with `spline_points` entries
636 E: result of integrating forces, on the same grid as `s`.
637 F: projected forces along MEP
638 """
639 # note we use standard Euclidean rather than preconditioned norm
640 # to compute the virtual work
641 fit = self.spline_fit(norm='euclidean')
642 forces = np.array([image.get_forces().reshape(-1)
643 for image in self.images])
644 f = CubicSpline(fit.s, forces, bc_type=bc_type)
646 s = np.linspace(0.0, 1.0, spline_points, endpoint=True)
647 dE = f(s) * fit.dx_ds(s)
648 F = dE.sum(axis=1)
649 E = -cumulative_trapezoid(F, s, initial=0.0)
650 return s, E, F
653class DyNEB(BaseNEB):
654 def __init__(self, images, k=0.1, fmax=0.05, climb=False, parallel=False,
655 remove_rotation_and_translation=False, world=None,
656 dynamic_relaxation=True, scale_fmax=0., method='aseneb',
657 allow_shared_calculator=False, precon=None):
658 """
659 Subclass of NEB that allows for scaled and dynamic optimizations of
660 images. This method, which only works in series, does not perform
661 force calls on images that are below the convergence criterion.
662 The convergence criteria can be scaled with a displacement metric
663 to focus the optimization on the saddle point region.
665 'Scaled and Dynamic Optimizations of Nudged Elastic Bands',
666 P. Lindgren, G. Kastlunger and A. A. Peterson,
667 J. Chem. Theory Comput. 15, 11, 5787-5793 (2019).
669 dynamic_relaxation: bool
670 True skips images with forces below the convergence criterion.
671 This is updated after each force call; if a previously converged
672 image goes out of tolerance (due to spring adjustments between
673 the image and its neighbors), it will be optimized again.
674 False reverts to the default NEB implementation.
676 fmax: float
677 Must be identical to the fmax of the optimizer.
679 scale_fmax: float
680 Scale convergence criteria along band based on the distance between
681 an image and the image with the highest potential energy. This
682 keyword determines how rapidly the convergence criteria are scaled.
683 """
684 super().__init__(
685 images, k=k, climb=climb, parallel=parallel,
686 remove_rotation_and_translation=remove_rotation_and_translation,
687 world=world, method=method,
688 allow_shared_calculator=allow_shared_calculator, precon=precon)
689 self.fmax = fmax
690 self.dynamic_relaxation = dynamic_relaxation
691 self.scale_fmax = scale_fmax
693 if not self.dynamic_relaxation and self.scale_fmax:
694 msg = ('Scaled convergence criteria only implemented in series '
695 'with dynamic relaxation.')
696 raise ValueError(msg)
698 def set_positions(self, positions):
699 if not self.dynamic_relaxation:
700 return super().set_positions(positions)
702 n1 = 0
703 for i, image in enumerate(self.images[1:-1]):
704 if self.parallel:
705 msg = ('Dynamic relaxation does not work efficiently '
706 'when parallelizing over images. Try AutoNEB '
707 'routine for freezing images in parallel.')
708 raise ValueError(msg)
709 else:
710 forces_dyn = self._fmax_all(self.images)
711 if forces_dyn[i] < self.fmax:
712 n1 += self.natoms
713 else:
714 n2 = n1 + self.natoms
715 image.set_positions(positions[n1:n2])
716 n1 = n2
718 def _fmax_all(self, images):
719 """Store maximum force acting on each image in list. This is used in
720 the dynamic optimization routine in the set_positions() function."""
721 n = self.natoms
722 forces = self.get_forces()
723 fmax_images = [
724 np.sqrt((forces[n * i:n + n * i] ** 2).sum(axis=1)).max()
725 for i in range(self.nimages - 2)]
726 return fmax_images
728 def get_forces(self):
729 forces = super().get_forces()
730 if not self.dynamic_relaxation:
731 return forces
733 """Get NEB forces and scale the convergence criteria to focus
734 optimization on saddle point region. The keyword scale_fmax
735 determines the rate of convergence scaling."""
736 n = self.natoms
737 for i in range(self.nimages - 2):
738 n1 = n * i
739 n2 = n1 + n
740 force = np.sqrt((forces[n1:n2] ** 2.).sum(axis=1)).max()
741 n_imax = (self.imax - 1) * n # Image with highest energy.
743 positions = self.get_positions()
744 pos_imax = positions[n_imax:n_imax + n]
746 """Scale convergence criteria based on distance between an
747 image and the image with the highest potential energy."""
748 rel_pos = np.sqrt(((positions[n1:n2] - pos_imax) ** 2).sum())
749 if force < self.fmax * (1 + rel_pos * self.scale_fmax):
750 if i == self.imax - 1:
751 # Keep forces at saddle point for the log file.
752 pass
753 else:
754 # Set forces to zero before they are sent to optimizer.
755 forces[n1:n2, :] = 0
756 return forces
759def _check_deprecation(keyword, kwargs):
760 if keyword in kwargs:
761 warnings.warn(f'Keyword {keyword} of NEB is deprecated. '
762 'Please use the DyNEB class instead for dynamic '
763 'relaxation', FutureWarning)
766class NEB(DyNEB):
767 def __init__(self, images, k=0.1, climb=False, parallel=False,
768 remove_rotation_and_translation=False, world=None,
769 method='aseneb', allow_shared_calculator=False,
770 precon=None, **kwargs):
771 """Nudged elastic band.
773 Paper I:
775 G. Henkelman and H. Jonsson, Chem. Phys, 113, 9978 (2000).
776 :doi:`10.1063/1.1323224`
778 Paper II:
780 G. Henkelman, B. P. Uberuaga, and H. Jonsson, Chem. Phys,
781 113, 9901 (2000).
782 :doi:`10.1063/1.1329672`
784 Paper III:
786 E. L. Kolsbjerg, M. N. Groves, and B. Hammer, J. Chem. Phys,
787 145, 094107 (2016)
788 :doi:`10.1063/1.4961868`
790 Paper IV:
792 S. Makri, C. Ortner and J. R. Kermode, J. Chem. Phys.
793 150, 094109 (2019)
794 https://dx.doi.org/10.1063/1.5064465
796 images: list of Atoms objects
797 Images defining path from initial to final state.
798 k: float or list of floats
799 Spring constant(s) in eV/Ang. One number or one for each spring.
800 climb: bool
801 Use a climbing image (default is no climbing image).
802 parallel: bool
803 Distribute images over processors.
804 remove_rotation_and_translation: bool
805 TRUE actives NEB-TR for removing translation and
806 rotation during NEB. By default applied non-periodic
807 systems
808 method: string of method
809 Choice betweeen five methods:
811 * aseneb: standard ase NEB implementation
812 * improvedtangent: Paper I NEB implementation
813 * eb: Paper III full spring force implementation
814 * spline: Paper IV spline interpolation (supports precon)
815 * string: Paper IV string method (supports precon)
816 allow_shared_calculator: bool
817 Allow images to share the same calculator between them.
818 Incompatible with parallelisation over images.
819 precon: string, :class:`ase.optimize.precon.Precon` instance or list of
820 instances. If present, enable preconditioing as in Paper IV. This is
821 possible using the 'spline' or 'string' methods only.
822 Default is no preconditioning (precon=None), which is converted to
823 a list of :class:`ase.precon.precon.IdentityPrecon` instances.
824 """
825 for keyword in 'dynamic_relaxation', 'fmax', 'scale_fmax':
826 _check_deprecation(keyword, kwargs)
827 defaults = dict(dynamic_relaxation=False,
828 fmax=0.05,
829 scale_fmax=0.0)
830 defaults.update(kwargs)
831 # Only reason for separating BaseNEB/NEB is that we are
832 # deprecating dynamic_relaxation.
833 #
834 # We can turn BaseNEB into NEB once we get rid of the
835 # deprecated variables.
836 #
837 # Then we can also move DyNEB into ase.dyneb without cyclic imports.
838 # We can do that in ase-3.22 or 3.23.
839 super().__init__(
840 images, k=k, climb=climb, parallel=parallel,
841 remove_rotation_and_translation=remove_rotation_and_translation,
842 world=world, method=method,
843 allow_shared_calculator=allow_shared_calculator,
844 precon=precon,
845 **defaults)
848class NEBOptimizer(Optimizer):
849 """
850 This optimizer applies an adaptive ODE solver to a NEB
852 Details of the adaptive ODE solver are described in paper IV
853 """
855 def __init__(self,
856 neb,
857 restart=None, logfile='-', trajectory=None,
858 master=None,
859 append_trajectory=False,
860 method='ODE',
861 alpha=0.01,
862 verbose=0,
863 rtol=0.1,
864 C1=1e-2,
865 C2=2.0):
867 super().__init__(atoms=neb, restart=restart,
868 logfile=logfile, trajectory=trajectory,
869 master=master,
870 append_trajectory=append_trajectory)
871 self.neb = neb
873 method = method.lower()
874 methods = ['ode', 'static', 'krylov']
875 if method not in methods:
876 raise ValueError(f'method must be one of {methods}')
877 self.method = method
879 self.alpha = alpha
880 self.verbose = verbose
881 self.rtol = rtol
882 self.C1 = C1
883 self.C2 = C2
885 def force_function(self, X):
886 positions = X.reshape((self.neb.nimages - 2) *
887 self.neb.natoms, 3)
888 self.neb.set_positions(positions)
889 forces = self.neb.get_forces().reshape(-1)
890 return forces
892 def get_residual(self, F=None, X=None):
893 return self.neb.get_residual()
895 def log(self):
896 fmax = self.get_residual()
897 T = time.localtime()
898 if self.logfile is not None:
899 name = f'{self.__class__.__name__}[{self.method}]'
900 if self.nsteps == 0:
901 args = (" " * len(name), "Step", "Time", "fmax")
902 msg = "%s %4s %8s %12s\n" % args
903 self.logfile.write(msg)
905 args = (name, self.nsteps, T[3], T[4], T[5], fmax)
906 msg = "%s: %3d %02d:%02d:%02d %12.4f\n" % args
907 self.logfile.write(msg)
908 self.logfile.flush()
910 def callback(self, X, F=None):
911 self.log()
912 self.call_observers()
913 self.nsteps += 1
915 def run_ode(self, fmax):
916 try:
917 ode12r(self.force_function,
918 self.neb.get_positions().reshape(-1),
919 fmax=fmax,
920 rtol=self.rtol,
921 C1=self.C1,
922 C2=self.C2,
923 steps=self.max_steps,
924 verbose=self.verbose,
925 callback=self.callback,
926 residual=self.get_residual)
927 return True
928 except OptimizerConvergenceError:
929 return False
931 def run_static(self, fmax):
932 X = self.neb.get_positions().reshape(-1)
933 for _ in range(self.max_steps):
934 F = self.force_function(X)
935 if self.neb.get_residual() <= fmax:
936 return True
937 X += self.alpha * F
938 self.callback(X)
939 return False
941 def run(self, fmax=0.05, steps=DEFAULT_MAX_STEPS, method=None):
942 """
943 Optimize images to obtain the minimum energy path
945 Parameters
946 ----------
947 fmax - desired force tolerance
948 steps - maximum number of steps
949 """
950 self.max_steps = steps
951 if method is None:
952 method = self.method
953 if method == 'ode':
954 return self.run_ode(fmax)
955 elif method == 'static':
956 return self.run_static(fmax)
957 else:
958 raise ValueError(f'unknown method: {self.method}')
961class IDPP(Calculator):
962 """Image dependent pair potential.
964 See:
965 Improved initial guess for minimum energy path calculations.
966 Søren Smidstrup, Andreas Pedersen, Kurt Stokbro and Hannes Jónsson
967 Chem. Phys. 140, 214106 (2014)
968 """
970 implemented_properties = ['energy', 'forces']
972 def __init__(self, target, mic):
973 Calculator.__init__(self)
974 self.target = target
975 self.mic = mic
977 def calculate(self, atoms, properties, system_changes):
978 Calculator.calculate(self, atoms, properties, system_changes)
980 P = atoms.get_positions()
981 d = []
982 D = []
983 for p in P:
984 Di = P - p
985 if self.mic:
986 Di, di = find_mic(Di, atoms.get_cell(), atoms.get_pbc())
987 else:
988 di = np.sqrt((Di ** 2).sum(1))
989 d.append(di)
990 D.append(Di)
991 d = np.array(d)
992 D = np.array(D)
994 dd = d - self.target
995 d.ravel()[::len(d) + 1] = 1 # avoid dividing by zero
996 d4 = d ** 4
997 e = 0.5 * (dd ** 2 / d4).sum()
998 f = -2 * ((dd * (1 - 2 * dd / d) / d ** 5)[..., np.newaxis] * D).sum(
999 0)
1000 self.results = {'energy': e, 'forces': f}
1003@deprecated("SingleCalculatorNEB is deprecated. "
1004 "Please use NEB(allow_shared_calculator=True) instead.")
1005class SingleCalculatorNEB(NEB):
1006 """
1007 .. deprecated:: 3.23.0
1008 Please use ``NEB(allow_shared_calculator=True)`` instead
1009 """
1011 def __init__(self, images, *args, **kwargs):
1012 kwargs["allow_shared_calculator"] = True
1013 super().__init__(images, *args, **kwargs)
1016def interpolate(images, mic=False, interpolate_cell=False,
1017 use_scaled_coord=False, apply_constraint=None):
1018 """Given a list of images, linearly interpolate the positions of the
1019 interior images.
1021 mic: bool
1022 Map movement into the unit cell by using the minimum image convention.
1023 interpolate_cell: bool
1024 Interpolate the three cell vectors linearly just like the atomic
1025 positions. Not implemented for NEB calculations!
1026 use_scaled_coord: bool
1027 Use scaled/internal/fractional coordinates instead of real ones for the
1028 interpolation. Not implemented for NEB calculations!
1029 apply_constraint: bool
1030 Controls if the constraints attached to the images
1031 are ignored or applied when setting the interpolated positions.
1032 Default value is None, in this case the resulting constrained positions
1033 (apply_constraint=True) are compared with unconstrained positions
1034 (apply_constraint=False), if the positions are not the same
1035 the user is required to specify the desired behaviour
1036 by setting up apply_constraint keyword argument to False or True.
1037 """
1038 if use_scaled_coord:
1039 pos1 = images[0].get_scaled_positions(wrap=mic)
1040 pos2 = images[-1].get_scaled_positions(wrap=mic)
1041 else:
1042 pos1 = images[0].get_positions()
1043 pos2 = images[-1].get_positions()
1044 d = pos2 - pos1
1045 if not use_scaled_coord and mic:
1046 d = find_mic(d, images[0].get_cell(), images[0].pbc)[0]
1047 d /= (len(images) - 1.0)
1048 if interpolate_cell:
1049 cell1 = images[0].get_cell()
1050 cell2 = images[-1].get_cell()
1051 cell_diff = cell2 - cell1
1052 cell_diff /= (len(images) - 1.0)
1053 for i in range(1, len(images) - 1):
1054 # first the new cell, otherwise scaled positions are wrong
1055 if interpolate_cell:
1056 images[i].set_cell(cell1 + i * cell_diff)
1057 new_pos = pos1 + i * d
1058 if use_scaled_coord:
1059 images[i].set_scaled_positions(new_pos)
1060 else:
1061 if apply_constraint is None:
1062 unconstrained_image = images[i].copy()
1063 unconstrained_image.set_positions(new_pos,
1064 apply_constraint=False)
1065 images[i].set_positions(new_pos, apply_constraint=True)
1066 if not np.allclose(unconstrained_image.positions,
1067 images[i].positions):
1068 raise RuntimeError(f"Constraints in image {i}\n"
1069 "affect the interpolation results.\n"
1070 "Please specify if you want to\n"
1071 "apply or ignore the constraints\n"
1072 "during the interpolation\n"
1073 "with the apply_constraint argument.")
1074 else:
1075 images[i].set_positions(new_pos,
1076 apply_constraint=apply_constraint)
1079def idpp_interpolate(images, traj='idpp.traj', log='idpp.log', fmax=0.1,
1080 optimizer=MDMin, mic=False, steps=100):
1081 """Interpolate using the IDPP method. 'images' can either be a plain
1082 list of images or an NEB object (containing a list of images)."""
1083 if hasattr(images, 'interpolate'):
1084 neb = images
1085 else:
1086 neb = NEB(images)
1088 d1 = neb.images[0].get_all_distances(mic=mic)
1089 d2 = neb.images[-1].get_all_distances(mic=mic)
1090 d = (d2 - d1) / (neb.nimages - 1)
1091 real_calcs = []
1092 for i, image in enumerate(neb.images):
1093 real_calcs.append(image.calc)
1094 image.calc = IDPP(d1 + i * d, mic=mic)
1096 with optimizer(neb, trajectory=traj, logfile=log) as opt:
1097 opt.run(fmax=fmax, steps=steps)
1099 for image, calc in zip(neb.images, real_calcs):
1100 image.calc = calc
1103class NEBTools:
1104 """Class to make many of the common tools for NEB analysis available to
1105 the user. Useful for scripting the output of many jobs. Initialize with
1106 list of images which make up one or more band of the NEB relaxation."""
1108 def __init__(self, images):
1109 self.images = images
1111 @deprecated('NEBTools.get_fit() is deprecated. '
1112 'Please use ase.utils.forcecurve.fit_images(images).')
1113 def get_fit(self):
1114 """
1115 .. deprecated:: 3.23.0
1116 Please use ``ase.utils.forcecurve.fit_images(images)``
1117 """
1118 return fit_images(self.images)
1120 def get_barrier(self, fit=True, raw=False):
1121 """Returns the barrier estimate from the NEB, along with the
1122 Delta E of the elementary reaction. If fit=True, the barrier is
1123 estimated based on the interpolated fit to the images; if
1124 fit=False, the barrier is taken as the maximum-energy image
1125 without interpolation. Set raw=True to get the raw energy of the
1126 transition state instead of the forward barrier."""
1127 forcefit = fit_images(self.images)
1128 energies = forcefit.energies
1129 fit_energies = forcefit.fit_energies
1130 dE = energies[-1] - energies[0]
1131 if fit:
1132 barrier = max(fit_energies)
1133 else:
1134 barrier = max(energies)
1135 if raw:
1136 barrier += self.images[0].get_potential_energy()
1137 return barrier, dE
1139 def get_fmax(self, **kwargs):
1140 """Returns fmax, as used by optimizers with NEB."""
1141 neb = NEB(self.images, **kwargs)
1142 forces = neb.get_forces()
1143 return np.sqrt((forces ** 2).sum(axis=1).max())
1145 def plot_band(self, ax=None):
1146 """Plots the NEB band on matplotlib axes object 'ax'. If ax=None
1147 returns a new figure object."""
1148 forcefit = fit_images(self.images)
1149 ax = forcefit.plot(ax=ax)
1150 return ax.figure
1152 def plot_bands(self, constant_x=False, constant_y=False,
1153 nimages=None, label='nebplots'):
1154 """Given a trajectory containing many steps of a NEB, makes
1155 plots of each band in the series in a single PDF.
1157 constant_x: bool
1158 Use the same x limits on all plots.
1159 constant_y: bool
1160 Use the same y limits on all plots.
1161 nimages: int
1162 Number of images per band. Guessed if not supplied.
1163 label: str
1164 Name for the output file. .pdf will be appended.
1165 """
1166 from matplotlib import pyplot
1167 from matplotlib.backends.backend_pdf import PdfPages
1168 if nimages is None:
1169 nimages = self._guess_nimages()
1170 nebsteps = len(self.images) // nimages
1171 if constant_x or constant_y:
1172 sys.stdout.write('Scaling axes.\n')
1173 sys.stdout.flush()
1174 # Plot all to one plot, then pull its x and y range.
1175 fig, ax = pyplot.subplots()
1176 for index in range(nebsteps):
1177 images = self.images[index * nimages:(index + 1) * nimages]
1178 NEBTools(images).plot_band(ax=ax)
1179 xlim = ax.get_xlim()
1180 ylim = ax.get_ylim()
1181 pyplot.close(fig) # Reference counting "bug" in pyplot.
1182 with PdfPages(label + '.pdf') as pdf:
1183 for index in range(nebsteps):
1184 sys.stdout.write('\rProcessing band {:10d} / {:10d}'
1185 .format(index, nebsteps))
1186 sys.stdout.flush()
1187 fig, ax = pyplot.subplots()
1188 images = self.images[index * nimages:(index + 1) * nimages]
1189 NEBTools(images).plot_band(ax=ax)
1190 if constant_x:
1191 ax.set_xlim(xlim)
1192 if constant_y:
1193 ax.set_ylim(ylim)
1194 pdf.savefig(fig)
1195 pyplot.close(fig) # Reference counting "bug" in pyplot.
1196 sys.stdout.write('\n')
1198 def _guess_nimages(self):
1199 """Attempts to guess the number of images per band from
1200 a trajectory, based solely on the repetition of the
1201 potential energy of images. This should also work for symmetric
1202 cases."""
1203 e_first = self.images[0].get_potential_energy()
1204 nimages = None
1205 for index, image in enumerate(self.images[1:], start=1):
1206 e = image.get_potential_energy()
1207 if e == e_first:
1208 # Need to check for symmetric case when e_first = e_last.
1209 try:
1210 e_next = self.images[index + 1].get_potential_energy()
1211 except IndexError:
1212 pass
1213 else:
1214 if e_next == e_first:
1215 nimages = index + 1 # Symmetric
1216 break
1217 nimages = index # Normal
1218 break
1219 if nimages is None:
1220 sys.stdout.write('Appears to be only one band in the images.\n')
1221 return len(self.images)
1222 # Sanity check that the energies of the last images line up too.
1223 e_last = self.images[nimages - 1].get_potential_energy()
1224 e_nextlast = self.images[2 * nimages - 1].get_potential_energy()
1225 if e_last != e_nextlast:
1226 raise RuntimeError('Could not guess number of images per band.')
1227 sys.stdout.write('Number of images per band guessed to be {:d}.\n'
1228 .format(nimages))
1229 return nimages
1232class NEBtools(NEBTools):
1233 @deprecated('NEBtools has been renamed; please use NEBTools.')
1234 def __init__(self, images):
1235 """
1236 .. deprecated:: 3.23.0
1237 Please use :class:`~ase.mep.NEBTools`.
1238 """
1239 NEBTools.__init__(self, images)
1242@deprecated('Please use NEBTools.plot_band_from_fit.')
1243def plot_band_from_fit(s, E, Sfit, Efit, lines, ax=None):
1244 """
1245 .. deprecated:: 3.23.0
1246 Please use :meth:`NEBTools.plot_band_from_fit`.
1247 """
1248 NEBTools.plot_band_from_fit(s, E, Sfit, Efit, lines, ax=None)