Coverage for ase / mep / neb.py: 83.15%
641 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-03-30 08:22 +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 norm = np.linalg.norm(tangent)
130 tangent /= norm if norm > 0 else 1.0
131 return tangent
133 def add_image_force(self, state, tangential_force, tangent, imgforce,
134 spring1, spring2, i):
135 imgforce -= tangential_force * tangent
136 # Improved parallel spring force (formula 12 of paper I)
137 imgforce += (spring2.nt * spring2.k - spring1.nt * spring1.k) * tangent
140class ASENEBMethod(NEBMethod):
141 """
142 Standard NEB implementation in ASE. The tangent of each image is
143 estimated from the spring closest to the saddle point in each
144 spring pair.
145 """
147 def get_tangent(self, state, spring1, spring2, i):
148 imax = self.neb.imax
149 if i < imax:
150 tangent = spring2.t
151 elif i > imax:
152 tangent = spring1.t
153 else:
154 tangent = spring1.t + spring2.t
155 return tangent
157 def add_image_force(self, state, tangential_force, tangent, imgforce,
158 spring1, spring2, i):
159 # Magnitude for normalizing. Ensure it is not 0
160 tangent_mag = np.vdot(tangent, tangent) or 1
161 factor = tangent / tangent_mag
162 imgforce -= tangential_force * factor
163 imgforce -= np.vdot(
164 spring1.t * spring1.k -
165 spring2.t * spring2.k, tangent) * factor
168class FullSpringMethod(NEBMethod):
169 """
170 Elastic band method. The full spring force is included.
171 """
173 def get_tangent(self, state, spring1, spring2, i):
174 # Tangents are bisections of spring-directions
175 # (formula C8 of paper III)
176 tangent = spring1.t / spring1.nt + spring2.t / spring2.nt
177 tangent /= np.linalg.norm(tangent)
178 return tangent
180 def add_image_force(self, state, tangential_force, tangent, imgforce,
181 spring1, spring2, i):
182 imgforce -= tangential_force * tangent
183 energies = state.energies
184 # Spring forces
185 # Eqs. C1, C5, C6 and C7 in paper III)
186 f1 = -(spring1.nt -
187 state.eqlength) * spring1.t / spring1.nt * spring1.k
188 f2 = (spring2.nt - state.eqlength) * spring2.t / spring2.nt * spring2.k
189 if self.neb.climb and abs(i - self.neb.imax) == 1:
190 deltavmax = max(abs(energies[i + 1] - energies[i]),
191 abs(energies[i - 1] - energies[i]))
192 deltavmin = min(abs(energies[i + 1] - energies[i]),
193 abs(energies[i - 1] - energies[i]))
194 imgforce += (f1 + f2) * deltavmin / deltavmax
195 else:
196 imgforce += f1 + f2
199class BaseSplineMethod(NEBMethod):
200 """
201 Base class for SplineNEB and String methods
203 Can optionally be preconditioned, as described in the following article:
205 S. Makri, C. Ortner and J. R. Kermode, J. Chem. Phys.
206 150, 094109 (2019)
207 https://dx.doi.org/10.1063/1.5064465
208 """
210 def __init__(self, neb):
211 NEBMethod.__init__(self, neb)
213 def get_tangent(self, state, spring1, spring2, i):
214 return state.precon.get_tangent(i)
216 def add_image_force(self, state, tangential_force, tangent, imgforce,
217 spring1, spring2, i):
218 # project out tangential component (Eqs 6 and 7 in Paper IV)
219 imgforce -= tangential_force * tangent
222class SplineMethod(BaseSplineMethod):
223 """
224 NEB using spline interpolation, plus optional preconditioning
225 """
227 def add_image_force(self, state, tangential_force, tangent, imgforce,
228 spring1, spring2, i):
229 super().add_image_force(state, tangential_force,
230 tangent, imgforce, spring1, spring2, i)
231 eta = state.precon.get_spring_force(i, spring1.k, spring2.k, tangent)
232 imgforce += eta
235class StringMethod(BaseSplineMethod):
236 """
237 String method using spline interpolation, plus optional preconditioning
238 """
240 def adjust_positions(self, positions):
241 # fit cubic spline to positions, reinterpolate to equispace images
242 # note this uses the preconditioned distance metric.
243 fit = self.neb.spline_fit(positions)
244 new_s = np.linspace(0.0, 1.0, self.neb.nimages)
245 new_positions = fit.x(new_s[1:-1]).reshape(-1, 3)
246 return new_positions
249def get_neb_method(neb, method):
250 if method == 'eb':
251 return FullSpringMethod(neb)
252 elif method == 'aseneb':
253 return ASENEBMethod(neb)
254 elif method == 'improvedtangent':
255 return ImprovedTangentMethod(neb)
256 elif method == 'spline':
257 return SplineMethod(neb)
258 elif method == 'string':
259 return StringMethod(neb)
260 else:
261 raise ValueError(f'Bad method: {method}')
264class NEBOptimizable(Optimizable):
265 def __init__(self, neb):
266 self.neb = neb
268 def get_gradient(self):
269 return -self.neb.get_forces().ravel()
271 def get_value(self):
272 return self.neb.get_potential_energy()
274 def get_x(self):
275 return self.neb.get_positions().ravel()
277 def set_x(self, x):
278 self.neb.set_positions(x.reshape(-1, 3))
280 def ndofs(self):
281 return 3 * len(self.neb)
283 def iterimages(self):
284 return self.neb.iterimages()
287class BaseNEB:
288 def __init__(
289 self,
290 images,
291 k=0.1,
292 climb=False,
293 parallel=False,
294 remove_rotation_and_translation=False,
295 world=None,
296 method=None,
297 allow_shared_calculator=False,
298 precon=None,
299 ):
301 self.images = images
302 self.climb = climb
303 self.parallel = parallel
304 self.allow_shared_calculator = allow_shared_calculator
306 for img in images:
307 if len(img) != self.natoms:
308 raise ValueError('Images have different numbers of atoms')
309 if np.any(img.pbc != images[0].pbc):
310 raise ValueError('Images have different boundary conditions')
311 if np.any(img.get_atomic_numbers() !=
312 images[0].get_atomic_numbers()):
313 raise ValueError('Images have atoms in different orders')
314 # check periodic cell directions
315 cell_ok = True
316 for pbc, vc, vc0 in zip(img.pbc, img.cell, images[0].cell):
317 if pbc and np.any(np.abs(vc - vc0) > 1e-8):
318 cell_ok = False
319 if not cell_ok:
320 raise NotImplementedError(
321 "Variable cell in periodic directions "
322 "is not implemented yet for NEB")
324 self.emax = np.nan
326 self.remove_rotation_and_translation = remove_rotation_and_translation
328 if method is None:
329 warnings.warn(
330 "The default method has changed from 'aseneb' to "
331 "'improvedtangent'. The 'aseneb' method is an unpublished, "
332 "custom implementation that is not recommended as it frequently "
333 "results in very poor bands. Please explicitly set "
334 "method='improvedtangent' to silence this warning, or set "
335 "method='aseneb' if you strictly require the old behavior "
336 "(results may vary). "
337 "See: https://gitlab.com/ase/ase/-/merge_requests/3952",
338 UserWarning,
339 )
340 method = 'improvedtangent'
342 if method in ['aseneb', 'eb', 'improvedtangent', 'spline', 'string']:
343 self.method = method
344 else:
345 raise NotImplementedError(method)
347 if precon is not None and method not in ['spline', 'string']:
348 raise NotImplementedError(f'no precon implemented: {method}')
349 self.precon = precon
351 self.neb_method = get_neb_method(self, method)
352 if isinstance(k, (float, int)):
353 k = [k] * (self.nimages - 1)
354 self.k = list(k)
356 if world is None:
357 world = ase.parallel.world
358 self.world = world
360 if parallel:
361 if self.allow_shared_calculator:
362 raise RuntimeError(
363 "Cannot use shared calculators in parallel in NEB.")
364 self.real_forces = None # ndarray of shape (nimages, natom, 3)
365 self.energies = None # ndarray of shape (nimages,)
366 self.residuals = None # ndarray of shape (nimages,)
368 def __ase_optimizable__(self):
369 return NEBOptimizable(self)
371 @property
372 def natoms(self):
373 return len(self.images[0])
375 @property
376 def nimages(self):
377 return len(self.images)
379 @staticmethod
380 def freeze_results_on_image(atoms: ase.Atoms,
381 **results_to_include):
382 atoms.calc = SinglePointCalculator(atoms=atoms, **results_to_include)
384 def interpolate(self, method='linear', mic=False, apply_constraint=None):
385 """Interpolate the positions of the interior images between the
386 initial state (image 0) and final state (image -1).
388 method: str
389 Method by which to interpolate: 'linear' or 'idpp'.
390 linear provides a standard straight-line interpolation, while
391 idpp uses an image-dependent pair potential.
392 mic: bool
393 Use the minimum-image convention when interpolating.
394 apply_constraint: bool
395 Controls if the constraints attached to the images
396 are ignored or applied when setting the interpolated positions.
397 Default value is None, in this case the resulting constrained
398 positions (apply_constraint=True) are compared with unconstrained
399 positions (apply_constraint=False),
400 if the positions are not the same
401 the user is required to specify the desired behaviour
402 by setting up apply_constraint keyword argument to False or True.
403 """
404 if self.remove_rotation_and_translation:
405 minimize_rotation_and_translation(self.images[0], self.images[-1])
407 interpolate(self.images, mic, apply_constraint=apply_constraint)
409 if method == 'idpp':
410 idpp_interpolate(images=self, traj=None, log=None, mic=mic)
412 @deprecated("Please use NEB's interpolate(method='idpp') method or "
413 "directly call the idpp_interpolate function from ase.mep")
414 def idpp_interpolate(self, traj='idpp.traj', log='idpp.log', fmax=0.1,
415 optimizer=MDMin, mic=False, steps=100):
416 """
417 .. deprecated:: 3.23.0
418 Please use :class:`~ase.mep.NEB`'s ``interpolate(method='idpp')``
419 method
420 """
421 idpp_interpolate(self, traj=traj, log=log, fmax=fmax,
422 optimizer=optimizer, mic=mic, steps=steps)
424 def get_positions(self):
425 positions = np.empty(((self.nimages - 2) * self.natoms, 3))
426 n1 = 0
427 for image in self.images[1:-1]:
428 n2 = n1 + self.natoms
429 positions[n1:n2] = image.get_positions()
430 n1 = n2
431 return positions
433 def set_positions(self, positions, adjust_positions=True):
434 if adjust_positions:
435 # optional reparameterisation step: some NEB methods need to adjust
436 # positions e.g. string method does this to equispace the images)
437 positions = self.neb_method.adjust_positions(positions)
438 n1 = 0
439 for image in self.images[1:-1]:
440 n2 = n1 + self.natoms
441 image.set_positions(positions[n1:n2])
442 n1 = n2
444 def get_forces(self):
445 """Evaluate and return the forces."""
446 images = self.images
448 if not self.allow_shared_calculator:
449 calculators = [image.calc for image in images
450 if image.calc is not None]
451 if len(set(calculators)) != len(calculators):
452 msg = ('One or more NEB images share the same calculator. '
453 'Each image must have its own calculator. '
454 'You may wish to use the ase.mep.SingleCalculatorNEB '
455 'class instead, although using separate calculators '
456 'is recommended.')
457 raise ValueError(msg)
459 forces = np.empty(((self.nimages - 2), self.natoms, 3))
460 energies = np.empty(self.nimages)
461 real_forces = np.zeros((self.nimages, self.natoms, 3))
463 if self.remove_rotation_and_translation:
464 for i in range(1, self.nimages):
465 minimize_rotation_and_translation(images[i - 1], images[i])
467 if self.method != 'aseneb':
468 energies[0] = images[0].get_potential_energy()
469 energies[-1] = images[-1].get_potential_energy()
471 if not self.parallel:
472 # Do all images - one at a time:
473 for i in range(1, self.nimages - 1):
474 forces[i - 1] = images[i].get_forces()
475 energies[i] = images[i].get_potential_energy()
476 real_forces[i] = images[i].get_forces(apply_constraint=False)
478 elif self.world.size == 1:
479 def run(image, energies, forces, real_forces):
480 forces[:] = image.get_forces()
481 energies[:] = image.get_potential_energy()
482 real_forces[:] = image.get_forces(apply_constraint=False)
484 threads = [threading.Thread(target=run,
485 args=(images[i],
486 energies[i:i + 1],
487 forces[i - 1:i],
488 real_forces[i:i + 1]))
489 for i in range(1, self.nimages - 1)]
490 for thread in threads:
491 thread.start()
492 for thread in threads:
493 thread.join()
494 else:
495 # Parallelize over images:
496 i = self.world.rank * (self.nimages - 2) // self.world.size + 1
497 try:
498 forces[i - 1] = images[i].get_forces()
499 energies[i] = images[i].get_potential_energy()
500 real_forces[i] = images[i].get_forces(apply_constraint=False)
501 except Exception:
502 # Make sure other images also fail:
503 error = self.world.sum(1.0)
504 raise
505 else:
506 error = self.world.sum(0.0)
507 if error:
508 raise RuntimeError('Parallel NEB failed!')
510 for i in range(1, self.nimages - 1):
511 root = (i - 1) * self.world.size // (self.nimages - 2)
512 self.world.broadcast(energies[i:i + 1], root)
513 self.world.broadcast(forces[i - 1], root)
514 self.world.broadcast(real_forces[i], root)
516 # if this is the first force call, we need to build the preconditioners
517 if self.precon is None or isinstance(self.precon, (str, Precon, list)):
518 self.precon = PreconImages(self.precon, images)
520 # apply preconditioners to transform forces
521 # for the default IdentityPrecon this does not change their values
522 precon_forces = self.precon.apply(forces, index=slice(1, -1))
524 # Save for later use in iterimages:
525 self.energies = energies
526 self.real_forces = real_forces
528 state = NEBState(self, images, energies)
530 # Can we get rid of self.energies, self.imax, self.emax etc.?
531 self.imax = state.imax
532 self.emax = state.emax
534 spring1 = state.spring(0)
536 self.residuals = []
537 for i in range(1, self.nimages - 1):
538 spring2 = state.spring(i)
539 tangent = self.neb_method.get_tangent(state, spring1, spring2, i)
541 # Get overlap between full PES-derived force and tangent
542 tangential_force = np.vdot(forces[i - 1], tangent)
544 # from now on we use the preconditioned forces (equal for precon=ID)
545 imgforce = precon_forces[i - 1]
547 if i == self.imax and self.climb:
548 """The climbing image, imax, is not affected by the spring
549 forces. This image feels the full PES-derived force,
550 but the tangential component is inverted:
551 see Eq. 5 in paper II."""
552 if self.method == 'aseneb':
553 tangent_mag = np.vdot(tangent, tangent) # For normalizing
554 imgforce -= 2 * tangential_force / tangent_mag * tangent
555 else:
556 imgforce -= 2 * tangential_force * tangent
557 else:
558 self.neb_method.add_image_force(state, tangential_force,
559 tangent, imgforce, spring1,
560 spring2, i)
561 # compute the residual - with ID precon, this is just max force
562 residual = self.precon.get_residual(i, imgforce)
563 self.residuals.append(residual)
565 spring1 = spring2
567 return precon_forces.reshape((-1, 3))
569 def get_residual(self):
570 """Return residual force along the band.
572 Typically this the maximum force component on any image. For
573 non-trivial preconditioners, the appropriate preconditioned norm
574 is used to compute the residual.
575 """
576 if self.residuals is None:
577 raise RuntimeError("get_residual() called before get_forces()")
578 return np.max(self.residuals)
580 def get_potential_energy(self):
581 """Return the maximum potential energy along the band."""
582 return self.emax
584 def set_calculators(self, calculators):
585 """Set new calculators to the images.
587 Parameters
588 ----------
589 calculators : Calculator / list(Calculator)
590 calculator(s) to attach to images
591 - single calculator, only if allow_shared_calculator=True
592 list of calculators if length:
593 - length nimages, set to all images
594 - length nimages-2, set to non-end images only
595 """
597 if not isinstance(calculators, list):
598 if self.allow_shared_calculator:
599 calculators = [calculators] * self.nimages
600 else:
601 raise RuntimeError("Cannot set shared calculator to NEB "
602 "with allow_shared_calculator=False")
604 n = len(calculators)
605 if n == self.nimages:
606 for i in range(self.nimages):
607 self.images[i].calc = calculators[i]
608 elif n == self.nimages - 2:
609 for i in range(1, self.nimages - 1):
610 self.images[i].calc = calculators[i - 1]
611 else:
612 raise RuntimeError(
613 'len(calculators)=%d does not fit to len(images)=%d'
614 % (n, self.nimages))
616 def __len__(self):
617 # Corresponds to number of optimizable degrees of freedom, i.e.
618 # virtual atom count for the optimization algorithm.
619 return (self.nimages - 2) * self.natoms
621 def iterimages(self):
622 # Allows trajectory to convert NEB into several images
623 for i, atoms in enumerate(self.images):
624 if i == 0 or i == self.nimages - 1:
625 yield atoms
626 else:
627 atoms = atoms.copy()
628 self.freeze_results_on_image(
629 atoms, energy=self.energies[i],
630 forces=self.real_forces[i])
632 yield atoms
634 def spline_fit(self, positions=None, norm='precon'):
635 """
636 Fit a cubic spline to this NEB
638 Args:
639 norm (str, optional): Norm to use: 'precon' (default) or 'euclidean'
641 Returns
642 -------
643 fit: ase.precon.precon.SplineFit instance
644 """
645 if norm == 'precon':
646 if self.precon is None or isinstance(self.precon, str):
647 self.precon = PreconImages(self.precon, self.images)
648 precon = self.precon
649 # if this is the first call, we need to build the preconditioners
650 elif norm == 'euclidean':
651 precon = PreconImages('ID', self.images)
652 else:
653 raise ValueError(f'unsupported norm {norm}')
654 return precon.spline_fit(positions)
656 def integrate_forces(self, spline_points=1000, bc_type='not-a-knot'):
657 """Use spline fit to integrate forces along MEP to approximate
658 energy differences using the virtual work approach.
660 Args:
661 spline_points (int, optional): Number of points. Defaults to 1000.
662 bc_type (str, optional): Boundary conditions, default 'not-a-knot'.
664 Returns
665 -------
666 s: reaction coordinate in range [0, 1], with `spline_points` entries
667 E: result of integrating forces, on the same grid as `s`.
668 F: projected forces along MEP
669 """
670 # note we use standard Euclidean rather than preconditioned norm
671 # to compute the virtual work
672 fit = self.spline_fit(norm='euclidean')
673 forces = np.array([image.get_forces().reshape(-1)
674 for image in self.images])
675 f = CubicSpline(fit.s, forces, bc_type=bc_type)
677 s = np.linspace(0.0, 1.0, spline_points, endpoint=True)
678 dE = f(s) * fit.dx_ds(s)
679 F = dE.sum(axis=1)
680 E = -cumulative_trapezoid(F, s, initial=0.0)
681 return s, E, F
684class DyNEB(BaseNEB):
685 def __init__(
686 self,
687 images,
688 k=0.1,
689 fmax=0.05,
690 climb=False,
691 parallel=False,
692 remove_rotation_and_translation=False,
693 world=None,
694 dynamic_relaxation=True,
695 scale_fmax=0.,
696 method=None,
697 allow_shared_calculator=False,
698 precon=None,
699 ):
700 """
701 Subclass of NEB that allows for scaled and dynamic optimizations of
702 images. This method, which only works in series, does not perform
703 force calls on images that are below the convergence criterion.
704 The convergence criteria can be scaled with a displacement metric
705 to focus the optimization on the saddle point region.
707 'Scaled and Dynamic Optimizations of Nudged Elastic Bands',
708 P. Lindgren, G. Kastlunger and A. A. Peterson,
709 J. Chem. Theory Comput. 15, 11, 5787-5793 (2019).
711 dynamic_relaxation: bool
712 True skips images with forces below the convergence criterion.
713 This is updated after each force call; if a previously converged
714 image goes out of tolerance (due to spring adjustments between
715 the image and its neighbors), it will be optimized again.
716 False reverts to the default NEB implementation.
718 fmax: float
719 Must be identical to the fmax of the optimizer.
721 scale_fmax: float
722 Scale convergence criteria along band based on the distance between
723 an image and the image with the highest potential energy. This
724 keyword determines how rapidly the convergence criteria are scaled.
725 """
726 super().__init__(
727 images, k=k, climb=climb, parallel=parallel,
728 remove_rotation_and_translation=remove_rotation_and_translation,
729 world=world, method=method,
730 allow_shared_calculator=allow_shared_calculator, precon=precon)
731 self.fmax = fmax
732 self.dynamic_relaxation = dynamic_relaxation
733 self.scale_fmax = scale_fmax
735 if not self.dynamic_relaxation and self.scale_fmax:
736 msg = ('Scaled convergence criteria only implemented in series '
737 'with dynamic relaxation.')
738 raise ValueError(msg)
740 def set_positions(self, positions):
741 if not self.dynamic_relaxation:
742 return super().set_positions(positions)
744 n1 = 0
745 for i, image in enumerate(self.images[1:-1]):
746 if self.parallel:
747 msg = ('Dynamic relaxation does not work efficiently '
748 'when parallelizing over images. Try AutoNEB '
749 'routine for freezing images in parallel.')
750 raise ValueError(msg)
751 else:
752 forces_dyn = self._fmax_all(self.images)
753 if forces_dyn[i] < self.fmax:
754 n1 += self.natoms
755 else:
756 n2 = n1 + self.natoms
757 image.set_positions(positions[n1:n2])
758 n1 = n2
760 def _fmax_all(self, images):
761 """Store maximum force acting on each image in list. This is used in
762 the dynamic optimization routine in the set_positions() function."""
763 n = self.natoms
764 forces = self.get_forces()
765 fmax_images = [
766 np.sqrt((forces[n * i:n + n * i] ** 2).sum(axis=1)).max()
767 for i in range(self.nimages - 2)]
768 return fmax_images
770 def get_forces(self):
771 forces = super().get_forces()
772 if not self.dynamic_relaxation:
773 return forces
775 """Get NEB forces and scale the convergence criteria to focus
776 optimization on saddle point region. The keyword scale_fmax
777 determines the rate of convergence scaling."""
778 n = self.natoms
779 for i in range(self.nimages - 2):
780 n1 = n * i
781 n2 = n1 + n
782 force = np.sqrt((forces[n1:n2] ** 2.).sum(axis=1)).max()
783 n_imax = (self.imax - 1) * n # Image with highest energy.
785 positions = self.get_positions()
786 pos_imax = positions[n_imax:n_imax + n]
788 """Scale convergence criteria based on distance between an
789 image and the image with the highest potential energy."""
790 rel_pos = np.sqrt(((positions[n1:n2] - pos_imax) ** 2).sum())
791 if force < self.fmax * (1 + rel_pos * self.scale_fmax):
792 if i == self.imax - 1:
793 # Keep forces at saddle point for the log file.
794 pass
795 else:
796 # Set forces to zero before they are sent to optimizer.
797 forces[n1:n2, :] = 0
798 return forces
801def _check_deprecation(keyword, kwargs):
802 if keyword in kwargs:
803 warnings.warn(f'Keyword {keyword} of NEB is deprecated. '
804 'Please use the DyNEB class instead for dynamic '
805 'relaxation', FutureWarning)
808class NEB(DyNEB):
809 def __init__(self, images, k=0.1, climb=False, parallel=False,
810 remove_rotation_and_translation=False, world=None,
811 method=None, allow_shared_calculator=False,
812 precon=None, **kwargs):
813 """Nudged elastic band.
815 Paper I:
817 G. Henkelman and H. Jonsson, Chem. Phys, 113, 9978 (2000).
818 :doi:`10.1063/1.1323224`
820 Paper II:
822 G. Henkelman, B. P. Uberuaga, and H. Jonsson, Chem. Phys,
823 113, 9901 (2000).
824 :doi:`10.1063/1.1329672`
826 Paper III:
828 E. L. Kolsbjerg, M. N. Groves, and B. Hammer, J. Chem. Phys,
829 145, 094107 (2016)
830 :doi:`10.1063/1.4961868`
832 Paper IV:
834 S. Makri, C. Ortner and J. R. Kermode, J. Chem. Phys.
835 150, 094109 (2019)
836 https://dx.doi.org/10.1063/1.5064465
838 images: list of Atoms objects
839 Images defining path from initial to final state.
840 k: float or list of floats
841 Spring constant(s) in eV/Ang. One number or one for each spring.
842 climb: bool
843 Use a climbing image (default is no climbing image).
844 parallel: bool
845 Distribute images over processors.
846 remove_rotation_and_translation: bool
847 TRUE actives NEB-TR for removing translation and
848 rotation during NEB. By default applied non-periodic
849 systems
850 method: string of method
851 Choice betweeen five methods:
853 * aseneb: legacy ase NEB implementation
854 * improvedtangent: Paper I NEB implementation (default)
855 * eb: Paper III full spring force implementation
856 * spline: Paper IV spline interpolation (supports precon)
857 * string: Paper IV string method (supports precon)
859 Defaults to 'improvedtangent' (with a warning if not specified).
861 .. versionchanged:: 3.27.1
862 The default changes from ``aseneb`` to ``improvedtangent``.
864 allow_shared_calculator: bool
865 Allow images to share the same calculator between them.
866 Incompatible with parallelisation over images.
867 precon: string, :class:`ase.optimize.precon.Precon` instance or list of
868 instances. If present, enable preconditioing as in Paper IV. This is
869 possible using the 'spline' or 'string' methods only.
870 Default is no preconditioning (precon=None), which is converted to
871 a list of :class:`ase.precon.precon.IdentityPrecon` instances.
872 """
873 for keyword in 'dynamic_relaxation', 'fmax', 'scale_fmax':
874 _check_deprecation(keyword, kwargs)
875 defaults = dict(dynamic_relaxation=False,
876 fmax=0.05,
877 scale_fmax=0.0)
878 defaults.update(kwargs)
879 # Only reason for separating BaseNEB/NEB is that we are
880 # deprecating dynamic_relaxation.
881 #
882 # We can turn BaseNEB into NEB once we get rid of the
883 # deprecated variables.
884 #
885 # Then we can also move DyNEB into ase.dyneb without cyclic imports.
886 # We can do that in ase-3.22 or 3.23.
887 super().__init__(
888 images, k=k, climb=climb, parallel=parallel,
889 remove_rotation_and_translation=remove_rotation_and_translation,
890 world=world, method=method,
891 allow_shared_calculator=allow_shared_calculator,
892 precon=precon,
893 **defaults)
896class NEBOptimizer(Optimizer):
897 """
898 This optimizer applies an adaptive ODE solver to a NEB
900 Details of the adaptive ODE solver are described in paper IV
901 """
903 def __init__(self,
904 neb,
905 restart=None, logfile='-', trajectory=None,
906 master=None,
907 append_trajectory=False,
908 method='ODE',
909 alpha=0.01,
910 verbose=0,
911 rtol=0.1,
912 C1=1e-2,
913 C2=2.0):
915 super().__init__(atoms=neb, restart=restart,
916 logfile=logfile, trajectory=trajectory,
917 master=master,
918 append_trajectory=append_trajectory)
919 self.neb = neb
921 method = method.lower()
922 methods = ['ode', 'static', 'krylov']
923 if method not in methods:
924 raise ValueError(f'method must be one of {methods}')
925 self.method = method
927 self.alpha = alpha
928 self.verbose = verbose
929 self.rtol = rtol
930 self.C1 = C1
931 self.C2 = C2
933 def force_function(self, X):
934 positions = X.reshape((self.neb.nimages - 2) *
935 self.neb.natoms, 3)
936 self.neb.set_positions(positions)
937 forces = self.neb.get_forces().reshape(-1)
938 return forces
940 def get_residual(self, F=None, X=None):
941 return self.neb.get_residual()
943 def log(self):
944 fmax = self.get_residual()
945 T = time.localtime()
946 if self.logfile is not None:
947 name = f'{self.__class__.__name__}[{self.method}]'
948 if self.nsteps == 0:
949 args = (" " * len(name), "Step", "Time", "fmax")
950 msg = "%s %4s %8s %12s\n" % args
951 self.logfile.write(msg)
953 args = (name, self.nsteps, T[3], T[4], T[5], fmax)
954 msg = "%s: %3d %02d:%02d:%02d %12.4f\n" % args
955 self.logfile.write(msg)
957 def callback(self, X, F=None):
958 self.log()
959 self.call_observers()
960 self.nsteps += 1
962 def run_ode(self, fmax):
963 try:
964 ode12r(self.force_function,
965 self.neb.get_positions().reshape(-1),
966 fmax=fmax,
967 rtol=self.rtol,
968 C1=self.C1,
969 C2=self.C2,
970 steps=self.max_steps,
971 verbose=self.verbose,
972 callback=self.callback,
973 residual=self.get_residual)
974 return True
975 except OptimizerConvergenceError:
976 return False
978 def run_static(self, fmax):
979 X = self.neb.get_positions().reshape(-1)
980 for _ in range(self.max_steps):
981 F = self.force_function(X)
982 if self.neb.get_residual() <= fmax:
983 return True
984 X += self.alpha * F
985 self.callback(X)
986 return False
988 def run(self, fmax=0.05, steps=DEFAULT_MAX_STEPS, method=None):
989 """
990 Optimize images to obtain the minimum energy path
992 Parameters
993 ----------
994 fmax - desired force tolerance
995 steps - maximum number of steps
996 """
997 self.max_steps = steps
998 if method is None:
999 method = self.method
1000 if method == 'ode':
1001 return self.run_ode(fmax)
1002 elif method == 'static':
1003 return self.run_static(fmax)
1004 else:
1005 raise ValueError(f'unknown method: {self.method}')
1008class IDPP(Calculator):
1009 """Image dependent pair potential.
1011 See:
1012 Improved initial guess for minimum energy path calculations.
1013 Søren Smidstrup, Andreas Pedersen, Kurt Stokbro and Hannes Jónsson
1014 Chem. Phys. 140, 214106 (2014)
1015 """
1017 implemented_properties = ['energy', 'forces']
1019 def __init__(self, target, mic):
1020 Calculator.__init__(self)
1021 self.target = target
1022 self.mic = mic
1024 def calculate(self, atoms, properties, system_changes):
1025 Calculator.calculate(self, atoms, properties, system_changes)
1027 P = atoms.get_positions()
1028 d = []
1029 D = []
1030 for p in P:
1031 Di = P - p
1032 if self.mic:
1033 Di, di = find_mic(Di, atoms.get_cell(), atoms.get_pbc())
1034 else:
1035 di = np.sqrt((Di ** 2).sum(1))
1036 d.append(di)
1037 D.append(Di)
1038 d = np.array(d)
1039 D = np.array(D)
1041 dd = d - self.target
1042 d.ravel()[::len(d) + 1] = 1 # avoid dividing by zero
1043 d4 = d ** 4
1044 e = 0.5 * (dd ** 2 / d4).sum()
1045 f = -2 * ((dd * (1 - 2 * dd / d) / d ** 5)[..., np.newaxis] * D).sum(
1046 0)
1047 self.results = {'energy': e, 'forces': f}
1050@deprecated("SingleCalculatorNEB is deprecated. "
1051 "Please use NEB(allow_shared_calculator=True) instead.")
1052class SingleCalculatorNEB(NEB):
1053 """
1054 .. deprecated:: 3.23.0
1055 Please use ``NEB(allow_shared_calculator=True)`` instead
1056 """
1058 def __init__(self, images, *args, **kwargs):
1059 kwargs["allow_shared_calculator"] = True
1060 super().__init__(images, *args, **kwargs)
1063def interpolate(images, mic=False, interpolate_cell=False,
1064 use_scaled_coord=False, apply_constraint=None):
1065 """Given a list of images, linearly interpolate the positions of the
1066 interior images.
1068 mic: bool
1069 Map movement into the unit cell by using the minimum image convention.
1070 interpolate_cell: bool
1071 Interpolate the three cell vectors linearly just like the atomic
1072 positions. Not implemented for NEB calculations!
1073 use_scaled_coord: bool
1074 Use scaled/internal/fractional coordinates instead of real ones for the
1075 interpolation. Not implemented for NEB calculations!
1076 apply_constraint: bool
1077 Controls if the constraints attached to the images
1078 are ignored or applied when setting the interpolated positions.
1079 Default value is None, in this case the resulting constrained positions
1080 (apply_constraint=True) are compared with unconstrained positions
1081 (apply_constraint=False), if the positions are not the same
1082 the user is required to specify the desired behaviour
1083 by setting up apply_constraint keyword argument to False or True.
1084 """
1085 if use_scaled_coord:
1086 pos1 = images[0].get_scaled_positions(wrap=mic)
1087 pos2 = images[-1].get_scaled_positions(wrap=mic)
1088 else:
1089 pos1 = images[0].get_positions()
1090 pos2 = images[-1].get_positions()
1091 d = pos2 - pos1
1092 if not use_scaled_coord and mic:
1093 d = find_mic(d, images[0].get_cell(), images[0].pbc)[0]
1094 d /= (len(images) - 1.0)
1095 if interpolate_cell:
1096 cell1 = images[0].get_cell()
1097 cell2 = images[-1].get_cell()
1098 cell_diff = cell2 - cell1
1099 cell_diff /= (len(images) - 1.0)
1100 for i in range(1, len(images) - 1):
1101 # first the new cell, otherwise scaled positions are wrong
1102 if interpolate_cell:
1103 images[i].set_cell(cell1 + i * cell_diff)
1104 new_pos = pos1 + i * d
1105 if use_scaled_coord:
1106 images[i].set_scaled_positions(new_pos)
1107 else:
1108 if apply_constraint is None:
1109 unconstrained_image = images[i].copy()
1110 unconstrained_image.set_positions(new_pos,
1111 apply_constraint=False)
1112 images[i].set_positions(new_pos, apply_constraint=True)
1113 if not np.allclose(unconstrained_image.positions,
1114 images[i].positions):
1115 raise RuntimeError(f"Constraints in image {i}\n"
1116 "affect the interpolation results.\n"
1117 "Please specify if you want to\n"
1118 "apply or ignore the constraints\n"
1119 "during the interpolation\n"
1120 "with the apply_constraint argument.")
1121 else:
1122 images[i].set_positions(new_pos,
1123 apply_constraint=apply_constraint)
1126def idpp_interpolate(images, traj='idpp.traj', log='idpp.log', fmax=0.1,
1127 optimizer=MDMin, mic=False, steps=100):
1128 """Interpolate using the IDPP method. 'images' can either be a plain
1129 list of images or an NEB object (containing a list of images)."""
1130 if hasattr(images, 'interpolate'):
1131 neb = images
1132 else:
1133 neb = NEB(images)
1135 d1 = neb.images[0].get_all_distances(mic=mic)
1136 d2 = neb.images[-1].get_all_distances(mic=mic)
1137 d = (d2 - d1) / (neb.nimages - 1)
1138 real_calcs = []
1139 for i, image in enumerate(neb.images):
1140 real_calcs.append(image.calc)
1141 image.calc = IDPP(d1 + i * d, mic=mic)
1143 with optimizer(neb, trajectory=traj, logfile=log) as opt:
1144 opt.run(fmax=fmax, steps=steps)
1146 for image, calc in zip(neb.images, real_calcs):
1147 image.calc = calc
1150class NEBTools:
1151 """Class to make many of the common tools for NEB analysis available to
1152 the user. Useful for scripting the output of many jobs. Initialize with
1153 list of images which make up one or more band of the NEB relaxation."""
1155 def __init__(self, images):
1156 self.images = images
1158 @deprecated('NEBTools.get_fit() is deprecated. '
1159 'Please use ase.utils.forcecurve.fit_images(images).')
1160 def get_fit(self):
1161 """
1162 .. deprecated:: 3.23.0
1163 Please use ``ase.utils.forcecurve.fit_images(images)``
1164 """
1165 return fit_images(self.images)
1167 def get_barrier(self, fit=True, raw=False):
1168 """Returns the barrier estimate from the NEB, along with the
1169 Delta E of the elementary reaction. If fit=True, the barrier is
1170 estimated based on the interpolated fit to the images; if
1171 fit=False, the barrier is taken as the maximum-energy image
1172 without interpolation. Set raw=True to get the raw energy of the
1173 transition state instead of the forward barrier."""
1174 forcefit = fit_images(self.images)
1175 energies = forcefit.energies
1176 fit_energies = forcefit.fit_energies
1177 dE = energies[-1] - energies[0]
1178 if fit:
1179 barrier = max(fit_energies)
1180 else:
1181 barrier = max(energies)
1182 if raw:
1183 barrier += self.images[0].get_potential_energy()
1184 return barrier, dE
1186 def get_fmax(self, **kwargs):
1187 """Returns fmax, as used by optimizers with NEB."""
1188 neb = NEB(self.images, **kwargs)
1189 forces = neb.get_forces()
1190 return np.sqrt((forces ** 2).sum(axis=1).max())
1192 def plot_band(self, ax=None):
1193 """Plots the NEB band on matplotlib axes object 'ax'. If ax=None
1194 returns a new figure object."""
1195 forcefit = fit_images(self.images)
1196 ax = forcefit.plot(ax=ax)
1197 return ax.figure
1199 def plot_bands(self, constant_x=False, constant_y=False,
1200 nimages=None, label='nebplots'):
1201 """Given a trajectory containing many steps of a NEB, makes
1202 plots of each band in the series in a single PDF.
1204 constant_x: bool
1205 Use the same x limits on all plots.
1206 constant_y: bool
1207 Use the same y limits on all plots.
1208 nimages: int
1209 Number of images per band. Guessed if not supplied.
1210 label: str
1211 Name for the output file. .pdf will be appended.
1212 """
1213 from matplotlib import pyplot
1214 from matplotlib.backends.backend_pdf import PdfPages
1215 if nimages is None:
1216 nimages = self._guess_nimages()
1217 nebsteps = len(self.images) // nimages
1218 if constant_x or constant_y:
1219 sys.stdout.write('Scaling axes.\n')
1220 sys.stdout.flush()
1221 # Plot all to one plot, then pull its x and y range.
1222 fig, ax = pyplot.subplots()
1223 for index in range(nebsteps):
1224 images = self.images[index * nimages:(index + 1) * nimages]
1225 NEBTools(images).plot_band(ax=ax)
1226 xlim = ax.get_xlim()
1227 ylim = ax.get_ylim()
1228 pyplot.close(fig) # Reference counting "bug" in pyplot.
1229 with PdfPages(label + '.pdf') as pdf:
1230 for index in range(nebsteps):
1231 sys.stdout.write('\rProcessing band {:10d} / {:10d}'
1232 .format(index, nebsteps))
1233 sys.stdout.flush()
1234 fig, ax = pyplot.subplots()
1235 images = self.images[index * nimages:(index + 1) * nimages]
1236 NEBTools(images).plot_band(ax=ax)
1237 if constant_x:
1238 ax.set_xlim(xlim)
1239 if constant_y:
1240 ax.set_ylim(ylim)
1241 pdf.savefig(fig)
1242 pyplot.close(fig) # Reference counting "bug" in pyplot.
1243 sys.stdout.write('\n')
1245 def _guess_nimages(self):
1246 """Attempts to guess the number of images per band from
1247 a trajectory, based solely on the repetition of the
1248 potential energy of images. This should also work for symmetric
1249 cases."""
1250 e_first = self.images[0].get_potential_energy()
1251 nimages = None
1252 for index, image in enumerate(self.images[1:], start=1):
1253 e = image.get_potential_energy()
1254 if e == e_first:
1255 # Need to check for symmetric case when e_first = e_last.
1256 try:
1257 e_next = self.images[index + 1].get_potential_energy()
1258 except IndexError:
1259 pass
1260 else:
1261 if e_next == e_first:
1262 nimages = index + 1 # Symmetric
1263 break
1264 nimages = index # Normal
1265 break
1266 if nimages is None:
1267 sys.stdout.write('Appears to be only one band in the images.\n')
1268 return len(self.images)
1269 # Sanity check that the energies of the last images line up too.
1270 e_last = self.images[nimages - 1].get_potential_energy()
1271 e_nextlast = self.images[2 * nimages - 1].get_potential_energy()
1272 if e_last != e_nextlast:
1273 raise RuntimeError('Could not guess number of images per band.')
1274 sys.stdout.write('Number of images per band guessed to be {:d}.\n'
1275 .format(nimages))
1276 return nimages
1279class NEBtools(NEBTools):
1280 @deprecated('NEBtools has been renamed; please use NEBTools.')
1281 def __init__(self, images):
1282 """
1283 .. deprecated:: 3.23.0
1284 Please use :class:`~ase.mep.NEBTools`.
1285 """
1286 NEBTools.__init__(self, images)
1289@deprecated('Please use NEBTools.plot_band_from_fit.')
1290def plot_band_from_fit(s, E, Sfit, Efit, lines, ax=None):
1291 """
1292 .. deprecated:: 3.23.0
1293 Please use :meth:`NEBTools.plot_band_from_fit`.
1294 """
1295 NEBTools.plot_band_from_fit(s, E, Sfit, Efit, lines, ax=None)