Coverage for /builds/ericyuan00000/ase/ase/optimize/climbfixinternals.py: 98.51%
67 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-06-18 01:20 +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 read(self):
126 (self.H, self.pos0, self.forces0, self.maxstep,
127 self.targetvalue) = self.load()
129 def step(self):
130 self.relax_remaining_dof() # optimization with optimizer 'B'
132 pos, dpos = self.pretend2climb() # with optimizer 'A'
133 self.update_positions_and_targetvalue(pos, dpos) # obey other constr.
135 self.dump((self.H, self.pos0, self.forces0, self.maxstep,
136 self.targetvalue))
138 def pretend2climb(self):
139 """Get directions for climbing and climb with optimizer 'A'."""
140 proj_forces = self.get_projected_forces()
141 pos = self.optimizable.get_positions()
142 dpos, steplengths = self.prepare_step(pos, proj_forces)
143 dpos = self.determine_step(dpos, steplengths)
144 return pos, dpos
146 def update_positions_and_targetvalue(self, pos, dpos):
147 """Adjust constrained targetvalue of constraint and update positions."""
148 self.constr2climb.adjust_positions(pos, pos + dpos) # update sigma
149 self.targetvalue += self.constr2climb.sigma # climb constraint
150 self.constr2climb.targetvalue = self.targetvalue # adjust positions
151 self.optimizable.set_positions(
152 self.optimizable.get_positions()) # to targetvalue
154 def relax_remaining_dof(self):
155 """Optimize remaining degrees of freedom with optimizer 'B'."""
156 if self.autolog:
157 self.optB_kwargs['logfile'] = f'optB_{self.targetvalue}.log'
158 if self.autotraj:
159 self.optB_kwargs['trajectory'] = f'optB_{self.targetvalue}.traj'
160 fmax = self.get_scaled_fmax()
161 with self.optB(self.optimizable.atoms, **self.optB_kwargs) as opt:
162 opt.run(fmax) # optimize with scaled fmax
163 if self.converged() and fmax > self.optB_fmax:
164 # (final) optimization with desired fmax
165 opt.run(self.optB_fmax)
167 def get_scaled_fmax(self):
168 """Return the adaptive 'fmax' based on the estimated distance to the
169 transition state."""
170 return (self.optB_fmax +
171 self.scaling * norm(self.constr2climb.projected_forces))
173 def get_projected_forces(self):
174 """Return the projected forces along the constrained coordinate in
175 uphill direction (negative sign)."""
176 forces = self.constr2climb.projected_forces
177 forces = -forces.reshape(self.optimizable.get_positions().shape)
178 return forces
180 def get_total_forces(self):
181 """Return forces obeying all constraints plus projected forces."""
182 return self.optimizable.get_forces() + self.get_projected_forces()
184 def converged(self, forces=None):
185 """Did the optimization converge based on the total forces?"""
186 forces = forces or self.get_total_forces()
187 return super().converged(forces=forces)
189 def log(self, forces=None):
190 forces = forces or self.get_total_forces()
191 super().log(forces=forces)
194def get_fixinternals(atoms):
195 """Get pointer to the FixInternals constraint on the atoms object."""
196 all_constr_types = list(map(type, atoms.constraints))
197 index = all_constr_types.index(FixInternals) # locate constraint
198 return atoms.constraints[index]
201def get_constr2climb(atoms, climb_coordinate):
202 """Get pointer to the subconstraint that is to be climbed.
203 Identification by its definition via indices (and coefficients)."""
204 constr = get_fixinternals(atoms)
205 return constr.get_subconstraint(atoms, climb_coordinate)