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