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