Source code for ase.optimize.climbfixinternals

from numpy.linalg import norm
from ase.optimize.bfgs import BFGS
from ase.constraints import FixInternals

[docs]class BFGSClimbFixInternals(BFGS): """Class for transition state search and optimization Climbs the 1D reaction coordinate defined as constrained internal coordinate via the :class:`~ase.constraints.FixInternals` class while minimizing all remaining degrees of freedom. Details: Two optimizers, 'A' and 'B', are applied orthogonal to each other. Optimizer 'A' climbs the constrained coordinate while optimizer 'B' optimizes the remaining degrees of freedom after each climbing step. Optimizer 'A' uses the BFGS algorithm to climb along the projected force of the selected constraint. Optimizer 'B' can be any ASE optimizer (default: BFGS). In combination with other constraints, the order of constraints matters. Generally, the FixInternals constraint should come first in the list of constraints. This has been tested with the :class:`~ase.constraints.FixAtoms` constraint. Inspired by concepts described by P. N. Plessow. [1]_ .. [1] Plessow, P. N., Efficient Transition State Optimization of Periodic Structures through Automated Relaxed Potential Energy Surface Scans. J. Chem. Theory Comput. 2018, 14 (2), 981–990. .. note:: Convergence is based on 'fmax' of the total forces which is the sum of the projected forces and the forces of the remaining degrees of freedom. This value is logged in the 'logfile'. Optimizer 'B' logs 'fmax' of the remaining degrees of freedom without the projected forces. The projected forces can be inspected using the :meth:`get_projected_forces` method: >>> for _ in dyn.irun(): ... projected_forces = dyn.get_projected_forces() Example ------- .. literalinclude:: ../../ase/test/optimize/ :end-before: # end example for documentation """ def __init__(self, atoms, restart=None, logfile='-', trajectory=None, maxstep=None, master=None, alpha=None, climb_coordinate=None, optB=BFGS, optB_kwargs=None, optB_fmax=0.05, optB_fmax_scaling=0.0): """Allowed parameters are similar to the parent class :class:`~ase.optimize.bfgs.BFGS` with the following additions: Parameters ---------- climb_coordinate: list Specifies which subconstraint of the :class:`~ase.constraints.FixInternals` constraint is to be climbed. Provide the corresponding nested list of indices (including coefficients in the case of Combo constraints). See the example above. optB: any ASE optimizer, optional Optimizer 'B' for optimization of the remaining degrees of freedom after each climbing step. optB_kwargs: dict, optional Specifies keyword arguments to be passed to optimizer 'B' at its initialization. By default, optimizer 'B' writes a logfile and trajectory (optB_{...}.log, optB_{...}.traj) where {...} is the current value of the ``climb_coordinate``. Set ``logfile`` to '-' for console output. Set ``trajectory`` to 'None' to suppress writing of the trajectory file. optB_fmax: float, optional Specifies the convergence criterion 'fmax' of optimizer 'B'. optB_fmax_scaling: float, optional Scaling factor to dynamically tighten 'fmax' of optimizer 'B' to the value of ``optB_fmax`` when close to convergence. Can speed up the climbing process. The scaling formula is 'fmax' = ``optB_fmax`` + ``optB_fmax_scaling`` :math:`\\cdot` norm_of_projected_forces The final optimization with optimizer 'B' is performed with ``optB_fmax`` independent of ``optB_fmax_scaling``. """ self.targetvalue = None # may be assigned during restart in super().__init__(atoms, restart=restart, logfile=logfile, trajectory=trajectory, maxstep=maxstep, master=master, alpha=alpha) self.constr2climb = self.get_constr2climb(self.atoms, climb_coordinate) self.targetvalue = self.targetvalue or self.constr2climb.targetvalue self.optB = optB self.optB_kwargs = optB_kwargs or {} self.optB_fmax = optB_fmax self.scaling = optB_fmax_scaling # log optimizer 'B' in logfiles named after current value of constraint self.autolog = 'logfile' not in self.optB_kwargs self.autotraj = 'trajectory' not in self.optB_kwargs def get_constr2climb(self, atoms, climb_coordinate): """Get pointer to the subconstraint that is to be climbed. Identification by its definition via indices (and coefficients).""" constr = self.get_fixinternals(atoms) return constr.get_subconstraint(atoms, climb_coordinate) def get_fixinternals(self, atoms): """Get pointer to the FixInternals constraint on the atoms object.""" all_constr_types = list(map(type, atoms.constraints)) index = all_constr_types.index(FixInternals) # locate constraint return atoms.constraints[index] def read(self): (self.H, self.pos0, self.forces0, self.maxstep, self.targetvalue) = self.load() def step(self): self.relax_remaining_dof() # optimization with optimizer 'B' pos, dpos = self.pretend2climb() # with optimizer 'A' self.update_positions_and_targetvalue(pos, dpos) # obey other constr. self.dump((self.H, self.pos0, self.forces0, self.maxstep, self.targetvalue)) def pretend2climb(self): """Get directions for climbing and climb with optimizer 'A'.""" proj_forces = self.get_projected_forces() pos = self.atoms.get_positions() dpos, steplengths = self.prepare_step(pos, proj_forces) dpos = self.determine_step(dpos, steplengths) return pos, dpos def update_positions_and_targetvalue(self, pos, dpos): """Adjust constrained targetvalue of constraint and update positions.""" self.constr2climb.adjust_positions(pos, pos + dpos) # update sigma self.targetvalue += self.constr2climb.sigma # climb constraint self.constr2climb.targetvalue = self.targetvalue # adjust positions self.atoms.set_positions(self.atoms.get_positions()) # to targetvalue def relax_remaining_dof(self): """Optimize remaining degrees of freedom with optimizer 'B'.""" if self.autolog: self.optB_kwargs['logfile'] = f'optB_{self.targetvalue}.log' if self.autotraj: self.optB_kwargs['trajectory'] = f'optB_{self.targetvalue}.traj' fmax = self.get_scaled_fmax() with self.optB(self.atoms, **self.optB_kwargs) as opt: # optimize with scaled fmax if self.converged() and fmax > self.optB_fmax: # (final) optimization with desired fmax
[docs] def get_scaled_fmax(self): """Return the adaptive 'fmax' based on the estimated distance to the transition state.""" return (self.optB_fmax + self.scaling * norm(self.constr2climb.projected_forces))
[docs] def get_projected_forces(self): """Return the projected forces along the constrained coordinate in uphill direction (negative sign).""" forces = self.constr2climb.projected_forces forces = -forces.reshape(self.atoms.positions.shape) return forces
def get_total_forces(self): """Return forces obeying all constraints plus projected forces.""" return self.atoms.get_forces() + self.get_projected_forces() def converged(self, forces=None): """Did the optimization converge based on the total forces?""" forces = forces or self.get_total_forces() return super().converged(forces=forces) def log(self, forces=None): forces = forces or self.get_total_forces() super().log(forces=forces)