Coverage for /builds/kinetik161/ase/ase/optimize/basin.py: 97.67%

1from typing import IO, Type, Union 

2 

3import numpy as np 

4 

5from ase import Atoms, units 

6from ase.io.trajectory import Trajectory 

7from ase.optimize.fire import FIRE 

8from ase.optimize.optimize import Dynamics, Optimizer 

9from ase.parallel import world 

10 

11 

12class BasinHopping(Dynamics): 

13 """Basin hopping algorithm. 

14 

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

16 

17 and 

18 

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

20 """ 

21 

22 def __init__( 

23 self, 

24 atoms: Atoms, 

25 temperature: float = 100 * units.kB, 

26 optimizer: Type[Optimizer] = FIRE, 

27 fmax: float = 0.1, 

28 dr: float = 0.1, 

29 logfile: Union[IO, str] = '-', 

30 trajectory: str = 'lowest.traj', 

31 optimizer_logfile: str = '-', 

32 local_minima_trajectory: str = 'local_minima.traj', 

33 adjust_cm: bool = True, 

34 ): 

35 """Parameters: 

36 

37 atoms: Atoms object 

38 The Atoms object to operate on. 

39 

40 trajectory: string 

41 Pickle file used to store trajectory of atomic movement. 

42 

43 logfile: file object or str 

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

45 Use '-' for stdout. 

46 """ 

47 self.kT = temperature 

48 self.optimizer = optimizer 

49 self.fmax = fmax 

50 self.dr = dr 

51 if adjust_cm: 

52 self.cm = atoms.get_center_of_mass() 

53 else: 

54 self.cm = None 

55 

56 self.optimizer_logfile = optimizer_logfile 

57 self.lm_trajectory = local_minima_trajectory 

58 if isinstance(local_minima_trajectory, str): 

59 self.lm_trajectory = self.closelater( 

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

61 

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

63 self.initialize() 

64 

65 def todict(self): 

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

67 'optimizer': self.__class__.__name__, 

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

69 'temperature': self.kT, 

70 'max-force': self.fmax, 

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

72 return d 

73 

74 def initialize(self): 

75 positions = self.optimizable.get_positions() 

76 self.positions = np.zeros_like(positions) 

77 self.Emin = self.get_energy(positions) or 1.e32 

78 self.rmin = self.optimizable.get_positions() 

79 self.positions = self.optimizable.get_positions() 

80 self.call_observers() 

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

82 

83 def run(self, steps): 

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

85 

86 ro = self.positions 

87 Eo = self.get_energy(ro) 

88 

89 for step in range(steps): 

90 En = None 

91 while En is None: 

92 rn = self.move(ro) 

93 En = self.get_energy(rn) 

94 

95 if En < self.Emin: 

96 # new minimum found 

97 self.Emin = En 

98 self.rmin = self.optimizable.get_positions() 

99 self.call_observers() 

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

101 

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

103 if accept: 

104 ro = rn.copy() 

105 Eo = En 

106 

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

108 if self.logfile is None: 

109 return 

110 name = self.__class__.__name__ 

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

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

113 self.logfile.flush() 

114 

115 def _atoms(self): 

116 from ase.optimize.optimize import OptimizableAtoms 

117 assert isinstance(self.optimizable, OptimizableAtoms) 

118 # Some parts of the basin code cannot work on Filter objects. 

119 # They evidently need an actual Atoms object - at least until 

120 # someone changes the code so it doesn't need that. 

121 return self.optimizable.atoms 

122 

123 def move(self, ro): 

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

125 atoms = self._atoms() 

126 # displace coordinates 

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

128 rn = ro + self.dr * disp 

129 atoms.set_positions(rn) 

130 if self.cm is not None: 

131 cm = atoms.get_center_of_mass() 

132 atoms.translate(self.cm - cm) 

133 rn = atoms.get_positions() 

134 world.broadcast(rn, 0) 

135 atoms.set_positions(rn) 

136 return atoms.get_positions() 

137 

138 def get_minimum(self): 

139 """Return minimal energy and configuration.""" 

140 atoms = self._atoms().copy() 

141 atoms.set_positions(self.rmin) 

142 return self.Emin, atoms 

143 

144 def get_energy(self, positions): 

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

146 if np.any(self.positions != positions): 

147 self.positions = positions 

148 self.optimizable.set_positions(positions) 

149 

150 with self.optimizer(self.optimizable, 

151 logfile=self.optimizer_logfile) as opt: 

152 opt.run(fmax=self.fmax) 

153 if self.lm_trajectory is not None: 

154 self.lm_trajectory.write(self.optimizable) 

155 

156 self.energy = self.optimizable.get_potential_energy() 

157 

158 return self.energy