1import warnings 

2 

3import numpy as np 

4from numpy.linalg import eigh 

5 

6from ase.optimize.optimize import Optimizer 

7 

8 

9class BFGS(Optimizer): 

10 # default parameters 

11 defaults = {**Optimizer.defaults, 'alpha': 70.0} 

12 

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

14 maxstep=None, master=None, alpha=None): 

15 """BFGS optimizer. 

16 

17 Parameters: 

18 

19 atoms: Atoms object 

20 The Atoms object to relax. 

21 

22 restart: string 

23 Pickle file used to store hessian matrix. If set, file with 

24 such a name will be searched and hessian matrix stored will 

25 be used, if the file exists. 

26 

27 trajectory: string 

28 Pickle file used to store trajectory of atomic movement. 

29 

30 logfile: file object or str 

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

32 Use '-' for stdout. 

33 

34 maxstep: float 

35 Used to set the maximum distance an atom can move per 

36 iteration (default value is 0.2 Å). 

37 

38 master: boolean 

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

40 set to true, this rank will save files. 

41 

42 alpha: float 

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

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

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

46 a lower value also means risk of instability. 

47 """ 

48 if maxstep is None: 

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

50 else: 

51 self.maxstep = maxstep 

52 

53 if self.maxstep > 1.0: 

54 warnings.warn('You are using a *very* large value for ' 

55 'the maximum step size: %.1f Å' % maxstep) 

56 

57 if alpha is None: 

58 self.alpha = self.defaults['alpha'] 

59 else: 

60 self.alpha = alpha 

61 

62 Optimizer.__init__(self, atoms, restart, logfile, trajectory, master) 

63 

64 def todict(self): 

65 d = Optimizer.todict(self) 

66 if hasattr(self, 'maxstep'): 

67 d.update(maxstep=self.maxstep) 

68 return d 

69 

70 def initialize(self): 

71 # initial hessian 

72 self.H0 = np.eye(3 * len(self.atoms)) * self.alpha 

73 

74 self.H = None 

75 self.r0 = None 

76 self.f0 = None 

77 

78 def read(self): 

79 self.H, self.r0, self.f0, self.maxstep = self.load() 

80 

81 def step(self, f=None): 

82 atoms = self.atoms 

83 

84 if f is None: 

85 f = atoms.get_forces() 

86 

87 r = atoms.get_positions() 

88 f = f.reshape(-1) 

89 self.update(r.flat, f, self.r0, self.f0) 

90 omega, V = eigh(self.H) 

91 

92 # FUTURE: Log this properly 

93 # # check for negative eigenvalues of the hessian 

94 # if any(omega < 0): 

95 # n_negative = len(omega[omega < 0]) 

96 # msg = '\n** BFGS Hessian has {} negative eigenvalues.'.format( 

97 # n_negative 

98 # ) 

99 # print(msg, flush=True) 

100 # if self.logfile is not None: 

101 # self.logfile.write(msg) 

102 # self.logfile.flush() 

103 

104 dr = np.dot(V, np.dot(f, V) / np.fabs(omega)).reshape((-1, 3)) 

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

106 dr = self.determine_step(dr, steplengths) 

107 atoms.set_positions(r + dr) 

108 self.r0 = r.flat.copy() 

109 self.f0 = f.copy() 

110 self.dump((self.H, self.r0, self.f0, self.maxstep)) 

111 

112 def determine_step(self, dr, steplengths): 

113 """Determine step to take according to maxstep 

114 

115 Normalize all steps as the largest step. This way 

116 we still move along the eigendirection. 

117 """ 

118 maxsteplength = np.max(steplengths) 

119 if maxsteplength >= self.maxstep: 

120 scale = self.maxstep / maxsteplength 

121 # FUTURE: Log this properly 

122 # msg = '\n** scale step by {:.3f} to be shorter than {}'.format( 

123 # scale, self.maxstep 

124 # ) 

125 # print(msg, flush=True) 

126 

127 dr *= scale 

128 

129 return dr 

130 

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

132 if self.H is None: 

133 self.H = self.H0 

134 return 

135 dr = r - r0 

136 

137 if np.abs(dr).max() < 1e-7: 

138 # Same configuration again (maybe a restart): 

139 return 

140 

141 df = f - f0 

142 a = np.dot(dr, df) 

143 dg = np.dot(self.H, dr) 

144 b = np.dot(dr, dg) 

145 self.H -= np.outer(df, df) / a + np.outer(dg, dg) / b 

146 

147 def replay_trajectory(self, traj): 

148 """Initialize hessian from old trajectory.""" 

149 if isinstance(traj, str): 

150 from ase.io.trajectory import Trajectory 

151 traj = Trajectory(traj, 'r') 

152 self.H = None 

153 atoms = traj[0] 

154 r0 = atoms.get_positions().ravel() 

155 f0 = atoms.get_forces().ravel() 

156 for atoms in traj: 

157 r = atoms.get_positions().ravel() 

158 f = atoms.get_forces().ravel() 

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

160 r0 = r 

161 f0 = f 

162 

163 self.r0 = r0 

164 self.f0 = f0 

165 

166 

167class oldBFGS(BFGS): 

168 def determine_step(self, dr, steplengths): 

169 """Old BFGS behaviour for scaling step lengths 

170 

171 This keeps the behaviour of truncating individual steps. Some might 

172 depend of this as some absurd kind of stimulated annealing to find the 

173 global minimum. 

174 """ 

175 dr /= np.maximum(steplengths / self.maxstep, 1.0).reshape(-1, 1) 

176 return dr