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