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