Source code for ase.mep.dimer

"""Minimum mode follower for finding saddle points in an unbiased way.

There is, currently, only one implemented method: The Dimer method.

import sys
import time
import warnings
from math import atan, cos, degrees, pi, sin, sqrt, tan
from typing import Any, Dict

import numpy as np

from ase.calculators.singlepoint import SinglePointCalculator
from ase.optimize.optimize import OptimizableAtoms, Optimizer
from ase.parallel import world
from ase.utils import IOContext

# Handy vector methods
norm = np.linalg.norm

class DimerOptimizable(OptimizableAtoms):
    def __init__(self, dimeratoms):
        self.dimeratoms = dimeratoms

    def converged(self, forces, fmax):
        forces_converged = super().converged(forces, fmax)
        return forces_converged and self.dimeratoms.get_curvature() < 0.0

def normalize(vector):
    """Create a unit vector along *vector*"""
    return vector / norm(vector)

def parallel_vector(vector, base):
    """Extract the components of *vector* that are parallel to *base*"""
    return np.vdot(vector, base) * base

def perpendicular_vector(vector, base):
    """Remove the components of *vector* that are parallel to *base*"""
    return vector - parallel_vector(vector, base)

def rotate_vectors(v1i, v2i, angle):
    """Rotate vectors *v1i* and *v2i* by *angle*"""
    cAng = cos(angle)
    sAng = sin(angle)
    v1o = v1i * cAng + v2i * sAng
    v2o = v2i * cAng - v1i * sAng
    # Ensure the length of the input and output vectors is equal
    return normalize(v1o) * norm(v1i), normalize(v2o) * norm(v2i)

[docs]class DimerEigenmodeSearch: """An implementation of the Dimer's minimum eigenvalue mode search. This class implements the rotational part of the dimer saddle point searching method. Parameters: atoms: MinModeAtoms object MinModeAtoms is an extension to the Atoms object, which includes information about the lowest eigenvalue mode. control: DimerControl object Contains the parameters necessary for the eigenmode search. If no control object is supplied a default DimerControl will be created and used. basis: list of xyz-values Eigenmode. Must be an ndarray of shape (n, 3). It is possible to constrain the eigenmodes to be orthogonal to this given eigenmode. Notes: The code is inspired, with permission, by code written by the Henkelman group, which can be found at References: * Henkelman and Jonsson, JCP 111, 7010 (1999) * Olsen, Kroes, Henkelman, Arnaldsson, and Jonsson, JCP 121, 9776 (2004). * Heyden, Bell, and Keil, JCP 123, 224101 (2005). * Kastner and Sherwood, JCP 128, 014106 (2008). """ def __init__(self, dimeratoms, control=None, eigenmode=None, basis=None, **kwargs): if hasattr(dimeratoms, 'get_eigenmode'): self.dimeratoms = dimeratoms else: e = 'The atoms object must be a MinModeAtoms object' raise TypeError(e) self.basis = basis if eigenmode is None: self.eigenmode = self.dimeratoms.get_eigenmode() else: self.eigenmode = eigenmode if control is None: self.control = DimerControl(**kwargs) w = 'Missing control object in ' + self.__class__.__name__ + \ '. Using default: DimerControl()' warnings.warn(w, UserWarning) if self.control.logfile is not None: self.control.logfile.write('DIM:WARN: ' + w + '\n') self.control.logfile.flush() else: self.control = control # kwargs must be empty if a control object is supplied for key in kwargs: e = f'__init__() got an unexpected keyword argument \'{(key)}\'' raise TypeError(e) self.dR = self.control.get_parameter('dimer_separation') self.logfile = self.control.get_logfile() def converge_to_eigenmode(self): """Perform an eigenmode search.""" self.set_up_for_eigenmode_search() stoprot = False # Load the relevant parameters from control f_rot_min = self.control.get_parameter('f_rot_min') f_rot_max = self.control.get_parameter('f_rot_max') trial_angle = self.control.get_parameter('trial_angle') max_num_rot = self.control.get_parameter('max_num_rot') extrapolate = self.control.get_parameter('extrapolate_forces') while not stoprot: if self.forces1E is None: self.update_virtual_forces() else: self.update_virtual_forces(extrapolated_forces=True) self.forces1A = self.forces1 self.update_curvature() f_rot_A = self.get_rotational_force() # Pre rotation stop criteria if norm(f_rot_A) <= f_rot_min: self.log(f_rot_A, None) stoprot = True else: n_A = self.eigenmode rot_unit_A = normalize(f_rot_A) # Get the curvature and its derivative c0 = self.get_curvature() c0d = np.vdot((self.forces2 - self.forces1), rot_unit_A) / \ self.dR # Trial rotation (no need to store the curvature) # NYI variable trial angles from [3] n_B, rot_unit_B = rotate_vectors(n_A, rot_unit_A, trial_angle) self.eigenmode = n_B self.update_virtual_forces() self.forces1B = self.forces1 # Get the curvature's derivative c1d = np.vdot((self.forces2 - self.forces1), rot_unit_B) / \ self.dR # Calculate the Fourier coefficients a1 = c0d * cos(2 * trial_angle) - c1d / \ (2 * sin(2 * trial_angle)) b1 = 0.5 * c0d a0 = 2 * (c0 - a1) # Estimate the rotational angle rotangle = atan(b1 / a1) / 2.0 # Make sure that you didn't find a maximum cmin = a0 / 2.0 + a1 * cos(2 * rotangle) + \ b1 * sin(2 * rotangle) if c0 < cmin: rotangle += pi / 2.0 # Rotate into the (hopefully) lowest eigenmode # NYI Conjugate gradient rotation n_min, dummy = rotate_vectors(n_A, rot_unit_A, rotangle) self.update_eigenmode(n_min) # Store the curvature estimate instead of the old curvature self.update_curvature(cmin) self.log(f_rot_A, rotangle) # Force extrapolation scheme from [4] if extrapolate: self.forces1E = sin(trial_angle - rotangle) / \ sin(trial_angle) * self.forces1A + sin(rotangle) / \ sin(trial_angle) * self.forces1B + \ (1 - cos(rotangle) - sin(rotangle) * tan(trial_angle / 2.0)) * self.forces0 else: self.forces1E = None # Post rotation stop criteria if not stoprot: if self.control.get_counter('rotcount') >= max_num_rot: stoprot = True elif norm(f_rot_A) <= f_rot_max: stoprot = True def log(self, f_rot_A, angle): """Log each rotational step.""" # NYI Log for the trial angle if self.logfile is not None: if angle: l = 'DIM:ROT: %7d %9d %9.4f %9.4f %9.4f\n' % \ (self.control.get_counter('optcount'), self.control.get_counter('rotcount'), self.get_curvature(), degrees(angle), norm(f_rot_A)) else: l = 'DIM:ROT: %7d %9d %9.4f %9s %9.4f\n' % \ (self.control.get_counter('optcount'), self.control.get_counter('rotcount'), self.get_curvature(), '---------', norm(f_rot_A)) self.logfile.write(l) self.logfile.flush() def get_rotational_force(self): """Calculate the rotational force that acts on the dimer.""" rot_force = perpendicular_vector((self.forces1 - self.forces2), self.eigenmode) / (2.0 * self.dR) if self.basis is not None: if ( len(self.basis) == len(self.dimeratoms) and len(self.basis[0]) == 3 and isinstance(self.basis[0][0], float)): rot_force = perpendicular_vector(rot_force, self.basis) else: for base in self.basis: rot_force = perpendicular_vector(rot_force, base) return rot_force def update_curvature(self, curv=None): """Update the curvature in the MinModeAtoms object.""" if curv: self.curvature = curv else: self.curvature = np.vdot((self.forces2 - self.forces1), self.eigenmode) / (2.0 * self.dR) def update_eigenmode(self, eigenmode): """Update the eigenmode in the MinModeAtoms object.""" self.eigenmode = eigenmode self.update_virtual_positions() self.control.increment_counter('rotcount') def get_eigenmode(self): """Returns the current eigenmode.""" return self.eigenmode def get_curvature(self): """Returns the curvature along the current eigenmode.""" return self.curvature def get_control(self): """Return the control object.""" return self.control def update_center_forces(self): """Get the forces at the center of the dimer.""" self.dimeratoms.set_positions(self.pos0) self.forces0 = self.dimeratoms.get_forces(real=True) self.energy0 = self.dimeratoms.get_potential_energy() def update_virtual_forces(self, extrapolated_forces=False): """Get the forces at the endpoints of the dimer.""" self.update_virtual_positions() # Estimate / Calculate the forces at pos1 if extrapolated_forces: self.forces1 = self.forces1E.copy() else: self.forces1 = self.dimeratoms.get_forces(real=True, pos=self.pos1) # Estimate / Calculate the forces at pos2 if self.control.get_parameter('use_central_forces'): self.forces2 = 2 * self.forces0 - self.forces1 else: self.forces2 = self.dimeratoms.get_forces(real=True, pos=self.pos2) def update_virtual_positions(self): """Update the end point positions.""" self.pos1 = self.pos0 + self.eigenmode * self.dR self.pos2 = self.pos0 - self.eigenmode * self.dR def set_up_for_eigenmode_search(self): """Before eigenmode search, prepare for rotation.""" self.pos0 = self.dimeratoms.get_positions() self.update_center_forces() self.update_virtual_positions() self.control.reset_counter('rotcount') self.forces1E = None def set_up_for_optimization_step(self): """At the end of rotation, prepare for displacement of the dimer.""" self.dimeratoms.set_positions(self.pos0) self.forces1E = None
class MinModeControl(IOContext): """A parent class for controlling minimum mode saddle point searches. Method specific control classes inherit this one. The only thing inheriting classes need to implement are the log() method and the *parameters* class variable with default values for ALL parameters needed by the method in question. When instantiating control classes default parameter values can be overwritten. """ parameters: Dict[str, Any] = {} def __init__(self, logfile='-', eigenmode_logfile=None, comm=world, **kwargs): # Overwrite the defaults with the input parameters given for key in kwargs: if key not in self.parameters: e = (f'Invalid parameter >>{key}<< with value >>' f'{kwargs[key]!s}<< in {self.__class__.__name__}') raise ValueError(e) else: self.set_parameter(key, kwargs[key], log=False) self.initialize_logfiles(comm=comm, logfile=logfile, eigenmode_logfile=eigenmode_logfile) self.counters = {'forcecalls': 0, 'rotcount': 0, 'optcount': 0} self.log() def initialize_logfiles(self, comm, logfile=None, eigenmode_logfile=None): self.logfile = self.openfile(file=logfile, comm=comm) self.eigenmode_logfile = self.openfile(file=eigenmode_logfile, comm=comm) def log(self, parameter=None): """Log the parameters of the eigenmode search.""" def set_parameter(self, parameter, value, log=True): """Change a parameter's value.""" if parameter not in self.parameters: e = f'Invalid parameter >>{parameter}<< with value >>{value!s}<<' raise ValueError(e) self.parameters[parameter] = value if log: self.log(parameter) def get_parameter(self, parameter): """Returns the value of a parameter.""" if parameter not in self.parameters: e = f'Invalid parameter >>{(parameter)}<<' raise ValueError(e) return self.parameters[parameter] def get_logfile(self): """Returns the log file.""" return self.logfile def get_eigenmode_logfile(self): """Returns the eigenmode log file.""" return self.eigenmode_logfile def get_counter(self, counter): """Returns a given counter.""" return self.counters[counter] def increment_counter(self, counter): """Increment a given counter.""" self.counters[counter] += 1 def reset_counter(self, counter): """Reset a given counter.""" self.counters[counter] = 0 def reset_all_counters(self): """Reset all counters.""" for key in self.counters: self.counters[key] = 0
[docs]class DimerControl(MinModeControl): """A class that takes care of the parameters needed for a Dimer search. Parameters: eigenmode_method: str The name of the eigenmode search method. f_rot_min: float Size of the rotational force under which no rotation will be performed. f_rot_max: float Size of the rotational force under which only one rotation will be performed. max_num_rot: int Maximum number of rotations per optimizer step. trial_angle: float Trial angle for the finite difference estimate of the rotational angle in radians. trial_trans_step: float Trial step size for the MinModeTranslate optimizer. maximum_translation: float Maximum step size and forced step size when the curvature is still positive for the MinModeTranslate optimizer. cg_translation: bool Conjugate Gradient for the MinModeTranslate optimizer. use_central_forces: bool Only calculate the forces at one end of the dimer and extrapolate the forces to the other. dimer_separation: float Separation of the dimer's images. initial_eigenmode_method: str How to construct the initial eigenmode of the dimer. If an eigenmode is given when creating the MinModeAtoms object, this will be ignored. Possible choices are: 'gauss' and 'displacement' extrapolate_forces: bool When more than one rotation is performed, an extrapolation scheme can be used to reduce the number of force evaluations. displacement_method: str How to displace the atoms. Possible choices are 'gauss' and 'vector'. gauss_std: float The standard deviation of the gauss curve used when doing random displacement. order: int How many lowest eigenmodes will be inverted. mask: list of bool Which atoms will be moved during displacement. displacement_center: int or [float, float, float] The center of displacement, nearby atoms will be displaced. displacement_radius: float When choosing which atoms to displace with the *displacement_center* keyword, this decides how many nearby atoms to displace. number_of_displacement_atoms: int The amount of atoms near *displacement_center* to displace. """ # Default parameters for the Dimer eigenmode search parameters = {'eigenmode_method': 'dimer', 'f_rot_min': 0.1, 'f_rot_max': 1.00, 'max_num_rot': 1, 'trial_angle': pi / 4.0, 'trial_trans_step': 0.001, 'maximum_translation': 0.1, 'cg_translation': True, 'use_central_forces': True, 'dimer_separation': 0.0001, 'initial_eigenmode_method': 'gauss', 'extrapolate_forces': False, 'displacement_method': 'gauss', 'gauss_std': 0.1, 'order': 1, 'mask': None, # NB mask should not be a "parameter" 'displacement_center': None, 'displacement_radius': None, 'number_of_displacement_atoms': None} # NB: Can maybe put this in EigenmodeSearch and MinModeControl def log(self, parameter=None): """Log the parameters of the eigenmode search.""" if self.logfile is not None: if parameter is not None: l = 'DIM:CONTROL: Updated Parameter: {} = {}\n'.format( parameter, str(self.get_parameter(parameter))) else: l = 'MINMODE:METHOD: Dimer\n' l += 'DIM:CONTROL: Search Parameters:\n' l += 'DIM:CONTROL: ------------------\n' for key in self.parameters: l += 'DIM:CONTROL: {} = {}\n'.format( key, str(self.get_parameter(key))) l += 'DIM:CONTROL: ------------------\n' l += 'DIM:ROT: OPT-STEP ROT-STEP CURVATURE ROT-ANGLE ' + \ 'ROT-FORCE\n' self.logfile.write(l) self.logfile.flush()
[docs]class MinModeAtoms: """Wrapper for Atoms with information related to minimum mode searching. Contains an Atoms object and pipes all unknown function calls to that object. Other information that is stored in this object are the estimate for the lowest eigenvalue, *curvature*, and its corresponding eigenmode, *eigenmode*. Furthermore, the original configuration of the Atoms object is stored for use in multiple minimum mode searches. The forces on the system are modified by inverting the component along the eigenmode estimate. This eventually brings the system to a saddle point. Parameters: atoms : Atoms object A regular Atoms object control : MinModeControl object Contains the parameters necessary for the eigenmode search. If no control object is supplied a default DimerControl will be created and used. mask: list of bool Determines which atoms will be moved when calling displace() random_seed: int The seed used for the random number generator. Defaults to modified version the current time. References: [1]_ [2]_ [3]_ [4]_ .. [1] Henkelman and Jonsson, JCP 111, 7010 (1999) .. [2] Olsen, Kroes, Henkelman, Arnaldsson, and Jonsson, JCP 121, 9776 (2004). .. [3] Heyden, Bell, and Keil, JCP 123, 224101 (2005). .. [4] Kastner and Sherwood, JCP 128, 014106 (2008). """ def __init__(self, atoms, control=None, eigenmodes=None, random_seed=None, comm=world, **kwargs): self.minmode_init = True self.atoms = atoms # Initialize to None to avoid strange behaviour due to __getattr__ self.eigenmodes = eigenmodes self.curvatures = None if control is None: self.control = DimerControl(**kwargs) w = 'Missing control object in ' + self.__class__.__name__ + \ '. Using default: DimerControl()' warnings.warn(w, UserWarning) if self.control.logfile is not None: self.control.logfile.write('DIM:WARN: ' + w + '\n') self.control.logfile.flush() else: self.control = control logfile = self.control.get_logfile() mlogfile = self.control.get_eigenmode_logfile() for key in kwargs: if key == 'logfile': logfile = kwargs[key] elif key == 'eigenmode_logfile': mlogfile = kwargs[key] else: self.control.set_parameter(key, kwargs[key]) self.control.initialize_logfiles(comm=comm, logfile=logfile, eigenmode_logfile=mlogfile) # Seed the randomness if random_seed is None: t = time.time() if world.size > 1: t = world.sum(t) / world.size # Harvest the latter part of the current time random_seed = int(('%30.9f' % t)[-9:]) self.random_state = np.random.RandomState(random_seed) # Check the order self.order = self.control.get_parameter('order') # Construct the curvatures list self.curvatures = [100.0] * self.order # Save the original state of the atoms. self.atoms0 = self.atoms.copy() self.save_original_forces() # Get a reference to the log files self.logfile = self.control.get_logfile() self.mlogfile = self.control.get_eigenmode_logfile() def __ase_optimizable__(self): return DimerOptimizable(self) def save_original_forces(self, force_calculation=False): """Store the forces (and energy) of the original state.""" # NB: Would be nice if atoms.copy() took care of this. if self.calc is not None: # Hack because some calculators do not have calculation_required if (hasattr(self.calc, 'calculation_required') and not self.calc.calculation_required(self.atoms, ['energy', 'forces'])) or force_calculation: calc = SinglePointCalculator( self.atoms0, energy=self.atoms.get_potential_energy(), forces=self.atoms.get_forces()) self.atoms0.calc = calc def initialize_eigenmodes(self, method=None, eigenmodes=None, gauss_std=None): """Make an initial guess for the eigenmode.""" if eigenmodes is None: pos = self.get_positions() old_pos = self.get_original_positions() if method is None: method = \ self.control.get_parameter('initial_eigenmode_method') if method.lower() == 'displacement' and (pos - old_pos).any(): eigenmode = normalize(pos - old_pos) elif method.lower() == 'gauss': self.displace(log=False, gauss_std=gauss_std, method=method) new_pos = self.get_positions() eigenmode = normalize(new_pos - pos) self.set_positions(pos) else: e = 'initial_eigenmode must use either \'gauss\' or ' + \ '\'displacement\', if the latter is used the atoms ' + \ 'must have moved away from the original positions.' + \ f'You have requested \'{method}\'.' raise NotImplementedError(e) # NYI eigenmodes = [eigenmode] # Create random higher order mode guesses if self.order > 1: if len(eigenmodes) == 1: for _ in range(1, self.order): pos = self.get_positions() self.displace(log=False, gauss_std=gauss_std, method=method) new_pos = self.get_positions() eigenmode = normalize(new_pos - pos) self.set_positions(pos) eigenmodes += [eigenmode] self.eigenmodes = eigenmodes # Ensure that the higher order mode guesses are all orthogonal if self.order > 1: for k in range(self.order): self.ensure_eigenmode_orthogonality(k) self.eigenmode_log() # NB maybe this name might be confusing in context to # calc.calculation_required() def calculation_required(self): """Check if a calculation is required.""" return self.minmode_init or self.check_atoms != self.atoms def calculate_real_forces_and_energies(self, **kwargs): """Calculate and store the potential energy and forces.""" if self.minmode_init: self.minmode_init = False self.initialize_eigenmodes(eigenmodes=self.eigenmodes) self.rotation_required = True self.forces0 = self.atoms.get_forces(**kwargs) self.energy0 = self.atoms.get_potential_energy() self.control.increment_counter('forcecalls') self.check_atoms = self.atoms.copy() def get_potential_energy(self): """Return the potential energy.""" if self.calculation_required(): self.calculate_real_forces_and_energies() return self.energy0 def get_forces(self, real=False, pos=None, **kwargs): """Return the forces, projected or real.""" if self.calculation_required() and pos is None: self.calculate_real_forces_and_energies(**kwargs) if real and pos is None: return self.forces0 elif real and pos is not None: old_pos = self.atoms.get_positions() self.atoms.set_positions(pos) forces = self.atoms.get_forces() self.control.increment_counter('forcecalls') self.atoms.set_positions(old_pos) return forces else: if self.rotation_required: self.find_eigenmodes(order=self.order) self.eigenmode_log() self.rotation_required = False self.control.increment_counter('optcount') return self.get_projected_forces() def ensure_eigenmode_orthogonality(self, order): mode = self.eigenmodes[order - 1].copy() for k in range(order - 1): mode = perpendicular_vector(mode, self.eigenmodes[k]) self.eigenmodes[order - 1] = normalize(mode) def find_eigenmodes(self, order=1): """Launch eigenmode searches.""" if self.control.get_parameter('eigenmode_method').lower() != 'dimer': e = 'Only the Dimer control object has been implemented.' raise NotImplementedError(e) # NYI for k in range(order): if k > 0: self.ensure_eigenmode_orthogonality(k + 1) search = DimerEigenmodeSearch( self, self.control, eigenmode=self.eigenmodes[k], basis=self.eigenmodes[:k]) search.converge_to_eigenmode() search.set_up_for_optimization_step() self.eigenmodes[k] = search.get_eigenmode() self.curvatures[k] = search.get_curvature() def get_projected_forces(self, pos=None): """Return the projected forces.""" if pos is not None: forces = self.get_forces(real=True, pos=pos).copy() else: forces = self.forces0.copy() # Loop through all the eigenmodes # NB: Can this be done with a linear combination, instead? for k, mode in enumerate(self.eigenmodes): # NYI This If statement needs to be overridable in the control if self.get_curvature(order=k) > 0.0 and self.order == 1: forces = -parallel_vector(forces, mode) else: forces -= 2 * parallel_vector(forces, mode) return forces def restore_original_positions(self): """Restore the MinModeAtoms object positions to the original state.""" self.atoms.set_positions(self.get_original_positions()) def get_barrier_energy(self): """The energy difference between the current and original states""" try: original_energy = self.get_original_potential_energy() dimer_energy = self.get_potential_energy() return dimer_energy - original_energy except RuntimeError: w = 'The potential energy is not available, without further ' + \ 'calculations, most likely at the original state.' warnings.warn(w, UserWarning) return np.nan def get_control(self): """Return the control object.""" return self.control def get_curvature(self, order='max'): """Return the eigenvalue estimate.""" if order == 'max': return max(self.curvatures) else: return self.curvatures[order - 1] def get_eigenmode(self, order=1): """Return the current eigenmode guess.""" return self.eigenmodes[order - 1] def get_atoms(self): """Return the unextended Atoms object.""" return self.atoms def set_atoms(self, atoms): """Set a new Atoms object""" self.atoms = atoms def set_eigenmode(self, eigenmode, order=1): """Set the eigenmode guess.""" self.eigenmodes[order - 1] = eigenmode def set_curvature(self, curvature, order=1): """Set the eigenvalue estimate.""" self.curvatures[order - 1] = curvature # Pipe all the stuff from Atoms that is not overwritten. # Pipe all requests for get_original_* to self.atoms0. def __getattr__(self, attr): """Return any value of the Atoms object""" if 'original' in attr.split('_'): attr = attr.replace('_original_', '_') return getattr(self.atoms0, attr) else: return getattr(self.atoms, attr) def __len__(self): return len(self.atoms) def displace(self, displacement_vector=None, mask=None, method=None, displacement_center=None, radius=None, number_of_atoms=None, gauss_std=None, mic=True, log=True): """Move the atoms away from their current position. This is one of the essential parts of minimum mode searches. The parameters can all be set in the control object and overwritten when this method is run, apart from *displacement_vector*. It is preferred to modify the control values rather than those here in order for the correct ones to show up in the log file. *method* can be either 'gauss' for random displacement or 'vector' to perform a predefined displacement. *gauss_std* is the standard deviation of the gauss curve that is used for random displacement. *displacement_center* can be either the number of an atom or a 3D position. It must be accompanied by a *radius* (all atoms within it will be displaced) or a *number_of_atoms* which decides how many of the closest atoms will be displaced. *mic* controls the usage of the Minimum Image Convention. If both *mask* and *displacement_center* are used, the atoms marked as False in the *mask* will not be affected even though they are within reach of the *displacement_center*. The parameters priority order: 1) displacement_vector 2) mask 3) displacement_center (with radius and/or number_of_atoms) If both *radius* and *number_of_atoms* are supplied with *displacement_center*, only atoms that fulfill both criteria will be displaced. """ # Fetch the default values from the control if mask is None: mask = self.control.get_parameter('mask') if method is None: method = self.control.get_parameter('displacement_method') if gauss_std is None: gauss_std = self.control.get_parameter('gauss_std') if displacement_center is None: displacement_center = \ self.control.get_parameter('displacement_center') if radius is None: radius = self.control.get_parameter('displacement_radius') if number_of_atoms is None: number_of_atoms = \ self.control.get_parameter('number_of_displacement_atoms') # Check for conflicts if displacement_vector is not None and method.lower() != 'vector': e = 'displacement_vector was supplied but a different method ' + \ f'(\'{method!s}\') was chosen.\n' raise ValueError(e) elif displacement_vector is None and method.lower() == 'vector': e = 'A displacement_vector must be supplied when using ' + \ f'method = \'{method!s}\'.\n' raise ValueError(e) elif displacement_center is not None and radius is None and \ number_of_atoms is None: e = 'When displacement_center is chosen, either radius or ' + \ 'number_of_atoms must be supplied.\n' raise ValueError(e) # Set up the center of displacement mask (c_mask) if displacement_center is not None: c = displacement_center # Construct a distance list # The center is an atom if isinstance(c, int): # Parse negative indexes c = displacement_center % len(self) d = [(k, self.get_distance(k, c, mic=mic)) for k in range(len(self))] # The center is a position in 3D space elif len(c) == 3 and [type(c_k) for c_k in c] == [float] * 3: # NB: MIC is not considered. d = [(k, norm(self.get_positions()[k] - c)) for k in range(len(self))] else: e = 'displacement_center must be either the number of an ' + \ 'atom in MinModeAtoms object or a 3D position ' + \ '(3-tuple of floats).' raise ValueError(e) # Set up the mask if radius is not None: r_mask = [dist[1] < radius for dist in d] else: r_mask = [True for _ in range(len(self))] if number_of_atoms is not None: d_sorted = [n[0] for n in sorted(d, key=lambda k: k[1])] n_nearest = d_sorted[:number_of_atoms] n_mask = [k in n_nearest for k in range(len(self))] else: n_mask = [True for _ in range(len(self))] # Resolve n_mask / r_mask conflicts c_mask = [n_mask[k] and r_mask[k] for k in range(len(self))] else: c_mask = None # Set up a True mask if there is no mask supplied if mask is None: mask = [True for _ in range(len(self))] if c_mask is None: w = 'It was not possible to figure out which atoms to ' + \ 'displace, Will try to displace all atoms.\n' warnings.warn(w, UserWarning) if self.logfile is not None: self.logfile.write('MINMODE:WARN: ' + w + '\n') self.logfile.flush() # Resolve mask / c_mask conflicts if c_mask is not None: mask = [mask[k] and c_mask[k] for k in range(len(self))] if displacement_vector is None: displacement_vector = [] for k in range(len(self)): if mask[k]: diff_line = [] for _ in range(3): if method.lower() == 'gauss': if not gauss_std: gauss_std = \ self.control.get_parameter('gauss_std') diff = self.random_state.normal(0.0, gauss_std) else: e = f'Invalid displacement method >>{method!s}<<' raise ValueError(e) diff_line.append(diff) displacement_vector.append(diff_line) else: displacement_vector.append([0.0] * 3) # Remove displacement of masked atoms for k in range(len(mask)): if not mask[k]: displacement_vector[k] = [0.0] * 3 # Perform the displacement and log it if log: pos0 = self.get_positions() self.set_positions(self.get_positions() + displacement_vector) if log: parameters = {'mask': mask, 'displacement_method': method, 'gauss_std': gauss_std, 'displacement_center': displacement_center, 'displacement_radius': radius, 'number_of_displacement_atoms': number_of_atoms} self.displacement_log(self.get_positions() - pos0, parameters) def eigenmode_log(self): """Log the eigenmodes (eigenmode estimates)""" if self.mlogfile is not None: l = 'MINMODE:MODE: Optimization Step: %i\n' % \ (self.control.get_counter('optcount')) for m_num, mode in enumerate(self.eigenmodes): l += 'MINMODE:MODE: Order: %i\n' % m_num for k in range(len(mode)): l += 'MINMODE:MODE: %7i %15.8f %15.8f %15.8f\n' % ( k, mode[k][0], mode[k][1], mode[k][2]) self.mlogfile.write(l) self.mlogfile.flush() def displacement_log(self, displacement_vector, parameters): """Log the displacement""" if self.logfile is not None: lp = 'MINMODE:DISP: Parameters, different from the control:\n' mod_para = False for key in parameters: if parameters[key] != self.control.get_parameter(key): lp += 'MINMODE:DISP: {} = {}\n'.format(str(key), str(parameters[key])) mod_para = True if mod_para: l = lp else: l = '' for k in range(len(displacement_vector)): l += 'MINMODE:DISP: %7i %15.8f %15.8f %15.8f\n' % ( k, displacement_vector[k][0], displacement_vector[k][1], displacement_vector[k][2]) self.logfile.write(l) self.logfile.flush() def summarize(self): """Summarize the Minimum mode search.""" if self.logfile is None: logfile = sys.stdout else: logfile = self.logfile c = self.control label = 'MINMODE:SUMMARY: ' l = label + '-------------------------\n' l += label + 'Barrier: %16.4f\n' % self.get_barrier_energy() l += label + 'Curvature: %14.4f\n' % self.get_curvature() l += label + 'Optimizer steps: %8i\n' % c.get_counter('optcount') l += label + 'Forcecalls: %13i\n' % c.get_counter('forcecalls') l += label + '-------------------------\n' logfile.write(l)
[docs]class MinModeTranslate(Optimizer): """An Optimizer specifically tailored to minimum mode following.""" def __init__(self, dimeratoms, logfile='-', trajectory=None): Optimizer.__init__(self, dimeratoms, None, logfile, trajectory) self.control = dimeratoms.get_control() self.dimeratoms = dimeratoms # Make a header for the log if self.logfile is not None: l = '' if isinstance(self.control, DimerControl): l = 'MinModeTranslate: STEP TIME ENERGY ' + \ 'MAX-FORCE STEPSIZE CURVATURE ROT-STEPS\n' self.logfile.write(l) self.logfile.flush() # Load the relevant parameters from control self.cg_on = self.control.get_parameter('cg_translation') self.trial_step = self.control.get_parameter('trial_trans_step') self.max_step = self.control.get_parameter('maximum_translation') # Start conjugate gradient if self.cg_on: self.cg_init = True def initialize(self): """Set initial values.""" self.r0 = None self.f0 = None def step(self, f=None): """Perform the optimization step.""" atoms = self.dimeratoms if f is None: f = atoms.get_forces() r = atoms.get_positions() curv = atoms.get_curvature() f0p = f.copy() r0 = r.copy() direction = f0p.copy() if self.cg_on: direction = self.get_cg_direction(direction) direction = normalize(direction) if curv > 0.0: step = direction * self.max_step else: r0t = r0 + direction * self.trial_step f0tp = self.dimeratoms.get_projected_forces(r0t) F = np.vdot((f0tp + f0p), direction) / 2.0 C = np.vdot((f0tp - f0p), direction) / self.trial_step step = (-F / C + self.trial_step / 2.0) * direction if norm(step) > self.max_step: step = direction * self.max_step self.log(f0p, norm(step)) atoms.set_positions(r + step) self.f0 = f.flat.copy() self.r0 = r.flat.copy() def get_cg_direction(self, direction): """Apply the Conjugate Gradient algorithm to the step direction.""" if self.cg_init: self.cg_init = False self.direction_old = direction.copy() self.cg_direction = direction.copy() old_norm = np.vdot(self.direction_old, self.direction_old) # Polak-Ribiere Conjugate Gradient if old_norm != 0.0: betaPR = np.vdot(direction, (direction - self.direction_old)) / \ old_norm else: betaPR = 0.0 if betaPR < 0.0: betaPR = 0.0 self.cg_direction = direction + self.cg_direction * betaPR self.direction_old = direction.copy() return self.cg_direction.copy() def log(self, f=None, stepsize=None): """Log each step of the optimization.""" if f is None: f = self.dimeratoms.get_forces() if self.logfile is not None: T = time.localtime() e = self.dimeratoms.get_potential_energy() fmax = sqrt((f**2).sum(axis=1).max()) rotsteps = self.dimeratoms.control.get_counter('rotcount') curvature = self.dimeratoms.get_curvature() l = '' if stepsize: if isinstance(self.control, DimerControl): l = '%s: %4d %02d:%02d:%02d %15.6f %12.4f %12.6f ' \ '%12.6f %10d\n' % ( 'MinModeTranslate', self.nsteps, T[3], T[4], T[5], e, fmax, stepsize, curvature, rotsteps) else: if isinstance(self.control, DimerControl): l = '%s: %4d %02d:%02d:%02d %15.6f %12.4f %s ' \ '%12.6f %10d\n' % ( 'MinModeTranslate', self.nsteps, T[3], T[4], T[5], e, fmax, ' --------', curvature, rotsteps) self.logfile.write(l) self.logfile.flush()
def read_eigenmode(mlog, index=-1): """Read an eigenmode. To access the pre optimization eigenmode set index = 'null'. """ mlog_is_str = isinstance(mlog, str) if mlog_is_str: fd = open(mlog) else: fd = mlog lines = fd.readlines() # Detect the amount of atoms and iterations k = 2 while lines[k].split()[1].lower() not in ['optimization', 'order']: k += 1 n = k - 2 n_itr = (len(lines) // (n + 1)) - 2 # Locate the correct image. if isinstance(index, str): if index.lower() == 'null': i = 0 else: i = int(index) + 1 else: if index >= 0: i = index + 1 else: if index < -n_itr - 1: raise IndexError('list index out of range') else: i = index mode = np.ndarray(shape=(n, 3), dtype=float) for k_atom, k in enumerate(range(1, n + 1)): line = lines[i * (n + 1) + k].split() for k_dim in range(3): mode[k_atom][k_dim] = float(line[k_dim + 2]) if mlog_is_str: fd.close() return mode # Aliases DimerAtoms = MinModeAtoms DimerTranslate = MinModeTranslate