1import numpy as np 

2 

3from ase.optimize.optimize import Optimizer 

4from ase.utils.linesearch import LineSearch 

5 

6 

7class LBFGS(Optimizer): 

8 """Limited memory BFGS optimizer. 

9 

10 A limited memory version of the bfgs algorithm. Unlike the bfgs algorithm 

11 used in bfgs.py, the inverse of Hessian matrix is updated. The inverse 

12 Hessian is represented only as a diagonal matrix to save memory 

13 

14 """ 

15 def __init__(self, atoms, restart=None, logfile='-', trajectory=None, 

16 maxstep=None, memory=100, damping=1.0, alpha=70.0, 

17 use_line_search=False, master=None, 

18 force_consistent=None): 

19 """Parameters: 

20 

21 atoms: Atoms object 

22 The Atoms object to relax. 

23 

24 restart: string 

25 Pickle file used to store vectors for updating the inverse of 

26 Hessian matrix. If set, file with such a name will be searched 

27 and information stored will be used, if the file exists. 

28 

29 logfile: file object or str 

30 If *logfile* is a string, a file with that name will be opened. 

31 Use '-' for stdout. 

32 

33 trajectory: string 

34 Pickle file used to store trajectory of atomic movement. 

35 

36 maxstep: float 

37 How far is a single atom allowed to move. This is useful for DFT 

38 calculations where wavefunctions can be reused if steps are small. 

39 Default is 0.2 Angstrom. 

40 

41 memory: int 

42 Number of steps to be stored. Default value is 100. Three numpy 

43 arrays of this length containing floats are stored. 

44 

45 damping: float 

46 The calculated step is multiplied with this number before added to 

47 the positions. 

48 

49 alpha: float 

50 Initial guess for the Hessian (curvature of energy surface). A 

51 conservative value of 70.0 is the default, but number of needed 

52 steps to converge might be less if a lower value is used. However, 

53 a lower value also means risk of instability. 

54 

55 master: boolean 

56 Defaults to None, which causes only rank 0 to save files. If 

57 set to true, this rank will save files. 

58 

59 force_consistent: boolean or None 

60 Use force-consistent energy calls (as opposed to the energy 

61 extrapolated to 0 K). By default (force_consistent=None) uses 

62 force-consistent energies if available in the calculator, but 

63 falls back to force_consistent=False if not. 

64 """ 

65 Optimizer.__init__(self, atoms, restart, logfile, trajectory, master, 

66 force_consistent=force_consistent) 

67 

68 if maxstep is not None: 

69 self.maxstep = maxstep 

70 else: 

71 self.maxstep = self.defaults['maxstep'] 

72 

73 if self.maxstep > 1.0: 

74 raise ValueError('You are using a much too large value for ' + 

75 'the maximum step size: %.1f Angstrom' % 

76 maxstep) 

77 

78 self.memory = memory 

79 # Initial approximation of inverse Hessian 1./70. is to emulate the 

80 # behaviour of BFGS. Note that this is never changed! 

81 self.H0 = 1. / alpha 

82 self.damping = damping 

83 self.use_line_search = use_line_search 

84 self.p = None 

85 self.function_calls = 0 

86 self.force_calls = 0 

87 

88 def initialize(self): 

89 """Initialize everything so no checks have to be done in step""" 

90 self.iteration = 0 

91 self.s = [] 

92 self.y = [] 

93 # Store also rho, to avoid calculating the dot product again and 

94 # again. 

95 self.rho = [] 

96 

97 self.r0 = None 

98 self.f0 = None 

99 self.e0 = None 

100 self.task = 'START' 

101 self.load_restart = False 

102 

103 def read(self): 

104 """Load saved arrays to reconstruct the Hessian""" 

105 self.iteration, self.s, self.y, self.rho, \ 

106 self.r0, self.f0, self.e0, self.task = self.load() 

107 self.load_restart = True 

108 

109 def step(self, f=None): 

110 """Take a single step 

111 

112 Use the given forces, update the history and calculate the next step -- 

113 then take it""" 

114 

115 if f is None: 

116 f = self.atoms.get_forces() 

117 

118 r = self.atoms.get_positions() 

119 

120 self.update(r, f, self.r0, self.f0) 

121 

122 s = self.s 

123 y = self.y 

124 rho = self.rho 

125 H0 = self.H0 

126 

127 loopmax = np.min([self.memory, self.iteration]) 

128 a = np.empty((loopmax,), dtype=np.float64) 

129 

130 # ## The algorithm itself: 

131 q = -f.reshape(-1) 

132 for i in range(loopmax - 1, -1, -1): 

133 a[i] = rho[i] * np.dot(s[i], q) 

134 q -= a[i] * y[i] 

135 z = H0 * q 

136 

137 for i in range(loopmax): 

138 b = rho[i] * np.dot(y[i], z) 

139 z += s[i] * (a[i] - b) 

140 

141 self.p = - z.reshape((-1, 3)) 

142 # ## 

143 

144 g = -f 

145 if self.use_line_search is True: 

146 e = self.func(r) 

147 self.line_search(r, g, e) 

148 dr = (self.alpha_k * self.p).reshape(len(self.atoms), -1) 

149 else: 

150 self.force_calls += 1 

151 self.function_calls += 1 

152 dr = self.determine_step(self.p) * self.damping 

153 self.atoms.set_positions(r + dr) 

154 

155 self.iteration += 1 

156 self.r0 = r 

157 self.f0 = -g 

158 self.dump((self.iteration, self.s, self.y, 

159 self.rho, self.r0, self.f0, self.e0, self.task)) 

160 

161 def determine_step(self, dr): 

162 """Determine step to take according to maxstep 

163 

164 Normalize all steps as the largest step. This way 

165 we still move along the eigendirection. 

166 """ 

167 steplengths = (dr**2).sum(1)**0.5 

168 longest_step = np.max(steplengths) 

169 if longest_step >= self.maxstep: 

170 dr *= self.maxstep / longest_step 

171 

172 return dr 

173 

174 def update(self, r, f, r0, f0): 

175 """Update everything that is kept in memory 

176 

177 This function is mostly here to allow for replay_trajectory. 

178 """ 

179 if self.iteration > 0: 

180 s0 = r.reshape(-1) - r0.reshape(-1) 

181 self.s.append(s0) 

182 

183 # We use the gradient which is minus the force! 

184 y0 = f0.reshape(-1) - f.reshape(-1) 

185 self.y.append(y0) 

186 

187 rho0 = 1.0 / np.dot(y0, s0) 

188 self.rho.append(rho0) 

189 

190 if self.iteration > self.memory: 

191 self.s.pop(0) 

192 self.y.pop(0) 

193 self.rho.pop(0) 

194 

195 def replay_trajectory(self, traj): 

196 """Initialize history from old trajectory.""" 

197 if isinstance(traj, str): 

198 from ase.io.trajectory import Trajectory 

199 traj = Trajectory(traj, 'r') 

200 r0 = None 

201 f0 = None 

202 # The last element is not added, as we get that for free when taking 

203 # the first qn-step after the replay 

204 for i in range(0, len(traj) - 1): 

205 r = traj[i].get_positions() 

206 f = traj[i].get_forces() 

207 self.update(r, f, r0, f0) 

208 r0 = r.copy() 

209 f0 = f.copy() 

210 self.iteration += 1 

211 self.r0 = r0 

212 self.f0 = f0 

213 

214 def func(self, x): 

215 """Objective function for use of the optimizers""" 

216 self.atoms.set_positions(x.reshape(-1, 3)) 

217 self.function_calls += 1 

218 return self.atoms.get_potential_energy( 

219 force_consistent=self.force_consistent) 

220 

221 def fprime(self, x): 

222 """Gradient of the objective function for use of the optimizers""" 

223 self.atoms.set_positions(x.reshape(-1, 3)) 

224 self.force_calls += 1 

225 # Remember that forces are minus the gradient! 

226 return - self.atoms.get_forces().reshape(-1) 

227 

228 def line_search(self, r, g, e): 

229 self.p = self.p.ravel() 

230 p_size = np.sqrt((self.p**2).sum()) 

231 if p_size <= np.sqrt(len(self.atoms) * 1e-10): 

232 self.p /= (p_size / np.sqrt(len(self.atoms) * 1e-10)) 

233 g = g.ravel() 

234 r = r.ravel() 

235 ls = LineSearch() 

236 self.alpha_k, e, self.e0, self.no_update = \ 

237 ls._line_search(self.func, self.fprime, r, self.p, g, e, self.e0, 

238 maxstep=self.maxstep, c1=.23, 

239 c2=.46, stpmax=50.) 

240 if self.alpha_k is None: 

241 raise RuntimeError('LineSearch failed!') 

242 

243 

244class LBFGSLineSearch(LBFGS): 

245 """This optimizer uses the LBFGS algorithm, but does a line search that 

246 fulfills the Wolff conditions. 

247 """ 

248 

249 def __init__(self, *args, **kwargs): 

250 kwargs['use_line_search'] = True 

251 LBFGS.__init__(self, *args, **kwargs) 

252 

253# """Modified version of LBFGS. 

254# 

255# This optimizer uses the LBFGS algorithm, but does a line search for the 

256# minimum along the search direction. This is done by issuing an additional 

257# force call for each step, thus doubling the number of calculations. 

258# 

259# Additionally the Hessian is reset if the new guess is not sufficiently 

260# better than the old one. 

261# """ 

262# def __init__(self, *args, **kwargs): 

263# self.dR = kwargs.pop('dR', 0.1) 

264# LBFGS.__init__(self, *args, **kwargs) 

265# 

266# def update(self, r, f, r0, f0): 

267# """Update everything that is kept in memory 

268# 

269# This function is mostly here to allow for replay_trajectory. 

270# """ 

271# if self.iteration > 0: 

272# a1 = abs(np.dot(f.reshape(-1), f0.reshape(-1))) 

273# a2 = np.dot(f0.reshape(-1), f0.reshape(-1)) 

274# if not (a1 <= 0.5 * a2 and a2 != 0): 

275# # Reset optimization 

276# self.initialize() 

277# 

278# # Note that the reset above will set self.iteration to 0 again 

279# # which is why we should check again 

280# if self.iteration > 0: 

281# s0 = r.reshape(-1) - r0.reshape(-1) 

282# self.s.append(s0) 

283# 

284# # We use the gradient which is minus the force! 

285# y0 = f0.reshape(-1) - f.reshape(-1) 

286# self.y.append(y0) 

287# 

288# rho0 = 1.0 / np.dot(y0, s0) 

289# self.rho.append(rho0) 

290# 

291# if self.iteration > self.memory: 

292# self.s.pop(0) 

293# self.y.pop(0) 

294# self.rho.pop(0) 

295# 

296# def determine_step(self, dr): 

297# f = self.atoms.get_forces() 

298# 

299# # Unit-vector along the search direction 

300# du = dr / np.sqrt(np.dot(dr.reshape(-1), dr.reshape(-1))) 

301# 

302# # We keep the old step determination before we figure 

303# # out what is the best to do. 

304# maxstep = self.maxstep * np.sqrt(3 * len(self.atoms)) 

305# 

306# # Finite difference step using temporary point 

307# self.atoms.positions += (du * self.dR) 

308# # Decide how much to move along the line du 

309# Fp1 = np.dot(f.reshape(-1), du.reshape(-1)) 

310# Fp2 = np.dot(self.atoms.get_forces().reshape(-1), du.reshape(-1)) 

311# CR = (Fp1 - Fp2) / self.dR 

312# #RdR = Fp1*0.1 

313# if CR < 0.0: 

314# #print "negcurve" 

315# RdR = maxstep 

316# #if(abs(RdR) > maxstep): 

317# # RdR = self.sign(RdR) * maxstep 

318# else: 

319# Fp = (Fp1 + Fp2) * 0.5 

320# RdR = Fp / CR 

321# if abs(RdR) > maxstep: 

322# RdR = np.sign(RdR) * maxstep 

323# else: 

324# RdR += self.dR * 0.5 

325# return du * RdR