/
example_solid_solution.py
361 lines (304 loc) · 14.9 KB
/
example_solid_solution.py
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
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit for the Earth and Planetary Sciences
# Copyright (C) 2012 - 2015 by the BurnMan team, released under the GNU
# GPL v2 or later.
"""
example_solid_solution
----------------------
This example shows how to create different solid solution models and output
thermodynamic and thermoelastic quantities.
There are four main types of solid solution currently implemented in
BurnMan:
1. Ideal solid solutions
2. Symmmetric solid solutions
3. Asymmetric solid solutions
4. Subregular solid solutions
These solid solutions can potentially deal with:
* Disordered endmembers (more than one element on a crystallographic site)
* Site vacancies
* More than one valence/spin state of the same element on a site
*Uses:*
* :doc:`mineral_database`
* :class:`burnman.solidsolution.SolidSolution`
* :class:`burnman.solutionmodel.SolutionModel`
*Demonstrates:*
* Different ways to define a solid solution
* How to set composition and state
* How to output thermodynamic and thermoelastic properties
"""
from __future__ import absolute_import
import numpy as np
import matplotlib.pyplot as plt
import burnman_path # adds the local burnman directory to the path
import burnman
from burnman import minerals
assert burnman_path # silence pyflakes warning
if __name__ == "__main__":
'''
First, let's pick a starting pressure and temperature
'''
P = 1.e5
T = 1000.
"""
We can create an instance of a solid solution in two ways.
We can create a single instance using the SolidSolution constructor
(here's an example of an ideal pyrope-almandine garnet) ...
"""
g1 = burnman.SolidSolution(name = 'Ideal pyrope-almandine garnet',
solution_type = 'ideal',
endmembers = [[minerals.HP_2011_ds62.py(),
'[Mg]3[Al]2Si3O12'],
[minerals.HP_2011_ds62.alm(),
'[Fe]3[Al]2Si3O12']],
molar_fractions = [0.5, 0.5])
g1.set_state(P, T)
gibbs = g1.gibbs
"""
... or we can create a class that derives from solid solution
(so that we can create multiple instances of the same solution)
"""
class mg_fe_garnet(burnman.SolidSolution):
def __init__(self, molar_fractions=None):
self.name = 'Ideal pyrope-almandine garnet'
self.solution_type = 'ideal'
self.endmembers = [[minerals.HP_2011_ds62.py(),
'[Mg]3[Al]2Si3O12'],
[minerals.HP_2011_ds62.alm(),
'[Fe]3[Al]2Si3O12']]
burnman.SolidSolution.__init__(self, molar_fractions=molar_fractions)
"""
Initialisation can optionally include setting the composition of the solid solution, i.e.
"""
g1 = mg_fe_garnet([0.0, 1.0])
g1.set_state(P, T)
alm_gibbs = g1.gibbs
"""
Alternatively, the composition can be set after initialisation
"""
g1.set_composition([1.0, 0.0])
g1.set_state(P, T)
py_gibbs = g1.gibbs
"""
Let's plot some interesting stuff...
"""
comp = np.linspace(0.001, 0.999, 100)
g1_gibbs = np.empty_like(comp)
g1_excess_gibbs = np.empty_like(comp)
g1_pyrope_activity = np.empty_like(comp)
g1_almandine_activity = np.empty_like(comp)
g1_pyrope_gamma = np.empty_like(comp)
g1_almandine_gamma = np.empty_like(comp)
for i, c in enumerate(comp):
molar_fractions = [1.0 - c, c]
g1.set_composition(molar_fractions)
g1.set_state(P, T)
g1_gibbs[i] = g1.gibbs
g1_excess_gibbs[i] = g1.excess_gibbs
g1_pyrope_activity[i] = g1.activities[0]
g1_almandine_activity[i] = g1.activities[1]
g1_pyrope_gamma[i] = g1.activity_coefficients[0]
g1_almandine_gamma[i] = g1.activity_coefficients[1]
plt.plot([0., 1.], [py_gibbs / 1000., alm_gibbs / 1000.],
'b-', linewidth=1., label='Mechanical mixing')
plt.plot(comp, g1_gibbs / 1000., 'r-',
linewidth=1., label='Ideal solid solution')
plt.title("Ideal pyrope-almandine join")
plt.ylabel("Gibbs free energy of solution (kJ/mol)")
plt.xlabel("Almandine fraction")
plt.legend(loc='lower right')
plt.show()
plt.plot(comp, g1_pyrope_activity, 'r-',
linewidth=1., label='Pyrope activity')
plt.plot(comp, g1_almandine_activity, 'b-',
linewidth=1., label='Almandine activity')
plt.plot(comp, g1_pyrope_gamma, 'r-.',
linewidth=1., label='Pyrope activity coefficient')
plt.plot(comp, g1_almandine_gamma, 'b--',
linewidth=1., label='Almandine activity coefficient')
plt.title("Ideal pyrope-almandine join")
plt.ylim(0.0, 1.01)
plt.ylabel("Activities")
plt.xlabel("Almandine fraction")
plt.legend(loc='lower right')
plt.show()
'''
Not included in this example document are ways to create solid solutions with spin transitions,
vacancies and mixed valence states. However, the formula parsing in BurnMan is comprehensive,
so any model should be reproducable. Some formatted formulae examples follow:
- [Fe]O: simple wuestite
- [Fef2/3Vac1/3]O: a theoretical ferric endmember of wuestite. Note the distinct element name for Fe3+.
- [Fe2/3Vac1/3]O: an alternative ferric endmember. Mixing between this and [Fe]O will produce a
different configurational entropy to the previous model because there is no distinction
between Fe2+ and Fe3+.
- [Fels]O: Low spin wuestite. Another example illustrating the free-form approach to element definition.
'''
'''
The solid solution corresponding to the pyrope-almandine join is actually not quite ideal.
It can be well-approximated with a symmetric regular solid solution model
'''
g2 = burnman.SolidSolution(name = 'Symmetric pyrope-almandine garnet',
solution_type = 'symmetric',
endmembers = [[minerals.HP_2011_ds62.py(),
'[Mg]3[Al]2Si3O12'],
[minerals.HP_2011_ds62.alm(),
'[Fe]3[Al]2Si3O12']],
energy_interaction = [[2.5e3]])
g2_excess_gibbs = np.empty_like(comp)
g2_pyrope_activity = np.empty_like(comp)
g2_almandine_activity = np.empty_like(comp)
g2_pyrope_gamma = np.empty_like(comp)
g2_almandine_gamma = np.empty_like(comp)
for i, c in enumerate(comp):
molar_fractions = [1. - c, c]
g2.set_composition(molar_fractions)
g2.set_state(P, T)
g2_excess_gibbs[i] = g2.excess_gibbs
g2_pyrope_activity[i] = g2.activities[0]
g2_almandine_activity[i] = g2.activities[1]
g2_pyrope_gamma[i] = g2.activity_coefficients[0]
g2_almandine_gamma[i] = g2.activity_coefficients[1]
plt.plot(comp, g1_excess_gibbs, 'r-', linewidth=1., label='Ideal solution')
plt.plot(comp, g2_excess_gibbs, 'g-', linewidth=1.,
label='Symmetric solution, 2.5 kJ/mol')
plt.title("Pyrope-almandine join (model comparison)")
plt.ylabel("Excess gibbs free energy of solution (J/mol)")
plt.xlabel("Almandine fraction")
plt.legend(loc='upper left')
plt.show()
plt.plot(comp, g2_pyrope_activity, 'r-',
linewidth=1., label='Pyrope activity')
plt.plot(comp, g2_almandine_activity, 'b-',
linewidth=1., label='Almandine activity')
plt.plot(comp, g2_pyrope_gamma, 'r-.',
linewidth=1., label='Pyrope activity coefficient')
plt.plot(comp, g2_almandine_gamma, 'b--',
linewidth=1., label='Almandine activity coefficient')
plt.title("Non-ideal pyrope-almandine join")
plt.ylabel("Activities")
plt.xlabel("Almandine fraction")
plt.legend(loc='lower right')
plt.show()
'''
Adding more endmembers is very straightforward.
Interaction terms must be added in the order [[12,13,14,...],[23,24,...],[34,...],...]
Here, the new endmember is majorite, illustrating the addition of endmembers with
multiple occupancy on a single site (here Mg and Si to replace Al)
Here we also illustrate the addition of excess entropy and volume terms.
V and S excess lists are optional parameters to the solid solution initialisation.
In the initialisation, the order of interaction terms is H, V, S
(S excesses are least commonly constrained).
An important note: In Holland and Powell datasets, the "DQF" terms
are equivalent to the excesses in burnman with the *exception* of the temperature correction,
which is the negative of the entropy correction. i.e.
DQF[0] = H_excess
DQF[1] = -S_excess
DQF[2] = V_excess
'''
g3 = burnman.SolidSolution(name = 'Symmetric pyrope-almandine-majorite garnet',
solution_type = 'symmetric',
endmembers = [[minerals.HP_2011_ds62.py(),
'[Mg]3[Al]2Si3O12'],
[minerals.HP_2011_ds62.alm(),
'[Fe]3[Al]2Si3O12'],
[minerals.HP_2011_ds62.maj(),
'[Mg]3[Mg1/2Si1/2]2Si3O12']],
energy_interaction = [[2.5e3, 0.0e3], [10.0e3]],
entropy_interaction = [[0.0e3, 0.0e3], [0.0e3]],
volume_interaction = [[0.0e3, 0.0e3], [0.0e3]])
g3_configurational_entropy = np.empty_like(comp)
g3_excess_entropy = np.empty_like(comp)
for i, c in enumerate(comp):
molar_fractions = [1.0 - c, 0., c]
g3.set_composition(molar_fractions)
g3.set_state(P, T)
g3_configurational_entropy[
i] = g3.solution_model._configurational_entropy(molar_fractions)
g3_excess_entropy[i] = -np.dot(
molar_fractions, g3.solution_model._ideal_excess_partial_gibbs(T, molar_fractions)) / T
plt.plot(comp, g3_configurational_entropy, 'g-',
linewidth=1., label='Configurational entropy')
plt.plot(comp, g3_excess_entropy, 'r-',
linewidth=1., label='Excess entropy')
plt.title("Pyrope-majorite join")
plt.ylabel("Entropy (J/K/mol)")
plt.xlabel("Majorite fraction")
plt.legend(loc='upper left')
plt.show()
'''
The addition of grossular garnet (Ca-Al garnet) illustrates the use of asymmetric solution models
This model is also found in the HP_2011_ds62 database
'''
g4 = burnman.SolidSolution(name = 'garnet',
solution_type = 'asymmetric',
endmembers = [[minerals.HP_2011_ds62.py(),
'[Mg]3[Al]2Si3O12'],
[minerals.HP_2011_ds62.alm(),
'[Fe]3[Al]2Si3O12'],
[minerals.HP_2011_ds62.gr(),
'[Ca]3[Al]2Si3O12'],
[minerals.HP_2011_ds62.andr(),
'[Ca]3[Fe]2Si3O12']],
alphas = [1.0, 1.0, 2.7, 2.7],
energy_interaction = [[2.5e3, 31.e3, 53.2e3],
[5.e3, 37.24e3],
[2.e3]])
g4_excess_gibbs_400 = np.empty_like(comp)
g4_excess_gibbs_800 = np.empty_like(comp)
g4_excess_gibbs_1200 = np.empty_like(comp)
for i, c in enumerate(comp):
molar_fractions = [1.0 - c, 0., c, 0.]
g4.set_composition(molar_fractions)
g4.set_state(P, 400.)
g4_excess_gibbs_400[i] = g4.excess_gibbs
g4.set_state(P, 800.)
g4_excess_gibbs_800[i] = g4.excess_gibbs
g4.set_state(P, 1200.)
g4_excess_gibbs_1200[i] = g4.excess_gibbs
plt.plot(comp, g4_excess_gibbs_400, 'r-', linewidth=1., label='400 K')
plt.plot(comp, g4_excess_gibbs_800, 'g-', linewidth=1., label='800 K')
plt.plot(comp, g4_excess_gibbs_1200, 'b-', linewidth=1., label='1200 K')
plt.title("Pyrope-grossular join (asymmetric model)")
plt.ylabel("Excess gibbs free energy of solution (J/mol)")
plt.xlabel("Grossular fraction")
plt.legend(loc='lower left')
plt.show()
'''
The subregular solution model (Helffrich and Wood, 1989)
provides a more flexible way of constructing an asymmetric
model. Here's a garnet model from Ganguly et al., 1996:
(N.B. Published excesses are on a 4-oxygen (1-cation) basis,
so we need to multiply each by 3.)
'''
mult = lambda x, n: [[[v*n for v in i] for i in j] for j in x]
g5 = burnman.SolidSolution(name = 'Subregular pyrope-almandine-grossular '
'garnet (Ganguly et al., 1996)',
solution_type = 'subregular',
endmembers = [[minerals.HP_2011_ds62.py(),
'[Mg]3[Al]2Si3O12'],
[minerals.HP_2011_ds62.alm(),
'[Fe]3[Al]2Si3O12'],
[minerals.HP_2011_ds62.gr(),
'[Ca]3[Al]2Si3O12'],
[minerals.HP_2011_ds62.spss(),
'[Mn]3[Al]2Si3O12']],
energy_interaction = mult([[[2117., 695.], [9834., 21627.], [12083., 12083.]],
[[6773., 873.], [539., 539.]],
[[0., 0.]]], 3.),
volume_interaction = mult([[[0.07e-5, 0.], [0.058e-5, 0.012e-5], [0.04e-5, 0.03e-5]],
[[0.03e-5, 0.], [0.04e-5, 0.01e-5]],
[[0., 0.]]], 3.),
entropy_interaction = mult([[[0., 0.], [5.78, 5.78], [7.67, 7.67]],
[[1.69, 1.69], [0., 0.]],
[[0., 0.]]], 3.))
g5_excess_gibbs = np.empty_like(comp)
for i, c in enumerate(comp):
molar_fractions = [1. - c, 0., c, 0.]
g5.set_composition(molar_fractions)
g5.set_state(1.e5, 298.15)
g5_excess_gibbs[i] = g5.excess_gibbs
plt.plot(comp, g5_excess_gibbs / 3., 'r-', linewidth=1.,
label='Py-Gr excess gibbs (J/cation-mole)')
plt.title("Asymmetric py-gr join (Ganguly et al., 1996; Figure 5)")
plt.ylabel("Excess gibbs free energy of solution (J/cation-mol)")
plt.xlabel("XCa")
plt.legend(loc='lower left')
plt.show()