 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"""

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

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>`_.

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

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)

43 Examples

44 --------

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

46 concentrations and energies::

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'])

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

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

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

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()

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

70 >>> low_energy_structures = hull.extract_low_energy_structures(

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

72 ... energy_tolerance=0.005)

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

75 <tutorial_enumerate_structures>`.

76 """

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) - 1

88 # Construct convex hull

89 hull = ConvexHullSciPy(points)

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])

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)

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)

109 self.energies = np.array(ces)

110 self.structures = np.array(ces)

111 else:

112 self.concentrations = concentrations

113 self.energies = energies

114 self.structures = structures

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

117 self._remove_points_above_tie_plane()

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.

124 Parameters

125 ----------

126 tol

127 Tolerance for what energy constitutes a lower one.

128 """

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

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

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')

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

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))

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.

192 Parameters

193 ----------

194 target_concentrations

195 concentrations at target points

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, Sized):

203 assert len(target_concentrations) == self.dimensions

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)

223 # Pick out the lowest energies found

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

225 return hull_energies

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.

235 Parameters

236 ----------

237 concentrations

238 concentrations of candidate structures

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)

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')

260 # Calculate energy at convex hull for specified concentrations

261 hull_energies = self.get_energy_at_convex_hull(concentrations)

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]

267 return close_to_hull