1import warnings 

2import numpy as np 

3 

4from ase.optimize.optimize import Optimizer 

5 

6 

7class FIRE(Optimizer): 

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

9 dt=0.1, maxstep=None, maxmove=None, dtmax=1.0, Nmin=5, 

10 finc=1.1, fdec=0.5, 

11 astart=0.1, fa=0.99, a=0.1, master=None, downhill_check=False, 

12 position_reset_callback=None, force_consistent=None): 

13 """Parameters: 

14 

15 atoms: Atoms object 

16 The Atoms object to relax. 

17 

18 restart: string 

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

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

21 be used, if the file exists. 

22 

23 trajectory: string 

24 Pickle file used to store trajectory of atomic movement. 

25 

26 logfile: file object or str 

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

28 Use '-' for stdout. 

29 

30 master: boolean 

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

32 set to true, this rank will save files. 

33 

34 downhill_check: boolean 

35 Downhill check directly compares potential energies of subsequent 

36 steps of the FIRE algorithm rather than relying on the current 

37 product v*f that is positive if the FIRE dynamics moves downhill. 

38 This can detect numerical issues where at large time steps the step 

39 is uphill in energy even though locally v*f is positive, i.e. the 

40 algorithm jumps over a valley because of a too large time step. 

41 

42 position_reset_callback: function(atoms, r, e, e_last) 

43 Function that takes current *atoms* object, an array of position 

44 *r* that the optimizer will revert to, current energy *e* and 

45 energy of last step *e_last*. This is only called if e > e_last. 

46 

47 force_consistent: boolean or None 

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

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

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

51 falls back to force_consistent=False if not. Only meaningful 

52 when downhill_check is True. 

53 """ 

54 Optimizer.__init__(self, atoms, restart, logfile, trajectory, 

55 master, force_consistent=force_consistent) 

56 

57 self.dt = dt 

58 

59 self.Nsteps = 0 

60 

61 if maxstep is not None: 

62 self.maxstep = maxstep 

63 elif maxmove is not None: 

64 self.maxstep = maxmove 

65 warnings.warn('maxmove is deprecated; please use maxstep', 

66 np.VisibleDeprecationWarning) 

67 else: 

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

69 

70 self.dtmax = dtmax 

71 self.Nmin = Nmin 

72 self.finc = finc 

73 self.fdec = fdec 

74 self.astart = astart 

75 self.fa = fa 

76 self.a = a 

77 self.downhill_check = downhill_check 

78 self.position_reset_callback = position_reset_callback 

79 

80 def initialize(self): 

81 self.v = None 

82 

83 def read(self): 

84 self.v, self.dt = self.load() 

85 

86 def step(self, f=None): 

87 atoms = self.atoms 

88 

89 if f is None: 

90 f = atoms.get_forces() 

91 

92 if self.v is None: 

93 self.v = np.zeros((len(atoms), 3)) 

94 if self.downhill_check: 

95 self.e_last = atoms.get_potential_energy( 

96 force_consistent=self.force_consistent) 

97 self.r_last = atoms.get_positions().copy() 

98 self.v_last = self.v.copy() 

99 else: 

100 is_uphill = False 

101 if self.downhill_check: 

102 e = atoms.get_potential_energy( 

103 force_consistent=self.force_consistent) 

104 # Check if the energy actually decreased 

105 if e > self.e_last: 

106 # If not, reset to old positions... 

107 if self.position_reset_callback is not None: 

108 self.position_reset_callback(atoms, self.r_last, e, 

109 self.e_last) 

110 atoms.set_positions(self.r_last) 

111 is_uphill = True 

112 self.e_last = atoms.get_potential_energy( 

113 force_consistent=self.force_consistent) 

114 self.r_last = atoms.get_positions().copy() 

115 self.v_last = self.v.copy() 

116 

117 vf = np.vdot(f, self.v) 

118 if vf > 0.0 and not is_uphill: 

119 self.v = (1.0 - self.a) * self.v + self.a * f / np.sqrt( 

120 np.vdot(f, f)) * np.sqrt(np.vdot(self.v, self.v)) 

121 if self.Nsteps > self.Nmin: 

122 self.dt = min(self.dt * self.finc, self.dtmax) 

123 self.a *= self.fa 

124 self.Nsteps += 1 

125 else: 

126 self.v[:] *= 0.0 

127 self.a = self.astart 

128 self.dt *= self.fdec 

129 self.Nsteps = 0 

130 

131 self.v += self.dt * f 

132 dr = self.dt * self.v 

133 normdr = np.sqrt(np.vdot(dr, dr)) 

134 if normdr > self.maxstep: 

135 dr = self.maxstep * dr / normdr 

136 r = atoms.get_positions() 

137 atoms.set_positions(r + dr) 

138 self.dump((self.v, self.dt))