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