Coverage for ase / mep / neb.py: 83.73%
633 statements
« prev ^ index » next coverage.py v7.13.3, created at 2026-02-04 10:20 +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)
909 def callback(self, X, F=None):
910 self.log()
911 self.call_observers()
912 self.nsteps += 1
914 def run_ode(self, fmax):
915 try:
916 ode12r(self.force_function,
917 self.neb.get_positions().reshape(-1),
918 fmax=fmax,
919 rtol=self.rtol,
920 C1=self.C1,
921 C2=self.C2,
922 steps=self.max_steps,
923 verbose=self.verbose,
924 callback=self.callback,
925 residual=self.get_residual)
926 return True
927 except OptimizerConvergenceError:
928 return False
930 def run_static(self, fmax):
931 X = self.neb.get_positions().reshape(-1)
932 for _ in range(self.max_steps):
933 F = self.force_function(X)
934 if self.neb.get_residual() <= fmax:
935 return True
936 X += self.alpha * F
937 self.callback(X)
938 return False
940 def run(self, fmax=0.05, steps=DEFAULT_MAX_STEPS, method=None):
941 """
942 Optimize images to obtain the minimum energy path
944 Parameters
945 ----------
946 fmax - desired force tolerance
947 steps - maximum number of steps
948 """
949 self.max_steps = steps
950 if method is None:
951 method = self.method
952 if method == 'ode':
953 return self.run_ode(fmax)
954 elif method == 'static':
955 return self.run_static(fmax)
956 else:
957 raise ValueError(f'unknown method: {self.method}')
960class IDPP(Calculator):
961 """Image dependent pair potential.
963 See:
964 Improved initial guess for minimum energy path calculations.
965 Søren Smidstrup, Andreas Pedersen, Kurt Stokbro and Hannes Jónsson
966 Chem. Phys. 140, 214106 (2014)
967 """
969 implemented_properties = ['energy', 'forces']
971 def __init__(self, target, mic):
972 Calculator.__init__(self)
973 self.target = target
974 self.mic = mic
976 def calculate(self, atoms, properties, system_changes):
977 Calculator.calculate(self, atoms, properties, system_changes)
979 P = atoms.get_positions()
980 d = []
981 D = []
982 for p in P:
983 Di = P - p
984 if self.mic:
985 Di, di = find_mic(Di, atoms.get_cell(), atoms.get_pbc())
986 else:
987 di = np.sqrt((Di ** 2).sum(1))
988 d.append(di)
989 D.append(Di)
990 d = np.array(d)
991 D = np.array(D)
993 dd = d - self.target
994 d.ravel()[::len(d) + 1] = 1 # avoid dividing by zero
995 d4 = d ** 4
996 e = 0.5 * (dd ** 2 / d4).sum()
997 f = -2 * ((dd * (1 - 2 * dd / d) / d ** 5)[..., np.newaxis] * D).sum(
998 0)
999 self.results = {'energy': e, 'forces': f}
1002@deprecated("SingleCalculatorNEB is deprecated. "
1003 "Please use NEB(allow_shared_calculator=True) instead.")
1004class SingleCalculatorNEB(NEB):
1005 """
1006 .. deprecated:: 3.23.0
1007 Please use ``NEB(allow_shared_calculator=True)`` instead
1008 """
1010 def __init__(self, images, *args, **kwargs):
1011 kwargs["allow_shared_calculator"] = True
1012 super().__init__(images, *args, **kwargs)
1015def interpolate(images, mic=False, interpolate_cell=False,
1016 use_scaled_coord=False, apply_constraint=None):
1017 """Given a list of images, linearly interpolate the positions of the
1018 interior images.
1020 mic: bool
1021 Map movement into the unit cell by using the minimum image convention.
1022 interpolate_cell: bool
1023 Interpolate the three cell vectors linearly just like the atomic
1024 positions. Not implemented for NEB calculations!
1025 use_scaled_coord: bool
1026 Use scaled/internal/fractional coordinates instead of real ones for the
1027 interpolation. Not implemented for NEB calculations!
1028 apply_constraint: bool
1029 Controls if the constraints attached to the images
1030 are ignored or applied when setting the interpolated positions.
1031 Default value is None, in this case the resulting constrained positions
1032 (apply_constraint=True) are compared with unconstrained positions
1033 (apply_constraint=False), if the positions are not the same
1034 the user is required to specify the desired behaviour
1035 by setting up apply_constraint keyword argument to False or True.
1036 """
1037 if use_scaled_coord:
1038 pos1 = images[0].get_scaled_positions(wrap=mic)
1039 pos2 = images[-1].get_scaled_positions(wrap=mic)
1040 else:
1041 pos1 = images[0].get_positions()
1042 pos2 = images[-1].get_positions()
1043 d = pos2 - pos1
1044 if not use_scaled_coord and mic:
1045 d = find_mic(d, images[0].get_cell(), images[0].pbc)[0]
1046 d /= (len(images) - 1.0)
1047 if interpolate_cell:
1048 cell1 = images[0].get_cell()
1049 cell2 = images[-1].get_cell()
1050 cell_diff = cell2 - cell1
1051 cell_diff /= (len(images) - 1.0)
1052 for i in range(1, len(images) - 1):
1053 # first the new cell, otherwise scaled positions are wrong
1054 if interpolate_cell:
1055 images[i].set_cell(cell1 + i * cell_diff)
1056 new_pos = pos1 + i * d
1057 if use_scaled_coord:
1058 images[i].set_scaled_positions(new_pos)
1059 else:
1060 if apply_constraint is None:
1061 unconstrained_image = images[i].copy()
1062 unconstrained_image.set_positions(new_pos,
1063 apply_constraint=False)
1064 images[i].set_positions(new_pos, apply_constraint=True)
1065 if not np.allclose(unconstrained_image.positions,
1066 images[i].positions):
1067 raise RuntimeError(f"Constraints in image {i}\n"
1068 "affect the interpolation results.\n"
1069 "Please specify if you want to\n"
1070 "apply or ignore the constraints\n"
1071 "during the interpolation\n"
1072 "with the apply_constraint argument.")
1073 else:
1074 images[i].set_positions(new_pos,
1075 apply_constraint=apply_constraint)
1078def idpp_interpolate(images, traj='idpp.traj', log='idpp.log', fmax=0.1,
1079 optimizer=MDMin, mic=False, steps=100):
1080 """Interpolate using the IDPP method. 'images' can either be a plain
1081 list of images or an NEB object (containing a list of images)."""
1082 if hasattr(images, 'interpolate'):
1083 neb = images
1084 else:
1085 neb = NEB(images)
1087 d1 = neb.images[0].get_all_distances(mic=mic)
1088 d2 = neb.images[-1].get_all_distances(mic=mic)
1089 d = (d2 - d1) / (neb.nimages - 1)
1090 real_calcs = []
1091 for i, image in enumerate(neb.images):
1092 real_calcs.append(image.calc)
1093 image.calc = IDPP(d1 + i * d, mic=mic)
1095 with optimizer(neb, trajectory=traj, logfile=log) as opt:
1096 opt.run(fmax=fmax, steps=steps)
1098 for image, calc in zip(neb.images, real_calcs):
1099 image.calc = calc
1102class NEBTools:
1103 """Class to make many of the common tools for NEB analysis available to
1104 the user. Useful for scripting the output of many jobs. Initialize with
1105 list of images which make up one or more band of the NEB relaxation."""
1107 def __init__(self, images):
1108 self.images = images
1110 @deprecated('NEBTools.get_fit() is deprecated. '
1111 'Please use ase.utils.forcecurve.fit_images(images).')
1112 def get_fit(self):
1113 """
1114 .. deprecated:: 3.23.0
1115 Please use ``ase.utils.forcecurve.fit_images(images)``
1116 """
1117 return fit_images(self.images)
1119 def get_barrier(self, fit=True, raw=False):
1120 """Returns the barrier estimate from the NEB, along with the
1121 Delta E of the elementary reaction. If fit=True, the barrier is
1122 estimated based on the interpolated fit to the images; if
1123 fit=False, the barrier is taken as the maximum-energy image
1124 without interpolation. Set raw=True to get the raw energy of the
1125 transition state instead of the forward barrier."""
1126 forcefit = fit_images(self.images)
1127 energies = forcefit.energies
1128 fit_energies = forcefit.fit_energies
1129 dE = energies[-1] - energies[0]
1130 if fit:
1131 barrier = max(fit_energies)
1132 else:
1133 barrier = max(energies)
1134 if raw:
1135 barrier += self.images[0].get_potential_energy()
1136 return barrier, dE
1138 def get_fmax(self, **kwargs):
1139 """Returns fmax, as used by optimizers with NEB."""
1140 neb = NEB(self.images, **kwargs)
1141 forces = neb.get_forces()
1142 return np.sqrt((forces ** 2).sum(axis=1).max())
1144 def plot_band(self, ax=None):
1145 """Plots the NEB band on matplotlib axes object 'ax'. If ax=None
1146 returns a new figure object."""
1147 forcefit = fit_images(self.images)
1148 ax = forcefit.plot(ax=ax)
1149 return ax.figure
1151 def plot_bands(self, constant_x=False, constant_y=False,
1152 nimages=None, label='nebplots'):
1153 """Given a trajectory containing many steps of a NEB, makes
1154 plots of each band in the series in a single PDF.
1156 constant_x: bool
1157 Use the same x limits on all plots.
1158 constant_y: bool
1159 Use the same y limits on all plots.
1160 nimages: int
1161 Number of images per band. Guessed if not supplied.
1162 label: str
1163 Name for the output file. .pdf will be appended.
1164 """
1165 from matplotlib import pyplot
1166 from matplotlib.backends.backend_pdf import PdfPages
1167 if nimages is None:
1168 nimages = self._guess_nimages()
1169 nebsteps = len(self.images) // nimages
1170 if constant_x or constant_y:
1171 sys.stdout.write('Scaling axes.\n')
1172 sys.stdout.flush()
1173 # Plot all to one plot, then pull its x and y range.
1174 fig, ax = pyplot.subplots()
1175 for index in range(nebsteps):
1176 images = self.images[index * nimages:(index + 1) * nimages]
1177 NEBTools(images).plot_band(ax=ax)
1178 xlim = ax.get_xlim()
1179 ylim = ax.get_ylim()
1180 pyplot.close(fig) # Reference counting "bug" in pyplot.
1181 with PdfPages(label + '.pdf') as pdf:
1182 for index in range(nebsteps):
1183 sys.stdout.write('\rProcessing band {:10d} / {:10d}'
1184 .format(index, nebsteps))
1185 sys.stdout.flush()
1186 fig, ax = pyplot.subplots()
1187 images = self.images[index * nimages:(index + 1) * nimages]
1188 NEBTools(images).plot_band(ax=ax)
1189 if constant_x:
1190 ax.set_xlim(xlim)
1191 if constant_y:
1192 ax.set_ylim(ylim)
1193 pdf.savefig(fig)
1194 pyplot.close(fig) # Reference counting "bug" in pyplot.
1195 sys.stdout.write('\n')
1197 def _guess_nimages(self):
1198 """Attempts to guess the number of images per band from
1199 a trajectory, based solely on the repetition of the
1200 potential energy of images. This should also work for symmetric
1201 cases."""
1202 e_first = self.images[0].get_potential_energy()
1203 nimages = None
1204 for index, image in enumerate(self.images[1:], start=1):
1205 e = image.get_potential_energy()
1206 if e == e_first:
1207 # Need to check for symmetric case when e_first = e_last.
1208 try:
1209 e_next = self.images[index + 1].get_potential_energy()
1210 except IndexError:
1211 pass
1212 else:
1213 if e_next == e_first:
1214 nimages = index + 1 # Symmetric
1215 break
1216 nimages = index # Normal
1217 break
1218 if nimages is None:
1219 sys.stdout.write('Appears to be only one band in the images.\n')
1220 return len(self.images)
1221 # Sanity check that the energies of the last images line up too.
1222 e_last = self.images[nimages - 1].get_potential_energy()
1223 e_nextlast = self.images[2 * nimages - 1].get_potential_energy()
1224 if e_last != e_nextlast:
1225 raise RuntimeError('Could not guess number of images per band.')
1226 sys.stdout.write('Number of images per band guessed to be {:d}.\n'
1227 .format(nimages))
1228 return nimages
1231class NEBtools(NEBTools):
1232 @deprecated('NEBtools has been renamed; please use NEBTools.')
1233 def __init__(self, images):
1234 """
1235 .. deprecated:: 3.23.0
1236 Please use :class:`~ase.mep.NEBTools`.
1237 """
1238 NEBTools.__init__(self, images)
1241@deprecated('Please use NEBTools.plot_band_from_fit.')
1242def plot_band_from_fit(s, E, Sfit, Efit, lines, ax=None):
1243 """
1244 .. deprecated:: 3.23.0
1245 Please use :meth:`NEBTools.plot_band_from_fit`.
1246 """
1247 NEBTools.plot_band_from_fit(s, E, Sfit, Efit, lines, ax=None)