/
pourbaix_diagram.py
1088 lines (910 loc) · 41.1 KB
/
pourbaix_diagram.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
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
"""
This module is intended to be used to compute Pourbaix diagrams of arbitrary compositions
and formation energies.
"""
from __future__ import annotations
import itertools
import logging
import re
import warnings
from copy import deepcopy
from functools import cmp_to_key, partial
from multiprocessing import Pool
from typing import TYPE_CHECKING, Any, no_type_check
import numpy as np
from monty.json import MontyDecoder, MSONable
from scipy.spatial import ConvexHull, HalfspaceIntersection
from scipy.special import comb
from pymatgen.analysis.phase_diagram import PDEntry, PhaseDiagram
from pymatgen.analysis.reaction_calculator import Reaction, ReactionError
from pymatgen.core import Composition, Element
from pymatgen.core.ion import Ion
from pymatgen.entries.compatibility import MU_H2O
from pymatgen.entries.computed_entries import ComputedEntry
from pymatgen.util.coord import Simplex
from pymatgen.util.due import Doi, due
from pymatgen.util.plotting import pretty_plot
from pymatgen.util.string import Stringify
if TYPE_CHECKING:
import matplotlib.pyplot as plt
__author__ = "Sai Jayaraman"
__copyright__ = "Copyright 2012, The Materials Project"
__version__ = "0.4"
__maintainer__ = "Joseph Montoya"
__credits__ = "Arunima Singh, Joseph Montoya, Anjli Patel"
__email__ = "joseph.montoya@tri.global"
__status__ = "Production"
__date__ = "Nov 1, 2012"
# If you use this module in your work, consider citing:
due.cite(
Doi("10.1103/PhysRevB.85.235438"),
description="Prediction of solid-aqueous equilibria: Scheme to combine first-principles "
"calculations of solids with experimental aqueous states",
)
due.cite(
Doi("10.1021/acs.chemmater.7b03980"),
description="Electrochemical Stability of Metastable Materials",
)
due.cite(
Doi("10.1021/acs.chemmater.7b03980"),
description="Fast computation of many-element Pourbaix diagrams",
)
logger = logging.getLogger(__name__)
PREFAC = 0.0591
# TODO: Revise to more closely reflect PDEntry, invoke from energy/composition
# TODO: PourbaixEntries depend implicitly on having entry energies be
# formation energies, should be a better way to get from raw energies
# TODO: uncorrected_energy is a bit of a misnomer, but not sure what to rename
class PourbaixEntry(MSONable, Stringify):
"""
An object encompassing all data relevant to a solid or ion
in a Pourbaix diagram. Each bulk solid/ion has an energy
g of the form: e = e0 + 0.0591 log10(conc) - nO mu_H2O
+ (nH - 2nO) pH + phi (-nH + 2nO + q).
Note that the energies corresponding to the input entries
should be formation energies with respect to hydrogen and
oxygen gas in order for the Pourbaix diagram formalism to
work. This may be changed to be more flexible in the future.
"""
def __init__(self, entry, entry_id=None, concentration=1e-6):
"""
Args:
entry (ComputedEntry/ComputedStructureEntry/PDEntry/IonEntry): An
entry object
entry_id ():
concentration ():
"""
self.entry = entry
if isinstance(entry, IonEntry):
self.concentration = concentration
self.phase_type = "Ion"
self.charge = entry.ion.charge
else:
self.concentration = 1.0
self.phase_type = "Solid"
self.charge = 0
self.uncorrected_energy = entry.energy
if entry_id is not None:
self.entry_id = entry_id
elif hasattr(entry, "entry_id") and entry.entry_id:
self.entry_id = entry.entry_id
else:
self.entry_id = None
@property
def npH(self):
"""Get the number of H."""
return self.entry.composition.get("H", 0) - 2 * self.entry.composition.get("O", 0)
@property
def nH2O(self):
"""Get the number of H2O."""
return self.entry.composition.get("O", 0)
@property
def nPhi(self):
"""Get the number of electrons."""
return self.npH - self.charge
@property
def name(self):
"""Get the name for entry."""
if self.phase_type == "Solid":
return f"{self.entry.reduced_formula}(s)"
return self.entry.name
@property
def energy(self):
"""
Returns (float): total energy of the Pourbaix
entry (at pH, V = 0 vs. SHE).
"""
# Note: this implicitly depends on formation energies as input
return self.uncorrected_energy + self.conc_term - (MU_H2O * self.nH2O)
@property
def energy_per_atom(self):
"""
energy per atom of the Pourbaix entry.
Returns (float): energy per atom
"""
return self.energy / self.composition.num_atoms
@property
def elements(self):
"""Returns elements in the entry."""
return self.entry.elements
def energy_at_conditions(self, pH, V):
"""
Get free energy for a given pH and V.
Args:
pH (float): pH at which to evaluate free energy
V (float): voltage at which to evaluate free energy
Returns:
free energy at conditions
"""
return self.energy + self.npH * PREFAC * pH + self.nPhi * V
def get_element_fraction(self, element):
"""
Gets the elemental fraction of a given non-OH element.
Args:
element (Element or str): string or element corresponding
to element to get from composition
Returns:
fraction of element / sum(all non-OH elements)
"""
return self.composition.get(element) * self.normalization_factor
@property
def normalized_energy(self):
"""
Returns:
energy normalized by number of non H or O atoms, e. g.
for Zn2O6, energy / 2 or for AgTe3(OH)3, energy / 4.
"""
return self.energy * self.normalization_factor
def normalized_energy_at_conditions(self, pH, V):
"""
Energy at an electrochemical condition, compatible with
numpy arrays for pH/V input.
Args:
pH (float): pH at condition
V (float): applied potential at condition
Returns:
energy normalized by number of non-O/H atoms at condition
"""
return self.energy_at_conditions(pH, V) * self.normalization_factor
@property
def conc_term(self):
"""
Returns the concentration contribution to the free energy,
and should only be present when there are ions in the entry.
"""
return PREFAC * np.log10(self.concentration)
# TODO: not sure if these are strictly necessary with refactor
def as_dict(self):
"""
Returns dict which contains Pourbaix Entry data.
Note that the pH, voltage, H2O factors are always calculated when
constructing a PourbaixEntry object.
"""
dct = {"@module": type(self).__module__, "@class": type(self).__name__}
if isinstance(self.entry, IonEntry):
dct["entry_type"] = "Ion"
else:
dct["entry_type"] = "Solid"
dct["entry"] = self.entry.as_dict()
dct["concentration"] = self.concentration
dct["entry_id"] = self.entry_id
return dct
@classmethod
def from_dict(cls, d):
"""Invokes a PourbaixEntry from a dictionary."""
entry_type = d["entry_type"]
entry = IonEntry.from_dict(d["entry"]) if entry_type == "Ion" else MontyDecoder().process_decoded(d["entry"])
entry_id = d["entry_id"]
concentration = d["concentration"]
return cls(entry, entry_id, concentration)
@property
def normalization_factor(self):
"""Sum of number of atoms minus the number of H and O in composition."""
return 1.0 / (self.num_atoms - self.composition.get("H", 0) - self.composition.get("O", 0))
@property
def composition(self):
"""Returns composition."""
return self.entry.composition
@property
def num_atoms(self):
"""Return number of atoms in current formula. Useful for normalization."""
return self.composition.num_atoms
def to_pretty_string(self) -> str:
"""A pretty string representation."""
if self.phase_type == "Solid":
return f"{self.entry.reduced_formula}(s)"
return self.entry.name
def __repr__(self):
energy, npH, nPhi, nH2O, entry_id = self.energy, self.npH, self.nPhi, self.nH2O, self.entry_id
return (
f"{type(self).__name__}({self.entry.composition} with {energy=:.4f}, {npH=}, "
f"{nPhi=}, {nH2O=}, {entry_id=})"
)
class MultiEntry(PourbaixEntry):
"""
PourbaixEntry-like object for constructing multi-elemental Pourbaix diagrams.
"""
def __init__(self, entry_list, weights=None):
"""
Initializes a MultiEntry.
Args:
entry_list ([PourbaixEntry]): List of component PourbaixEntries
weights ([float]): Weights associated with each entry. Default is None
"""
self.weights = weights or [1.0] * len(entry_list)
self.entry_list = entry_list
def __getattr__(self, attr):
"""
Because most of the attributes here are just weighted averages of the entry_list,
we save some space by having a set of conditionals to define the attributes.
"""
# Attributes that are weighted averages of entry attributes
if attr in ["energy", "npH", "nH2O", "nPhi", "conc_term", "composition", "uncorrected_energy", "elements"]:
# TODO: Composition could be changed for compat with sum
start = Composition({}) if attr == "composition" else 0
weighted_values = (getattr(entry, attr) * weight for entry, weight in zip(self.entry_list, self.weights))
return sum(weighted_values, start)
# Attributes that are just lists of entry attributes
if attr in ["entry_id", "phase_type"]:
return [getattr(e, attr) for e in self.entry_list]
# normalization_factor, num_atoms should work from superclass
return self.__getattribute__(attr)
@property
def name(self):
"""MultiEntry name, i. e. the name of each entry joined by ' + '."""
return " + ".join(e.name for e in self.entry_list)
def __repr__(self):
energy, npH, nPhi, nH2O, entry_id = self.energy, self.npH, self.nPhi, self.nH2O, self.entry_id
cls_name, species = type(self).__name__, self.name
return f"Pourbaix{cls_name}({energy=:.4f}, {npH=}, {nPhi=}, {nH2O=}, {entry_id=}, {species=})"
def as_dict(self):
"""Returns: MSONable dict."""
return {
"@module": type(self).__module__,
"@class": type(self).__name__,
"entry_list": [e.as_dict() for e in self.entry_list],
"weights": self.weights,
}
@classmethod
def from_dict(cls, dct):
"""
Args:
dct (dict): Dict representation.
Returns:
MultiEntry
"""
entry_list = [PourbaixEntry.from_dict(entry) for entry in dct.get("entry_list")]
return cls(entry_list, dct.get("weights"))
# TODO: this class isn't particularly useful in its current form, could be
# refactored to include information about the reference solid
class IonEntry(PDEntry):
"""
Object similar to PDEntry, but contains an Ion object instead of a
Composition object.
Attributes:
name (str): A name for the entry. This is the string shown in the phase diagrams.
By default, this is the reduced formula for the composition, but can be
set to some other string for display purposes.
"""
def __init__(self, ion: Ion, energy: float, name: str | None = None, attribute=None):
"""
Args:
ion: Ion object
energy: Energy for composition.
name: Optional parameter to name the entry. Defaults to the
chemical formula.
attribute: Optional attribute of the entry, e.g., band gap.
"""
self.ion = ion
# Auto-assign name
name = name or self.ion.reduced_formula
super().__init__(composition=ion.composition, energy=energy, name=name, attribute=attribute)
@classmethod
def from_dict(cls, d):
"""Returns an IonEntry object from a dict."""
return cls(Ion.from_dict(d["ion"]), d["energy"], d.get("name"), d.get("attribute"))
def as_dict(self):
"""Creates a dict of composition, energy, and ion name."""
return {"ion": self.ion.as_dict(), "energy": self.energy, "name": self.name}
def __repr__(self):
return f"IonEntry : {self.composition} with energy = {self.energy:.4f}"
def ion_or_solid_comp_object(formula):
"""
Returns either an ion object or composition object given
a formula.
Args:
formula: String formula. Eg. of ion: NaOH(aq), Na[+];
Eg. of solid: Fe2O3(s), Fe(s), Na2O
Returns:
Composition/Ion object
"""
match = re.search(r"\[([^\[\]]+)\]|\(aq\)", formula)
if match:
comp_obj = Ion.from_formula(formula)
elif re.search(r"\(s\)", formula):
comp_obj = Composition(formula[:-3])
else:
comp_obj = Composition(formula)
return comp_obj
ELEMENTS_HO = {Element("H"), Element("O")}
# TODO: the solids filter breaks some of the functionality of the
# heatmap plotter, because the reference states for decomposition
# don't include oxygen/hydrogen in the OER/HER regions
# TODO: create a from_phase_diagram class method for non-formation energy invocation
# TODO: invocation from a MultiEntry entry list could be a bit more robust
# TODO: serialization is still a bit rough around the edges
class PourbaixDiagram(MSONable):
"""Class to create a Pourbaix diagram from entries."""
def __init__(
self,
entries: list[PourbaixEntry] | list[MultiEntry],
comp_dict: dict[str, float] | None = None,
conc_dict: dict[str, float] | None = None,
filter_solids: bool = True,
nproc: int | None = None,
):
"""
Args:
entries ([PourbaixEntry] or [MultiEntry]): Entries list
containing Solids and Ions or a list of MultiEntries
comp_dict (dict[str, float]): Dictionary of compositions,
defaults to equal parts of each elements
conc_dict (dict[str, float]): Dictionary of ion concentrations,
defaults to 1e-6 for each element
filter_solids (bool): applying this filter to a Pourbaix
diagram ensures all included solid phases are filtered by
stability on the compositional phase diagram. Defaults to True.
The practical consequence of this is that highly oxidized or reduced
phases that might show up in experiments due to kinetic limitations
on oxygen/hydrogen evolution won't appear in the diagram, but they are
not actually "stable" (and are frequently overstabilized from DFT errors).
Hence, including only the stable solid phases generally leads to the
most accurate Pourbaix diagrams.
nproc (int): number of processes to generate multi-entries with
in parallel. Defaults to None (serial processing).
"""
entries = deepcopy(entries)
self.filter_solids = filter_solids
# Get non-OH elements
self.pbx_elts = list(
set(itertools.chain.from_iterable([entry.composition.elements for entry in entries])) - ELEMENTS_HO
)
self.dim = len(self.pbx_elts) - 1
# Process multi-entry inputs
if isinstance(entries[0], MultiEntry):
self._processed_entries = entries
# Extract individual entries
single_entries = list(set(itertools.chain.from_iterable([e.entry_list for e in entries])))
self._unprocessed_entries = single_entries
self._filtered_entries = single_entries
self._conc_dict = None
self._elt_comp = {k: v for k, v in entries[0].composition.items() if k not in ELEMENTS_HO}
self._multi_element = True
# Process single entry inputs
else:
# Set default conc/comp dicts
if not comp_dict:
comp_dict = {elt.symbol: 1 / len(self.pbx_elts) for elt in self.pbx_elts}
if not conc_dict:
conc_dict = {elt.symbol: 1e-6 for elt in self.pbx_elts}
self._conc_dict = conc_dict
self._elt_comp = comp_dict
self.pourbaix_elements = self.pbx_elts
solid_entries = [entry for entry in entries if entry.phase_type == "Solid"]
ion_entries = [entry for entry in entries if entry.phase_type == "Ion"]
# If a conc_dict is specified, override individual entry concentrations
for entry in ion_entries:
ion_elts = list(set(entry.elements) - ELEMENTS_HO)
# TODO: the logic here for ion concentration setting is in two
# places, in PourbaixEntry and here, should be consolidated
if len(ion_elts) == 1:
entry.concentration = conc_dict[ion_elts[0].symbol] * entry.normalization_factor
elif len(ion_elts) > 1 and not entry.concentration:
raise ValueError("Elemental concentration not compatible with multi-element ions")
self._unprocessed_entries = solid_entries + ion_entries
if len(solid_entries + ion_entries) != len(entries):
raise ValueError('All supplied entries must have a phase type of either "Solid" or "Ion"')
if self.filter_solids:
# O is 2.46 b/c pbx entry finds energies referenced to H2O
entries_HO = [ComputedEntry("H", 0), ComputedEntry("O", 2.46)]
solid_pd = PhaseDiagram(solid_entries + entries_HO)
solid_entries = list(set(solid_pd.stable_entries) - set(entries_HO))
self._filtered_entries = solid_entries + ion_entries
if len(comp_dict) > 1:
self._multi_element = True
self._processed_entries = self._preprocess_pourbaix_entries(self._filtered_entries, nproc=nproc)
else:
self._processed_entries = self._filtered_entries
self._multi_element = False
self._stable_domains, self._stable_domain_vertices = self.get_pourbaix_domains(self._processed_entries)
def _convert_entries_to_points(self, pourbaix_entries):
"""
Args:
pourbaix_entries ([PourbaixEntry]): list of Pourbaix entries
to process into vectors in nph-nphi-composition space.
Returns:
list of vectors, [[nph, nphi, e0, x1, x2, ..., xn-1]]
corresponding to each entry in nph-nphi-composition space
"""
vecs = [
[entry.npH, entry.nPhi, entry.energy] + [entry.composition.get(elt) for elt in self.pbx_elts[:-1]]
for entry in pourbaix_entries
]
vecs = np.array(vecs)
norms = np.transpose([[entry.normalization_factor for entry in pourbaix_entries]])
vecs *= norms
return vecs
def _get_hull_in_nph_nphi_space(self, entries) -> tuple[list[PourbaixEntry], list[Simplex]]:
"""
Generates convex hull of Pourbaix diagram entries in composition,
npH, and nphi space. This enables filtering of multi-entries
such that only compositionally stable combinations of entries
are included.
Args:
entries ([PourbaixEntry]): list of PourbaixEntries to construct
the convex hull
Returns:
tuple[list[PourbaixEntry], list[Simplex]]: PourbaixEntry list and stable
facets corresponding to that list
"""
ion_entries = [entry for entry in entries if entry.phase_type == "Ion"]
solid_entries = [entry for entry in entries if entry.phase_type == "Solid"]
# Pre-filter solids based on min at each composition
logger.debug("Pre-filtering solids by min energy at each composition")
sorted_entries = sorted(
solid_entries,
key=lambda x: (x.composition.reduced_composition, x.entry.energy_per_atom),
)
grouped_by_composition = itertools.groupby(sorted_entries, key=lambda x: x.composition.reduced_composition)
min_entries = [next(iter(grouped_entries)) for comp, grouped_entries in grouped_by_composition]
min_entries += ion_entries
logger.debug("Constructing nph-nphi-composition points for qhull")
vecs = self._convert_entries_to_points(min_entries)
maxes = np.max(vecs[:, :3], axis=0)
extra_point = np.concatenate([maxes, np.ones(self.dim) / self.dim], axis=0)
# Add padding for extra point
pad = 1000
extra_point[2] += pad
points = np.concatenate([vecs, np.array([extra_point])], axis=0)
logger.debug("Constructing convex hull in nph-nphi-composition space")
hull = ConvexHull(points, qhull_options="QJ i")
# Create facets and remove top
facets = [facet for facet in hull.simplices if len(points) - 1 not in facet]
if self.dim > 1:
logger.debug("Filtering facets by Pourbaix composition")
valid_facets = []
for facet in facets:
comps = vecs[facet][:, 3:]
full_comps = np.concatenate([comps, 1 - np.sum(comps, axis=1).reshape(len(comps), 1)], axis=1)
# Ensure an compositional interior point exists in the simplex
if np.linalg.matrix_rank(full_comps) > self.dim:
valid_facets.append(facet)
else:
valid_facets = facets
return min_entries, valid_facets
def _preprocess_pourbaix_entries(self, entries, nproc=None):
"""
Generates multi-entries for Pourbaix diagram.
Args:
entries ([PourbaixEntry]): list of PourbaixEntries to preprocess
into MultiEntries
nproc (int): number of processes to be used in parallel
treatment of entry combos
Returns:
list[MultiEntry]: stable MultiEntry candidates
"""
# Get composition
tot_comp = Composition(self._elt_comp)
min_entries, valid_facets = self._get_hull_in_nph_nphi_space(entries)
combos = []
for facet in valid_facets:
for idx in range(1, self.dim + 2):
these_combos = []
for combo in itertools.combinations(facet, idx):
these_entries = [min_entries[i] for i in combo]
these_combos.append(frozenset(these_entries))
combos.append(these_combos)
all_combos = set(itertools.chain.from_iterable(combos))
list_combos = []
for idx in all_combos:
list_combos.append(list(idx))
all_combos = list_combos
multi_entries = []
# Parallel processing of multi-entry generation
if nproc is not None:
func = partial(self.process_multientry, prod_comp=tot_comp)
with Pool(nproc) as proc_pool:
multi_entries = list(proc_pool.imap(func, all_combos))
multi_entries = list(filter(bool, multi_entries))
else:
# Serial processing of multi-entry generation
for combo in all_combos:
multi_entry = self.process_multientry(combo, prod_comp=tot_comp)
if multi_entry:
multi_entries.append(multi_entry)
return multi_entries
def _generate_multielement_entries(self, entries, nproc=None):
"""
Create entries for multi-element Pourbaix construction.
This works by finding all possible linear combinations
of entries that can result in the specified composition
from the initialized comp_dict.
Args:
entries ([PourbaixEntries]): list of Pourbaix entries
to process into MultiEntries
nproc (int): number of processes to be used in parallel
treatment of entry combos
"""
n_elems = len(self._elt_comp) # No. of elements
total_comp = Composition(self._elt_comp)
# generate all combinations of compounds that have all elements
entry_combos = [itertools.combinations(entries, idx + 1) for idx in range(n_elems)]
entry_combos = itertools.chain.from_iterable(entry_combos)
entry_combos = filter(lambda x: total_comp < MultiEntry(x).composition, entry_combos)
# Generate and filter entries
processed_entries = []
total = sum(comb(len(entries), idx + 1) for idx in range(n_elems))
if total > 1e6:
warnings.warn(f"Your Pourbaix diagram includes {total} entries and may take a long time to generate.")
# Parallel processing of multi-entry generation
if nproc is not None:
func = partial(self.process_multientry, prod_comp=total_comp)
with Pool(nproc) as proc_pool:
processed_entries = list(proc_pool.imap(func, entry_combos))
processed_entries = list(filter(bool, processed_entries))
# Serial processing of multi-entry generation
else:
for entry_combo in entry_combos:
processed_entry = self.process_multientry(entry_combo, total_comp)
if processed_entry is not None:
processed_entries.append(processed_entry)
return processed_entries
@staticmethod
def process_multientry(entry_list, prod_comp, coeff_threshold=1e-4):
"""
Static method for finding a multientry based on
a list of entries and a product composition.
Essentially checks to see if a valid aqueous
reaction exists between the entries and the
product composition and returns a MultiEntry
with weights according to the coefficients if so.
Args:
entry_list ([Entry]): list of entries from which to
create a MultiEntry
prod_comp (Composition): composition constraint for setting
weights of MultiEntry
coeff_threshold (float): threshold of stoichiometric
coefficients to filter, if weights are lower than
this value, the entry is not returned
"""
dummy_oh = [Composition("H"), Composition("O")]
try:
# Get balanced reaction coeffs, ensuring all < 0 or conc thresh
# Note that we get reduced compositions for solids and non-reduced
# compositions for ions because ions aren't normalized due to
# their charge state.
entry_comps = [e.composition for e in entry_list]
rxn = Reaction(entry_comps + dummy_oh, [prod_comp])
react_coeffs = [-coeff for coeff in rxn.coeffs[: len(entry_list)]]
all_coeffs = [*react_coeffs, rxn.get_coeff(prod_comp)]
# Check if reaction coeff threshold met for Pourbaix compounds
# All reactant/product coefficients must be positive nonzero
if all(coeff > coeff_threshold for coeff in all_coeffs):
return MultiEntry(entry_list, weights=react_coeffs)
return None
except ReactionError:
return None
@staticmethod
def get_pourbaix_domains(pourbaix_entries, limits=None):
"""
Returns a set of Pourbaix stable domains (i. e. polygons) in
pH-V space from a list of pourbaix_entries.
This function works by using scipy's HalfspaceIntersection
function to construct all of the 2-D polygons that form the
boundaries of the planes corresponding to individual entry
gibbs free energies as a function of pH and V. Hyperplanes
of the form a*pH + b*V + 1 - g(0, 0) are constructed and
supplied to HalfspaceIntersection, which then finds the
boundaries of each Pourbaix region using the intersection
points.
Args:
pourbaix_entries ([PourbaixEntry]): Pourbaix entries
with which to construct stable Pourbaix domains
limits ([[float]]): limits in which to do the pourbaix
analysis
Returns:
Returns a dict of the form {entry: [boundary_points]}.
The list of boundary points are the sides of the N-1
dim polytope bounding the allowable ph-V range of each entry.
"""
if limits is None:
limits = [[-2, 16], [-4, 4]]
# Get hyperplanes
hyperplanes = [
np.array([-PREFAC * entry.npH, -entry.nPhi, 0, -entry.energy]) * entry.normalization_factor
for entry in pourbaix_entries
]
hyperplanes = np.array(hyperplanes)
hyperplanes[:, 2] = 1
max_contribs = np.max(np.abs(hyperplanes), axis=0)
g_max = np.dot(-max_contribs, [limits[0][1], limits[1][1], 0, 1])
# Add border hyperplanes and generate HalfspaceIntersection
border_hyperplanes = [
[-1, 0, 0, limits[0][0]],
[1, 0, 0, -limits[0][1]],
[0, -1, 0, limits[1][0]],
[0, 1, 0, -limits[1][1]],
[0, 0, -1, 2 * g_max],
]
hs_hyperplanes = np.vstack([hyperplanes, border_hyperplanes])
interior_point = [*np.average(limits, axis=1).tolist(), g_max]
hs_int = HalfspaceIntersection(hs_hyperplanes, np.array(interior_point))
# organize the boundary points by entry
pourbaix_domains = {entry: [] for entry in pourbaix_entries}
for intersection, facet in zip(hs_int.intersections, hs_int.dual_facets):
for v in facet:
if v < len(pourbaix_entries):
this_entry = pourbaix_entries[v]
pourbaix_domains[this_entry].append(intersection)
# Remove entries with no Pourbaix region
pourbaix_domains = {k: v for k, v in pourbaix_domains.items() if v}
pourbaix_domain_vertices = {}
for entry, points in pourbaix_domains.items():
points = np.array(points)[:, :2]
# Initial sort to ensure consistency
points = points[np.lexsort(np.transpose(points))]
center = np.average(points, axis=0)
points_centered = points - center
# Sort points by cross product of centered points,
# isn't strictly necessary but useful for plotting tools
points_centered = sorted(points_centered, key=cmp_to_key(lambda x, y: x[0] * y[1] - x[1] * y[0]))
points = points_centered + center
# Create simplices corresponding to Pourbaix boundary
simplices = [Simplex(points[indices]) for indices in ConvexHull(points).simplices]
pourbaix_domains[entry] = simplices
pourbaix_domain_vertices[entry] = points
return pourbaix_domains, pourbaix_domain_vertices
def find_stable_entry(self, pH, V):
"""
Finds stable entry at a pH,V condition
Args:
pH (float): pH to find stable entry
V (float): V to find stable entry.
Returns:
PourbaixEntry: stable entry at pH, V
"""
energies_at_conditions = [e.normalized_energy_at_conditions(pH, V) for e in self.stable_entries]
return self.stable_entries[np.argmin(energies_at_conditions)]
def get_decomposition_energy(self, entry, pH, V):
"""
Finds decomposition to most stable entries in eV/atom,
supports vectorized inputs for pH and V.
Args:
entry (PourbaixEntry): PourbaixEntry corresponding to
compound to find the decomposition for
pH (float, [float]): pH at which to find the decomposition
V (float, [float]): voltage at which to find the decomposition
Returns:
Decomposition energy for the entry, i. e. the energy above
the "Pourbaix hull" in eV/atom at the given conditions
"""
# Check composition consistency between entry and Pourbaix diagram:
pbx_comp = Composition(self._elt_comp).fractional_composition
entry_pbx_comp = Composition(
{elt: coeff for elt, coeff in entry.composition.items() if elt not in ELEMENTS_HO}
).fractional_composition
if entry_pbx_comp != pbx_comp:
raise ValueError("Composition of stability entry does not match Pourbaix Diagram")
entry_normalized_energy = entry.normalized_energy_at_conditions(pH, V)
hull_energy = self.get_hull_energy(pH, V)
decomposition_energy = entry_normalized_energy - hull_energy
# Convert to eV/atom instead of eV/normalized formula unit
decomposition_energy /= entry.normalization_factor
decomposition_energy /= entry.composition.num_atoms
return decomposition_energy
def get_hull_energy(self, pH, V):
"""
Gets the minimum energy of the Pourbaix "basin" that is formed
from the stable Pourbaix planes. Vectorized.
Args:
pH (float or [float]): pH at which to find the hull energy
V (float or [float]): V at which to find the hull energy
Returns:
np.array: minimum Pourbaix energy at conditions
"""
all_gs = np.array([e.normalized_energy_at_conditions(pH, V) for e in self.stable_entries])
return np.min(all_gs, axis=0)
def get_stable_entry(self, pH, V):
"""
Gets the stable entry at a given pH, V condition.
Args:
pH (float): pH at a given condition
V (float): V at a given condition
Returns:
PourbaixEntry | MultiEntry: Pourbaix or multi-entry
corresponding to the minimum energy entry at a given pH, V condition
"""
all_gs = np.array([e.normalized_energy_at_conditions(pH, V) for e in self.stable_entries])
return self.stable_entries[np.argmin(all_gs)]
@property
def stable_entries(self):
"""Returns the stable entries in the Pourbaix diagram."""
return list(self._stable_domains)
@property
def unstable_entries(self):
"""Returns all unstable entries in the Pourbaix diagram."""
return [e for e in self.all_entries if e not in self.stable_entries]
@property
def all_entries(self):
"""Return all entries used to generate the Pourbaix diagram."""
return self._processed_entries
@property
def unprocessed_entries(self):
"""Return unprocessed entries."""
return self._unprocessed_entries
def as_dict(self):
"""
Returns:
MSONable dict.
"""
return {
"@module": type(self).__module__,
"@class": type(self).__name__,
"entries": [e.as_dict() for e in self._unprocessed_entries],
"comp_dict": self._elt_comp,
"conc_dict": self._conc_dict,
"filter_solids": self.filter_solids,
}
@classmethod
def from_dict(cls, d):
"""
Args:
d (): Dict representation.
Returns:
PourbaixDiagram
"""
decoded_entries = MontyDecoder().process_decoded(d["entries"])
return cls(decoded_entries, d.get("comp_dict"), d.get("conc_dict"), d.get("filter_solids"))
class PourbaixPlotter:
"""A plotter class for phase diagrams."""
def __init__(self, pourbaix_diagram):
"""
Args:
pourbaix_diagram (PourbaixDiagram): A PourbaixDiagram object.
"""
self._pbx = pourbaix_diagram
def show(self, *args, **kwargs):
"""
Shows the Pourbaix plot.
Args:
*args: args to get_pourbaix_plot
**kwargs: kwargs to get_pourbaix_plot
"""
plt = self.get_pourbaix_plot(*args, **kwargs)
plt.show()
@no_type_check
def get_pourbaix_plot(
self,
limits: tuple[float, float] | None = None,
title: str = "",
label_domains: bool = True,
label_fontsize: int = 20,
show_water_lines: bool = True,
show_neutral_axes: bool = True,
ax: plt.Axes = None,
) -> plt.Axes:
"""
Plot Pourbaix diagram.
Args:
limits: 2D list containing limits of the Pourbaix diagram
of the form [[xlo, xhi], [ylo, yhi]]
title (str): Title to display on plot
label_domains (bool): whether to label Pourbaix domains
label_fontsize: font size for domain labels
show_water_lines: whether to show dashed lines indicating the region
of water stability.
show_neutral_axes; whether to show dashed horizontal and vertical lines
at 0 V and pH 7, respectively.
ax (Axes): Matplotlib Axes instance for plotting
Returns:
Axes: matplotlib Axes object with Pourbaix diagram
"""
if limits is None:
limits = [[-2, 16], [-3, 3]]
ax = ax or pretty_plot(16)
xlim, ylim = limits
lw = 3
if show_water_lines:
h_line = np.transpose([[xlim[0], -xlim[0] * PREFAC], [xlim[1], -xlim[1] * PREFAC]])
o_line = np.transpose([[xlim[0], -xlim[0] * PREFAC + 1.23], [xlim[1], -xlim[1] * PREFAC + 1.23]])
ax.plot(h_line[0], h_line[1], "r--", linewidth=lw)
ax.plot(o_line[0], o_line[1], "r--", linewidth=lw)
if show_neutral_axes:
neutral_line = np.transpose([[7, ylim[0]], [7, ylim[1]]])
V0_line = np.transpose([[xlim[0], 0], [xlim[1], 0]])
ax.plot(neutral_line[0], neutral_line[1], "k-.", linewidth=lw)
ax.plot(V0_line[0], V0_line[1], "k-.", linewidth=lw)
for entry, vertices in self._pbx._stable_domain_vertices.items():
center = np.average(vertices, axis=0)
x, y = np.transpose(np.vstack([vertices, vertices[0]]))
ax.plot(x, y, "k-", linewidth=lw)
if label_domains:
ax.annotate(
generate_entry_label(entry),
center,
ha="center",
va="center",
fontsize=label_fontsize,