Coverage for /builds/kinetik161/ase/ase/optimize/climbfixinternals.py: 98.51%
67 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +0000
1from typing import IO, Any, Dict, List, Optional, Type, Union
3from numpy.linalg import norm
5from ase import Atoms
6from ase.constraints import FixInternals
7from ase.optimize.bfgs import BFGS
8from ase.optimize.optimize import Optimizer
11class BFGSClimbFixInternals(BFGS):
12 """Class for transition state search and optimization
14 Climbs the 1D reaction coordinate defined as constrained internal coordinate
15 via the :class:`~ase.constraints.FixInternals` class while minimizing all
16 remaining degrees of freedom.
18 Details: Two optimizers, 'A' and 'B', are applied orthogonal to each other.
19 Optimizer 'A' climbs the constrained coordinate while optimizer 'B'
20 optimizes the remaining degrees of freedom after each climbing step.
21 Optimizer 'A' uses the BFGS algorithm to climb along the projected force of
22 the selected constraint. Optimizer 'B' can be any ASE optimizer
23 (default: BFGS).
25 In combination with other constraints, the order of constraints matters.
26 Generally, the FixInternals constraint should come first in the list of
27 constraints.
28 This has been tested with the :class:`~ase.constraints.FixAtoms` constraint.
30 Inspired by concepts described by P. N. Plessow. [1]_
32 .. [1] Plessow, P. N., Efficient Transition State Optimization of Periodic
33 Structures through Automated Relaxed Potential Energy Surface Scans.
34 J. Chem. Theory Comput. 2018, 14 (2), 981–990.
35 https://doi.org/10.1021/acs.jctc.7b01070.
37 .. note::
38 Convergence is based on 'fmax' of the total forces which is the sum of
39 the projected forces and the forces of the remaining degrees of freedom.
40 This value is logged in the 'logfile'. Optimizer 'B' logs 'fmax' of the
41 remaining degrees of freedom without the projected forces. The projected
42 forces can be inspected using the :meth:`get_projected_forces` method:
44 >>> for _ in dyn.irun():
45 ... projected_forces = dyn.get_projected_forces()
47 Example
48 -------
49 .. literalinclude:: ../../ase/test/optimize/test_climb_fix_internals.py
50 :end-before: # end example for documentation
51 """
53 def __init__(
54 self,
55 atoms: Atoms,
56 restart: Optional[str] = None,
57 logfile: Union[IO, str] = '-',
58 trajectory: Optional[str] = None,
59 maxstep: Optional[float] = None,
60 master: Optional[bool] = None,
61 alpha: Optional[float] = None,
62 climb_coordinate: Optional[List[FixInternals]] = None,
63 optB: Type[Optimizer] = BFGS,
64 optB_kwargs: Optional[Dict[str, Any]] = None,
65 optB_fmax: float = 0.05,
66 optB_fmax_scaling: float = 0.0,
67 ):
68 """Allowed parameters are similar to the parent class
69 :class:`~ase.optimize.bfgs.BFGS` with the following additions:
71 Parameters
72 ----------
73 climb_coordinate: list
74 Specifies which subconstraint of the
75 :class:`~ase.constraints.FixInternals` constraint is to be climbed.
76 Provide the corresponding nested list of indices
77 (including coefficients in the case of Combo constraints).
78 See the example above.
80 optB: any ASE optimizer, optional
81 Optimizer 'B' for optimization of the remaining degrees of freedom
82 after each climbing step.
84 optB_kwargs: dict, optional
85 Specifies keyword arguments to be passed to optimizer 'B' at its
86 initialization. By default, optimizer 'B' writes a logfile and
87 trajectory (optB_{...}.log, optB_{...}.traj) where {...} is the
88 current value of the ``climb_coordinate``. Set ``logfile`` to '-'
89 for console output. Set ``trajectory`` to 'None' to suppress
90 writing of the trajectory file.
92 optB_fmax: float, optional
93 Specifies the convergence criterion 'fmax' of optimizer 'B'.
95 optB_fmax_scaling: float, optional
96 Scaling factor to dynamically tighten 'fmax' of optimizer 'B' to
97 the value of ``optB_fmax`` when close to convergence.
98 Can speed up the climbing process. The scaling formula is
100 'fmax' = ``optB_fmax`` + ``optB_fmax_scaling``
101 :math:`\\cdot` norm_of_projected_forces
103 The final optimization with optimizer 'B' is
104 performed with ``optB_fmax`` independent of ``optB_fmax_scaling``.
105 """
106 self.targetvalue = None # may be assigned during restart in self.read()
107 super().__init__(atoms, restart=restart, logfile=logfile,
108 trajectory=trajectory, maxstep=maxstep, master=master,
109 alpha=alpha)
111 self.constr2climb = get_constr2climb(
112 self.optimizable.atoms, climb_coordinate)
113 self.targetvalue = self.targetvalue or self.constr2climb.targetvalue
115 self.optB = optB
116 self.optB_kwargs = optB_kwargs or {}
117 self.optB_fmax = optB_fmax
118 self.scaling = optB_fmax_scaling
119 # log optimizer 'B' in logfiles named after current value of constraint
120 self.autolog = 'logfile' not in self.optB_kwargs
121 self.autotraj = 'trajectory' not in self.optB_kwargs
123 def read(self):
124 (self.H, self.pos0, self.forces0, self.maxstep,
125 self.targetvalue) = self.load()
127 def step(self):
128 self.relax_remaining_dof() # optimization with optimizer 'B'
130 pos, dpos = self.pretend2climb() # with optimizer 'A'
131 self.update_positions_and_targetvalue(pos, dpos) # obey other constr.
133 self.dump((self.H, self.pos0, self.forces0, self.maxstep,
134 self.targetvalue))
136 def pretend2climb(self):
137 """Get directions for climbing and climb with optimizer 'A'."""
138 proj_forces = self.get_projected_forces()
139 pos = self.optimizable.get_positions()
140 dpos, steplengths = self.prepare_step(pos, proj_forces)
141 dpos = self.determine_step(dpos, steplengths)
142 return pos, dpos
144 def update_positions_and_targetvalue(self, pos, dpos):
145 """Adjust constrained targetvalue of constraint and update positions."""
146 self.constr2climb.adjust_positions(pos, pos + dpos) # update sigma
147 self.targetvalue += self.constr2climb.sigma # climb constraint
148 self.constr2climb.targetvalue = self.targetvalue # adjust positions
149 self.optimizable.set_positions(
150 self.optimizable.get_positions()) # to targetvalue
152 def relax_remaining_dof(self):
153 """Optimize remaining degrees of freedom with optimizer 'B'."""
154 if self.autolog:
155 self.optB_kwargs['logfile'] = f'optB_{self.targetvalue}.log'
156 if self.autotraj:
157 self.optB_kwargs['trajectory'] = f'optB_{self.targetvalue}.traj'
158 fmax = self.get_scaled_fmax()
159 with self.optB(self.optimizable.atoms, **self.optB_kwargs) as opt:
160 opt.run(fmax) # optimize with scaled fmax
161 if self.converged() and fmax > self.optB_fmax:
162 # (final) optimization with desired fmax
163 opt.run(self.optB_fmax)
165 def get_scaled_fmax(self):
166 """Return the adaptive 'fmax' based on the estimated distance to the
167 transition state."""
168 return (self.optB_fmax +
169 self.scaling * norm(self.constr2climb.projected_forces))
171 def get_projected_forces(self):
172 """Return the projected forces along the constrained coordinate in
173 uphill direction (negative sign)."""
174 forces = self.constr2climb.projected_forces
175 forces = -forces.reshape(self.optimizable.get_positions().shape)
176 return forces
178 def get_total_forces(self):
179 """Return forces obeying all constraints plus projected forces."""
180 return self.optimizable.get_forces() + self.get_projected_forces()
182 def converged(self, forces=None):
183 """Did the optimization converge based on the total forces?"""
184 forces = forces or self.get_total_forces()
185 return super().converged(forces=forces)
187 def log(self, forces=None):
188 forces = forces or self.get_total_forces()
189 super().log(forces=forces)
192def get_fixinternals(atoms):
193 """Get pointer to the FixInternals constraint on the atoms object."""
194 all_constr_types = list(map(type, atoms.constraints))
195 index = all_constr_types.index(FixInternals) # locate constraint
196 return atoms.constraints[index]
199def get_constr2climb(atoms, climb_coordinate):
200 """Get pointer to the subconstraint that is to be climbed.
201 Identification by its definition via indices (and coefficients)."""
202 constr = get_fixinternals(atoms)
203 return constr.get_subconstraint(atoms, climb_coordinate)