-
Notifications
You must be signed in to change notification settings - Fork 854
/
interface_reactions.py
745 lines (618 loc) · 28.1 KB
/
interface_reactions.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
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
# Copyright (c) Pymatgen Development Team.
# Distributed under the terms of the MIT License.
"""
This module provides a class to predict and analyze interfacial reactions between two
solids, with or without an open element (e.g., flowing O2).
Please consider citing one or both of the following papers if you use this
code in your own work.
References:
(1) Richards, W. D., Miara, L. J., Wang, Y., Kim, J. C., & Ceder, G. (2015).
Interface stability in solid-state batteries. Chemistry of Materials, 28(1),
266–273. https://doi.org/10.1021/acs.chemmater.5b04082
(2) Xiao, Y., Wang, Y., Bo, S.-H., Kim, J. C., Miara, L. J., & Ceder, G. (2019).
Understanding interface stability in solid-state batteries.
Nature Reviews Materials, 5(2), 105–126. https://doi.org/10.1038/s41578-019-0157-5
"""
from __future__ import annotations
import json
import os
import warnings
from typing import Literal
import matplotlib.pyplot as plt
import numpy as np
from monty.json import MSONable
from pandas import DataFrame
from plotly.graph_objects import Figure, Scatter
from pymatgen.analysis.phase_diagram import GrandPotentialPhaseDiagram, PhaseDiagram
from pymatgen.analysis.reaction_calculator import Reaction
from pymatgen.core.composition import Composition
from pymatgen.util.plotting import pretty_plot
from pymatgen.util.string import htmlify, latexify
__author__ = "Yihan Xiao, Matthew McDermott"
__maintainer__ = "Matthew McDermott"
__email__ = "mcdermott@lbl.gov"
__date__ = "Sep 1, 2021"
with open(os.path.join(os.path.dirname(__file__), "..", "util", "plotly_interface_rxn_layouts.json")) as f:
plotly_layouts = json.load(f)
class InterfacialReactivity(MSONable):
"""
Class for modeling an interface between two solids and its possible reactions.
The two reactants are provided as Composition objects (c1 and c2), along with the
relevant compositional PhaseDiagram object. Possible reactions are calculated by
finding all points along a tie-line between c1 and c2 where there is a "kink" in
the phase diagram; i.e. a point or facet of the phase diagram.
Please consider citing one or both of the following papers if you use this code
in your own work.
References:
Richards, W. D., Miara, L. J., Wang, Y., Kim, J. C., & Ceder, G. (2015).
Interface stability in solid-state batteries. Chemistry of Materials, 28(1),
266–273. https://doi.org/10.1021/acs.chemmater.5b04082
Xiao, Y., Wang, Y., Bo, S.-H., Kim, J. C., Miara, L. J., & Ceder, G. (2019).
Understanding interface stability in solid-state batteries.
Nature Reviews Materials, 5(2), 105–126.
https://doi.org/10.1038/s41578-019-0157-5
"""
EV_TO_KJ_PER_MOL = 96.4853
def __init__(
self,
c1: Composition,
c2: Composition,
pd: PhaseDiagram,
norm: bool = True,
use_hull_energy: bool = False,
**kwargs,
):
"""
Args:
c1: Reactant 1 composition
c2: Reactant 2 composition
pd: Phase diagram object built from all elements in composition c1 and c2.
norm: Whether or not the total number of atoms in composition
of reactant will be normalized to 1.
use_hull_energy: Whether or not use the convex hull energy for
a given composition for reaction energy calculation. If false,
the energy of ground state structure will be used instead.
Note that in case when ground state can not be found for a
composition, convex hull energy will be used associated with a
warning message.
"""
bypass_grand_warning = kwargs.get("bypass_grand_warning", False)
if isinstance(pd, GrandPotentialPhaseDiagram) and not bypass_grand_warning:
raise ValueError(
"Please use the GrandPotentialInterfacialReactivity "
"class for interfacial reactions with open elements!"
)
self.c1 = c1
self.c2 = c2
self.pd = pd
self.norm = norm
self.use_hull_energy = use_hull_energy
self.c1_original = c1
self.c2_original = c2
self.comp1 = c1
self.comp2 = c2
self.grand = False
# Factor is the compositional ratio between composition self.c1 and
# processed composition self.comp1. For example, the factor for
# Composition('SiO2') and Composition('O') is 2.0. This factor will be used
# to convert mixing ratio in self.comp1 - self.comp2 tie line to that in
# self.c1 - self.c2 tie line.
self.factor1 = 1.0
self.factor2 = 1.0
if self.norm:
self.c1 = c1.fractional_composition
self.c2 = c2.fractional_composition
self.comp1 = self.comp1.fractional_composition
self.comp2 = self.comp2.fractional_composition
if not bypass_grand_warning:
# Computes energies for reactants in different scenarios.
if self.use_hull_energy:
self.e1 = self.pd.get_hull_energy(self.comp1)
self.e2 = self.pd.get_hull_energy(self.comp2)
else:
self.e1 = self._get_entry_energy(self.pd, self.comp1)
self.e2 = self._get_entry_energy(self.pd, self.comp2)
def get_kinks(self) -> list[tuple[int, float, float, Reaction, float]]:
"""
Finds all the kinks in mixing ratio where reaction products changes
along the tie-line of composition self.c1 and composition self.c2.
Returns:
List object of tuples, each of which contains 5 elements:
(index, mixing ratio, reaction energy in eV/atom, Reaction object, reaction
energy per mol of formula in kJ/mol).
"""
c1_coord = self.pd.pd_coords(self.comp1)
c2_coord = self.pd.pd_coords(self.comp2)
n1 = self.comp1.num_atoms
n2 = self.comp2.num_atoms
critical_comp = self.pd.get_critical_compositions(self.comp1, self.comp2)
x_kink, energy_kink, react_kink, energy_per_rxt_formula = [], [], [], []
if (c1_coord == c2_coord).all():
x_kink = [0, 1]
energy_kink = [self._get_energy(x) for x in x_kink]
react_kink = [self._get_reaction(x) for x in x_kink]
num_atoms = [(x * self.comp1.num_atoms + (1 - x) * self.comp2.num_atoms) for x in x_kink]
energy_per_rxt_formula = [
energy_kink[i]
* self._get_elmt_amt_in_rxn(react_kink[i])
/ num_atoms[i]
* InterfacialReactivity.EV_TO_KJ_PER_MOL
for i in range(2)
]
else:
for i in reversed(critical_comp):
# Gets mixing ratio x at kinks.
c = self.pd.pd_coords(i)
x = np.linalg.norm(c - c2_coord) / np.linalg.norm(c1_coord - c2_coord)
# Modifies mixing ratio in case compositions self.comp1 and
# self.comp2 are not normalized.
x = x * n2 / (n1 + x * (n2 - n1))
n_atoms = x * self.comp1.num_atoms + (1 - x) * self.comp2.num_atoms
# Converts mixing ratio in comp1 - comp2 tie line to that in
# c1 - c2 tie line.
x_converted = self._convert(x, self.factor1, self.factor2)
x_kink.append(x_converted)
# Gets reaction energy at kinks
normalized_energy = self._get_energy(x)
energy_kink.append(normalized_energy)
# Gets balanced reaction at kinks
rxt = self._get_reaction(x)
react_kink.append(rxt)
rxt_energy = normalized_energy * self._get_elmt_amt_in_rxn(rxt) / n_atoms
energy_per_rxt_formula.append(rxt_energy * self.EV_TO_KJ_PER_MOL)
index_kink = range(1, len(critical_comp) + 1)
return list(zip(index_kink, x_kink, energy_kink, react_kink, energy_per_rxt_formula))
def plot(self, backend: Literal["plotly", "matplotlib"] = "plotly") -> Figure | plt.Figure:
"""
Plots reaction energy as a function of mixing ratio x in self.c1 - self.c2
tie line.
Args:
backend ("plotly" | "matplotlib"): Plotting library used to create the plot. Defaults to
"plotly" but can also be "matplotlib".
Returns:
Plot of reaction energies as a function of mixing ratio
"""
if backend.lower() == "plotly":
fig = self._get_plotly_figure()
elif backend.lower() in ["matplotlib", "mpl", "plt"]:
fig = self._get_matplotlib_figure()
else:
raise ValueError("The provided backend is not a valid option!")
return fig
def get_dataframe(self) -> DataFrame:
"""
Returns a pandas DataFrame representation of the data produced by the
get_kinks() method.
"""
rxns = [
{
"Atomic fraction": round(ratio, 3),
"Reaction": rxn,
"E$_{\textrm{rxn}}$ (kJ/mol)": round(rxn_energy, 1),
"E$_{\textrm{rxn}}$ (eV/atom)": round(reactivity, 3),
}
for _, ratio, reactivity, rxn, rxn_energy in self.get_kinks()
]
df = DataFrame(rxns)
return df
def get_critical_original_kink_ratio(self):
"""
Returns a list of molar mixing ratio for each kink between ORIGINAL
(instead of processed) reactant compositions. This is the
same list as mixing ratio obtained from get_kinks method
if self.norm = False.
Returns:
A list of floats representing molar mixing ratios between
the original reactant compositions for each kink.
"""
ratios = []
if self.c1_original == self.c2_original:
return [0, 1]
reaction_kink = [k[3] for k in self.get_kinks()]
for rxt in reaction_kink:
ratios.append(abs(self._get_original_composition_ratio(rxt)))
return ratios
def _get_original_composition_ratio(self, reaction):
"""
Returns the molar mixing ratio between the reactants with ORIGINAL (
instead of processed) compositions for a reaction.
Args:
reaction (Reaction): Reaction object that contains the original
reactant compositions.
Returns:
The molar mixing ratio between the original reactant
compositions for a reaction.
"""
if self.c1_original == self.c2_original:
return 1
c1_coeff = reaction.get_coeff(self.c1_original) if self.c1_original in reaction.reactants else 0
c2_coeff = reaction.get_coeff(self.c2_original) if self.c2_original in reaction.reactants else 0
return c1_coeff * 1.0 / (c1_coeff + c2_coeff)
def _get_energy(self, x):
"""
Computes reaction energy in eV/atom at mixing ratio x : (1-x) for
self.comp1 : self.comp2.
Args:
x (float): Mixing ratio x of reactants, a float between 0 and 1.
Returns:
Reaction energy.
"""
return self.pd.get_hull_energy(self.comp1 * x + self.comp2 * (1 - x)) - self.e1 * x - self.e2 * (1 - x)
def _get_reactants(self, x: float) -> list[Composition]:
"""Returns a list of relevant reactant compositions given an x coordinate"""
# Uses original composition for reactants.
if np.isclose(x, 0):
reactants = [self.c2_original]
elif np.isclose(x, 1):
reactants = [self.c1_original]
else:
reactants = list({self.c1_original, self.c2_original})
return reactants
def _get_reaction(self, x: float) -> Reaction:
"""
Generates balanced reaction at mixing ratio x : (1-x) for
self.comp1 : self.comp2.
Args:
x (float): Mixing ratio x of reactants, a float between 0 and 1.
Returns:
Reaction object.
"""
mix_comp = self.comp1 * x + self.comp2 * (1 - x)
decomp = self.pd.get_decomposition(mix_comp)
reactants = self._get_reactants(x)
product = [Composition(k.name) for k, v in decomp.items()]
reaction = Reaction(reactants, product)
x_original = self._get_original_composition_ratio(reaction)
if np.isclose(x_original, 1):
reaction.normalize_to(self.c1_original, x_original)
else:
reaction.normalize_to(self.c2_original, 1 - x_original)
return reaction
def _get_elmt_amt_in_rxn(self, rxn: Reaction) -> int:
"""
Computes total number of atoms in a reaction formula for elements
not in external reservoir. This method is used in the calculation
of reaction energy per mol of reaction formula.
Args:
rxn: a Reaction object.
Returns:
Total number of atoms for non_reservoir elements.
"""
return sum(rxn.get_el_amount(e) for e in self.pd.elements)
def _get_plotly_figure(self) -> Figure:
"""Returns a Plotly figure of reaction kinks diagram"""
kinks = map(list, zip(*self.get_kinks()))
_, x, energy, reactions, _ = kinks
lines = Scatter(
x=x,
y=energy,
mode="lines",
name="Lines",
line=dict(color="navy", dash="solid", width=5.0),
hoverinfo="none",
)
annotations = self._get_plotly_annotations(x, energy, reactions) # type: ignore
min_idx = energy.index(min(energy)) # type: ignore
x_min = x.pop(min_idx)
e_min = energy.pop(min_idx)
rxn_min = reactions.pop(min_idx)
labels = [
rf"{htmlify(str(r))} <br>" + "\u0394" + f"E<sub>rxn</sub> = {round(e, 3)} eV/atom" # type: ignore
for r, e in zip(reactions, energy)
]
markers = Scatter(
x=x,
y=energy,
mode="markers",
name="Reactions",
hoverinfo="text",
hovertext=labels,
marker=dict(
color="black",
size=12,
opacity=0.8,
line=dict(color="black", width=3),
),
hoverlabel=dict(bgcolor="navy"),
)
min_label = (
rf"{htmlify(str(rxn_min))} <br>" + "\u0394" + f"E<sub>rxn</sub> = {round(e_min, 3)} eV/atom" # type: ignore
)
minimum = Scatter(
x=[x_min],
y=[e_min],
mode="markers",
hoverinfo="text",
hovertext=[min_label],
marker=dict(color="darkred", size=24, symbol="star"),
name="Suggested reaction",
)
data = [lines, markers, minimum]
layout = plotly_layouts["default_interface_rxn_layout"]
layout["xaxis"]["title"] = self._get_xaxis_title(latex=False)
layout["annotations"] = annotations
fig = Figure(data=data, layout=layout)
return fig
def _get_matplotlib_figure(self) -> plt.Figure:
"""Returns a matplotlib figure of reaction kinks diagram"""
pretty_plot(8, 5)
plt.xlim([-0.05, 1.05]) # plot boundary is 5% wider on each side
kinks = list(zip(*self.get_kinks()))
_, x, energy, reactions, _ = kinks
plt.plot(x, energy, "o-", markersize=8, c="navy", zorder=1)
plt.scatter(self.minimum[0], self.minimum[1], marker="*", c="red", s=400, zorder=2)
for x_coord, y_coord, rxn in zip(x, energy, reactions):
products = ", ".join(
[latexify(p.reduced_formula) for p in rxn.products if not np.isclose(rxn.get_coeff(p), 0)]
)
plt.annotate(
products,
xy=(x_coord, y_coord),
xytext=(10, -30),
textcoords="offset points",
ha="right",
va="bottom",
arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=0"),
)
if self.norm:
plt.ylabel("Energy (eV/atom)")
else:
plt.ylabel("Energy (eV/f.u.)")
plt.xlabel(self._get_xaxis_title())
plt.ylim(self.minimum[1] + 0.05 * self.minimum[1]) # plot boundary is 5% lower
fig = plt.gcf()
plt.close(fig)
return fig
def _get_xaxis_title(self, latex: bool = True) -> str:
"""Returns the formatted title of the x axis (using either html/latex)"""
if latex:
f1 = latexify(self.c1.reduced_formula)
f2 = latexify(self.c2.reduced_formula)
title = f"$x$ in $x${f1} + $(1-x)${f2}"
else:
f1 = htmlify(self.c1.reduced_formula)
f2 = htmlify(self.c2.reduced_formula)
title = f"<i>x</i> in <i>x</i>{f1} + (1-<i>x</i>){f2}"
return title
@staticmethod
def _get_plotly_annotations(x: list[float], y: list[float], reactions: list[Reaction]):
"""Returns dictionary of annotations for the Plotly figure layout"""
annotations = []
for x_coord, y_coord, rxn in zip(x, y, reactions):
products = ", ".join(
[htmlify(p.reduced_formula) for p in rxn.products if not np.isclose(rxn.get_coeff(p), 0)]
)
annotation = dict(x=x_coord, y=y_coord, text=products, font=dict(size=18), ax=-25, ay=55)
annotations.append(annotation)
return annotations
@staticmethod
def _get_entry_energy(pd: PhaseDiagram, composition: Composition):
"""
Finds the lowest entry energy for entries matching the composition.
Entries with non-negative formation energies are excluded. If no
entry is found, use the convex hull energy for the composition.
Args:
pd: Phase diagram object
composition: Composition object that the target entry should match
Returns:
The lowest entry energy among entries matching the composition.
"""
candidate = [
i.energy_per_atom
for i in pd.qhull_entries
if i.composition.fractional_composition == composition.fractional_composition
]
if not candidate:
warnings.warn(
"The reactant " + composition.reduced_formula + " has no matching entry with negative formation"
" energy, instead convex hull energy for this"
" composition will be used for reaction energy "
"calculation. "
)
return pd.get_hull_energy(composition)
min_entry_energy = min(candidate)
return min_entry_energy * composition.num_atoms
@staticmethod
def _convert(x: float, factor1: float, factor2: float):
"""
Converts mixing ratio x in comp1 - comp2 tie line to that in
c1 - c2 tie line.
Args:
x: Mixing ratio x in comp1 - comp2 tie line, a float
between 0 and 1.
factor1: Compositional ratio between composition c1 and
processed composition comp1. E.g., factor for
Composition('SiO2') and Composition('O') is 2.0.
factor2: Compositional ratio between composition c2 and
processed composition comp2.
Returns:
Mixing ratio in c1 - c2 tie line, a float between 0 and 1.
"""
return x * factor2 / ((1 - x) * factor1 + x * factor2)
@staticmethod
def _reverse_convert(x: float, factor1: float, factor2: float):
"""
Converts mixing ratio x in c1 - c2 tie line to that in
comp1 - comp2 tie line.
Args:
x: Mixing ratio x in c1 - c2 tie line, a float between
0 and 1.
factor1: Compositional ratio between composition c1 and
processed composition comp1. E.g., factor for
Composition('SiO2') and Composition('O') is 2.
factor2: Compositional ratio between composition c2 and
processed composition comp2.
Returns:
Mixing ratio in comp1 - comp2 tie line, a float between 0 and 1.
"""
return x * factor1 / ((1 - x) * factor2 + x * factor1)
@classmethod
def get_chempot_correction(cls, element: str, temp: float, pres: float):
"""
Get the normalized correction term Δμ for chemical potential of a gas
phase consisting of element at given temperature and pressure,
referenced to that in the standard state (T_std = 298.15 K,
T_std = 1 bar). The gas phase is limited to be one of O2, N2, Cl2,
F2, H2. Calculation formula can be found in the documentation of
Materials Project website.
Args:
element: The string representing the element.
temp: The temperature of the gas phase in Kelvin.
pres: The pressure of the gas phase in Pa.
Returns:
The correction of chemical potential in eV/atom of the gas
phase at given temperature and pressure.
"""
if element not in ["O", "N", "Cl", "F", "H"]:
warnings.warn(f"Element {element} not one of valid options: ['O', 'N', 'Cl', 'F', 'H']")
return 0
std_temp = 298.15
std_pres = 1e5
ideal_gas_const = 8.3144598
# Cp and S at standard state in J/(K.mol). Data from NIST-JANAF tables
# Tables: O-029, N-O23, Cl-073, F-054, H-050
cp_dict = {"O": 29.376, "N": 29.124, "Cl": 33.949, "F": 31.302, "H": 28.836}
s_dict = {"O": 205.147, "N": 191.609, "Cl": 223.079, "F": 202.789, "H": 130.680}
cp_std = cp_dict[element]
s_std = s_dict[element]
pv_correction = ideal_gas_const * temp * np.log(pres / std_pres)
ts_correction = (
-cp_std * (temp * np.log(temp) - std_temp * np.log(std_temp))
+ cp_std * (temp - std_temp) * (1 + np.log(std_temp))
- s_std * (temp - std_temp)
)
dg = pv_correction + ts_correction
dg /= 1000 * cls.EV_TO_KJ_PER_MOL # convert to eV/f.u.
dg /= 2 # convert from eV/f.u. to eV/atom (2 atoms in a diatomic gas)
return dg
@property
def labels(self):
"""
Returns a dictionary containing kink information:
{index: 'x= mixing_ratio energy= reaction_energy reaction_equation'}.
E.g., {1: 'x= 0.0 energy = 0.0 Mn -> Mn',
2: 'x= 0.5 energy = -15.0 O2 + Mn -> MnO2',
3: 'x= 1.0 energy = 0.0 O2 -> O2'}.
"""
return {
j: "x= " + str(round(x, 4)) + " energy in eV/atom = " + str(round(energy, 4)) + " " + str(reaction)
for j, x, energy, reaction, _ in self.get_kinks()
}
@property
def minimum(self):
"""
Finds the minimum reaction energy E_min and corresponding
mixing ratio x_min.
Returns:
Tuple (x_min, E_min).
"""
return min(((x, energy) for _, x, energy, _, _ in self.get_kinks()), key=lambda i: i[1])
@property
def products(self):
"""
List of formulas of potential products. E.g., ['Li','O2','Mn'].
"""
products = set()
for _, _, _, react, _ in self.get_kinks():
products = products.union({k.reduced_formula for k in react.products})
return list(products)
class GrandPotentialInterfacialReactivity(InterfacialReactivity):
"""
Extends upon InterfacialReactivity to allow for modelling possible reactions
at the interface between two solids in the presence of an open element. The
thermodynamics of the open system are provided by the user via the
GrandPotentialPhaseDiagram class.
"""
def __init__(
self,
c1: Composition,
c2: Composition,
grand_pd: GrandPotentialPhaseDiagram,
pd_non_grand: PhaseDiagram,
include_no_mixing_energy: bool = False,
norm: bool = True,
use_hull_energy: bool = True,
):
"""
Args:
c1: Reactant 1 composition
c2: Reactant 2 composition
grand_pd: Grand potential phase diagram object built from all elements in
composition c1 and c2.
include_no_mixing_energy: No_mixing_energy for a reactant is the
opposite number of its energy above grand potential convex hull. In
cases where reactions involve elements reservoir, this param
determines whether no_mixing_energy of reactants will be included
in the final reaction energy calculation. By definition, if pd is
not a GrandPotentialPhaseDiagram object, this param is False.
pd_non_grand: PhaseDiagram object but not
GrandPotentialPhaseDiagram object built from elements in c1 and c2.
norm: Whether or not the total number of atoms in composition
of reactant will be normalized to 1.
use_hull_energy: Whether or not use the convex hull energy for
a given composition for reaction energy calculation. If false,
the energy of ground state structure will be used instead.
Note that in case when ground state can not be found for a
composition, convex hull energy will be used associated with a
warning message.
"""
if not isinstance(grand_pd, GrandPotentialPhaseDiagram):
raise ValueError("Please use the InterfacialReactivity class if using a regular phase diagram!")
super().__init__(
c1=c1, c2=c2, pd=grand_pd, norm=norm, use_hull_energy=use_hull_energy, bypass_grand_warning=True
)
self.pd_non_grand = pd_non_grand
self.grand = True
self.comp1 = Composition({k: v for k, v in c1.items() if k not in grand_pd.chempots})
self.comp2 = Composition({k: v for k, v in c2.items() if k not in grand_pd.chempots})
if self.norm:
self.factor1 = self.comp1.num_atoms / c1.num_atoms
self.factor2 = self.comp2.num_atoms / c2.num_atoms
self.comp1 = self.comp1.fractional_composition
self.comp2 = self.comp2.fractional_composition
if include_no_mixing_energy:
self.e1 = self._get_grand_potential(self.c1)
self.e2 = self._get_grand_potential(self.c2)
else:
self.e1 = self.pd.get_hull_energy(self.comp1)
self.e2 = self.pd.get_hull_energy(self.comp2)
def get_no_mixing_energy(self):
"""
Generates the opposite number of energy above grand potential
convex hull for both reactants.
Returns:
[(reactant1, no_mixing_energy1),(reactant2,no_mixing_energy2)].
"""
energy1 = self.pd.get_hull_energy(self.comp1) - self._get_grand_potential(self.c1)
energy2 = self.pd.get_hull_energy(self.comp2) - self._get_grand_potential(self.c2)
unit = "eV/f.u."
if self.norm:
unit = "eV/atom"
return [
(f"{self.c1_original.reduced_formula} ({unit})", energy1),
(f"{self.c2_original.reduced_formula} ({unit})", energy2),
]
def _get_reactants(self, x: float) -> list[Composition]:
"""Returns a list of relevant reactant compositions given an x coordinate"""
reactants = super()._get_reactants(x)
reactants += [Composition(e.symbol) for e, v in self.pd.chempots.items()]
return reactants
def _get_grand_potential(self, composition: Composition) -> float:
"""
Computes the grand potential Phi at a given composition and
chemical potential(s).
Args:
composition: Composition object.
Returns:
Grand potential at a given composition at chemical potential(s).
"""
if self.use_hull_energy:
grand_potential = self.pd_non_grand.get_hull_energy(composition)
else:
grand_potential = self._get_entry_energy(self.pd_non_grand, composition)
grand_potential -= sum(composition[e] * mu for e, mu in self.pd.chempots.items())
if self.norm:
# Normalizes energy to the composition excluding element(s)
# from reservoir.
grand_potential /= sum(composition[el] for el in composition if el not in self.pd.chempots)
return grand_potential