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 self.maxstep = maxstep 

49 if self.maxstep is None: 

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

51 

52 if self.maxstep > 1.0: 

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

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

55 

56 self.alpha = alpha 

57 if self.alpha is None: 

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

59 

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

61 

62 def initialize(self): 

63 # initial hessian 

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

65 

66 self.H = None 

67 self.pos0 = None 

68 self.forces0 = None 

69 

70 def read(self): 

71 self.H, self.pos0, self.forces0, self.maxstep = self.load() 

72 

73 def step(self, forces=None): 

74 atoms = self.atoms 

75 

76 if forces is None: 

77 forces = atoms.get_forces() 

78 

79 pos = atoms.get_positions() 

80 dpos, steplengths = self.prepare_step(pos, forces) 

81 dpos = self.determine_step(dpos, steplengths) 

82 atoms.set_positions(pos + dpos) 

83 self.dump((self.H, self.pos0, self.forces0, self.maxstep)) 

84 

85 def prepare_step(self, pos, forces): 

86 forces = forces.reshape(-1) 

87 self.update(pos.flat, forces, self.pos0, self.forces0) 

88 omega, V = eigh(self.H) 

89 

90 # FUTURE: Log this properly 

91 # # check for negative eigenvalues of the hessian 

92 # if any(omega < 0): 

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

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

95 # n_negative 

96 # ) 

97 # print(msg, flush=True) 

98 # if self.logfile is not None: 

99 # self.logfile.write(msg) 

100 # self.logfile.flush() 

101 

102 dpos = np.dot(V, np.dot(forces, V) / np.fabs(omega)).reshape((-1, 3)) 

103 steplengths = (dpos**2).sum(1)**0.5 

104 self.pos0 = pos.flat.copy() 

105 self.forces0 = forces.copy() 

106 return dpos, steplengths 

107 

108 def determine_step(self, dpos, steplengths): 

109 """Determine step to take according to maxstep 

110 

111 Normalize all steps as the largest step. This way 

112 we still move along the eigendirection. 

113 """ 

114 maxsteplength = np.max(steplengths) 

115 if maxsteplength >= self.maxstep: 

116 scale = self.maxstep / maxsteplength 

117 # FUTURE: Log this properly 

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

119 # scale, self.maxstep 

120 # ) 

121 # print(msg, flush=True) 

122 

123 dpos *= scale 

124 return dpos 

125 

126 def update(self, pos, forces, pos0, forces0): 

127 if self.H is None: 

128 self.H = self.H0 

129 return 

130 dpos = pos - pos0 

131 

132 if np.abs(dpos).max() < 1e-7: 

133 # Same configuration again (maybe a restart): 

134 return 

135 

136 dforces = forces - forces0 

137 a = np.dot(dpos, dforces) 

138 dg = np.dot(self.H, dpos) 

139 b = np.dot(dpos, dg) 

140 self.H -= np.outer(dforces, dforces) / a + np.outer(dg, dg) / b 

141 

142 def replay_trajectory(self, traj): 

143 """Initialize hessian from old trajectory.""" 

144 if isinstance(traj, str): 

145 from ase.io.trajectory import Trajectory 

146 traj = Trajectory(traj, 'r') 

147 self.H = None 

148 atoms = traj[0] 

149 pos0 = atoms.get_positions().ravel() 

150 forces0 = atoms.get_forces().ravel() 

151 for atoms in traj: 

152 pos = atoms.get_positions().ravel() 

153 forces = atoms.get_forces().ravel() 

154 self.update(pos, forces, pos0, forces0) 

155 pos0 = pos 

156 forces0 = forces 

157 

158 self.pos0 = pos0 

159 self.forces0 = forces0 

160 

161 

162class oldBFGS(BFGS): 

163 def determine_step(self, dpos, steplengths): 

164 """Old BFGS behaviour for scaling step lengths 

165 

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

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

168 global minimum. 

169 """ 

170 dpos /= np.maximum(steplengths / self.maxstep, 1.0).reshape(-1, 1) 

171 return dpos