Coverage for ase / optimize / bfgslinesearch.py: 83.97%
131 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-03-30 08:22 +0000
1# fmt: off
3# ******NOTICE***************
4# optimize.py module by Travis E. Oliphant
5#
6# You may copy and use this module as you see fit with no
7# guarantee implied provided you keep this notice in all copies.
8# *****END NOTICE************
10import time
11from typing import IO
13import numpy as np
14from numpy import absolute, eye, isinf
16from ase import Atoms
17from ase.optimize.optimize import Optimizer
18from ase.utils.linesearch import LineSearch
20# These have been copied from Numeric's MLab.py
21# I don't think they made the transition to scipy_core
23# Modified from scipy_optimize
24abs = absolute
25pymin = min
26pymax = max
27__version__ = '0.1'
30class BFGSLineSearch(Optimizer):
31 def __init__(
32 self,
33 atoms: Atoms,
34 restart: str | None = None,
35 logfile: IO | str = '-',
36 maxstep: float | None = None,
37 trajectory: str | None = None,
38 c1: float = 0.23,
39 c2: float = 0.46,
40 alpha: float = 10.0,
41 stpmax: float = 50.0,
42 **kwargs,
43 ):
44 """Optimize atomic positions in the BFGSLineSearch algorithm, which
45 uses both forces and potential energy information.
47 Parameters
48 ----------
49 atoms: :class:`~ase.Atoms`
50 The Atoms object to relax.
52 restart: str
53 JSON file used to store hessian matrix. If set, file with
54 such a name will be searched and hessian matrix stored will
55 be used, if the file exists.
57 trajectory: str
58 Trajectory file used to store optimisation path.
60 maxstep: float
61 Used to set the maximum distance an atom can move per
62 iteration (default value is 0.2 Angstroms).
64 logfile: file object or str
65 If *logfile* is a string, a file with that name will be opened.
66 Use '-' for stdout.
68 kwargs : dict, optional
69 Extra arguments passed to
70 :class:`~ase.optimize.optimize.Optimizer`.
72 """
73 if maxstep is None:
74 self.maxstep = self.defaults['maxstep']
75 else:
76 self.maxstep = maxstep
77 self.stpmax = stpmax
78 self.alpha = alpha
79 self.H = None
80 self.c1 = c1
81 self.c2 = c2
82 self.force_calls = 0
83 self.function_calls = 0
84 self.r0 = None
85 self.g0 = None
86 self.e0 = None
87 self.load_restart = False
88 self.task = 'START'
89 self.rep_count = 0
90 self.p = None
91 self.alpha_k = None
92 self.no_update = False
93 self.replay = False
95 super().__init__(atoms, restart, logfile, trajectory, **kwargs)
97 def read(self):
98 self.r0, self.g0, self.e0, self.task, self.H = self.load()
99 self.load_restart = True
101 def reset(self):
102 self.H = None
103 self.r0 = None
104 self.g0 = None
105 self.e0 = None
106 self.rep_count = 0
108 def step(self, gradient=None):
109 optimizable = self.optimizable
110 gradient = self._get_gradient(gradient)
111 r = optimizable.get_x()
112 g = gradient / self.alpha
113 p0 = self.p
114 self.update(r, g, self.r0, self.g0, p0)
115 # o,v = np.linalg.eigh(self.B)
116 e = self.func(r)
118 self.p = -np.dot(self.H, g)
119 p_size = np.sqrt((self.p**2).sum())
120 if p_size <= np.sqrt(optimizable.ndofs() / 3 * 1e-10):
121 self.p /= (p_size / np.sqrt(optimizable.ndofs() / 3 * 1e-10))
122 ls = LineSearch(get_gradient_norm=self.optimizable.gradient_norm)
123 self.alpha_k, e, self.e0, self.no_update = \
124 ls._line_search(self.func, self.fprime, r, self.p, g, e, self.e0,
125 maxstep=self.maxstep, c1=self.c1,
126 c2=self.c2, stpmax=self.stpmax)
127 if self.alpha_k is None:
128 raise RuntimeError("LineSearch failed!")
130 dr = self.alpha_k * self.p
131 optimizable.set_x(r + dr)
132 self.r0 = r
133 self.g0 = g
134 self.dump((self.r0, self.g0, self.e0, self.task, self.H))
136 def update(self, r, g, r0, g0, p0):
137 self.I = eye(self.optimizable.ndofs(), dtype=int)
138 if self.H is None:
139 self.H = eye(self.optimizable.ndofs())
140 # self.B = np.linalg.inv(self.H)
141 return
142 else:
143 dr = r - r0
144 dg = g - g0
145 # self.alpha_k can be None!!!
146 if not (((self.alpha_k or 0) > 0 and
147 abs(np.dot(g, p0)) - abs(np.dot(g0, p0)) < 0) or
148 self.replay):
149 return
150 if self.no_update is True:
151 print('skip update')
152 return
154 try: # this was handled in numeric, let it remain for more safety
155 rhok = 1.0 / (np.dot(dg, dr))
156 except ZeroDivisionError:
157 rhok = 1000.0
158 print("Divide-by-zero encountered: rhok assumed large")
159 if isinf(rhok): # this is patch for np
160 rhok = 1000.0
161 print("Divide-by-zero encountered: rhok assumed large")
162 A1 = self.I - dr[:, np.newaxis] * dg[np.newaxis, :] * rhok
163 A2 = self.I - dg[:, np.newaxis] * dr[np.newaxis, :] * rhok
164 self.H = (np.dot(A1, np.dot(self.H, A2)) +
165 rhok * dr[:, np.newaxis] * dr[np.newaxis, :])
166 # self.B = np.linalg.inv(self.H)
168 def func(self, x):
169 """Objective function for use of the optimizers"""
170 self.optimizable.set_x(x)
171 self.function_calls += 1
172 # Scale the problem as SciPy uses I as initial Hessian.
173 return self.optimizable.get_value() / self.alpha
175 def fprime(self, x):
176 """Gradient of the objective function for use of the optimizers"""
177 self.optimizable.set_x(x)
178 self.force_calls += 1
179 # Scale the problem as SciPy uses I as initial Hessian.
180 gradient = self.optimizable.get_gradient()
181 return gradient / self.alpha
183 def replay_trajectory(self, traj):
184 """Initialize hessian from old trajectory."""
185 self.replay = True
186 from ase.utils import IOContext
188 with IOContext() as files:
189 if isinstance(traj, str):
190 from ase.io.trajectory import Trajectory
191 traj = files.closelater(Trajectory(traj, mode='r'))
193 r0 = None
194 g0 = None
195 for i in range(len(traj) - 1):
196 r = traj[i].get_positions().ravel()
197 g = - traj[i].get_forces().ravel() / self.alpha
198 self.update(r, g, r0, g0, self.p)
199 self.p = -np.dot(self.H, g)
200 r0 = r.copy()
201 g0 = g.copy()
202 self.r0 = r0
203 self.g0 = g0
205 def log(self, gradient):
206 if self.logfile is None:
207 return
208 fmax = self.optimizable.gradient_norm(gradient)
209 e = self.optimizable.get_value()
210 T = time.localtime()
211 name = self.__class__.__name__
212 w = self.logfile.write
213 if self.nsteps == 0:
214 w('%s %4s[%3s] %8s %15s %12s\n' %
215 (' ' * len(name), 'Step', 'FC', 'Time', 'Energy', 'fmax'))
216 w('%s: %3d[%3d] %02d:%02d:%02d %15.6f %12.4f\n'
217 % (name, self.nsteps, self.force_calls, T[3], T[4], T[5], e,
218 fmax))
221def wrap_function(function, args):
222 ncalls = [0]
224 def function_wrapper(x):
225 ncalls[0] += 1
226 return function(x, *args)
227 return ncalls, function_wrapper