Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1""" 

2This module provides the ConvexHull class. 

3""" 

4 

5import itertools 

6import numpy as np 

7from typing import List, Sized, Union 

8from scipy.interpolate import griddata 

9from scipy.spatial import ConvexHull as ConvexHullSciPy 

10from scipy.spatial.qhull import QhullError 

11 

12 

13class ConvexHull: 

14 """This class provides functionality for extracting the convex hull 

15 of the (free) energy of mixing. It is based on the `convex hull 

16 calculator in SciPy 

17 <http://docs.scipy.org/doc/scipy-dev/reference/\ 

18generated/scipy.spatial.ConvexHull.html>`_. 

19 

20 Parameters 

21 ---------- 

22 concentrations : list(float) or list(list(float)) 

23 concentrations for each structure listed as ``[[c1, c2], [c1, c2], 

24 ...]``; for binaries, in which case there is only one independent 

25 concentration, the format ``[c1, c2, c3, ...]`` works as well. 

26 energies : list(float) 

27 energy (or energy of mixing) for each structure 

28 

29 Attributes 

30 ---------- 

31 concentrations : np.ndarray 

32 concentrations of the `N` structures on the convex hull 

33 energies : np.ndarray 

34 energies of the `N` structures on the convex hull 

35 dimensions : int 

36 number of independent concentrations needed to specify a point in 

37 concentration space (1 for binaries, 2 for ternaries etc.) 

38 structures : list(int) 

39 indices of structures that constitute the convex hull (indices are 

40 defined by the order of their concentrations and energies are fed when 

41 initializing the ConvexHull object) 

42 

43 Examples 

44 -------- 

45 A `ConvexHull` object is easily initialized by providing lists of 

46 concentrations and energies:: 

47 

48 >>> data = {'concentration': [0, 0.2, 0.2, 0.3, 0.4, 0.5, 0.8, 1.0], 

49 ... 'mixing_energy': [0.1, -0.2, -0.1, -0.2, 0.2, -0.4, -0.2, -0.1]} 

50 >>> hull = ConvexHull(data['concentration'], data['mixing_energy']) 

51 

52 Now one can for example access the points along the convex hull directly:: 

53 

54 >>> for c, e in zip(hull.concentrations, hull.energies): 

55 ... print(c, e) 

56 0.0 0.1 

57 0.2 -0.2 

58 0.5 -0.4 

59 1.0 -0.1 

60 

61 or plot the convex hull along with the original data using e.g., matplotlib:: 

62 

63 >>> import matplotlib.pyplot as plt 

64 >>> plt.scatter(data['concentration'], data['mixing_energy'], color='darkred') 

65 >>> plt.plot(hull.concentrations, hull.energies) 

66 >>> plt.show() 

67 

68 It is also possible to extract structures at or close to the convex hull:: 

69 

70 >>> low_energy_structures = hull.extract_low_energy_structures( 

71 ... data['concentration'], data['mixing_energy'], 

72 ... energy_tolerance=0.005) 

73 

74 A complete example can be found in the :ref:`basic tutorial 

75 <tutorial_enumerate_structures>`. 

76 """ 

77 

78 def __init__(self, 

79 concentrations: Union[List[float], List[List[float]]], 

80 energies: List[float]) -> None: 

81 assert len(concentrations) == len(energies) 

82 # Prepare data in format suitable for SciPy-ConvexHull 

83 concentrations = np.array(concentrations) 

84 energies = np.array(energies) 

85 points = np.column_stack((concentrations, energies)) 

86 self.dimensions = len(points[0]) - 1 

87 

88 # Construct convex hull 

89 hull = ConvexHullSciPy(points) 

90 

91 # Collect convex hull points in handy arrays 

92 concentrations = [] # type: ignore 

93 energies = [] # type: ignore 

94 for vertex in hull.vertices: 

95 if self.dimensions == 1: 

96 concentrations.append(points[vertex][0]) 

97 else: 

98 concentrations.append(points[vertex][0:-1]) 

99 energies.append(points[vertex][-1]) 

100 concentrations = np.array(concentrations) 

101 energies = np.array(energies) 

102 

103 structures = hull.vertices 

104 # If there is just one independent concentration, we'd better sort 

105 # according to it 

106 if self.dimensions == 1: 

107 ces = list(zip(*sorted(zip(concentrations, energies, structures)))) 

108 self.concentrations = np.array(ces[0]) 

109 self.energies = np.array(ces[1]) 

110 self.structures = np.array(ces[2]) 

111 else: 

112 self.concentrations = concentrations 

113 self.energies = energies 

114 self.structures = structures 

115 

116 # Remove points that are above the "pure components plane" 

117 self._remove_points_above_tie_plane() 

118 

119 def _remove_points_above_tie_plane(self, tol: float = 1e-6) -> None: 

120 """ 

121 Remove all points on the convex hull that correspond to maximum rather 

122 than minimum energy. 

123 

124 Parameters 

125 ---------- 

126 tol 

127 Tolerance for what energy constitutes a lower one. 

128 """ 

129 

130 # Identify the "complex concentration hull", i.e. the extremal 

131 # concentrations. In the simplest case, these should simply be the 

132 # pure components. 

133 if self.dimensions == 1: 

134 # Then the ConvexHullScipy function doesn't work, so we just pick 

135 # the indices of the lowest and highest concentrations. 

136 vertices = [] 

137 vertices.append(np.argmin(self.concentrations)) 

138 vertices.append(np.argmax(self.concentrations)) 

139 vertices = np.array(vertices) 

140 else: 

141 concentration_hull = ConvexHullSciPy(self.concentrations) 

142 vertices = concentration_hull.vertices 

143 

144 # Remove all points of the convex energy hull that have an energy that 

145 # is higher than what would be gotten with pure components at the same 

146 # concentration. These points are mathematically on the convex hull, 

147 # but in the physically uninteresting upper part, i.e. they maximize 

148 # rather than minimize energy. 

149 to_delete = [] 

150 for i, concentration in enumerate(self.concentrations): 

151 # The points on the convex concentration hull should always be 

152 # included, so skip them. 

153 if i in vertices: 

154 continue 

155 

156 # The energy obtained as a linear combination of concentrations on 

157 # the convex hull is the "z coordinate" of the position on a 

158 # (hyper)plane in the (number of independent concentrations + 

159 # 1)-dimensional (N-D) space. This plane is spanned by N points. 

160 # If there are more vertices on the convex hull, we need to loop 

161 # over all combinations of N vertices. 

162 for plane in itertools.combinations(vertices, 

163 min(len(vertices), 

164 self.dimensions + 1)): 

165 # Calculate energy that would be gotten with pure components 

166 # with ascribed concentration. 

167 energy_pure = griddata(self.concentrations[np.array(plane)], 

168 self.energies[np.array(plane)], 

169 concentration, 

170 method='linear') 

171 

172 # Prepare to delete if the energy was lowered. `griddata` gives 

173 # NaN if the concentration is outside the triangle formed by 

174 # the three vertices. The result of the below comparison is 

175 # then False, which is what we want. 

176 if energy_pure < self.energies[i] - tol: 

177 to_delete.append(i) 

178 break 

179 

180 # Finally remove all points 

181 self.concentrations = np.delete(self.concentrations, to_delete, 0) 

182 self.energies = np.delete(self.energies, to_delete, 0) 

183 self.structures = list(np.delete(self.structures, to_delete, 0)) 

184 

185 def get_energy_at_convex_hull(self, target_concentrations: 

186 Union[List[float], 

187 List[List[float]]]) -> np.ndarray: 

188 """Returns the energy of the convex hull at specified concentrations. 

189 If any concentration is outside the allowed range, NaN is 

190 returned. 

191 

192 Parameters 

193 ---------- 

194 target_concentrations 

195 concentrations at target points 

196 

197 If there is one independent concentration, a list of 

198 floats is sufficient. Otherwise, the concentrations ought 

199 to be provided as a list of lists, such as ``[[0.1, 0.2], 

200 [0.3, 0.1], ...]``. 

201 """ 

202 if self.dimensions > 1 and isinstance(target_concentrations[0], Sized): 

203 assert len(target_concentrations[0]) == self.dimensions 

204 

205 # Loop over all complexes of N+1 points to make sure that the lowest 

206 # energy plane is used in the end. This is needed in two dimensions 

207 # but in higher. 

208 hull_candidate_energies = [] 

209 for plane in itertools.combinations(range(len(self.energies)), 

210 min(len(self.energies), 

211 self.dimensions + 1)): 

212 try: 

213 plane_energies = griddata(self.concentrations[list(plane)], 

214 self.energies[list(plane)], 

215 np.array(target_concentrations), 

216 method='linear') 

217 except QhullError: 

218 # If the points lie on a line, the convex hull will fail, but 

219 # we do not need to care about these "planes" anyway 

220 continue 

221 hull_candidate_energies.append(plane_energies) 

222 

223 # Pick out the lowest energies found 

224 hull_energies = np.nanmin(hull_candidate_energies, axis=0) 

225 return hull_energies 

226 

227 def extract_low_energy_structures(self, concentrations: 

228 Union[List[float], 

229 List[List[float]]], 

230 energies: List[float], 

231 energy_tolerance: float) -> List[int]: 

232 """Returns the indices of energies that lie within a certain 

233 tolerance of the convex hull. 

234 

235 Parameters 

236 ---------- 

237 concentrations 

238 concentrations of candidate structures 

239 

240 If there is one independent concentration, a list of 

241 floats is sufficient. Otherwise, the concentrations must 

242 be provided as a list of lists, such as ``[[0.1, 0.2], 

243 [0.3, 0.1], ...]``. 

244 energies 

245 energies of candidate structures 

246 energy_tolerance 

247 include structures with an energy that is at most this far 

248 from the convex hull 

249 """ 

250 # Convert to numpy arrays, can be necessary if, for example, 

251 # they are Pandas Series with "gaps" 

252 concentrations = np.array(concentrations) 

253 energies = np.array(energies) 

254 

255 n_points = len(concentrations) 

256 if len(energies) != n_points: 256 ↛ 257line 256 didn't jump to line 257, because the condition on line 256 was never true

257 raise ValueError('concentrations and energies must have ' 

258 'the same length') 

259 

260 # Calculate energy at convex hull for specified concentrations 

261 hull_energies = self.get_energy_at_convex_hull(concentrations) 

262 

263 # Extract those that are close enough 

264 close_to_hull = [i for i in range(n_points) 

265 if energies[i] <= hull_energies[i] + energy_tolerance] 

266 

267 return close_to_hull