Coverage for ase / mep / dimer.py: 71.79%
585 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-03-30 08:22 +0000
1# fmt: off
3"""Minimum mode follower for finding saddle points in an unbiased way.
5There is, currently, only one implemented method: The Dimer method.
6"""
8import sys
9import time
10import warnings
11from math import atan, cos, degrees, pi, sin, sqrt, tan
12from typing import Any
14import numpy as np
16from ase.calculators.singlepoint import SinglePointCalculator
17from ase.optimize.optimize import OptimizableAtoms, Optimizer
18from ase.parallel import world
19from ase.utils import IOContext
21# Handy vector methods
22norm = np.linalg.norm
25class DimerOptimizable(OptimizableAtoms):
26 def __init__(self, dimeratoms):
27 self.dimeratoms = dimeratoms
28 super().__init__(dimeratoms)
30 def converged(self, forces, fmax):
31 forces_converged = super().converged(forces, fmax)
32 return forces_converged and self.dimeratoms.get_curvature() < 0.0
35def normalize(vector):
36 """Create a unit vector along *vector*"""
37 return vector / norm(vector)
40def parallel_vector(vector, base):
41 """Extract the components of *vector* that are parallel to *base*"""
42 return np.vdot(vector, base) * base
45def perpendicular_vector(vector, base):
46 """Remove the components of *vector* that are parallel to *base*"""
47 return vector - parallel_vector(vector, base)
50def rotate_vectors(v1i, v2i, angle):
51 """Rotate vectors *v1i* and *v2i* by *angle*"""
52 cAng = cos(angle)
53 sAng = sin(angle)
54 v1o = v1i * cAng + v2i * sAng
55 v2o = v2i * cAng - v1i * sAng
56 # Ensure the length of the input and output vectors is equal
57 return normalize(v1o) * norm(v1i), normalize(v2o) * norm(v2i)
60class DimerEigenmodeSearch:
61 """An implementation of the Dimer's minimum eigenvalue mode search.
63 This class implements the rotational part of the dimer saddle point
64 searching method.
66 Parameters
67 ----------
69 atoms: MinModeAtoms object
70 MinModeAtoms is an extension to the Atoms object, which includes
71 information about the lowest eigenvalue mode.
72 control: DimerControl object
73 Contains the parameters necessary for the eigenmode search.
74 If no control object is supplied a default DimerControl
75 will be created and used.
76 basis: list of xyz-values
77 Eigenmode. Must be an ndarray of shape (n, 3).
78 It is possible to constrain the eigenmodes to be orthogonal
79 to this given eigenmode.
81 Notes
82 -----
84 The code is inspired, with permission, by code written by the Henkelman
85 group, which can be found at http://theory.cm.utexas.edu/vtsttools/code/
87 References
88 ----------
90 * Henkelman and Jonsson, JCP 111, 7010 (1999)
91 * Olsen, Kroes, Henkelman, Arnaldsson, and Jonsson, JCP 121,
92 9776 (2004).
93 * Heyden, Bell, and Keil, JCP 123, 224101 (2005).
94 * Kastner and Sherwood, JCP 128, 014106 (2008).
96 """
98 def __init__(self, dimeratoms, control=None, eigenmode=None, basis=None,
99 **kwargs):
100 if hasattr(dimeratoms, 'get_eigenmode'):
101 self.dimeratoms = dimeratoms
102 else:
103 e = 'The atoms object must be a MinModeAtoms object'
104 raise TypeError(e)
105 self.basis = basis
107 if eigenmode is None:
108 self.eigenmode = self.dimeratoms.get_eigenmode()
109 else:
110 self.eigenmode = eigenmode
112 if control is None:
113 self.control = DimerControl(**kwargs)
114 w = 'Missing control object in ' + self.__class__.__name__ + \
115 '. Using default: DimerControl()'
116 warnings.warn(w, UserWarning)
117 if self.control.logfile is not None:
118 self.control.logfile.write('DIM:WARN: ' + w + '\n')
119 self.control.logfile.flush()
120 else:
121 self.control = control
122 # kwargs must be empty if a control object is supplied
123 for key in kwargs:
124 e = f'__init__() got an unexpected keyword argument \'{(key)}\''
125 raise TypeError(e)
127 self.dR = self.control.get_parameter('dimer_separation')
128 self.logfile = self.control.get_logfile()
130 def converge_to_eigenmode(self):
131 """Perform an eigenmode search."""
132 self.set_up_for_eigenmode_search()
133 stoprot = False
135 # Load the relevant parameters from control
136 f_rot_min = self.control.get_parameter('f_rot_min')
137 f_rot_max = self.control.get_parameter('f_rot_max')
138 trial_angle = self.control.get_parameter('trial_angle')
139 max_num_rot = self.control.get_parameter('max_num_rot')
140 extrapolate = self.control.get_parameter('extrapolate_forces')
142 while not stoprot:
143 if self.forces1E is None:
144 self.update_virtual_forces()
145 else:
146 self.update_virtual_forces(extrapolated_forces=True)
147 self.forces1A = self.forces1
148 self.update_curvature()
149 f_rot_A = self.get_rotational_force()
151 # Pre rotation stop criteria
152 if norm(f_rot_A) <= f_rot_min:
153 self.log(f_rot_A, None)
154 stoprot = True
155 else:
156 n_A = self.eigenmode
157 rot_unit_A = normalize(f_rot_A)
159 # Get the curvature and its derivative
160 c0 = self.get_curvature()
161 c0d = np.vdot((self.forces2 - self.forces1), rot_unit_A) / \
162 self.dR
164 # Trial rotation (no need to store the curvature)
165 # NYI variable trial angles from [3]
166 n_B, rot_unit_B = rotate_vectors(n_A, rot_unit_A, trial_angle)
167 self.eigenmode = n_B
168 self.update_virtual_forces()
169 self.forces1B = self.forces1
171 # Get the curvature's derivative
172 c1d = np.vdot((self.forces2 - self.forces1), rot_unit_B) / \
173 self.dR
175 # Calculate the Fourier coefficients
176 a1 = c0d * cos(2 * trial_angle) - c1d / \
177 (2 * sin(2 * trial_angle))
178 b1 = 0.5 * c0d
179 a0 = 2 * (c0 - a1)
181 # Estimate the rotational angle
182 rotangle = atan(b1 / a1) / 2.0
184 # Make sure that you didn't find a maximum
185 cmin = a0 / 2.0 + a1 * cos(2 * rotangle) + \
186 b1 * sin(2 * rotangle)
187 if c0 < cmin:
188 rotangle += pi / 2.0
190 # Rotate into the (hopefully) lowest eigenmode
191 # NYI Conjugate gradient rotation
192 n_min, _dummy = rotate_vectors(n_A, rot_unit_A, rotangle)
193 self.update_eigenmode(n_min)
195 # Store the curvature estimate instead of the old curvature
196 self.update_curvature(cmin)
198 self.log(f_rot_A, rotangle)
200 # Force extrapolation scheme from [4]
201 if extrapolate:
202 self.forces1E = sin(trial_angle - rotangle) / \
203 sin(trial_angle) * self.forces1A + sin(rotangle) / \
204 sin(trial_angle) * self.forces1B + \
205 (1 - cos(rotangle) - sin(rotangle) *
206 tan(trial_angle / 2.0)) * self.forces0
207 else:
208 self.forces1E = None
210 # Post rotation stop criteria
211 if not stoprot:
212 if self.control.get_counter('rotcount') >= max_num_rot:
213 stoprot = True
214 elif norm(f_rot_A) <= f_rot_max:
215 stoprot = True
217 def log(self, f_rot_A, angle):
218 """Log each rotational step."""
219 # NYI Log for the trial angle
220 if self.logfile is not None:
221 if angle:
222 l = 'DIM:ROT: %7d %9d %9.4f %9.4f %9.4f\n' % \
223 (self.control.get_counter('optcount'),
224 self.control.get_counter('rotcount'),
225 self.get_curvature(), degrees(angle), norm(f_rot_A))
226 else:
227 l = 'DIM:ROT: %7d %9d %9.4f %9s %9.4f\n' % \
228 (self.control.get_counter('optcount'),
229 self.control.get_counter('rotcount'),
230 self.get_curvature(), '---------', norm(f_rot_A))
231 self.logfile.write(l)
233 def get_rotational_force(self):
234 """Calculate the rotational force that acts on the dimer."""
235 rot_force = perpendicular_vector((self.forces1 - self.forces2),
236 self.eigenmode) / (2.0 * self.dR)
237 if self.basis is not None:
238 if (
239 len(self.basis) == len(self.dimeratoms)
240 and len(self.basis[0]) == 3
241 and isinstance(self.basis[0][0], float)):
242 rot_force = perpendicular_vector(rot_force, self.basis)
243 else:
244 for base in self.basis:
245 rot_force = perpendicular_vector(rot_force, base)
246 return rot_force
248 def update_curvature(self, curv=None):
249 """Update the curvature in the MinModeAtoms object."""
250 if curv:
251 self.curvature = curv
252 else:
253 self.curvature = np.vdot((self.forces2 - self.forces1),
254 self.eigenmode) / (2.0 * self.dR)
256 def update_eigenmode(self, eigenmode):
257 """Update the eigenmode in the MinModeAtoms object."""
258 self.eigenmode = eigenmode
259 self.update_virtual_positions()
260 self.control.increment_counter('rotcount')
262 def get_eigenmode(self):
263 """Returns the current eigenmode."""
264 return self.eigenmode
266 def get_curvature(self):
267 """Returns the curvature along the current eigenmode."""
268 return self.curvature
270 def get_control(self):
271 """Return the control object."""
272 return self.control
274 def update_center_forces(self):
275 """Get the forces at the center of the dimer."""
276 self.dimeratoms.set_positions(self.pos0)
277 self.forces0 = self.dimeratoms.get_forces(real=True)
278 self.energy0 = self.dimeratoms.get_potential_energy()
280 def update_virtual_forces(self, extrapolated_forces=False):
281 """Get the forces at the endpoints of the dimer."""
282 self.update_virtual_positions()
284 # Estimate / Calculate the forces at pos1
285 if extrapolated_forces:
286 self.forces1 = self.forces1E.copy()
287 else:
288 self.forces1 = self.dimeratoms.get_forces(real=True, pos=self.pos1)
290 # Estimate / Calculate the forces at pos2
291 if self.control.get_parameter('use_central_forces'):
292 self.forces2 = 2 * self.forces0 - self.forces1
293 else:
294 self.forces2 = self.dimeratoms.get_forces(real=True, pos=self.pos2)
296 def update_virtual_positions(self):
297 """Update the end point positions."""
298 self.pos1 = self.pos0 + self.eigenmode * self.dR
299 self.pos2 = self.pos0 - self.eigenmode * self.dR
301 def set_up_for_eigenmode_search(self):
302 """Before eigenmode search, prepare for rotation."""
303 self.pos0 = self.dimeratoms.get_positions()
304 self.update_center_forces()
305 self.update_virtual_positions()
306 self.control.reset_counter('rotcount')
307 self.forces1E = None
309 def set_up_for_optimization_step(self):
310 """At the end of rotation, prepare for displacement of the dimer."""
311 self.dimeratoms.set_positions(self.pos0)
312 self.forces1E = None
315class MinModeControl(IOContext):
316 """A parent class for controlling minimum mode saddle point searches.
318 Method specific control classes inherit this one. The only thing
319 inheriting classes need to implement are the log() method and
320 the *parameters* class variable with default values for ALL
321 parameters needed by the method in question.
322 When instantiating control classes default parameter values can
323 be overwritten.
325 """
326 parameters: dict[str, Any] = {}
328 def __init__(self, logfile='-', eigenmode_logfile=None, comm=world,
329 **kwargs):
330 # Overwrite the defaults with the input parameters given
331 for key in kwargs:
332 if key not in self.parameters:
333 e = (f'Invalid parameter >>{key}<< with value >>'
334 f'{kwargs[key]!s}<< in {self.__class__.__name__}')
335 raise ValueError(e)
336 else:
337 self.set_parameter(key, kwargs[key], log=False)
339 self.initialize_logfiles(comm=comm, logfile=logfile,
340 eigenmode_logfile=eigenmode_logfile)
341 self.counters = {'forcecalls': 0, 'rotcount': 0, 'optcount': 0}
342 self.log()
344 def initialize_logfiles(self, comm, logfile=None, eigenmode_logfile=None):
345 self.logfile = self.openfile(file=logfile, comm=comm)
346 self.eigenmode_logfile = self.openfile(file=eigenmode_logfile,
347 comm=comm)
349 def log(self, parameter=None):
350 """Log the parameters of the eigenmode search."""
352 def set_parameter(self, parameter, value, log=True):
353 """Change a parameter's value."""
354 if parameter not in self.parameters:
355 e = f'Invalid parameter >>{parameter}<< with value >>{value!s}<<'
356 raise ValueError(e)
357 self.parameters[parameter] = value
358 if log:
359 self.log(parameter)
361 def get_parameter(self, parameter):
362 """Returns the value of a parameter."""
363 if parameter not in self.parameters:
364 e = f'Invalid parameter >>{(parameter)}<<'
365 raise ValueError(e)
366 return self.parameters[parameter]
368 def get_logfile(self):
369 """Returns the log file."""
370 return self.logfile
372 def get_eigenmode_logfile(self):
373 """Returns the eigenmode log file."""
374 return self.eigenmode_logfile
376 def get_counter(self, counter):
377 """Returns a given counter."""
378 return self.counters[counter]
380 def increment_counter(self, counter):
381 """Increment a given counter."""
382 self.counters[counter] += 1
384 def reset_counter(self, counter):
385 """Reset a given counter."""
386 self.counters[counter] = 0
388 def reset_all_counters(self):
389 """Reset all counters."""
390 for key in self.counters:
391 self.counters[key] = 0
394class DimerControl(MinModeControl):
395 """A class that takes care of the parameters needed for a Dimer search.
397 Parameters
398 ----------
400 eigenmode_method: str
401 The name of the eigenmode search method.
402 f_rot_min: float
403 Size of the rotational force under which no rotation will be
404 performed.
405 f_rot_max: float
406 Size of the rotational force under which only one rotation will be
407 performed.
408 max_num_rot: int
409 Maximum number of rotations per optimizer step.
410 trial_angle: float
411 Trial angle for the finite difference estimate of the rotational
412 angle in radians.
413 trial_trans_step: float
414 Trial step size for the MinModeTranslate optimizer.
415 maximum_translation: float
416 Maximum step size and forced step size when the curvature is still
417 positive for the MinModeTranslate optimizer.
418 cg_translation: bool
419 Conjugate Gradient for the MinModeTranslate optimizer.
420 use_central_forces: bool
421 Only calculate the forces at one end of the dimer and extrapolate
422 the forces to the other.
423 dimer_separation: float
424 Separation of the dimer's images.
425 initial_eigenmode_method: str
426 How to construct the initial eigenmode of the dimer. If an eigenmode
427 is given when creating the MinModeAtoms object, this will be ignored.
428 Possible choices are: 'gauss' and 'displacement'
429 extrapolate_forces: bool
430 When more than one rotation is performed, an extrapolation scheme can
431 be used to reduce the number of force evaluations.
432 displacement_method: str
433 How to displace the atoms. Possible choices are 'gauss' and 'vector'.
434 gauss_std: float
435 The standard deviation of the gauss curve used when doing random
436 displacement.
437 order: int
438 How many lowest eigenmodes will be inverted.
439 mask: list of bool
440 Which atoms will be moved during displacement.
441 displacement_center: int or [float, float, float]
442 The center of displacement, nearby atoms will be displaced.
443 displacement_radius: float
444 When choosing which atoms to displace with the *displacement_center*
445 keyword, this decides how many nearby atoms to displace.
446 number_of_displacement_atoms: int
447 The amount of atoms near *displacement_center* to displace.
449 """
450 # Default parameters for the Dimer eigenmode search
451 parameters = {'eigenmode_method': 'dimer',
452 'f_rot_min': 0.1,
453 'f_rot_max': 1.00,
454 'max_num_rot': 1,
455 'trial_angle': pi / 4.0,
456 'trial_trans_step': 0.001,
457 'maximum_translation': 0.1,
458 'cg_translation': True,
459 'use_central_forces': True,
460 'dimer_separation': 0.0001,
461 'initial_eigenmode_method': 'gauss',
462 'extrapolate_forces': False,
463 'displacement_method': 'gauss',
464 'gauss_std': 0.1,
465 'order': 1,
466 'mask': None, # NB mask should not be a "parameter"
467 'displacement_center': None,
468 'displacement_radius': None,
469 'number_of_displacement_atoms': None}
471 # NB: Can maybe put this in EigenmodeSearch and MinModeControl
472 def log(self, parameter=None):
473 """Log the parameters of the eigenmode search."""
474 if self.logfile is not None:
475 if parameter is not None:
476 l = 'DIM:CONTROL: Updated Parameter: {} = {}\n'.format(
477 parameter, str(self.get_parameter(parameter)))
478 else:
479 l = 'MINMODE:METHOD: Dimer\n'
480 l += 'DIM:CONTROL: Search Parameters:\n'
481 l += 'DIM:CONTROL: ------------------\n'
482 for key in self.parameters:
483 l += 'DIM:CONTROL: {} = {}\n'.format(
484 key, str(self.get_parameter(key)))
485 l += 'DIM:CONTROL: ------------------\n'
486 l += 'DIM:ROT: OPT-STEP ROT-STEP CURVATURE ROT-ANGLE ' + \
487 'ROT-FORCE\n'
488 self.logfile.write(l)
489 self.logfile.flush()
492class MinModeAtoms:
493 """Wrapper for Atoms with information related to minimum mode searching.
495 Contains an Atoms object and pipes all unknown function calls to that
496 object.
497 Other information that is stored in this object are the estimate for
498 the lowest eigenvalue, *curvature*, and its corresponding eigenmode,
499 *eigenmode*. Furthermore, the original configuration of the Atoms
500 object is stored for use in multiple minimum mode searches.
501 The forces on the system are modified by inverting the component
502 along the eigenmode estimate. This eventually brings the system to
503 a saddle point.
505 [1]_ [2]_ [3]_ [4]_
507 Parameters
508 ----------
510 atoms : Atoms object
511 A regular Atoms object
512 control : MinModeControl object
513 Contains the parameters necessary for the eigenmode search.
514 If no control object is supplied a default DimerControl
515 will be created and used.
516 mask: list of bool
517 Determines which atoms will be moved when calling displace()
518 random_seed: int
519 The seed used for the random number generator. Defaults to
520 modified version the current time.
522 References
523 ----------
524 .. [1] Henkelman and Jonsson, JCP 111, 7010 (1999)
525 .. [2] Olsen, Kroes, Henkelman, Arnaldsson, and Jonsson, JCP 121,
526 9776 (2004).
527 .. [3] Heyden, Bell, and Keil, JCP 123, 224101 (2005).
528 .. [4] Kastner and Sherwood, JCP 128, 014106 (2008).
530 """
532 def __init__(self, atoms, control=None, eigenmodes=None,
533 random_seed=None, comm=world, **kwargs):
534 self.minmode_init = True
535 self.atoms = atoms
537 # Initialize to None to avoid strange behaviour due to __getattr__
538 self.eigenmodes = eigenmodes
539 self.curvatures = None
541 if control is None:
542 self.control = DimerControl(**kwargs)
543 w = 'Missing control object in ' + self.__class__.__name__ + \
544 '. Using default: DimerControl()'
545 warnings.warn(w, UserWarning)
546 if self.control.logfile is not None:
547 self.control.logfile.write('DIM:WARN: ' + w + '\n')
548 self.control.logfile.flush()
549 else:
550 self.control = control
551 logfile = self.control.get_logfile()
552 mlogfile = self.control.get_eigenmode_logfile()
553 for key in kwargs:
554 if key == 'logfile':
555 logfile = kwargs[key]
556 elif key == 'eigenmode_logfile':
557 mlogfile = kwargs[key]
558 else:
559 self.control.set_parameter(key, kwargs[key])
560 self.control.initialize_logfiles(comm=comm, logfile=logfile,
561 eigenmode_logfile=mlogfile)
563 # Seed the randomness
564 if random_seed is None:
565 t = time.time()
566 if world.size > 1:
567 t = world.sum(t) / world.size
568 # Harvest the latter part of the current time
569 random_seed = int(('%30.9f' % t)[-9:])
570 self.random_state = np.random.RandomState(random_seed)
572 # Check the order
573 self.order = self.control.get_parameter('order')
575 # Construct the curvatures list
576 self.curvatures = [100.0] * self.order
578 # Save the original state of the atoms.
579 self.atoms0 = self.atoms.copy()
580 self.save_original_forces()
582 # Get a reference to the log files
583 self.logfile = self.control.get_logfile()
584 self.mlogfile = self.control.get_eigenmode_logfile()
586 def __ase_optimizable__(self):
587 return DimerOptimizable(self)
589 def save_original_forces(self, force_calculation=False):
590 """Store the forces (and energy) of the original state."""
591 # NB: Would be nice if atoms.copy() took care of this.
592 if self.calc is not None:
593 # Hack because some calculators do not have calculation_required
594 if (hasattr(self.calc, 'calculation_required')
595 and not self.calc.calculation_required(self.atoms,
596 ['energy', 'forces'])) or force_calculation:
597 calc = SinglePointCalculator(
598 self.atoms0,
599 energy=self.atoms.get_potential_energy(),
600 forces=self.atoms.get_forces())
601 self.atoms0.calc = calc
603 def initialize_eigenmodes(self, method=None, eigenmodes=None,
604 gauss_std=None):
605 """Make an initial guess for the eigenmode."""
606 if eigenmodes is None:
607 pos = self.get_positions()
608 old_pos = self.get_original_positions()
609 if method is None:
610 method = \
611 self.control.get_parameter('initial_eigenmode_method')
612 if method.lower() == 'displacement' and (pos - old_pos).any():
613 eigenmode = normalize(pos - old_pos)
614 elif method.lower() == 'gauss':
615 self.displace(log=False, gauss_std=gauss_std,
616 method=method)
617 new_pos = self.get_positions()
618 eigenmode = normalize(new_pos - pos)
619 self.set_positions(pos)
620 else:
621 e = 'initial_eigenmode must use either \'gauss\' or ' + \
622 '\'displacement\', if the latter is used the atoms ' + \
623 'must have moved away from the original positions.' + \
624 f'You have requested \'{method}\'.'
625 raise NotImplementedError(e) # NYI
626 eigenmodes = [eigenmode]
628 # Create random higher order mode guesses
629 if self.order > 1:
630 if len(eigenmodes) == 1:
631 for _ in range(1, self.order):
632 pos = self.get_positions()
633 self.displace(log=False, gauss_std=gauss_std,
634 method=method)
635 new_pos = self.get_positions()
636 eigenmode = normalize(new_pos - pos)
637 self.set_positions(pos)
638 eigenmodes += [eigenmode]
640 self.eigenmodes = eigenmodes
641 # Ensure that the higher order mode guesses are all orthogonal
642 if self.order > 1:
643 for k in range(self.order):
644 self.ensure_eigenmode_orthogonality(k)
645 self.eigenmode_log()
647 # NB maybe this name might be confusing in context to
648 # calc.calculation_required()
649 def calculation_required(self):
650 """Check if a calculation is required."""
651 return self.minmode_init or self.check_atoms != self.atoms
653 def calculate_real_forces_and_energies(self, **kwargs):
654 """Calculate and store the potential energy and forces."""
655 if self.minmode_init:
656 self.minmode_init = False
657 self.initialize_eigenmodes(eigenmodes=self.eigenmodes)
658 self.rotation_required = True
659 self.forces0 = self.atoms.get_forces(**kwargs)
660 self.energy0 = self.atoms.get_potential_energy()
661 self.control.increment_counter('forcecalls')
662 self.check_atoms = self.atoms.copy()
664 def get_potential_energy(self):
665 """Return the potential energy."""
666 if self.calculation_required():
667 self.calculate_real_forces_and_energies()
668 return self.energy0
670 def get_forces(self, real=False, pos=None, **kwargs):
671 """Return the forces, projected or real."""
672 if self.calculation_required() and pos is None:
673 self.calculate_real_forces_and_energies(**kwargs)
674 if real and pos is None:
675 return self.forces0
676 elif real and pos is not None:
677 old_pos = self.atoms.get_positions()
678 self.atoms.set_positions(pos)
679 forces = self.atoms.get_forces()
680 self.control.increment_counter('forcecalls')
681 self.atoms.set_positions(old_pos)
682 return forces
683 else:
684 if self.rotation_required:
685 self.find_eigenmodes(order=self.order)
686 self.eigenmode_log()
687 self.rotation_required = False
688 self.control.increment_counter('optcount')
689 return self.get_projected_forces()
691 def ensure_eigenmode_orthogonality(self, order):
692 mode = self.eigenmodes[order - 1].copy()
693 for k in range(order - 1):
694 mode = perpendicular_vector(mode, self.eigenmodes[k])
695 self.eigenmodes[order - 1] = normalize(mode)
697 def find_eigenmodes(self, order=1):
698 """Launch eigenmode searches."""
699 if self.control.get_parameter('eigenmode_method').lower() != 'dimer':
700 e = 'Only the Dimer control object has been implemented.'
701 raise NotImplementedError(e) # NYI
702 for k in range(order):
703 if k > 0:
704 self.ensure_eigenmode_orthogonality(k + 1)
705 search = DimerEigenmodeSearch(
706 self, self.control,
707 eigenmode=self.eigenmodes[k], basis=self.eigenmodes[:k])
708 search.converge_to_eigenmode()
709 search.set_up_for_optimization_step()
710 self.eigenmodes[k] = search.get_eigenmode()
711 self.curvatures[k] = search.get_curvature()
713 def get_projected_forces(self, pos=None):
714 """Return the projected forces."""
715 if pos is not None:
716 forces = self.get_forces(real=True, pos=pos).copy()
717 else:
718 forces = self.forces0.copy()
720 # Loop through all the eigenmodes
721 # NB: Can this be done with a linear combination, instead?
722 for k, mode in enumerate(self.eigenmodes):
723 # NYI This If statement needs to be overridable in the control
724 if self.get_curvature(order=k) > 0.0 and self.order == 1:
725 forces = -parallel_vector(forces, mode)
726 else:
727 forces -= 2 * parallel_vector(forces, mode)
728 return forces
730 def restore_original_positions(self):
731 """Restore the MinModeAtoms object positions to the original state."""
732 self.atoms.set_positions(self.get_original_positions())
734 def get_barrier_energy(self):
735 """The energy difference between the current and original states"""
736 try:
737 original_energy = self.get_original_potential_energy()
738 dimer_energy = self.get_potential_energy()
739 return dimer_energy - original_energy
740 except RuntimeError:
741 w = 'The potential energy is not available, without further ' + \
742 'calculations, most likely at the original state.'
743 warnings.warn(w, UserWarning)
744 return np.nan
746 def get_control(self):
747 """Return the control object."""
748 return self.control
750 def get_curvature(self, order='max'):
751 """Return the eigenvalue estimate."""
752 if order == 'max':
753 return max(self.curvatures)
754 else:
755 return self.curvatures[order - 1]
757 def get_eigenmode(self, order=1):
758 """Return the current eigenmode guess."""
759 return self.eigenmodes[order - 1]
761 def get_atoms(self):
762 """Return the unextended Atoms object."""
763 return self.atoms
765 def set_atoms(self, atoms):
766 """Set a new Atoms object"""
767 self.atoms = atoms
769 def set_eigenmode(self, eigenmode, order=1):
770 """Set the eigenmode guess."""
771 self.eigenmodes[order - 1] = eigenmode
773 def set_curvature(self, curvature, order=1):
774 """Set the eigenvalue estimate."""
775 self.curvatures[order - 1] = curvature
777 # Pipe all the stuff from Atoms that is not overwritten.
778 # Pipe all requests for get_original_* to self.atoms0.
779 def __getattr__(self, attr):
780 """Return any value of the Atoms object"""
781 if 'original' in attr.split('_'):
782 attr = attr.replace('_original_', '_')
783 return getattr(self.atoms0, attr)
784 else:
785 return getattr(self.atoms, attr)
787 def __len__(self):
788 return len(self.atoms)
790 def displace(self, displacement_vector=None, mask=None, method=None,
791 displacement_center=None, radius=None, number_of_atoms=None,
792 gauss_std=None, mic=True, log=True):
793 """Move the atoms away from their current position.
795 This is one of the essential parts of minimum mode searches.
796 The parameters can all be set in the control object and overwritten
797 when this method is run, apart from *displacement_vector*.
798 It is preferred to modify the control values rather than those here
799 in order for the correct ones to show up in the log file.
801 *method* can be either 'gauss' for random displacement or 'vector'
802 to perform a predefined displacement.
804 *gauss_std* is the standard deviation of the gauss curve that is
805 used for random displacement.
807 *displacement_center* can be either the number of an atom or a 3D
808 position. It must be accompanied by a *radius* (all atoms within it
809 will be displaced) or a *number_of_atoms* which decides how many of
810 the closest atoms will be displaced.
812 *mic* controls the usage of the Minimum Image Convention.
814 If both *mask* and *displacement_center* are used, the atoms marked
815 as False in the *mask* will not be affected even though they are
816 within reach of the *displacement_center*.
818 The parameters priority order:
819 1) displacement_vector
820 2) mask
821 3) displacement_center (with radius and/or number_of_atoms)
823 If both *radius* and *number_of_atoms* are supplied with
824 *displacement_center*, only atoms that fulfill both criteria will
825 be displaced.
827 """
829 # Fetch the default values from the control
830 if mask is None:
831 mask = self.control.get_parameter('mask')
832 if method is None:
833 method = self.control.get_parameter('displacement_method')
834 if gauss_std is None:
835 gauss_std = self.control.get_parameter('gauss_std')
836 if displacement_center is None:
837 displacement_center = \
838 self.control.get_parameter('displacement_center')
839 if radius is None:
840 radius = self.control.get_parameter('displacement_radius')
841 if number_of_atoms is None:
842 number_of_atoms = \
843 self.control.get_parameter('number_of_displacement_atoms')
845 # Check for conflicts
846 if displacement_vector is not None and method.lower() != 'vector':
847 e = 'displacement_vector was supplied but a different method ' + \
848 f'(\'{method!s}\') was chosen.\n'
849 raise ValueError(e)
850 elif displacement_vector is None and method.lower() == 'vector':
851 e = 'A displacement_vector must be supplied when using ' + \
852 f'method = \'{method!s}\'.\n'
853 raise ValueError(e)
854 elif displacement_center is not None and radius is None and \
855 number_of_atoms is None:
856 e = 'When displacement_center is chosen, either radius or ' + \
857 'number_of_atoms must be supplied.\n'
858 raise ValueError(e)
860 # Set up the center of displacement mask (c_mask)
861 if displacement_center is not None:
862 c = displacement_center
863 # Construct a distance list
864 # The center is an atom
865 if isinstance(c, int):
866 # Parse negative indexes
867 c = displacement_center % len(self)
868 d = [(k, self.get_distance(k, c, mic=mic)) for k in
869 range(len(self))]
870 # The center is a position in 3D space
871 elif len(c) == 3 and [type(c_k) for c_k in c] == [float] * 3:
872 # NB: MIC is not considered.
873 d = [(k, norm(self.get_positions()[k] - c))
874 for k in range(len(self))]
875 else:
876 e = 'displacement_center must be either the number of an ' + \
877 'atom in MinModeAtoms object or a 3D position ' + \
878 '(3-tuple of floats).'
879 raise ValueError(e)
881 # Set up the mask
882 if radius is not None:
883 r_mask = [dist[1] < radius for dist in d]
884 else:
885 r_mask = [True for _ in range(len(self))]
887 if number_of_atoms is not None:
888 d_sorted = [n[0] for n in sorted(d, key=lambda k: k[1])]
889 n_nearest = d_sorted[:number_of_atoms]
890 n_mask = [k in n_nearest for k in range(len(self))]
891 else:
892 n_mask = [True for _ in range(len(self))]
894 # Resolve n_mask / r_mask conflicts
895 c_mask = [n_mask[k] and r_mask[k] for k in range(len(self))]
896 else:
897 c_mask = None
899 # Set up a True mask if there is no mask supplied
900 if mask is None:
901 mask = [True for _ in range(len(self))]
902 if c_mask is None:
903 w = 'It was not possible to figure out which atoms to ' + \
904 'displace, Will try to displace all atoms.\n'
905 warnings.warn(w, UserWarning)
906 if self.logfile is not None:
907 self.logfile.write('MINMODE:WARN: ' + w + '\n')
908 self.logfile.flush()
910 # Resolve mask / c_mask conflicts
911 if c_mask is not None:
912 mask = [mask[k] and c_mask[k] for k in range(len(self))]
914 if displacement_vector is None:
915 displacement_vector = []
916 for k in range(len(self)):
917 if mask[k]:
918 diff_line = []
919 for _ in range(3):
920 if method.lower() == 'gauss':
921 if not gauss_std:
922 gauss_std = \
923 self.control.get_parameter('gauss_std')
924 diff = self.random_state.normal(0.0, gauss_std)
925 else:
926 e = f'Invalid displacement method >>{method!s}<<'
927 raise ValueError(e)
928 diff_line.append(diff)
929 displacement_vector.append(diff_line)
930 else:
931 displacement_vector.append([0.0] * 3)
933 # Remove displacement of masked atoms
934 for k in range(len(mask)):
935 if not mask[k]:
936 displacement_vector[k] = [0.0] * 3
938 # Perform the displacement and log it
939 if log:
940 pos0 = self.get_positions()
941 self.set_positions(self.get_positions() + displacement_vector)
942 if log:
943 parameters = {'mask': mask,
944 'displacement_method': method,
945 'gauss_std': gauss_std,
946 'displacement_center': displacement_center,
947 'displacement_radius': radius,
948 'number_of_displacement_atoms': number_of_atoms}
949 self.displacement_log(self.get_positions() - pos0, parameters)
951 def eigenmode_log(self):
952 """Log the eigenmodes (eigenmode estimates)"""
953 if self.mlogfile is not None:
954 l = 'MINMODE:MODE: Optimization Step: %i\n' % \
955 (self.control.get_counter('optcount'))
956 for m_num, mode in enumerate(self.eigenmodes):
957 l += 'MINMODE:MODE: Order: %i\n' % m_num
958 for k in range(len(mode)):
959 l += 'MINMODE:MODE: %7i %15.8f %15.8f %15.8f\n' % (
960 k, mode[k][0], mode[k][1], mode[k][2])
961 self.mlogfile.write(l)
962 self.mlogfile.flush()
964 def displacement_log(self, displacement_vector, parameters):
965 """Log the displacement"""
966 if self.logfile is not None:
967 lp = 'MINMODE:DISP: Parameters, different from the control:\n'
968 mod_para = False
969 for key in parameters:
970 if parameters[key] != self.control.get_parameter(key):
971 lp += 'MINMODE:DISP: {} = {}\n'.format(str(key),
972 str(parameters[key]))
973 mod_para = True
974 if mod_para:
975 l = lp
976 else:
977 l = ''
978 for k in range(len(displacement_vector)):
979 l += 'MINMODE:DISP: %7i %15.8f %15.8f %15.8f\n' % (
980 k,
981 displacement_vector[k][0], displacement_vector[k][1],
982 displacement_vector[k][2])
983 self.logfile.write(l)
984 self.logfile.flush()
986 def summarize(self):
987 """Summarize the Minimum mode search."""
988 if self.logfile is None:
989 logfile = sys.stdout
990 else:
991 logfile = self.logfile
993 c = self.control
994 label = 'MINMODE:SUMMARY: '
996 l = label + '-------------------------\n'
997 l += label + 'Barrier: %16.4f\n' % self.get_barrier_energy()
998 l += label + 'Curvature: %14.4f\n' % self.get_curvature()
999 l += label + 'Optimizer steps: %8i\n' % c.get_counter('optcount')
1000 l += label + 'Forcecalls: %13i\n' % c.get_counter('forcecalls')
1001 l += label + '-------------------------\n'
1003 logfile.write(l)
1006class MinModeTranslate(Optimizer):
1007 """An Optimizer specifically tailored to minimum mode following."""
1009 def __init__(self, dimeratoms, logfile='-', trajectory=None):
1010 Optimizer.__init__(self, dimeratoms, None, logfile, trajectory)
1012 self.control = dimeratoms.get_control()
1013 self.dimeratoms = dimeratoms
1015 # Make a header for the log
1016 if self.logfile is not None:
1017 l = ''
1018 if isinstance(self.control, DimerControl):
1019 l = 'MinModeTranslate: STEP TIME ENERGY ' + \
1020 'MAX-FORCE STEPSIZE CURVATURE ROT-STEPS\n'
1021 self.logfile.write(l)
1023 # Load the relevant parameters from control
1024 self.cg_on = self.control.get_parameter('cg_translation')
1025 self.trial_step = self.control.get_parameter('trial_trans_step')
1026 self.max_step = self.control.get_parameter('maximum_translation')
1028 # Start conjugate gradient
1029 if self.cg_on:
1030 self.cg_init = True
1032 def initialize(self):
1033 """Set initial values."""
1034 self.r0 = None
1035 self.f0 = None
1037 def step(self, f=None):
1038 """Perform the optimization step."""
1040 atoms = self.dimeratoms
1041 if f is None:
1042 f = atoms.get_forces()
1043 r = atoms.get_positions()
1044 curv = atoms.get_curvature()
1045 f0p = f.copy()
1046 r0 = r.copy()
1047 direction = f0p.copy()
1048 if self.cg_on:
1049 direction = self.get_cg_direction(direction)
1050 direction = normalize(direction)
1051 if curv > 0.0:
1052 step = direction * self.max_step
1053 else:
1054 r0t = r0 + direction * self.trial_step
1055 f0tp = self.dimeratoms.get_projected_forces(r0t)
1056 F = np.vdot((f0tp + f0p), direction) / 2.0
1057 C = np.vdot((f0tp - f0p), direction) / self.trial_step
1058 step = (-F / C + self.trial_step / 2.0) * direction
1059 if norm(step) > self.max_step:
1060 step = direction * self.max_step
1061 self.log(f0p, norm(step))
1063 atoms.set_positions(r + step)
1065 self.f0 = f.flat.copy()
1066 self.r0 = r.flat.copy()
1068 def get_cg_direction(self, direction):
1069 """Apply the Conjugate Gradient algorithm to the step direction."""
1070 if self.cg_init:
1071 self.cg_init = False
1072 self.direction_old = direction.copy()
1073 self.cg_direction = direction.copy()
1074 old_norm = np.vdot(self.direction_old, self.direction_old)
1075 # Polak-Ribiere Conjugate Gradient
1076 if old_norm != 0.0:
1077 betaPR = np.vdot(direction, (direction - self.direction_old)) / \
1078 old_norm
1079 else:
1080 betaPR = 0.0
1081 if betaPR < 0.0:
1082 betaPR = 0.0
1083 self.cg_direction = direction + self.cg_direction * betaPR
1084 self.direction_old = direction.copy()
1085 return self.cg_direction.copy()
1087 def log(self, gradient, stepsize=None):
1088 """Log each step of the optimization."""
1089 f = self.dimeratoms.get_forces()
1090 if self.logfile is not None:
1091 T = time.localtime()
1092 e = self.dimeratoms.get_potential_energy()
1093 fmax = sqrt((f**2).sum(axis=1).max())
1094 rotsteps = self.dimeratoms.control.get_counter('rotcount')
1095 curvature = self.dimeratoms.get_curvature()
1096 l = ''
1097 if stepsize:
1098 if isinstance(self.control, DimerControl):
1099 l = '%s: %4d %02d:%02d:%02d %15.6f %12.4f %12.6f ' \
1100 '%12.6f %10d\n' % (
1101 'MinModeTranslate', self.nsteps,
1102 T[3], T[4], T[5], e, fmax, stepsize, curvature,
1103 rotsteps)
1104 else:
1105 if isinstance(self.control, DimerControl):
1106 l = '%s: %4d %02d:%02d:%02d %15.6f %12.4f %s ' \
1107 '%12.6f %10d\n' % (
1108 'MinModeTranslate', self.nsteps,
1109 T[3], T[4], T[5], e, fmax, ' --------',
1110 curvature, rotsteps)
1111 self.logfile.write(l)
1114def read_eigenmode(mlog, index=-1):
1115 """Read an eigenmode.
1116 To access the pre optimization eigenmode set index = 'null'.
1118 """
1119 mlog_is_str = isinstance(mlog, str)
1120 if mlog_is_str:
1121 fd = open(mlog)
1122 else:
1123 fd = mlog
1125 lines = fd.readlines()
1127 # Detect the amount of atoms and iterations
1128 k = 2
1129 while lines[k].split()[1].lower() not in ['optimization', 'order']:
1130 k += 1
1131 n = k - 2
1132 n_itr = (len(lines) // (n + 1)) - 2
1134 # Locate the correct image.
1135 if isinstance(index, str):
1136 if index.lower() == 'null':
1137 i = 0
1138 else:
1139 i = int(index) + 1
1140 else:
1141 if index >= 0:
1142 i = index + 1
1143 else:
1144 if index < -n_itr - 1:
1145 raise IndexError('list index out of range')
1146 else:
1147 i = index
1149 mode = np.ndarray(shape=(n, 3), dtype=float)
1150 for k_atom, k in enumerate(range(1, n + 1)):
1151 line = lines[i * (n + 1) + k].split()
1152 for k_dim in range(3):
1153 mode[k_atom][k_dim] = float(line[k_dim + 2])
1154 if mlog_is_str:
1155 fd.close()
1157 return mode
1160# Aliases
1161DimerAtoms = MinModeAtoms
1162DimerTranslate = MinModeTranslate