Coverage for /builds/debichem-team/python-ase/ase/mep/dimer.py: 71.99%
589 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-03-06 04:00 +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, comm=world,
325 **kwargs):
326 # Overwrite the defaults with the input parameters given
327 for key in kwargs:
328 if key not in self.parameters:
329 e = (f'Invalid parameter >>{key}<< with value >>'
330 f'{kwargs[key]!s}<< in {self.__class__.__name__}')
331 raise ValueError(e)
332 else:
333 self.set_parameter(key, kwargs[key], log=False)
335 self.initialize_logfiles(comm=comm, logfile=logfile,
336 eigenmode_logfile=eigenmode_logfile)
337 self.counters = {'forcecalls': 0, 'rotcount': 0, 'optcount': 0}
338 self.log()
340 def initialize_logfiles(self, comm, logfile=None, eigenmode_logfile=None):
341 self.logfile = self.openfile(file=logfile, comm=comm)
342 self.eigenmode_logfile = self.openfile(file=eigenmode_logfile,
343 comm=comm)
345 def log(self, parameter=None):
346 """Log the parameters of the eigenmode search."""
348 def set_parameter(self, parameter, value, log=True):
349 """Change a parameter's value."""
350 if parameter not in self.parameters:
351 e = f'Invalid parameter >>{parameter}<< with value >>{value!s}<<'
352 raise ValueError(e)
353 self.parameters[parameter] = value
354 if log:
355 self.log(parameter)
357 def get_parameter(self, parameter):
358 """Returns the value of a parameter."""
359 if parameter not in self.parameters:
360 e = f'Invalid parameter >>{(parameter)}<<'
361 raise ValueError(e)
362 return self.parameters[parameter]
364 def get_logfile(self):
365 """Returns the log file."""
366 return self.logfile
368 def get_eigenmode_logfile(self):
369 """Returns the eigenmode log file."""
370 return self.eigenmode_logfile
372 def get_counter(self, counter):
373 """Returns a given counter."""
374 return self.counters[counter]
376 def increment_counter(self, counter):
377 """Increment a given counter."""
378 self.counters[counter] += 1
380 def reset_counter(self, counter):
381 """Reset a given counter."""
382 self.counters[counter] = 0
384 def reset_all_counters(self):
385 """Reset all counters."""
386 for key in self.counters:
387 self.counters[key] = 0
390class DimerControl(MinModeControl):
391 """A class that takes care of the parameters needed for a Dimer search.
393 Parameters:
395 eigenmode_method: str
396 The name of the eigenmode search method.
397 f_rot_min: float
398 Size of the rotational force under which no rotation will be
399 performed.
400 f_rot_max: float
401 Size of the rotational force under which only one rotation will be
402 performed.
403 max_num_rot: int
404 Maximum number of rotations per optimizer step.
405 trial_angle: float
406 Trial angle for the finite difference estimate of the rotational
407 angle in radians.
408 trial_trans_step: float
409 Trial step size for the MinModeTranslate optimizer.
410 maximum_translation: float
411 Maximum step size and forced step size when the curvature is still
412 positive for the MinModeTranslate optimizer.
413 cg_translation: bool
414 Conjugate Gradient for the MinModeTranslate optimizer.
415 use_central_forces: bool
416 Only calculate the forces at one end of the dimer and extrapolate
417 the forces to the other.
418 dimer_separation: float
419 Separation of the dimer's images.
420 initial_eigenmode_method: str
421 How to construct the initial eigenmode of the dimer. If an eigenmode
422 is given when creating the MinModeAtoms object, this will be ignored.
423 Possible choices are: 'gauss' and 'displacement'
424 extrapolate_forces: bool
425 When more than one rotation is performed, an extrapolation scheme can
426 be used to reduce the number of force evaluations.
427 displacement_method: str
428 How to displace the atoms. Possible choices are 'gauss' and 'vector'.
429 gauss_std: float
430 The standard deviation of the gauss curve used when doing random
431 displacement.
432 order: int
433 How many lowest eigenmodes will be inverted.
434 mask: list of bool
435 Which atoms will be moved during displacement.
436 displacement_center: int or [float, float, float]
437 The center of displacement, nearby atoms will be displaced.
438 displacement_radius: float
439 When choosing which atoms to displace with the *displacement_center*
440 keyword, this decides how many nearby atoms to displace.
441 number_of_displacement_atoms: int
442 The amount of atoms near *displacement_center* to displace.
444 """
445 # Default parameters for the Dimer eigenmode search
446 parameters = {'eigenmode_method': 'dimer',
447 'f_rot_min': 0.1,
448 'f_rot_max': 1.00,
449 'max_num_rot': 1,
450 'trial_angle': pi / 4.0,
451 'trial_trans_step': 0.001,
452 'maximum_translation': 0.1,
453 'cg_translation': True,
454 'use_central_forces': True,
455 'dimer_separation': 0.0001,
456 'initial_eigenmode_method': 'gauss',
457 'extrapolate_forces': False,
458 'displacement_method': 'gauss',
459 'gauss_std': 0.1,
460 'order': 1,
461 'mask': None, # NB mask should not be a "parameter"
462 'displacement_center': None,
463 'displacement_radius': None,
464 'number_of_displacement_atoms': None}
466 # NB: Can maybe put this in EigenmodeSearch and MinModeControl
467 def log(self, parameter=None):
468 """Log the parameters of the eigenmode search."""
469 if self.logfile is not None:
470 if parameter is not None:
471 l = 'DIM:CONTROL: Updated Parameter: {} = {}\n'.format(
472 parameter, str(self.get_parameter(parameter)))
473 else:
474 l = 'MINMODE:METHOD: Dimer\n'
475 l += 'DIM:CONTROL: Search Parameters:\n'
476 l += 'DIM:CONTROL: ------------------\n'
477 for key in self.parameters:
478 l += 'DIM:CONTROL: {} = {}\n'.format(
479 key, str(self.get_parameter(key)))
480 l += 'DIM:CONTROL: ------------------\n'
481 l += 'DIM:ROT: OPT-STEP ROT-STEP CURVATURE ROT-ANGLE ' + \
482 'ROT-FORCE\n'
483 self.logfile.write(l)
484 self.logfile.flush()
487class MinModeAtoms:
488 """Wrapper for Atoms with information related to minimum mode searching.
490 Contains an Atoms object and pipes all unknown function calls to that
491 object.
492 Other information that is stored in this object are the estimate for
493 the lowest eigenvalue, *curvature*, and its corresponding eigenmode,
494 *eigenmode*. Furthermore, the original configuration of the Atoms
495 object is stored for use in multiple minimum mode searches.
496 The forces on the system are modified by inverting the component
497 along the eigenmode estimate. This eventually brings the system to
498 a saddle point.
500 Parameters:
502 atoms : Atoms object
503 A regular Atoms object
504 control : MinModeControl object
505 Contains the parameters necessary for the eigenmode search.
506 If no control object is supplied a default DimerControl
507 will be created and used.
508 mask: list of bool
509 Determines which atoms will be moved when calling displace()
510 random_seed: int
511 The seed used for the random number generator. Defaults to
512 modified version the current time.
514 References: [1]_ [2]_ [3]_ [4]_
516 .. [1] Henkelman and Jonsson, JCP 111, 7010 (1999)
517 .. [2] Olsen, Kroes, Henkelman, Arnaldsson, and Jonsson, JCP 121,
518 9776 (2004).
519 .. [3] Heyden, Bell, and Keil, JCP 123, 224101 (2005).
520 .. [4] Kastner and Sherwood, JCP 128, 014106 (2008).
522 """
524 def __init__(self, atoms, control=None, eigenmodes=None,
525 random_seed=None, comm=world, **kwargs):
526 self.minmode_init = True
527 self.atoms = atoms
529 # Initialize to None to avoid strange behaviour due to __getattr__
530 self.eigenmodes = eigenmodes
531 self.curvatures = None
533 if control is None:
534 self.control = DimerControl(**kwargs)
535 w = 'Missing control object in ' + self.__class__.__name__ + \
536 '. Using default: DimerControl()'
537 warnings.warn(w, UserWarning)
538 if self.control.logfile is not None:
539 self.control.logfile.write('DIM:WARN: ' + w + '\n')
540 self.control.logfile.flush()
541 else:
542 self.control = control
543 logfile = self.control.get_logfile()
544 mlogfile = self.control.get_eigenmode_logfile()
545 for key in kwargs:
546 if key == 'logfile':
547 logfile = kwargs[key]
548 elif key == 'eigenmode_logfile':
549 mlogfile = kwargs[key]
550 else:
551 self.control.set_parameter(key, kwargs[key])
552 self.control.initialize_logfiles(comm=comm, logfile=logfile,
553 eigenmode_logfile=mlogfile)
555 # Seed the randomness
556 if random_seed is None:
557 t = time.time()
558 if world.size > 1:
559 t = world.sum(t) / world.size
560 # Harvest the latter part of the current time
561 random_seed = int(('%30.9f' % t)[-9:])
562 self.random_state = np.random.RandomState(random_seed)
564 # Check the order
565 self.order = self.control.get_parameter('order')
567 # Construct the curvatures list
568 self.curvatures = [100.0] * self.order
570 # Save the original state of the atoms.
571 self.atoms0 = self.atoms.copy()
572 self.save_original_forces()
574 # Get a reference to the log files
575 self.logfile = self.control.get_logfile()
576 self.mlogfile = self.control.get_eigenmode_logfile()
578 def __ase_optimizable__(self):
579 return DimerOptimizable(self)
581 def save_original_forces(self, force_calculation=False):
582 """Store the forces (and energy) of the original state."""
583 # NB: Would be nice if atoms.copy() took care of this.
584 if self.calc is not None:
585 # Hack because some calculators do not have calculation_required
586 if (hasattr(self.calc, 'calculation_required')
587 and not self.calc.calculation_required(self.atoms,
588 ['energy', 'forces'])) or force_calculation:
589 calc = SinglePointCalculator(
590 self.atoms0,
591 energy=self.atoms.get_potential_energy(),
592 forces=self.atoms.get_forces())
593 self.atoms0.calc = calc
595 def initialize_eigenmodes(self, method=None, eigenmodes=None,
596 gauss_std=None):
597 """Make an initial guess for the eigenmode."""
598 if eigenmodes is None:
599 pos = self.get_positions()
600 old_pos = self.get_original_positions()
601 if method is None:
602 method = \
603 self.control.get_parameter('initial_eigenmode_method')
604 if method.lower() == 'displacement' and (pos - old_pos).any():
605 eigenmode = normalize(pos - old_pos)
606 elif method.lower() == 'gauss':
607 self.displace(log=False, gauss_std=gauss_std,
608 method=method)
609 new_pos = self.get_positions()
610 eigenmode = normalize(new_pos - pos)
611 self.set_positions(pos)
612 else:
613 e = 'initial_eigenmode must use either \'gauss\' or ' + \
614 '\'displacement\', if the latter is used the atoms ' + \
615 'must have moved away from the original positions.' + \
616 f'You have requested \'{method}\'.'
617 raise NotImplementedError(e) # NYI
618 eigenmodes = [eigenmode]
620 # Create random higher order mode guesses
621 if self.order > 1:
622 if len(eigenmodes) == 1:
623 for _ in range(1, self.order):
624 pos = self.get_positions()
625 self.displace(log=False, gauss_std=gauss_std,
626 method=method)
627 new_pos = self.get_positions()
628 eigenmode = normalize(new_pos - pos)
629 self.set_positions(pos)
630 eigenmodes += [eigenmode]
632 self.eigenmodes = eigenmodes
633 # Ensure that the higher order mode guesses are all orthogonal
634 if self.order > 1:
635 for k in range(self.order):
636 self.ensure_eigenmode_orthogonality(k)
637 self.eigenmode_log()
639 # NB maybe this name might be confusing in context to
640 # calc.calculation_required()
641 def calculation_required(self):
642 """Check if a calculation is required."""
643 return self.minmode_init or self.check_atoms != self.atoms
645 def calculate_real_forces_and_energies(self, **kwargs):
646 """Calculate and store the potential energy and forces."""
647 if self.minmode_init:
648 self.minmode_init = False
649 self.initialize_eigenmodes(eigenmodes=self.eigenmodes)
650 self.rotation_required = True
651 self.forces0 = self.atoms.get_forces(**kwargs)
652 self.energy0 = self.atoms.get_potential_energy()
653 self.control.increment_counter('forcecalls')
654 self.check_atoms = self.atoms.copy()
656 def get_potential_energy(self):
657 """Return the potential energy."""
658 if self.calculation_required():
659 self.calculate_real_forces_and_energies()
660 return self.energy0
662 def get_forces(self, real=False, pos=None, **kwargs):
663 """Return the forces, projected or real."""
664 if self.calculation_required() and pos is None:
665 self.calculate_real_forces_and_energies(**kwargs)
666 if real and pos is None:
667 return self.forces0
668 elif real and pos is not None:
669 old_pos = self.atoms.get_positions()
670 self.atoms.set_positions(pos)
671 forces = self.atoms.get_forces()
672 self.control.increment_counter('forcecalls')
673 self.atoms.set_positions(old_pos)
674 return forces
675 else:
676 if self.rotation_required:
677 self.find_eigenmodes(order=self.order)
678 self.eigenmode_log()
679 self.rotation_required = False
680 self.control.increment_counter('optcount')
681 return self.get_projected_forces()
683 def ensure_eigenmode_orthogonality(self, order):
684 mode = self.eigenmodes[order - 1].copy()
685 for k in range(order - 1):
686 mode = perpendicular_vector(mode, self.eigenmodes[k])
687 self.eigenmodes[order - 1] = normalize(mode)
689 def find_eigenmodes(self, order=1):
690 """Launch eigenmode searches."""
691 if self.control.get_parameter('eigenmode_method').lower() != 'dimer':
692 e = 'Only the Dimer control object has been implemented.'
693 raise NotImplementedError(e) # NYI
694 for k in range(order):
695 if k > 0:
696 self.ensure_eigenmode_orthogonality(k + 1)
697 search = DimerEigenmodeSearch(
698 self, self.control,
699 eigenmode=self.eigenmodes[k], basis=self.eigenmodes[:k])
700 search.converge_to_eigenmode()
701 search.set_up_for_optimization_step()
702 self.eigenmodes[k] = search.get_eigenmode()
703 self.curvatures[k] = search.get_curvature()
705 def get_projected_forces(self, pos=None):
706 """Return the projected forces."""
707 if pos is not None:
708 forces = self.get_forces(real=True, pos=pos).copy()
709 else:
710 forces = self.forces0.copy()
712 # Loop through all the eigenmodes
713 # NB: Can this be done with a linear combination, instead?
714 for k, mode in enumerate(self.eigenmodes):
715 # NYI This If statement needs to be overridable in the control
716 if self.get_curvature(order=k) > 0.0 and self.order == 1:
717 forces = -parallel_vector(forces, mode)
718 else:
719 forces -= 2 * parallel_vector(forces, mode)
720 return forces
722 def restore_original_positions(self):
723 """Restore the MinModeAtoms object positions to the original state."""
724 self.atoms.set_positions(self.get_original_positions())
726 def get_barrier_energy(self):
727 """The energy difference between the current and original states"""
728 try:
729 original_energy = self.get_original_potential_energy()
730 dimer_energy = self.get_potential_energy()
731 return dimer_energy - original_energy
732 except RuntimeError:
733 w = 'The potential energy is not available, without further ' + \
734 'calculations, most likely at the original state.'
735 warnings.warn(w, UserWarning)
736 return np.nan
738 def get_control(self):
739 """Return the control object."""
740 return self.control
742 def get_curvature(self, order='max'):
743 """Return the eigenvalue estimate."""
744 if order == 'max':
745 return max(self.curvatures)
746 else:
747 return self.curvatures[order - 1]
749 def get_eigenmode(self, order=1):
750 """Return the current eigenmode guess."""
751 return self.eigenmodes[order - 1]
753 def get_atoms(self):
754 """Return the unextended Atoms object."""
755 return self.atoms
757 def set_atoms(self, atoms):
758 """Set a new Atoms object"""
759 self.atoms = atoms
761 def set_eigenmode(self, eigenmode, order=1):
762 """Set the eigenmode guess."""
763 self.eigenmodes[order - 1] = eigenmode
765 def set_curvature(self, curvature, order=1):
766 """Set the eigenvalue estimate."""
767 self.curvatures[order - 1] = curvature
769 # Pipe all the stuff from Atoms that is not overwritten.
770 # Pipe all requests for get_original_* to self.atoms0.
771 def __getattr__(self, attr):
772 """Return any value of the Atoms object"""
773 if 'original' in attr.split('_'):
774 attr = attr.replace('_original_', '_')
775 return getattr(self.atoms0, attr)
776 else:
777 return getattr(self.atoms, attr)
779 def __len__(self):
780 return len(self.atoms)
782 def displace(self, displacement_vector=None, mask=None, method=None,
783 displacement_center=None, radius=None, number_of_atoms=None,
784 gauss_std=None, mic=True, log=True):
785 """Move the atoms away from their current position.
787 This is one of the essential parts of minimum mode searches.
788 The parameters can all be set in the control object and overwritten
789 when this method is run, apart from *displacement_vector*.
790 It is preferred to modify the control values rather than those here
791 in order for the correct ones to show up in the log file.
793 *method* can be either 'gauss' for random displacement or 'vector'
794 to perform a predefined displacement.
796 *gauss_std* is the standard deviation of the gauss curve that is
797 used for random displacement.
799 *displacement_center* can be either the number of an atom or a 3D
800 position. It must be accompanied by a *radius* (all atoms within it
801 will be displaced) or a *number_of_atoms* which decides how many of
802 the closest atoms will be displaced.
804 *mic* controls the usage of the Minimum Image Convention.
806 If both *mask* and *displacement_center* are used, the atoms marked
807 as False in the *mask* will not be affected even though they are
808 within reach of the *displacement_center*.
810 The parameters priority order:
811 1) displacement_vector
812 2) mask
813 3) displacement_center (with radius and/or number_of_atoms)
815 If both *radius* and *number_of_atoms* are supplied with
816 *displacement_center*, only atoms that fulfill both criteria will
817 be displaced.
819 """
821 # Fetch the default values from the control
822 if mask is None:
823 mask = self.control.get_parameter('mask')
824 if method is None:
825 method = self.control.get_parameter('displacement_method')
826 if gauss_std is None:
827 gauss_std = self.control.get_parameter('gauss_std')
828 if displacement_center is None:
829 displacement_center = \
830 self.control.get_parameter('displacement_center')
831 if radius is None:
832 radius = self.control.get_parameter('displacement_radius')
833 if number_of_atoms is None:
834 number_of_atoms = \
835 self.control.get_parameter('number_of_displacement_atoms')
837 # Check for conflicts
838 if displacement_vector is not None and method.lower() != 'vector':
839 e = 'displacement_vector was supplied but a different method ' + \
840 f'(\'{method!s}\') was chosen.\n'
841 raise ValueError(e)
842 elif displacement_vector is None and method.lower() == 'vector':
843 e = 'A displacement_vector must be supplied when using ' + \
844 f'method = \'{method!s}\'.\n'
845 raise ValueError(e)
846 elif displacement_center is not None and radius is None and \
847 number_of_atoms is None:
848 e = 'When displacement_center is chosen, either radius or ' + \
849 'number_of_atoms must be supplied.\n'
850 raise ValueError(e)
852 # Set up the center of displacement mask (c_mask)
853 if displacement_center is not None:
854 c = displacement_center
855 # Construct a distance list
856 # The center is an atom
857 if isinstance(c, int):
858 # Parse negative indexes
859 c = displacement_center % len(self)
860 d = [(k, self.get_distance(k, c, mic=mic)) for k in
861 range(len(self))]
862 # The center is a position in 3D space
863 elif len(c) == 3 and [type(c_k) for c_k in c] == [float] * 3:
864 # NB: MIC is not considered.
865 d = [(k, norm(self.get_positions()[k] - c))
866 for k in range(len(self))]
867 else:
868 e = 'displacement_center must be either the number of an ' + \
869 'atom in MinModeAtoms object or a 3D position ' + \
870 '(3-tuple of floats).'
871 raise ValueError(e)
873 # Set up the mask
874 if radius is not None:
875 r_mask = [dist[1] < radius for dist in d]
876 else:
877 r_mask = [True for _ in range(len(self))]
879 if number_of_atoms is not None:
880 d_sorted = [n[0] for n in sorted(d, key=lambda k: k[1])]
881 n_nearest = d_sorted[:number_of_atoms]
882 n_mask = [k in n_nearest for k in range(len(self))]
883 else:
884 n_mask = [True for _ in range(len(self))]
886 # Resolve n_mask / r_mask conflicts
887 c_mask = [n_mask[k] and r_mask[k] for k in range(len(self))]
888 else:
889 c_mask = None
891 # Set up a True mask if there is no mask supplied
892 if mask is None:
893 mask = [True for _ in range(len(self))]
894 if c_mask is None:
895 w = 'It was not possible to figure out which atoms to ' + \
896 'displace, Will try to displace all atoms.\n'
897 warnings.warn(w, UserWarning)
898 if self.logfile is not None:
899 self.logfile.write('MINMODE:WARN: ' + w + '\n')
900 self.logfile.flush()
902 # Resolve mask / c_mask conflicts
903 if c_mask is not None:
904 mask = [mask[k] and c_mask[k] for k in range(len(self))]
906 if displacement_vector is None:
907 displacement_vector = []
908 for k in range(len(self)):
909 if mask[k]:
910 diff_line = []
911 for _ in range(3):
912 if method.lower() == 'gauss':
913 if not gauss_std:
914 gauss_std = \
915 self.control.get_parameter('gauss_std')
916 diff = self.random_state.normal(0.0, gauss_std)
917 else:
918 e = f'Invalid displacement method >>{method!s}<<'
919 raise ValueError(e)
920 diff_line.append(diff)
921 displacement_vector.append(diff_line)
922 else:
923 displacement_vector.append([0.0] * 3)
925 # Remove displacement of masked atoms
926 for k in range(len(mask)):
927 if not mask[k]:
928 displacement_vector[k] = [0.0] * 3
930 # Perform the displacement and log it
931 if log:
932 pos0 = self.get_positions()
933 self.set_positions(self.get_positions() + displacement_vector)
934 if log:
935 parameters = {'mask': mask,
936 'displacement_method': method,
937 'gauss_std': gauss_std,
938 'displacement_center': displacement_center,
939 'displacement_radius': radius,
940 'number_of_displacement_atoms': number_of_atoms}
941 self.displacement_log(self.get_positions() - pos0, parameters)
943 def eigenmode_log(self):
944 """Log the eigenmodes (eigenmode estimates)"""
945 if self.mlogfile is not None:
946 l = 'MINMODE:MODE: Optimization Step: %i\n' % \
947 (self.control.get_counter('optcount'))
948 for m_num, mode in enumerate(self.eigenmodes):
949 l += 'MINMODE:MODE: Order: %i\n' % m_num
950 for k in range(len(mode)):
951 l += 'MINMODE:MODE: %7i %15.8f %15.8f %15.8f\n' % (
952 k, mode[k][0], mode[k][1], mode[k][2])
953 self.mlogfile.write(l)
954 self.mlogfile.flush()
956 def displacement_log(self, displacement_vector, parameters):
957 """Log the displacement"""
958 if self.logfile is not None:
959 lp = 'MINMODE:DISP: Parameters, different from the control:\n'
960 mod_para = False
961 for key in parameters:
962 if parameters[key] != self.control.get_parameter(key):
963 lp += 'MINMODE:DISP: {} = {}\n'.format(str(key),
964 str(parameters[key]))
965 mod_para = True
966 if mod_para:
967 l = lp
968 else:
969 l = ''
970 for k in range(len(displacement_vector)):
971 l += 'MINMODE:DISP: %7i %15.8f %15.8f %15.8f\n' % (
972 k,
973 displacement_vector[k][0], displacement_vector[k][1],
974 displacement_vector[k][2])
975 self.logfile.write(l)
976 self.logfile.flush()
978 def summarize(self):
979 """Summarize the Minimum mode search."""
980 if self.logfile is None:
981 logfile = sys.stdout
982 else:
983 logfile = self.logfile
985 c = self.control
986 label = 'MINMODE:SUMMARY: '
988 l = label + '-------------------------\n'
989 l += label + 'Barrier: %16.4f\n' % self.get_barrier_energy()
990 l += label + 'Curvature: %14.4f\n' % self.get_curvature()
991 l += label + 'Optimizer steps: %8i\n' % c.get_counter('optcount')
992 l += label + 'Forcecalls: %13i\n' % c.get_counter('forcecalls')
993 l += label + '-------------------------\n'
995 logfile.write(l)
998class MinModeTranslate(Optimizer):
999 """An Optimizer specifically tailored to minimum mode following."""
1001 def __init__(self, dimeratoms, logfile='-', trajectory=None):
1002 Optimizer.__init__(self, dimeratoms, None, logfile, trajectory)
1004 self.control = dimeratoms.get_control()
1005 self.dimeratoms = dimeratoms
1007 # Make a header for the log
1008 if self.logfile is not None:
1009 l = ''
1010 if isinstance(self.control, DimerControl):
1011 l = 'MinModeTranslate: STEP TIME ENERGY ' + \
1012 'MAX-FORCE STEPSIZE CURVATURE ROT-STEPS\n'
1013 self.logfile.write(l)
1014 self.logfile.flush()
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."""
1032 atoms = self.dimeratoms
1033 if f is None:
1034 f = atoms.get_forces()
1035 r = atoms.get_positions()
1036 curv = atoms.get_curvature()
1037 f0p = f.copy()
1038 r0 = r.copy()
1039 direction = f0p.copy()
1040 if self.cg_on:
1041 direction = self.get_cg_direction(direction)
1042 direction = normalize(direction)
1043 if curv > 0.0:
1044 step = direction * self.max_step
1045 else:
1046 r0t = r0 + direction * self.trial_step
1047 f0tp = self.dimeratoms.get_projected_forces(r0t)
1048 F = np.vdot((f0tp + f0p), direction) / 2.0
1049 C = np.vdot((f0tp - f0p), direction) / self.trial_step
1050 step = (-F / C + self.trial_step / 2.0) * direction
1051 if norm(step) > self.max_step:
1052 step = direction * self.max_step
1053 self.log(f0p, norm(step))
1055 atoms.set_positions(r + step)
1057 self.f0 = f.flat.copy()
1058 self.r0 = r.flat.copy()
1060 def get_cg_direction(self, direction):
1061 """Apply the Conjugate Gradient algorithm to the step direction."""
1062 if self.cg_init:
1063 self.cg_init = False
1064 self.direction_old = direction.copy()
1065 self.cg_direction = direction.copy()
1066 old_norm = np.vdot(self.direction_old, self.direction_old)
1067 # Polak-Ribiere Conjugate Gradient
1068 if old_norm != 0.0:
1069 betaPR = np.vdot(direction, (direction - self.direction_old)) / \
1070 old_norm
1071 else:
1072 betaPR = 0.0
1073 if betaPR < 0.0:
1074 betaPR = 0.0
1075 self.cg_direction = direction + self.cg_direction * betaPR
1076 self.direction_old = direction.copy()
1077 return self.cg_direction.copy()
1079 def log(self, f=None, stepsize=None):
1080 """Log each step of the optimization."""
1081 if f is None:
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)
1105 self.logfile.flush()
1108def read_eigenmode(mlog, index=-1):
1109 """Read an eigenmode.
1110 To access the pre optimization eigenmode set index = 'null'.
1112 """
1113 mlog_is_str = isinstance(mlog, str)
1114 if mlog_is_str:
1115 fd = open(mlog)
1116 else:
1117 fd = mlog
1119 lines = fd.readlines()
1121 # Detect the amount of atoms and iterations
1122 k = 2
1123 while lines[k].split()[1].lower() not in ['optimization', 'order']:
1124 k += 1
1125 n = k - 2
1126 n_itr = (len(lines) // (n + 1)) - 2
1128 # Locate the correct image.
1129 if isinstance(index, str):
1130 if index.lower() == 'null':
1131 i = 0
1132 else:
1133 i = int(index) + 1
1134 else:
1135 if index >= 0:
1136 i = index + 1
1137 else:
1138 if index < -n_itr - 1:
1139 raise IndexError('list index out of range')
1140 else:
1141 i = index
1143 mode = np.ndarray(shape=(n, 3), dtype=float)
1144 for k_atom, k in enumerate(range(1, n + 1)):
1145 line = lines[i * (n + 1) + k].split()
1146 for k_dim in range(3):
1147 mode[k_atom][k_dim] = float(line[k_dim + 2])
1148 if mlog_is_str:
1149 fd.close()
1151 return mode
1154# Aliases
1155DimerAtoms = MinModeAtoms
1156DimerTranslate = MinModeTranslate