1import numpy as np 

2 

3from ase.optimize.optimize import Dynamics 

4from ase.optimize.fire import FIRE 

5from ase.units import kB 

6from ase.parallel import world 

7from ase.io.trajectory import Trajectory 

8 

9 

10class BasinHopping(Dynamics): 

11 """Basin hopping algorithm. 

12 

13 After Wales and Doye, J. Phys. Chem. A, vol 101 (1997) 5111-5116 

14 

15 and 

16 

17 David J. Wales and Harold A. Scheraga, Science, Vol. 285, 1368 (1999) 

18 """ 

19 

20 def __init__(self, atoms, 

21 temperature=100 * kB, 

22 optimizer=FIRE, 

23 fmax=0.1, 

24 dr=0.1, 

25 logfile='-', 

26 trajectory='lowest.traj', 

27 optimizer_logfile='-', 

28 local_minima_trajectory='local_minima.traj', 

29 adjust_cm=True): 

30 """Parameters: 

31 

32 atoms: Atoms object 

33 The Atoms object to operate on. 

34 

35 trajectory: string 

36 Pickle file used to store trajectory of atomic movement. 

37 

38 logfile: file object or str 

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

40 Use '-' for stdout. 

41 """ 

42 self.kT = temperature 

43 self.optimizer = optimizer 

44 self.fmax = fmax 

45 self.dr = dr 

46 if adjust_cm: 

47 self.cm = atoms.get_center_of_mass() 

48 else: 

49 self.cm = None 

50 

51 self.optimizer_logfile = optimizer_logfile 

52 self.lm_trajectory = local_minima_trajectory 

53 if isinstance(local_minima_trajectory, str): 

54 self.lm_trajectory = self.closelater( 

55 Trajectory(local_minima_trajectory, 'w', atoms)) 

56 

57 Dynamics.__init__(self, atoms, logfile, trajectory) 

58 self.initialize() 

59 

60 def todict(self): 

61 d = {'type': 'optimization', 

62 'optimizer': self.__class__.__name__, 

63 'local-minima-optimizer': self.optimizer.__name__, 

64 'temperature': self.kT, 

65 'max-force': self.fmax, 

66 'maximal-step-width': self.dr} 

67 return d 

68 

69 def initialize(self): 

70 self.positions = 0.0 * self.atoms.get_positions() 

71 self.Emin = self.get_energy(self.atoms.get_positions()) or 1.e32 

72 self.rmin = self.atoms.get_positions() 

73 self.positions = self.atoms.get_positions() 

74 self.call_observers() 

75 self.log(-1, self.Emin, self.Emin) 

76 

77 def run(self, steps): 

78 """Hop the basins for defined number of steps.""" 

79 

80 ro = self.positions 

81 Eo = self.get_energy(ro) 

82 

83 for step in range(steps): 

84 En = None 

85 while En is None: 

86 rn = self.move(ro) 

87 En = self.get_energy(rn) 

88 

89 if En < self.Emin: 

90 # new minimum found 

91 self.Emin = En 

92 self.rmin = self.atoms.get_positions() 

93 self.call_observers() 

94 self.log(step, En, self.Emin) 

95 

96 accept = np.exp((Eo - En) / self.kT) > np.random.uniform() 

97 if accept: 

98 ro = rn.copy() 

99 Eo = En 

100 

101 def log(self, step, En, Emin): 

102 if self.logfile is None: 

103 return 

104 name = self.__class__.__name__ 

105 self.logfile.write('%s: step %d, energy %15.6f, emin %15.6f\n' 

106 % (name, step, En, Emin)) 

107 self.logfile.flush() 

108 

109 def move(self, ro): 

110 """Move atoms by a random step.""" 

111 atoms = self.atoms 

112 # displace coordinates 

113 disp = np.random.uniform(-1., 1., (len(atoms), 3)) 

114 rn = ro + self.dr * disp 

115 atoms.set_positions(rn) 

116 if self.cm is not None: 

117 cm = atoms.get_center_of_mass() 

118 atoms.translate(self.cm - cm) 

119 rn = atoms.get_positions() 

120 world.broadcast(rn, 0) 

121 atoms.set_positions(rn) 

122 return atoms.get_positions() 

123 

124 def get_minimum(self): 

125 """Return minimal energy and configuration.""" 

126 atoms = self.atoms.copy() 

127 atoms.set_positions(self.rmin) 

128 return self.Emin, atoms 

129 

130 def get_energy(self, positions): 

131 """Return the energy of the nearest local minimum.""" 

132 if np.sometrue(self.positions != positions): 

133 self.positions = positions 

134 self.atoms.set_positions(positions) 

135 

136 with self.optimizer(self.atoms, 

137 logfile=self.optimizer_logfile) as opt: 

138 opt.run(fmax=self.fmax) 

139 if self.lm_trajectory is not None: 

140 self.lm_trajectory.write(self.atoms) 

141 

142 self.energy = self.atoms.get_potential_energy() 

143 

144 return self.energy