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