Coverage for /builds/kinetik161/ase/ase/mep/neb.py: 83.44%
640 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +0000
1import sys
2import threading
3import time
4import warnings
5from abc import ABC, abstractmethod
7import numpy as np
8from scipy.integrate import cumtrapz
9from scipy.interpolate import CubicSpline
11import ase.parallel
12from ase.build import minimize_rotation_and_translation
13from ase.calculators.calculator import Calculator
14from ase.calculators.singlepoint import SinglePointCalculator
15from ase.geometry import find_mic
16from ase.optimize import MDMin
17from ase.optimize.ode import ode12r
18from ase.optimize.optimize import DEFAULT_MAX_STEPS, Optimizer
19from ase.optimize.precon import Precon, PreconImages
20from ase.optimize.sciopt import OptimizerConvergenceError
21from ase.utils import deprecated, lazyproperty
22from ase.utils.abc import Optimizable
23from ase.utils.forcecurve import fit_images
26class Spring:
27 def __init__(self, atoms1, atoms2, energy1, energy2, k):
28 self.atoms1 = atoms1
29 self.atoms2 = atoms2
30 self.energy1 = energy1
31 self.energy2 = energy2
32 self.k = k
34 def _find_mic(self):
35 pos1 = self.atoms1.get_positions()
36 pos2 = self.atoms2.get_positions()
37 # XXX If we want variable cells we will need to edit this.
38 mic, _ = find_mic(pos2 - pos1, self.atoms1.cell, self.atoms1.pbc)
39 return mic
41 @lazyproperty
42 def t(self):
43 return self._find_mic()
45 @lazyproperty
46 def nt(self):
47 return np.linalg.norm(self.t)
50class NEBState:
51 def __init__(self, neb, images, energies):
52 self.neb = neb
53 self.images = images
54 self.energies = energies
56 def spring(self, i):
57 return Spring(self.images[i], self.images[i + 1],
58 self.energies[i], self.energies[i + 1],
59 self.neb.k[i])
61 @lazyproperty
62 def imax(self):
63 return 1 + np.argsort(self.energies[1:-1])[-1]
65 @property
66 def emax(self):
67 return self.energies[self.imax]
69 @lazyproperty
70 def eqlength(self):
71 images = self.images
72 beeline = (images[self.neb.nimages - 1].get_positions() -
73 images[0].get_positions())
74 beelinelength = np.linalg.norm(beeline)
75 return beelinelength / (self.neb.nimages - 1)
77 @lazyproperty
78 def nimages(self):
79 return len(self.images)
81 @property
82 def precon(self):
83 return self.neb.precon
86class NEBMethod(ABC):
87 def __init__(self, neb):
88 self.neb = neb
90 @abstractmethod
91 def get_tangent(self, state, spring1, spring2, i):
92 ...
94 @abstractmethod
95 def add_image_force(self, state, tangential_force, tangent, imgforce,
96 spring1, spring2, i):
97 ...
99 def adjust_positions(self, positions):
100 return positions
103class ImprovedTangentMethod(NEBMethod):
104 """
105 Tangent estimates are improved according to Eqs. 8-11 in paper I.
106 Tangents are weighted at extrema to ensure smooth transitions between
107 the positive and negative tangents.
108 """
110 def get_tangent(self, state, spring1, spring2, i):
111 energies = state.energies
112 if energies[i + 1] > energies[i] > energies[i - 1]:
113 tangent = spring2.t.copy()
114 elif energies[i + 1] < energies[i] < energies[i - 1]:
115 tangent = spring1.t.copy()
116 else:
117 deltavmax = max(abs(energies[i + 1] - energies[i]),
118 abs(energies[i - 1] - energies[i]))
119 deltavmin = min(abs(energies[i + 1] - energies[i]),
120 abs(energies[i - 1] - energies[i]))
121 if energies[i + 1] > energies[i - 1]:
122 tangent = spring2.t * deltavmax + spring1.t * deltavmin
123 else:
124 tangent = spring2.t * deltavmin + spring1.t * deltavmax
125 # Normalize the tangent vector
126 tangent /= np.linalg.norm(tangent)
127 return tangent
129 def add_image_force(self, state, tangential_force, tangent, imgforce,
130 spring1, spring2, i):
131 imgforce -= tangential_force * tangent
132 # Improved parallel spring force (formula 12 of paper I)
133 imgforce += (spring2.nt * spring2.k - spring1.nt * spring1.k) * tangent
136class ASENEBMethod(NEBMethod):
137 """
138 Standard NEB implementation in ASE. The tangent of each image is
139 estimated from the spring closest to the saddle point in each
140 spring pair.
141 """
143 def get_tangent(self, state, spring1, spring2, i):
144 imax = self.neb.imax
145 if i < imax:
146 tangent = spring2.t
147 elif i > imax:
148 tangent = spring1.t
149 else:
150 tangent = spring1.t + spring2.t
151 return tangent
153 def add_image_force(self, state, tangential_force, tangent, imgforce,
154 spring1, spring2, i):
155 # Magnitude for normalizing. Ensure it is not 0
156 tangent_mag = np.vdot(tangent, tangent) or 1
157 factor = tangent / tangent_mag
158 imgforce -= tangential_force * factor
159 imgforce -= np.vdot(
160 spring1.t * spring1.k -
161 spring2.t * spring2.k, tangent) * factor
164class FullSpringMethod(NEBMethod):
165 """
166 Elastic band method. The full spring force is included.
167 """
169 def get_tangent(self, state, spring1, spring2, i):
170 # Tangents are bisections of spring-directions
171 # (formula C8 of paper III)
172 tangent = spring1.t / spring1.nt + spring2.t / spring2.nt
173 tangent /= np.linalg.norm(tangent)
174 return tangent
176 def add_image_force(self, state, tangential_force, tangent, imgforce,
177 spring1, spring2, i):
178 imgforce -= tangential_force * tangent
179 energies = state.energies
180 # Spring forces
181 # Eqs. C1, C5, C6 and C7 in paper III)
182 f1 = -(spring1.nt -
183 state.eqlength) * spring1.t / spring1.nt * spring1.k
184 f2 = (spring2.nt - state.eqlength) * spring2.t / spring2.nt * spring2.k
185 if self.neb.climb and abs(i - self.neb.imax) == 1:
186 deltavmax = max(abs(energies[i + 1] - energies[i]),
187 abs(energies[i - 1] - energies[i]))
188 deltavmin = min(abs(energies[i + 1] - energies[i]),
189 abs(energies[i - 1] - energies[i]))
190 imgforce += (f1 + f2) * deltavmin / deltavmax
191 else:
192 imgforce += f1 + f2
195class BaseSplineMethod(NEBMethod):
196 """
197 Base class for SplineNEB and String methods
199 Can optionally be preconditioned, as described in the following article:
201 S. Makri, C. Ortner and J. R. Kermode, J. Chem. Phys.
202 150, 094109 (2019)
203 https://dx.doi.org/10.1063/1.5064465
204 """
206 def __init__(self, neb):
207 NEBMethod.__init__(self, neb)
209 def get_tangent(self, state, spring1, spring2, i):
210 return state.precon.get_tangent(i)
212 def add_image_force(self, state, tangential_force, tangent, imgforce,
213 spring1, spring2, i):
214 # project out tangential component (Eqs 6 and 7 in Paper IV)
215 imgforce -= tangential_force * tangent
218class SplineMethod(BaseSplineMethod):
219 """
220 NEB using spline interpolation, plus optional preconditioning
221 """
223 def add_image_force(self, state, tangential_force, tangent, imgforce,
224 spring1, spring2, i):
225 super().add_image_force(state, tangential_force,
226 tangent, imgforce, spring1, spring2, i)
227 eta = state.precon.get_spring_force(i, spring1.k, spring2.k, tangent)
228 imgforce += eta
231class StringMethod(BaseSplineMethod):
232 """
233 String method using spline interpolation, plus optional preconditioning
234 """
236 def adjust_positions(self, positions):
237 # fit cubic spline to positions, reinterpolate to equispace images
238 # note this uses the preconditioned distance metric.
239 fit = self.neb.spline_fit(positions)
240 new_s = np.linspace(0.0, 1.0, self.neb.nimages)
241 new_positions = fit.x(new_s[1:-1]).reshape(-1, 3)
242 return new_positions
245def get_neb_method(neb, method):
246 if method == 'eb':
247 return FullSpringMethod(neb)
248 elif method == 'aseneb':
249 return ASENEBMethod(neb)
250 elif method == 'improvedtangent':
251 return ImprovedTangentMethod(neb)
252 elif method == 'spline':
253 return SplineMethod(neb)
254 elif method == 'string':
255 return StringMethod(neb)
256 else:
257 raise ValueError(f'Bad method: {method}')
260class NEBOptimizable(Optimizable):
261 def __init__(self, neb):
262 self.neb = neb
264 def get_forces(self):
265 return self.neb.get_forces()
267 def get_potential_energy(self):
268 return self.neb.get_potential_energy()
270 def is_neb(self):
271 return True
273 def get_positions(self):
274 return self.neb.get_positions()
276 def set_positions(self, positions):
277 self.neb.set_positions(positions)
279 def __len__(self):
280 return 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 = -cumtrapz(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 step 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 """
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 try:
1066 np.testing.assert_allclose(unconstrained_image.positions,
1067 images[i].positions)
1068 except AssertionError:
1069 raise RuntimeError(f"Constraint(s) in image number {i} \n"
1070 "affect the interpolation results.\n"
1071 "Please specify if you want to \n"
1072 "apply or ignore the constraints \n"
1073 "during the interpolation \n"
1074 "with apply_constraint argument.")
1075 else:
1076 images[i].set_positions(new_pos,
1077 apply_constraint=apply_constraint)
1080def idpp_interpolate(images, traj='idpp.traj', log='idpp.log', fmax=0.1,
1081 optimizer=MDMin, mic=False, steps=100):
1082 """Interpolate using the IDPP method. 'images' can either be a plain
1083 list of images or an NEB object (containing a list of images)."""
1084 if hasattr(images, 'interpolate'):
1085 neb = images
1086 else:
1087 neb = NEB(images)
1089 d1 = neb.images[0].get_all_distances(mic=mic)
1090 d2 = neb.images[-1].get_all_distances(mic=mic)
1091 d = (d2 - d1) / (neb.nimages - 1)
1092 real_calcs = []
1093 for i, image in enumerate(neb.images):
1094 real_calcs.append(image.calc)
1095 image.calc = IDPP(d1 + i * d, mic=mic)
1097 with optimizer(neb, trajectory=traj, logfile=log) as opt:
1098 opt.run(fmax=fmax, steps=steps)
1100 for image, calc in zip(neb.images, real_calcs):
1101 image.calc = calc
1104class NEBTools:
1105 """Class to make many of the common tools for NEB analysis available to
1106 the user. Useful for scripting the output of many jobs. Initialize with
1107 list of images which make up one or more band of the NEB relaxation."""
1109 def __init__(self, images):
1110 self.images = images
1112 @deprecated('NEBTools.get_fit() is deprecated. '
1113 'Please use ase.utils.forcecurve.fit_images(images).')
1114 def get_fit(self):
1115 """
1116 .. deprecated:: 3.23.0
1117 Please use ``ase.utils.forcecurve.fit_images(images)``
1118 """
1119 return fit_images(self.images)
1121 def get_barrier(self, fit=True, raw=False):
1122 """Returns the barrier estimate from the NEB, along with the
1123 Delta E of the elementary reaction. If fit=True, the barrier is
1124 estimated based on the interpolated fit to the images; if
1125 fit=False, the barrier is taken as the maximum-energy image
1126 without interpolation. Set raw=True to get the raw energy of the
1127 transition state instead of the forward barrier."""
1128 forcefit = fit_images(self.images)
1129 energies = forcefit.energies
1130 fit_energies = forcefit.fit_energies
1131 dE = energies[-1] - energies[0]
1132 if fit:
1133 barrier = max(fit_energies)
1134 else:
1135 barrier = max(energies)
1136 if raw:
1137 barrier += self.images[0].get_potential_energy()
1138 return barrier, dE
1140 def get_fmax(self, **kwargs):
1141 """Returns fmax, as used by optimizers with NEB."""
1142 neb = NEB(self.images, **kwargs)
1143 forces = neb.get_forces()
1144 return np.sqrt((forces ** 2).sum(axis=1).max())
1146 def plot_band(self, ax=None):
1147 """Plots the NEB band on matplotlib axes object 'ax'. If ax=None
1148 returns a new figure object."""
1149 forcefit = fit_images(self.images)
1150 ax = forcefit.plot(ax=ax)
1151 return ax.figure
1153 def plot_bands(self, constant_x=False, constant_y=False,
1154 nimages=None, label='nebplots'):
1155 """Given a trajectory containing many steps of a NEB, makes
1156 plots of each band in the series in a single PDF.
1158 constant_x: bool
1159 Use the same x limits on all plots.
1160 constant_y: bool
1161 Use the same y limits on all plots.
1162 nimages: int
1163 Number of images per band. Guessed if not supplied.
1164 label: str
1165 Name for the output file. .pdf will be appended.
1166 """
1167 from matplotlib import pyplot
1168 from matplotlib.backends.backend_pdf import PdfPages
1169 if nimages is None:
1170 nimages = self._guess_nimages()
1171 nebsteps = len(self.images) // nimages
1172 if constant_x or constant_y:
1173 sys.stdout.write('Scaling axes.\n')
1174 sys.stdout.flush()
1175 # Plot all to one plot, then pull its x and y range.
1176 fig, ax = pyplot.subplots()
1177 for index in range(nebsteps):
1178 images = self.images[index * nimages:(index + 1) * nimages]
1179 NEBTools(images).plot_band(ax=ax)
1180 xlim = ax.get_xlim()
1181 ylim = ax.get_ylim()
1182 pyplot.close(fig) # Reference counting "bug" in pyplot.
1183 with PdfPages(label + '.pdf') as pdf:
1184 for index in range(nebsteps):
1185 sys.stdout.write('\rProcessing band {:10d} / {:10d}'
1186 .format(index, nebsteps))
1187 sys.stdout.flush()
1188 fig, ax = pyplot.subplots()
1189 images = self.images[index * nimages:(index + 1) * nimages]
1190 NEBTools(images).plot_band(ax=ax)
1191 if constant_x:
1192 ax.set_xlim(xlim)
1193 if constant_y:
1194 ax.set_ylim(ylim)
1195 pdf.savefig(fig)
1196 pyplot.close(fig) # Reference counting "bug" in pyplot.
1197 sys.stdout.write('\n')
1199 def _guess_nimages(self):
1200 """Attempts to guess the number of images per band from
1201 a trajectory, based solely on the repetition of the
1202 potential energy of images. This should also work for symmetric
1203 cases."""
1204 e_first = self.images[0].get_potential_energy()
1205 nimages = None
1206 for index, image in enumerate(self.images[1:], start=1):
1207 e = image.get_potential_energy()
1208 if e == e_first:
1209 # Need to check for symmetric case when e_first = e_last.
1210 try:
1211 e_next = self.images[index + 1].get_potential_energy()
1212 except IndexError:
1213 pass
1214 else:
1215 if e_next == e_first:
1216 nimages = index + 1 # Symmetric
1217 break
1218 nimages = index # Normal
1219 break
1220 if nimages is None:
1221 sys.stdout.write('Appears to be only one band in the images.\n')
1222 return len(self.images)
1223 # Sanity check that the energies of the last images line up too.
1224 e_last = self.images[nimages - 1].get_potential_energy()
1225 e_nextlast = self.images[2 * nimages - 1].get_potential_energy()
1226 if not (e_last == e_nextlast):
1227 raise RuntimeError('Could not guess number of images per band.')
1228 sys.stdout.write('Number of images per band guessed to be {:d}.\n'
1229 .format(nimages))
1230 return nimages
1233class NEBtools(NEBTools):
1234 @deprecated('NEBtools has been renamed; please use NEBTools.')
1235 def __init__(self, images):
1236 """
1237 .. deprecated:: 3.23.0
1238 Please use :class:`~ase.mep.NEBTools`.
1239 """
1240 NEBTools.__init__(self, images)
1243@deprecated('Please use NEBTools.plot_band_from_fit.')
1244def plot_band_from_fit(s, E, Sfit, Efit, lines, ax=None):
1245 """
1246 .. deprecated:: 3.23.0
1247 Please use :meth:`NEBTools.plot_band_from_fit`.
1248 """
1249 NEBTools.plot_band_from_fit(s, E, Sfit, Efit, lines, ax=None)