Coverage for /builds/ase/ase/ase/optimize/lbfgs.py: 80.33%
122 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +0000
1# fmt: off
3from typing import IO, Optional, Union
5import numpy as np
7from ase import Atoms
8from ase.optimize.optimize import Optimizer
9from ase.utils.linesearch import LineSearch
12class LBFGS(Optimizer):
13 """Limited memory BFGS optimizer.
15 A limited memory version of the bfgs algorithm. Unlike the bfgs algorithm
16 used in bfgs.py, the inverse of Hessian matrix is updated. The inverse
17 Hessian is represented only as a diagonal matrix to save memory
19 """
21 def __init__(
22 self,
23 atoms: Atoms,
24 restart: Optional[str] = None,
25 logfile: Union[IO, str] = '-',
26 trajectory: Optional[str] = None,
27 maxstep: Optional[float] = None,
28 memory: int = 100,
29 damping: float = 1.0,
30 alpha: float = 70.0,
31 use_line_search: bool = False,
32 **kwargs,
33 ):
34 """
36 Parameters
37 ----------
38 atoms: :class:`~ase.Atoms`
39 The Atoms object to relax.
41 restart: str
42 JSON file used to store vectors for updating the inverse of
43 Hessian matrix. If set, file with such a name will be searched
44 and information stored will be used, if the file exists.
46 logfile: file object or str
47 If *logfile* is a string, a file with that name will be opened.
48 Use '-' for stdout.
50 trajectory: string
51 Trajectory file used to store optimisation path.
53 maxstep: float
54 How far is a single atom allowed to move. This is useful for DFT
55 calculations where wavefunctions can be reused if steps are small.
56 Default is 0.2 Angstrom.
58 memory: int
59 Number of steps to be stored. Default value is 100. Three numpy
60 arrays of this length containing floats are stored.
62 damping: float
63 The calculated step is multiplied with this number before added to
64 the positions.
66 alpha: float
67 Initial guess for the Hessian (curvature of energy surface). A
68 conservative value of 70.0 is the default, but number of needed
69 steps to converge might be less if a lower value is used. However,
70 a lower value also means risk of instability.
72 kwargs : dict, optional
73 Extra arguments passed to
74 :class:`~ase.optimize.optimize.Optimizer`.
76 """
77 Optimizer.__init__(self, atoms, restart, logfile, trajectory, **kwargs)
79 if maxstep is not None:
80 self.maxstep = maxstep
81 else:
82 self.maxstep = self.defaults['maxstep']
84 if self.maxstep > 1.0:
85 raise ValueError('You are using a much too large value for ' +
86 'the maximum step size: %.1f Angstrom' %
87 self.maxstep)
89 self.memory = memory
90 # Initial approximation of inverse Hessian 1./70. is to emulate the
91 # behaviour of BFGS. Note that this is never changed!
92 self.H0 = 1. / alpha
93 self.damping = damping
94 self.use_line_search = use_line_search
95 self.p = None
96 self.function_calls = 0
97 self.force_calls = 0
99 def initialize(self):
100 """Initialize everything so no checks have to be done in step"""
101 self.iteration = 0
102 self.s = []
103 self.y = []
104 # Store also rho, to avoid calculating the dot product again and
105 # again.
106 self.rho = []
108 self.r0 = None
109 self.f0 = None
110 self.e0 = None
111 self.task = 'START'
112 self.load_restart = False
114 def read(self):
115 """Load saved arrays to reconstruct the Hessian"""
116 self.iteration, self.s, self.y, self.rho, \
117 self.r0, self.f0, self.e0, self.task = self.load()
118 self.load_restart = True
120 def step(self, forces=None):
121 """Take a single step
123 Use the given forces, update the history and calculate the next step --
124 then take it"""
126 if forces is None:
127 forces = self.optimizable.get_gradient().reshape(-1, 3)
129 pos = self.optimizable.get_x().reshape(-1, 3)
131 self.update(pos, forces, self.r0, self.f0)
133 s = self.s
134 y = self.y
135 rho = self.rho
136 H0 = self.H0
138 loopmax = np.min([self.memory, self.iteration])
139 a = np.empty((loopmax,), dtype=np.float64)
141 # ## The algorithm itself:
142 q = -forces.reshape(-1)
143 for i in range(loopmax - 1, -1, -1):
144 a[i] = rho[i] * np.dot(s[i], q)
145 q -= a[i] * y[i]
146 z = H0 * q
148 for i in range(loopmax):
149 b = rho[i] * np.dot(y[i], z)
150 z += s[i] * (a[i] - b)
152 self.p = - z.reshape((-1, 3))
153 # ##
155 g = -forces
156 if self.use_line_search is True:
157 e = self.func(pos)
158 self.line_search(pos, g, e)
159 dr = (self.alpha_k * self.p).reshape(-1, 3)
160 else:
161 self.force_calls += 1
162 self.function_calls += 1
163 dr = self.determine_step(self.p) * self.damping
164 self.optimizable.set_x((pos + dr).ravel())
166 self.iteration += 1
167 self.r0 = pos
168 self.f0 = -g
169 self.dump((self.iteration, self.s, self.y,
170 self.rho, self.r0, self.f0, self.e0, self.task))
172 def determine_step(self, dr):
173 """Determine step to take according to maxstep
175 Normalize all steps as the largest step. This way
176 we still move along the eigendirection.
177 """
178 steplengths = (dr**2).sum(1)**0.5
179 longest_step = np.max(steplengths)
180 if longest_step >= self.maxstep:
181 dr *= self.maxstep / longest_step
183 return dr
185 def update(self, pos, forces, r0, f0):
186 """Update everything that is kept in memory
188 This function is mostly here to allow for replay_trajectory.
189 """
190 if self.iteration > 0:
191 s0 = pos.reshape(-1) - r0.reshape(-1)
192 self.s.append(s0)
194 # We use the gradient which is minus the force!
195 y0 = f0.reshape(-1) - forces.reshape(-1)
196 self.y.append(y0)
198 rho0 = 1.0 / np.dot(y0, s0)
199 self.rho.append(rho0)
201 if self.iteration > self.memory:
202 self.s.pop(0)
203 self.y.pop(0)
204 self.rho.pop(0)
206 def replay_trajectory(self, traj):
207 """Initialize history from old trajectory."""
208 if isinstance(traj, str):
209 from ase.io.trajectory import Trajectory
210 traj = Trajectory(traj, 'r')
211 r0 = None
212 f0 = None
213 # The last element is not added, as we get that for free when taking
214 # the first qn-step after the replay
215 for i in range(len(traj) - 1):
216 pos = traj[i].get_positions()
217 forces = traj[i].get_forces()
218 self.update(pos, forces, r0, f0)
219 r0 = pos.copy()
220 f0 = forces.copy()
221 self.iteration += 1
222 self.r0 = r0
223 self.f0 = f0
225 def func(self, x):
226 """Objective function for use of the optimizers"""
227 self.optimizable.set_x(x)
228 self.function_calls += 1
229 return self.optimizable.get_value()
231 def fprime(self, x):
232 """Gradient of the objective function for use of the optimizers"""
233 self.optimizable.set_x(x)
234 self.force_calls += 1
235 # Remember that forces are minus the gradient!
236 return -self.optimizable.get_gradient()
238 def line_search(self, r, g, e):
239 self.p = self.p.ravel()
240 p_size = np.sqrt((self.p**2).sum())
241 if p_size <= np.sqrt(self.optimizable.ndofs() / 3 * 1e-10):
242 self.p /= (p_size / np.sqrt(self.optimizable.ndofs() / 3 * 1e-10))
243 g = g.ravel()
244 r = r.ravel()
245 ls = LineSearch()
246 self.alpha_k, e, self.e0, self.no_update = \
247 ls._line_search(self.func, self.fprime, r, self.p, g, e, self.e0,
248 maxstep=self.maxstep, c1=.23,
249 c2=.46, stpmax=50.)
250 if self.alpha_k is None:
251 raise RuntimeError('LineSearch failed!')
254class LBFGSLineSearch(LBFGS):
255 """This optimizer uses the LBFGS algorithm, but does a line search that
256 fulfills the Wolff conditions.
257 """
259 def __init__(self, *args, **kwargs):
260 kwargs['use_line_search'] = True
261 LBFGS.__init__(self, *args, **kwargs)
263# """Modified version of LBFGS.
264#
265# This optimizer uses the LBFGS algorithm, but does a line search for the
266# minimum along the search direction. This is done by issuing an additional
267# force call for each step, thus doubling the number of calculations.
268#
269# Additionally the Hessian is reset if the new guess is not sufficiently
270# better than the old one.
271# """
272# def __init__(self, *args, **kwargs):
273# self.dR = kwargs.pop('dR', 0.1)
274# LBFGS.__init__(self, *args, **kwargs)
275#
276# def update(self, r, f, r0, f0):
277# """Update everything that is kept in memory
278#
279# This function is mostly here to allow for replay_trajectory.
280# """
281# if self.iteration > 0:
282# a1 = abs(np.dot(f.reshape(-1), f0.reshape(-1)))
283# a2 = np.dot(f0.reshape(-1), f0.reshape(-1))
284# if not (a1 <= 0.5 * a2 and a2 != 0):
285# # Reset optimization
286# self.initialize()
287#
288# # Note that the reset above will set self.iteration to 0 again
289# # which is why we should check again
290# if self.iteration > 0:
291# s0 = r.reshape(-1) - r0.reshape(-1)
292# self.s.append(s0)
293#
294# # We use the gradient which is minus the force!
295# y0 = f0.reshape(-1) - f.reshape(-1)
296# self.y.append(y0)
297#
298# rho0 = 1.0 / np.dot(y0, s0)
299# self.rho.append(rho0)
300#
301# if self.iteration > self.memory:
302# self.s.pop(0)
303# self.y.pop(0)
304# self.rho.pop(0)
305#
306# def determine_step(self, dr):
307# f = self.atoms.get_forces()
308#
309# # Unit-vector along the search direction
310# du = dr / np.sqrt(np.dot(dr.reshape(-1), dr.reshape(-1)))
311#
312# # We keep the old step determination before we figure
313# # out what is the best to do.
314# maxstep = self.maxstep * np.sqrt(3 * len(self.atoms))
315#
316# # Finite difference step using temporary point
317# self.atoms.positions += (du * self.dR)
318# # Decide how much to move along the line du
319# Fp1 = np.dot(f.reshape(-1), du.reshape(-1))
320# Fp2 = np.dot(self.atoms.get_forces().reshape(-1), du.reshape(-1))
321# CR = (Fp1 - Fp2) / self.dR
322# #RdR = Fp1*0.1
323# if CR < 0.0:
324# #print "negcurve"
325# RdR = maxstep
326# #if(abs(RdR) > maxstep):
327# # RdR = self.sign(RdR) * maxstep
328# else:
329# Fp = (Fp1 + Fp2) * 0.5
330# RdR = Fp / CR
331# if abs(RdR) > maxstep:
332# RdR = np.sign(RdR) * maxstep
333# else:
334# RdR += self.dR * 0.5
335# return du * RdR