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

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

""" 

This module provides the ConvexHull class. 

""" 

 

import itertools 

import numpy as np 

from typing import List, Union 

from scipy.interpolate import griddata 

from scipy.spatial import ConvexHull as ConvexHullSciPy 

from scipy.spatial.qhull import QhullError 

 

 

class ConvexHull: 

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

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

calculator in SciPy 

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

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

 

Parameters 

---------- 

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

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

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

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

energies : list(float) 

energy (or energy of mixing) for each structure 

 

Attributes 

---------- 

concentrations : np.ndarray 

concentrations of the `N` structures on the convex hull 

energies : np.ndarray 

energies of the `N` structures on the convex hull 

dimensions : int 

number of independent concentrations needed to specify a point in 

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

structures : list(int) 

indices of structures that constitute the convex hull (indices are 

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

initializing the ConvexHull object) 

 

Examples 

-------- 

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

concentrations and energies:: 

 

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

 

after which one can for example plot the data (assuming a 

matplotlib axis object ``ax``):: 

 

ax.plot(hull.concentrations, hull.energies) 

 

or extract structures at or close to the convex hull:: 

 

low_energy_structures = hull.extract_low_energy_structures( 

data['concentration'], data['mixing_energy'], 

energy_tolerance=0.005) 

 

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

<tutorial_enumerate_structures>`. 

 

""" 

 

def __init__(self, 

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

energies: List[float]) -> None: 

assert len(concentrations) == len(energies) 

# Prepare data in format suitable for SciPy-ConvexHull 

concentrations = np.array(concentrations) 

energies = np.array(energies) 

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

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

 

# Construct convex hull 

hull = ConvexHullSciPy(points) 

 

# Collect convex hull points in handy arrays 

concentrations = [] 

energies = [] 

for vertex in hull.vertices: 

if self.dimensions == 1: 

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

else: 

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

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

concentrations = np.array(concentrations) 

energies = np.array(energies) 

 

structures = hull.vertices 

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

# according to it 

if self.dimensions == 1: 

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

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

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

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

else: 

self.concentrations = concentrations 

self.energies = energies 

self.structures = structures 

 

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

self._remove_points_above_tie_plane() 

 

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

""" 

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

than minimum energy. 

 

Parameters 

---------- 

tol 

Tolerance for what energy constitutes a lower one. 

""" 

 

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

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

# pure components. 

if self.dimensions == 1: 

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

# the indices of the lowest and highest concentrations. 

vertices = [] 

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

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

vertices = np.array(vertices) 

else: 

concentration_hull = ConvexHullSciPy(self.concentrations) 

vertices = concentration_hull.vertices 

 

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

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

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

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

# rather than minimize energy. 

to_delete = [] 

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

# The points on the convex concentration hull should always be 

# included, so skip them. 

if i in vertices: 

continue 

 

# The energy obtained as a linear combination of concentrations on 

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

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

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

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

# over all combinations of N vertices. 

for plane in itertools.combinations(vertices, 

min(len(vertices), 

self.dimensions + 1)): 

# Calculate energy that would be gotten with pure components 

# with ascribed concentration. 

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

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

concentration, 

method='linear') 

 

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

# NaN if the concentration is outside the triangle formed by 

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

# then False, which is what we want. 

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

to_delete.append(i) 

break 

 

# Finally remove all points 

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

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

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

 

def get_energy_at_convex_hull(self, target_concentrations: 

Union[List[float], 

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

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

If any concentration is outside the allowed range, NaN is 

returned. 

 

Parameters 

---------- 

target_concentrations 

concentrations at target points 

 

If there is one independent concentration, a list of 

floats is sufficient. Otherwise, the concentrations ought 

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

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

""" 

if self.dimensions > 1: 

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

 

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

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

# but in higher. 

hull_candidate_energies = [] 

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

min(len(self.energies), 

self.dimensions + 1)): 

try: 

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

self.energies[list(plane)], 

np.array(target_concentrations), 

method='linear') 

except QhullError: 

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

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

continue 

hull_candidate_energies.append(plane_energies) 

 

# Pick out the lowest energies found 

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

return hull_energies 

 

def extract_low_energy_structures(self, concentrations: 

Union[List[float], 

List[List[float]]], 

energies: List[float], 

energy_tolerance: float) -> List[int]: 

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

tolerance of the convex hull. 

 

Parameters 

---------- 

concentrations 

concentrations of candidate structures 

 

If there is one independent concentration, a list of 

floats is sufficient. Otherwise, the concentrations must 

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

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

energies 

energies of candidate structures 

energy_tolerance 

include structures with an energy that is at most this far 

from the convex hull 

""" 

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

# they are Pandas Series with "gaps" 

concentrations = np.array(concentrations) 

energies = np.array(energies) 

 

n_points = len(concentrations) 

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

raise ValueError('concentrations and energies must have ' 

'the same length') 

 

# Calculate energy at convex hull for specified concentrations 

hull_energies = self.get_energy_at_convex_hull(concentrations) 

 

# Extract those that are close enough 

close_to_hull = [i for i in range(n_points) 

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

 

return close_to_hull