Coverage for /builds/ase/ase/ase/mep/dimer.py: 71.94%
588 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +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)
229 self.logfile.flush()
231 def get_rotational_force(self):
232 """Calculate the rotational force that acts on the dimer."""
233 rot_force = perpendicular_vector((self.forces1 - self.forces2),
234 self.eigenmode) / (2.0 * self.dR)
235 if self.basis is not None:
236 if (
237 len(self.basis) == len(self.dimeratoms)
238 and len(self.basis[0]) == 3
239 and isinstance(self.basis[0][0], float)):
240 rot_force = perpendicular_vector(rot_force, self.basis)
241 else:
242 for base in self.basis:
243 rot_force = perpendicular_vector(rot_force, base)
244 return rot_force
246 def update_curvature(self, curv=None):
247 """Update the curvature in the MinModeAtoms object."""
248 if curv:
249 self.curvature = curv
250 else:
251 self.curvature = np.vdot((self.forces2 - self.forces1),
252 self.eigenmode) / (2.0 * self.dR)
254 def update_eigenmode(self, eigenmode):
255 """Update the eigenmode in the MinModeAtoms object."""
256 self.eigenmode = eigenmode
257 self.update_virtual_positions()
258 self.control.increment_counter('rotcount')
260 def get_eigenmode(self):
261 """Returns the current eigenmode."""
262 return self.eigenmode
264 def get_curvature(self):
265 """Returns the curvature along the current eigenmode."""
266 return self.curvature
268 def get_control(self):
269 """Return the control object."""
270 return self.control
272 def update_center_forces(self):
273 """Get the forces at the center of the dimer."""
274 self.dimeratoms.set_positions(self.pos0)
275 self.forces0 = self.dimeratoms.get_forces(real=True)
276 self.energy0 = self.dimeratoms.get_potential_energy()
278 def update_virtual_forces(self, extrapolated_forces=False):
279 """Get the forces at the endpoints of the dimer."""
280 self.update_virtual_positions()
282 # Estimate / Calculate the forces at pos1
283 if extrapolated_forces:
284 self.forces1 = self.forces1E.copy()
285 else:
286 self.forces1 = self.dimeratoms.get_forces(real=True, pos=self.pos1)
288 # Estimate / Calculate the forces at pos2
289 if self.control.get_parameter('use_central_forces'):
290 self.forces2 = 2 * self.forces0 - self.forces1
291 else:
292 self.forces2 = self.dimeratoms.get_forces(real=True, pos=self.pos2)
294 def update_virtual_positions(self):
295 """Update the end point positions."""
296 self.pos1 = self.pos0 + self.eigenmode * self.dR
297 self.pos2 = self.pos0 - self.eigenmode * self.dR
299 def set_up_for_eigenmode_search(self):
300 """Before eigenmode search, prepare for rotation."""
301 self.pos0 = self.dimeratoms.get_positions()
302 self.update_center_forces()
303 self.update_virtual_positions()
304 self.control.reset_counter('rotcount')
305 self.forces1E = None
307 def set_up_for_optimization_step(self):
308 """At the end of rotation, prepare for displacement of the dimer."""
309 self.dimeratoms.set_positions(self.pos0)
310 self.forces1E = None
313class MinModeControl(IOContext):
314 """A parent class for controlling minimum mode saddle point searches.
316 Method specific control classes inherit this one. The only thing
317 inheriting classes need to implement are the log() method and
318 the *parameters* class variable with default values for ALL
319 parameters needed by the method in question.
320 When instantiating control classes default parameter values can
321 be overwritten.
323 """
324 parameters: Dict[str, Any] = {}
326 def __init__(self, logfile='-', eigenmode_logfile=None, comm=world,
327 **kwargs):
328 # Overwrite the defaults with the input parameters given
329 for key in kwargs:
330 if key not in self.parameters:
331 e = (f'Invalid parameter >>{key}<< with value >>'
332 f'{kwargs[key]!s}<< in {self.__class__.__name__}')
333 raise ValueError(e)
334 else:
335 self.set_parameter(key, kwargs[key], log=False)
337 self.initialize_logfiles(comm=comm, logfile=logfile,
338 eigenmode_logfile=eigenmode_logfile)
339 self.counters = {'forcecalls': 0, 'rotcount': 0, 'optcount': 0}
340 self.log()
342 def initialize_logfiles(self, comm, logfile=None, eigenmode_logfile=None):
343 self.logfile = self.openfile(file=logfile, comm=comm)
344 self.eigenmode_logfile = self.openfile(file=eigenmode_logfile,
345 comm=comm)
347 def log(self, parameter=None):
348 """Log the parameters of the eigenmode search."""
350 def set_parameter(self, parameter, value, log=True):
351 """Change a parameter's value."""
352 if parameter not in self.parameters:
353 e = f'Invalid parameter >>{parameter}<< with value >>{value!s}<<'
354 raise ValueError(e)
355 self.parameters[parameter] = value
356 if log:
357 self.log(parameter)
359 def get_parameter(self, parameter):
360 """Returns the value of a parameter."""
361 if parameter not in self.parameters:
362 e = f'Invalid parameter >>{(parameter)}<<'
363 raise ValueError(e)
364 return self.parameters[parameter]
366 def get_logfile(self):
367 """Returns the log file."""
368 return self.logfile
370 def get_eigenmode_logfile(self):
371 """Returns the eigenmode log file."""
372 return self.eigenmode_logfile
374 def get_counter(self, counter):
375 """Returns a given counter."""
376 return self.counters[counter]
378 def increment_counter(self, counter):
379 """Increment a given counter."""
380 self.counters[counter] += 1
382 def reset_counter(self, counter):
383 """Reset a given counter."""
384 self.counters[counter] = 0
386 def reset_all_counters(self):
387 """Reset all counters."""
388 for key in self.counters:
389 self.counters[key] = 0
392class DimerControl(MinModeControl):
393 """A class that takes care of the parameters needed for a Dimer search.
395 Parameters:
397 eigenmode_method: str
398 The name of the eigenmode search method.
399 f_rot_min: float
400 Size of the rotational force under which no rotation will be
401 performed.
402 f_rot_max: float
403 Size of the rotational force under which only one rotation will be
404 performed.
405 max_num_rot: int
406 Maximum number of rotations per optimizer step.
407 trial_angle: float
408 Trial angle for the finite difference estimate of the rotational
409 angle in radians.
410 trial_trans_step: float
411 Trial step size for the MinModeTranslate optimizer.
412 maximum_translation: float
413 Maximum step size and forced step size when the curvature is still
414 positive for the MinModeTranslate optimizer.
415 cg_translation: bool
416 Conjugate Gradient for the MinModeTranslate optimizer.
417 use_central_forces: bool
418 Only calculate the forces at one end of the dimer and extrapolate
419 the forces to the other.
420 dimer_separation: float
421 Separation of the dimer's images.
422 initial_eigenmode_method: str
423 How to construct the initial eigenmode of the dimer. If an eigenmode
424 is given when creating the MinModeAtoms object, this will be ignored.
425 Possible choices are: 'gauss' and 'displacement'
426 extrapolate_forces: bool
427 When more than one rotation is performed, an extrapolation scheme can
428 be used to reduce the number of force evaluations.
429 displacement_method: str
430 How to displace the atoms. Possible choices are 'gauss' and 'vector'.
431 gauss_std: float
432 The standard deviation of the gauss curve used when doing random
433 displacement.
434 order: int
435 How many lowest eigenmodes will be inverted.
436 mask: list of bool
437 Which atoms will be moved during displacement.
438 displacement_center: int or [float, float, float]
439 The center of displacement, nearby atoms will be displaced.
440 displacement_radius: float
441 When choosing which atoms to displace with the *displacement_center*
442 keyword, this decides how many nearby atoms to displace.
443 number_of_displacement_atoms: int
444 The amount of atoms near *displacement_center* to displace.
446 """
447 # Default parameters for the Dimer eigenmode search
448 parameters = {'eigenmode_method': 'dimer',
449 'f_rot_min': 0.1,
450 'f_rot_max': 1.00,
451 'max_num_rot': 1,
452 'trial_angle': pi / 4.0,
453 'trial_trans_step': 0.001,
454 'maximum_translation': 0.1,
455 'cg_translation': True,
456 'use_central_forces': True,
457 'dimer_separation': 0.0001,
458 'initial_eigenmode_method': 'gauss',
459 'extrapolate_forces': False,
460 'displacement_method': 'gauss',
461 'gauss_std': 0.1,
462 'order': 1,
463 'mask': None, # NB mask should not be a "parameter"
464 'displacement_center': None,
465 'displacement_radius': None,
466 'number_of_displacement_atoms': None}
468 # NB: Can maybe put this in EigenmodeSearch and MinModeControl
469 def log(self, parameter=None):
470 """Log the parameters of the eigenmode search."""
471 if self.logfile is not None:
472 if parameter is not None:
473 l = 'DIM:CONTROL: Updated Parameter: {} = {}\n'.format(
474 parameter, str(self.get_parameter(parameter)))
475 else:
476 l = 'MINMODE:METHOD: Dimer\n'
477 l += 'DIM:CONTROL: Search Parameters:\n'
478 l += 'DIM:CONTROL: ------------------\n'
479 for key in self.parameters:
480 l += 'DIM:CONTROL: {} = {}\n'.format(
481 key, str(self.get_parameter(key)))
482 l += 'DIM:CONTROL: ------------------\n'
483 l += 'DIM:ROT: OPT-STEP ROT-STEP CURVATURE ROT-ANGLE ' + \
484 'ROT-FORCE\n'
485 self.logfile.write(l)
486 self.logfile.flush()
489class MinModeAtoms:
490 """Wrapper for Atoms with information related to minimum mode searching.
492 Contains an Atoms object and pipes all unknown function calls to that
493 object.
494 Other information that is stored in this object are the estimate for
495 the lowest eigenvalue, *curvature*, and its corresponding eigenmode,
496 *eigenmode*. Furthermore, the original configuration of the Atoms
497 object is stored for use in multiple minimum mode searches.
498 The forces on the system are modified by inverting the component
499 along the eigenmode estimate. This eventually brings the system to
500 a saddle point.
502 Parameters:
504 atoms : Atoms object
505 A regular Atoms object
506 control : MinModeControl object
507 Contains the parameters necessary for the eigenmode search.
508 If no control object is supplied a default DimerControl
509 will be created and used.
510 mask: list of bool
511 Determines which atoms will be moved when calling displace()
512 random_seed: int
513 The seed used for the random number generator. Defaults to
514 modified version the current time.
516 References: [1]_ [2]_ [3]_ [4]_
518 .. [1] Henkelman and Jonsson, JCP 111, 7010 (1999)
519 .. [2] Olsen, Kroes, Henkelman, Arnaldsson, and Jonsson, JCP 121,
520 9776 (2004).
521 .. [3] Heyden, Bell, and Keil, JCP 123, 224101 (2005).
522 .. [4] Kastner and Sherwood, JCP 128, 014106 (2008).
524 """
526 def __init__(self, atoms, control=None, eigenmodes=None,
527 random_seed=None, comm=world, **kwargs):
528 self.minmode_init = True
529 self.atoms = atoms
531 # Initialize to None to avoid strange behaviour due to __getattr__
532 self.eigenmodes = eigenmodes
533 self.curvatures = None
535 if control is None:
536 self.control = DimerControl(**kwargs)
537 w = 'Missing control object in ' + self.__class__.__name__ + \
538 '. Using default: DimerControl()'
539 warnings.warn(w, UserWarning)
540 if self.control.logfile is not None:
541 self.control.logfile.write('DIM:WARN: ' + w + '\n')
542 self.control.logfile.flush()
543 else:
544 self.control = control
545 logfile = self.control.get_logfile()
546 mlogfile = self.control.get_eigenmode_logfile()
547 for key in kwargs:
548 if key == 'logfile':
549 logfile = kwargs[key]
550 elif key == 'eigenmode_logfile':
551 mlogfile = kwargs[key]
552 else:
553 self.control.set_parameter(key, kwargs[key])
554 self.control.initialize_logfiles(comm=comm, logfile=logfile,
555 eigenmode_logfile=mlogfile)
557 # Seed the randomness
558 if random_seed is None:
559 t = time.time()
560 if world.size > 1:
561 t = world.sum(t) / world.size
562 # Harvest the latter part of the current time
563 random_seed = int(('%30.9f' % t)[-9:])
564 self.random_state = np.random.RandomState(random_seed)
566 # Check the order
567 self.order = self.control.get_parameter('order')
569 # Construct the curvatures list
570 self.curvatures = [100.0] * self.order
572 # Save the original state of the atoms.
573 self.atoms0 = self.atoms.copy()
574 self.save_original_forces()
576 # Get a reference to the log files
577 self.logfile = self.control.get_logfile()
578 self.mlogfile = self.control.get_eigenmode_logfile()
580 def __ase_optimizable__(self):
581 return DimerOptimizable(self)
583 def save_original_forces(self, force_calculation=False):
584 """Store the forces (and energy) of the original state."""
585 # NB: Would be nice if atoms.copy() took care of this.
586 if self.calc is not None:
587 # Hack because some calculators do not have calculation_required
588 if (hasattr(self.calc, 'calculation_required')
589 and not self.calc.calculation_required(self.atoms,
590 ['energy', 'forces'])) or force_calculation:
591 calc = SinglePointCalculator(
592 self.atoms0,
593 energy=self.atoms.get_potential_energy(),
594 forces=self.atoms.get_forces())
595 self.atoms0.calc = calc
597 def initialize_eigenmodes(self, method=None, eigenmodes=None,
598 gauss_std=None):
599 """Make an initial guess for the eigenmode."""
600 if eigenmodes is None:
601 pos = self.get_positions()
602 old_pos = self.get_original_positions()
603 if method is None:
604 method = \
605 self.control.get_parameter('initial_eigenmode_method')
606 if method.lower() == 'displacement' and (pos - old_pos).any():
607 eigenmode = normalize(pos - old_pos)
608 elif method.lower() == 'gauss':
609 self.displace(log=False, gauss_std=gauss_std,
610 method=method)
611 new_pos = self.get_positions()
612 eigenmode = normalize(new_pos - pos)
613 self.set_positions(pos)
614 else:
615 e = 'initial_eigenmode must use either \'gauss\' or ' + \
616 '\'displacement\', if the latter is used the atoms ' + \
617 'must have moved away from the original positions.' + \
618 f'You have requested \'{method}\'.'
619 raise NotImplementedError(e) # NYI
620 eigenmodes = [eigenmode]
622 # Create random higher order mode guesses
623 if self.order > 1:
624 if len(eigenmodes) == 1:
625 for _ in range(1, self.order):
626 pos = self.get_positions()
627 self.displace(log=False, gauss_std=gauss_std,
628 method=method)
629 new_pos = self.get_positions()
630 eigenmode = normalize(new_pos - pos)
631 self.set_positions(pos)
632 eigenmodes += [eigenmode]
634 self.eigenmodes = eigenmodes
635 # Ensure that the higher order mode guesses are all orthogonal
636 if self.order > 1:
637 for k in range(self.order):
638 self.ensure_eigenmode_orthogonality(k)
639 self.eigenmode_log()
641 # NB maybe this name might be confusing in context to
642 # calc.calculation_required()
643 def calculation_required(self):
644 """Check if a calculation is required."""
645 return self.minmode_init or self.check_atoms != self.atoms
647 def calculate_real_forces_and_energies(self, **kwargs):
648 """Calculate and store the potential energy and forces."""
649 if self.minmode_init:
650 self.minmode_init = False
651 self.initialize_eigenmodes(eigenmodes=self.eigenmodes)
652 self.rotation_required = True
653 self.forces0 = self.atoms.get_forces(**kwargs)
654 self.energy0 = self.atoms.get_potential_energy()
655 self.control.increment_counter('forcecalls')
656 self.check_atoms = self.atoms.copy()
658 def get_potential_energy(self):
659 """Return the potential energy."""
660 if self.calculation_required():
661 self.calculate_real_forces_and_energies()
662 return self.energy0
664 def get_forces(self, real=False, pos=None, **kwargs):
665 """Return the forces, projected or real."""
666 if self.calculation_required() and pos is None:
667 self.calculate_real_forces_and_energies(**kwargs)
668 if real and pos is None:
669 return self.forces0
670 elif real and pos is not None:
671 old_pos = self.atoms.get_positions()
672 self.atoms.set_positions(pos)
673 forces = self.atoms.get_forces()
674 self.control.increment_counter('forcecalls')
675 self.atoms.set_positions(old_pos)
676 return forces
677 else:
678 if self.rotation_required:
679 self.find_eigenmodes(order=self.order)
680 self.eigenmode_log()
681 self.rotation_required = False
682 self.control.increment_counter('optcount')
683 return self.get_projected_forces()
685 def ensure_eigenmode_orthogonality(self, order):
686 mode = self.eigenmodes[order - 1].copy()
687 for k in range(order - 1):
688 mode = perpendicular_vector(mode, self.eigenmodes[k])
689 self.eigenmodes[order - 1] = normalize(mode)
691 def find_eigenmodes(self, order=1):
692 """Launch eigenmode searches."""
693 if self.control.get_parameter('eigenmode_method').lower() != 'dimer':
694 e = 'Only the Dimer control object has been implemented.'
695 raise NotImplementedError(e) # NYI
696 for k in range(order):
697 if k > 0:
698 self.ensure_eigenmode_orthogonality(k + 1)
699 search = DimerEigenmodeSearch(
700 self, self.control,
701 eigenmode=self.eigenmodes[k], basis=self.eigenmodes[:k])
702 search.converge_to_eigenmode()
703 search.set_up_for_optimization_step()
704 self.eigenmodes[k] = search.get_eigenmode()
705 self.curvatures[k] = search.get_curvature()
707 def get_projected_forces(self, pos=None):
708 """Return the projected forces."""
709 if pos is not None:
710 forces = self.get_forces(real=True, pos=pos).copy()
711 else:
712 forces = self.forces0.copy()
714 # Loop through all the eigenmodes
715 # NB: Can this be done with a linear combination, instead?
716 for k, mode in enumerate(self.eigenmodes):
717 # NYI This If statement needs to be overridable in the control
718 if self.get_curvature(order=k) > 0.0 and self.order == 1:
719 forces = -parallel_vector(forces, mode)
720 else:
721 forces -= 2 * parallel_vector(forces, mode)
722 return forces
724 def restore_original_positions(self):
725 """Restore the MinModeAtoms object positions to the original state."""
726 self.atoms.set_positions(self.get_original_positions())
728 def get_barrier_energy(self):
729 """The energy difference between the current and original states"""
730 try:
731 original_energy = self.get_original_potential_energy()
732 dimer_energy = self.get_potential_energy()
733 return dimer_energy - original_energy
734 except RuntimeError:
735 w = 'The potential energy is not available, without further ' + \
736 'calculations, most likely at the original state.'
737 warnings.warn(w, UserWarning)
738 return np.nan
740 def get_control(self):
741 """Return the control object."""
742 return self.control
744 def get_curvature(self, order='max'):
745 """Return the eigenvalue estimate."""
746 if order == 'max':
747 return max(self.curvatures)
748 else:
749 return self.curvatures[order - 1]
751 def get_eigenmode(self, order=1):
752 """Return the current eigenmode guess."""
753 return self.eigenmodes[order - 1]
755 def get_atoms(self):
756 """Return the unextended Atoms object."""
757 return self.atoms
759 def set_atoms(self, atoms):
760 """Set a new Atoms object"""
761 self.atoms = atoms
763 def set_eigenmode(self, eigenmode, order=1):
764 """Set the eigenmode guess."""
765 self.eigenmodes[order - 1] = eigenmode
767 def set_curvature(self, curvature, order=1):
768 """Set the eigenvalue estimate."""
769 self.curvatures[order - 1] = curvature
771 # Pipe all the stuff from Atoms that is not overwritten.
772 # Pipe all requests for get_original_* to self.atoms0.
773 def __getattr__(self, attr):
774 """Return any value of the Atoms object"""
775 if 'original' in attr.split('_'):
776 attr = attr.replace('_original_', '_')
777 return getattr(self.atoms0, attr)
778 else:
779 return getattr(self.atoms, attr)
781 def __len__(self):
782 return len(self.atoms)
784 def displace(self, displacement_vector=None, mask=None, method=None,
785 displacement_center=None, radius=None, number_of_atoms=None,
786 gauss_std=None, mic=True, log=True):
787 """Move the atoms away from their current position.
789 This is one of the essential parts of minimum mode searches.
790 The parameters can all be set in the control object and overwritten
791 when this method is run, apart from *displacement_vector*.
792 It is preferred to modify the control values rather than those here
793 in order for the correct ones to show up in the log file.
795 *method* can be either 'gauss' for random displacement or 'vector'
796 to perform a predefined displacement.
798 *gauss_std* is the standard deviation of the gauss curve that is
799 used for random displacement.
801 *displacement_center* can be either the number of an atom or a 3D
802 position. It must be accompanied by a *radius* (all atoms within it
803 will be displaced) or a *number_of_atoms* which decides how many of
804 the closest atoms will be displaced.
806 *mic* controls the usage of the Minimum Image Convention.
808 If both *mask* and *displacement_center* are used, the atoms marked
809 as False in the *mask* will not be affected even though they are
810 within reach of the *displacement_center*.
812 The parameters priority order:
813 1) displacement_vector
814 2) mask
815 3) displacement_center (with radius and/or number_of_atoms)
817 If both *radius* and *number_of_atoms* are supplied with
818 *displacement_center*, only atoms that fulfill both criteria will
819 be displaced.
821 """
823 # Fetch the default values from the control
824 if mask is None:
825 mask = self.control.get_parameter('mask')
826 if method is None:
827 method = self.control.get_parameter('displacement_method')
828 if gauss_std is None:
829 gauss_std = self.control.get_parameter('gauss_std')
830 if displacement_center is None:
831 displacement_center = \
832 self.control.get_parameter('displacement_center')
833 if radius is None:
834 radius = self.control.get_parameter('displacement_radius')
835 if number_of_atoms is None:
836 number_of_atoms = \
837 self.control.get_parameter('number_of_displacement_atoms')
839 # Check for conflicts
840 if displacement_vector is not None and method.lower() != 'vector':
841 e = 'displacement_vector was supplied but a different method ' + \
842 f'(\'{method!s}\') was chosen.\n'
843 raise ValueError(e)
844 elif displacement_vector is None and method.lower() == 'vector':
845 e = 'A displacement_vector must be supplied when using ' + \
846 f'method = \'{method!s}\'.\n'
847 raise ValueError(e)
848 elif displacement_center is not None and radius is None and \
849 number_of_atoms is None:
850 e = 'When displacement_center is chosen, either radius or ' + \
851 'number_of_atoms must be supplied.\n'
852 raise ValueError(e)
854 # Set up the center of displacement mask (c_mask)
855 if displacement_center is not None:
856 c = displacement_center
857 # Construct a distance list
858 # The center is an atom
859 if isinstance(c, int):
860 # Parse negative indexes
861 c = displacement_center % len(self)
862 d = [(k, self.get_distance(k, c, mic=mic)) for k in
863 range(len(self))]
864 # The center is a position in 3D space
865 elif len(c) == 3 and [type(c_k) for c_k in c] == [float] * 3:
866 # NB: MIC is not considered.
867 d = [(k, norm(self.get_positions()[k] - c))
868 for k in range(len(self))]
869 else:
870 e = 'displacement_center must be either the number of an ' + \
871 'atom in MinModeAtoms object or a 3D position ' + \
872 '(3-tuple of floats).'
873 raise ValueError(e)
875 # Set up the mask
876 if radius is not None:
877 r_mask = [dist[1] < radius for dist in d]
878 else:
879 r_mask = [True for _ in range(len(self))]
881 if number_of_atoms is not None:
882 d_sorted = [n[0] for n in sorted(d, key=lambda k: k[1])]
883 n_nearest = d_sorted[:number_of_atoms]
884 n_mask = [k in n_nearest for k in range(len(self))]
885 else:
886 n_mask = [True for _ in range(len(self))]
888 # Resolve n_mask / r_mask conflicts
889 c_mask = [n_mask[k] and r_mask[k] for k in range(len(self))]
890 else:
891 c_mask = None
893 # Set up a True mask if there is no mask supplied
894 if mask is None:
895 mask = [True for _ in range(len(self))]
896 if c_mask is None:
897 w = 'It was not possible to figure out which atoms to ' + \
898 'displace, Will try to displace all atoms.\n'
899 warnings.warn(w, UserWarning)
900 if self.logfile is not None:
901 self.logfile.write('MINMODE:WARN: ' + w + '\n')
902 self.logfile.flush()
904 # Resolve mask / c_mask conflicts
905 if c_mask is not None:
906 mask = [mask[k] and c_mask[k] for k in range(len(self))]
908 if displacement_vector is None:
909 displacement_vector = []
910 for k in range(len(self)):
911 if mask[k]:
912 diff_line = []
913 for _ in range(3):
914 if method.lower() == 'gauss':
915 if not gauss_std:
916 gauss_std = \
917 self.control.get_parameter('gauss_std')
918 diff = self.random_state.normal(0.0, gauss_std)
919 else:
920 e = f'Invalid displacement method >>{method!s}<<'
921 raise ValueError(e)
922 diff_line.append(diff)
923 displacement_vector.append(diff_line)
924 else:
925 displacement_vector.append([0.0] * 3)
927 # Remove displacement of masked atoms
928 for k in range(len(mask)):
929 if not mask[k]:
930 displacement_vector[k] = [0.0] * 3
932 # Perform the displacement and log it
933 if log:
934 pos0 = self.get_positions()
935 self.set_positions(self.get_positions() + displacement_vector)
936 if log:
937 parameters = {'mask': mask,
938 'displacement_method': method,
939 'gauss_std': gauss_std,
940 'displacement_center': displacement_center,
941 'displacement_radius': radius,
942 'number_of_displacement_atoms': number_of_atoms}
943 self.displacement_log(self.get_positions() - pos0, parameters)
945 def eigenmode_log(self):
946 """Log the eigenmodes (eigenmode estimates)"""
947 if self.mlogfile is not None:
948 l = 'MINMODE:MODE: Optimization Step: %i\n' % \
949 (self.control.get_counter('optcount'))
950 for m_num, mode in enumerate(self.eigenmodes):
951 l += 'MINMODE:MODE: Order: %i\n' % m_num
952 for k in range(len(mode)):
953 l += 'MINMODE:MODE: %7i %15.8f %15.8f %15.8f\n' % (
954 k, mode[k][0], mode[k][1], mode[k][2])
955 self.mlogfile.write(l)
956 self.mlogfile.flush()
958 def displacement_log(self, displacement_vector, parameters):
959 """Log the displacement"""
960 if self.logfile is not None:
961 lp = 'MINMODE:DISP: Parameters, different from the control:\n'
962 mod_para = False
963 for key in parameters:
964 if parameters[key] != self.control.get_parameter(key):
965 lp += 'MINMODE:DISP: {} = {}\n'.format(str(key),
966 str(parameters[key]))
967 mod_para = True
968 if mod_para:
969 l = lp
970 else:
971 l = ''
972 for k in range(len(displacement_vector)):
973 l += 'MINMODE:DISP: %7i %15.8f %15.8f %15.8f\n' % (
974 k,
975 displacement_vector[k][0], displacement_vector[k][1],
976 displacement_vector[k][2])
977 self.logfile.write(l)
978 self.logfile.flush()
980 def summarize(self):
981 """Summarize the Minimum mode search."""
982 if self.logfile is None:
983 logfile = sys.stdout
984 else:
985 logfile = self.logfile
987 c = self.control
988 label = 'MINMODE:SUMMARY: '
990 l = label + '-------------------------\n'
991 l += label + 'Barrier: %16.4f\n' % self.get_barrier_energy()
992 l += label + 'Curvature: %14.4f\n' % self.get_curvature()
993 l += label + 'Optimizer steps: %8i\n' % c.get_counter('optcount')
994 l += label + 'Forcecalls: %13i\n' % c.get_counter('forcecalls')
995 l += label + '-------------------------\n'
997 logfile.write(l)
1000class MinModeTranslate(Optimizer):
1001 """An Optimizer specifically tailored to minimum mode following."""
1003 def __init__(self, dimeratoms, logfile='-', trajectory=None):
1004 Optimizer.__init__(self, dimeratoms, None, logfile, trajectory)
1006 self.control = dimeratoms.get_control()
1007 self.dimeratoms = dimeratoms
1009 # Make a header for the log
1010 if self.logfile is not None:
1011 l = ''
1012 if isinstance(self.control, DimerControl):
1013 l = 'MinModeTranslate: STEP TIME ENERGY ' + \
1014 'MAX-FORCE STEPSIZE CURVATURE ROT-STEPS\n'
1015 self.logfile.write(l)
1016 self.logfile.flush()
1018 # Load the relevant parameters from control
1019 self.cg_on = self.control.get_parameter('cg_translation')
1020 self.trial_step = self.control.get_parameter('trial_trans_step')
1021 self.max_step = self.control.get_parameter('maximum_translation')
1023 # Start conjugate gradient
1024 if self.cg_on:
1025 self.cg_init = True
1027 def initialize(self):
1028 """Set initial values."""
1029 self.r0 = None
1030 self.f0 = None
1032 def step(self, f=None):
1033 """Perform the optimization step."""
1034 atoms = self.dimeratoms
1035 if f is None:
1036 f = atoms.get_forces()
1037 r = atoms.get_positions()
1038 curv = atoms.get_curvature()
1039 f0p = f.copy()
1040 r0 = r.copy()
1041 direction = f0p.copy()
1042 if self.cg_on:
1043 direction = self.get_cg_direction(direction)
1044 direction = normalize(direction)
1045 if curv > 0.0:
1046 step = direction * self.max_step
1047 else:
1048 r0t = r0 + direction * self.trial_step
1049 f0tp = self.dimeratoms.get_projected_forces(r0t)
1050 F = np.vdot((f0tp + f0p), direction) / 2.0
1051 C = np.vdot((f0tp - f0p), direction) / self.trial_step
1052 step = (-F / C + self.trial_step / 2.0) * direction
1053 if norm(step) > self.max_step:
1054 step = direction * self.max_step
1055 self.log(f0p, norm(step))
1057 atoms.set_positions(r + step)
1059 self.f0 = f.flat.copy()
1060 self.r0 = r.flat.copy()
1062 def get_cg_direction(self, direction):
1063 """Apply the Conjugate Gradient algorithm to the step direction."""
1064 if self.cg_init:
1065 self.cg_init = False
1066 self.direction_old = direction.copy()
1067 self.cg_direction = direction.copy()
1068 old_norm = np.vdot(self.direction_old, self.direction_old)
1069 # Polak-Ribiere Conjugate Gradient
1070 if old_norm != 0.0:
1071 betaPR = np.vdot(direction, (direction - self.direction_old)) / \
1072 old_norm
1073 else:
1074 betaPR = 0.0
1075 if betaPR < 0.0:
1076 betaPR = 0.0
1077 self.cg_direction = direction + self.cg_direction * betaPR
1078 self.direction_old = direction.copy()
1079 return self.cg_direction.copy()
1081 def log(self, gradient, stepsize=None):
1082 """Log each step of the optimization."""
1083 f = self.dimeratoms.get_forces()
1084 if self.logfile is not None:
1085 T = time.localtime()
1086 e = self.dimeratoms.get_potential_energy()
1087 fmax = sqrt((f**2).sum(axis=1).max())
1088 rotsteps = self.dimeratoms.control.get_counter('rotcount')
1089 curvature = self.dimeratoms.get_curvature()
1090 l = ''
1091 if stepsize:
1092 if isinstance(self.control, DimerControl):
1093 l = '%s: %4d %02d:%02d:%02d %15.6f %12.4f %12.6f ' \
1094 '%12.6f %10d\n' % (
1095 'MinModeTranslate', self.nsteps,
1096 T[3], T[4], T[5], e, fmax, stepsize, curvature,
1097 rotsteps)
1098 else:
1099 if isinstance(self.control, DimerControl):
1100 l = '%s: %4d %02d:%02d:%02d %15.6f %12.4f %s ' \
1101 '%12.6f %10d\n' % (
1102 'MinModeTranslate', self.nsteps,
1103 T[3], T[4], T[5], e, fmax, ' --------',
1104 curvature, rotsteps)
1105 self.logfile.write(l)
1106 self.logfile.flush()
1109def read_eigenmode(mlog, index=-1):
1110 """Read an eigenmode.
1111 To access the pre optimization eigenmode set index = 'null'.
1113 """
1114 mlog_is_str = isinstance(mlog, str)
1115 if mlog_is_str:
1116 fd = open(mlog)
1117 else:
1118 fd = mlog
1120 lines = fd.readlines()
1122 # Detect the amount of atoms and iterations
1123 k = 2
1124 while lines[k].split()[1].lower() not in ['optimization', 'order']:
1125 k += 1
1126 n = k - 2
1127 n_itr = (len(lines) // (n + 1)) - 2
1129 # Locate the correct image.
1130 if isinstance(index, str):
1131 if index.lower() == 'null':
1132 i = 0
1133 else:
1134 i = int(index) + 1
1135 else:
1136 if index >= 0:
1137 i = index + 1
1138 else:
1139 if index < -n_itr - 1:
1140 raise IndexError('list index out of range')
1141 else:
1142 i = index
1144 mode = np.ndarray(shape=(n, 3), dtype=float)
1145 for k_atom, k in enumerate(range(1, n + 1)):
1146 line = lines[i * (n + 1) + k].split()
1147 for k_dim in range(3):
1148 mode[k_atom][k_dim] = float(line[k_dim + 2])
1149 if mlog_is_str:
1150 fd.close()
1152 return mode
1155# Aliases
1156DimerAtoms = MinModeAtoms
1157DimerTranslate = MinModeTranslate