Coverage for /builds/ase/ase/ase/pourbaix.py: 94.24%

1# fmt: off 

2 

3import functools 

4import re 

5from collections import Counter 

6from dataclasses import dataclass 

7from fractions import Fraction 

8from itertools import chain, combinations, product 

9from typing import List, Tuple, Union 

10 

11import numpy as np 

12from scipy.linalg import null_space 

13 

14from ase.formula import Formula 

15from ase.units import kB 

16 

17CONST = kB * np.log(10) # Nernst constant 

18 

19PREDEF_ENERGIES = { # Default chemical potentials 

20 'H+': 0.0, # for water, protons and electrons 

21 'e-': 0.0, 

22 'H2O': -2.4583 

23} 

24 

25U_STD_AGCL = 0.222 # Standard redox potential of AgCl electrode 

26U_STD_SCE = 0.244 # Standard redox potential of SCE electrode 

27 

28 

29def parse_formula(formula: str, fmt: str = 'metal'): 

30 aq = formula.endswith('(aq)') 

31 charge = formula.count('+') - formula.count('-') 

32 formula_strip = formula.replace('(aq)', '').rstrip('+-') 

33 formula_obj = Formula(formula_strip, format=fmt) 

34 return formula_obj, charge, aq 

35 

36 

37def initialize_refs(refs_dct, reduce=False, fmt='metal'): 

38 """Convert dictionary entries to Species instances""" 

39 refs = {} 

40 for label, energy in refs_dct.items(): 

41 spec = Species.from_string(label, energy, reduce, fmt) 

42 refs[label] = spec 

43 return refs 

44 

45 

46def get_product_combos(reactant, refs): 

47 """Obtain all possible combinations of products. 

48 

49 Obtain - from the available references and based 

50 on the stoichiometry of the target material (reactant) - 

51 different combinations of products to be inserted 

52 in (electro)chemical reactions. 

53 """ 

54 ref_elem = reactant._main_elements 

55 allcombos = [[] for _ in range(len(ref_elem))] 

56 for ref in refs.values(): 

57 contained = ref.contains(ref_elem) 

58 for w in np.argwhere(contained).flatten(): 

59 allcombos[w].append(ref) 

60 return [np.unique(combo) for combo in product(*allcombos)] 

61 

62 

63def get_phases(reactant, refs, T, conc, reference, tol=1e-3): 

64 """Obtain all the possible decomposition pathways 

65 for a given reactant as a collection of RedOx objects. 

66 """ 

67 phases = [] 

68 phase_matrix = [] 

69 

70 for products in get_product_combos(reactant, refs): 

71 phase = RedOx.from_species( 

72 reactant, products, T, conc, reference, tol 

73 ) 

74 if phase is not None: 

75 phases.append(phase) 

76 phase_matrix.append(phase._vector) 

77 

78 if len(phase_matrix) == 0: 

79 raise ValueError( 

80 "No valid decomposition pathways have been found" + 

81 " given this set of references." 

82 ) 

83 

84 return phases, np.array(phase_matrix).astype('float64') 

85 

86 

87def get_main_products(count): 

88 """Obtain the reaction products excluded protons, 

89 water and electrons. 

90 """ 

91 return [spec for spec, coef in count.items() 

92 if coef > 0 and spec not in ['H+', 'H2O', 'e-']] 

93 

94 

95def format_label(count) -> str: 

96 """Obtain phase labels formatted in LaTeX math style.""" 

97 formatted = [] 

98 for prod in get_main_products(count): 

99 label = re.sub(r'(\S)([+-]+)', r'\1$^{\2}$', prod) 

100 label = re.sub(r'(\d+)', r'$_{\1}$', label) 

101 for symbol in ['+', '-']: 

102 count = label.count(symbol) 

103 if count > 1: 

104 label = label.replace(count * symbol, f'{count}{symbol}') 

105 if count == 1: 

106 label = label.replace(count * symbol, symbol) 

107 formatted.append(label) 

108 return ', '.join(f for f in formatted) 

109 

110 

111def make_coeff_nice(coeff, max_denom) -> str: 

112 """Convert a fraction into a string while limiting the denominator""" 

113 frac = abs(Fraction(coeff).limit_denominator(max_denom)) 

114 if frac.numerator == frac.denominator: 

115 return '' 

116 return str(frac) 

117 

118 

119def add_numbers(ax, text) -> None: 

120 """Add number identifiers to the different domains of a Pourbaix diagram.""" 

121 import matplotlib.patheffects as pfx 

122 for i, (x, y, _) in enumerate(text): 

123 txt = ax.text( 

124 y, x, f'{i}', 

125 fontsize=20, 

126 horizontalalignment='center' 

127 ) 

128 txt.set_path_effects([pfx.withStroke(linewidth=2.0, foreground='w')]) 

129 

130 

131def add_labels(ax, text) -> None: 

132 """Add phase labels to the different domains of a Pourbaix diagram.""" 

133 import matplotlib.patheffects as pfx 

134 for i, (x, y, species) in enumerate(text): 

135 label = format_label(species) 

136 annotation = ax.annotate( 

137 label, xy=(y, x), color='w', 

138 fontsize=16, horizontalalignment='center' 

139 ) 

140 annotation.set_path_effects( 

141 [pfx.withStroke(linewidth=2.0, foreground='k')] 

142 ) 

143 annotation.draggable() 

144 ax.add_artist(annotation) 

145 

146 

147def wrap_text(text) -> str: 

148 import textwrap 

149 

150 textlines = [] 

151 for i, (_, _, species) in enumerate(text): 

152 label = format_label(species) 

153 textlines.append( 

154 textwrap.fill( 

155 f'({i}) {label}', 

156 width=40, 

157 subsequent_indent=' ' 

158 ) 

159 ) 

160 return '\n'.join(textlines) 

161 

162 

163def add_phase_labels(fig, text, offset=0.0): 

164 """Add phase labels to the right of the diagram""" 

165 

166 fig.text( 

167 0.75 + offset, 0.5, 

168 wrap_text(text), 

169 fontsize=16, 

170 va='center', 

171 ha='left') 

172 

173 

174def add_redox_lines(axes, pH, reference, color='k') -> None: 

175 """Add water redox potentials to a Pourbaix diagram""" 

176 const = -0.5 * PREDEF_ENERGIES['H2O'] 

177 corr = { 

178 'SHE': 0, 

179 'RHE': 0, 

180 'Pt': 0, 

181 'AgCl': -U_STD_AGCL, 

182 'SCE': -U_STD_SCE, 

183 } 

184 kwargs = { 

185 'c': color, 

186 'ls': '--', 

187 'zorder': 2 

188 } 

189 if reference in ['SHE', 'AgCl', 'SCE']: 

190 slope = -59.2e-3 

191 axes.plot(pH, slope * pH + corr[reference], **kwargs) 

192 axes.plot(pH, slope * pH + const + corr[reference], **kwargs) 

193 elif reference in ['Pt', 'RHE']: 

194 axes.axhline(0 + corr[reference], **kwargs) 

195 axes.axhline(const + corr[reference], **kwargs) 

196 else: 

197 raise ValueError('The specified reference electrode doesnt exist') 

198 

199 

200@functools.total_ordering 

201class Species: 

202 """Class representing an individual chemical species, 

203 grouping relevant properties. 

204 

205 Initialization 

206 -------------- 

207 name: str 

208 A label representing the species 

209 

210 formula: Formula 

211 

212 charge: int 

213 the electric charge of the species, if ionic 

214 

215 aq: bool 

216 whether the species is solid (False) 

217 or acqueous (True) 

218 

219 energy: float 

220 the chemical potential of the species 

221 """ 

222 def __init__(self, 

223 name: str, 

224 formula: Formula, 

225 charge: int, 

226 aq: bool, 

227 energy: float): 

228 

229 self.name = name 

230 self.formula = formula 

231 self.energy = energy 

232 self.charge = charge 

233 self.aq = aq 

234 

235 self.count = formula.count() 

236 self.natoms = sum(self.count.values()) 

237 self._main_elements = [ 

238 e for e in self.count.keys() if e not in ['H', 'O'] 

239 ] 

240 

241 def __eq__(self, other): 

242 if not isinstance(other, Species): 

243 return NotImplemented 

244 return self.name == other.name 

245 

246 def __lt__(self, other): 

247 if not isinstance(other, Species): 

248 return NotImplemented 

249 return self.name < other.name 

250 

251 @classmethod 

252 def from_string(cls, string: str, energy: float, 

253 reduce: bool = True, fmt: str = 'metal'): 

254 """Initialize the class provided a formula and an energy. 

255 

256 string: str 

257 The chemical formula of the species (e.g. ``ZnO``). 

258 For solid species, the formula is automatically reduced to the 

259 unit formula, and the chemical potential normalized accordingly. 

260 Acqueous species are specified by expliciting 

261 all the positive or negative charges and by appending ``(aq)``. 

262 Parentheses for grouping functional groups are acceptable. 

263 e.g. 

264 Be3(OH)3[3+] ➜ wrong 

265 Be3(OH)3+++ ➜ wrong 

266 Be3(OH)3+++(aq) ➜ correct 

267 Be3O3H3+++(aq) ➜ correct 

268 

269 energy: float 

270 the energy (chemical potential) associated with the species. 

271 

272 reduce: bool 

273 reduce to the unit formula and normalize the energy accordingly. 

274 Formulae and energies of acqueous species are never reduced. 

275 

276 fmt: str 

277 Formula formatting according to the available options in 

278 ase.formula.Formula 

279 """ 

280 formula, charge, aq = parse_formula(string, fmt=fmt) 

281 if not aq: 

282 if reduce: 

283 formula, n_fu = formula.reduce() 

284 energy /= n_fu 

285 name = str(formula) 

286 else: 

287 name = string 

288 return cls(name, formula, charge, aq, energy) 

289 

290 def get_chemsys(self): 

291 """Get the possible combinations of elements based 

292 on the stoichiometry. Useful for database queries. 

293 """ 

294 elements = set(self.count.keys()) 

295 elements.update(['H', 'O']) 

296 chemsys = list( 

297 chain.from_iterable( 

298 [combinations(elements, i + 1) 

299 for i in range(len(elements))] 

300 ) 

301 ) 

302 return chemsys 

303 

304 def balance_electrochemistry(self): 

305 """Obtain number of H2O, H+, e- "carried" by the species 

306 in electrochemical reactions. 

307 """ 

308 n_H2O = -self.count.get('O', 0) 

309 n_H = -2 * n_H2O - self.count.get('H', 0) 

310 n_e = n_H + self.charge 

311 return n_H2O, n_H, n_e 

312 

313 def _count_array(self, elements): 

314 return np.array([self.count.get(e, 0) for e in elements]) 

315 

316 def contains(self, elements): 

317 return [elem in self._main_elements for elem in elements] 

318 

319 def get_fractional_composition(self, elements): 

320 """Obtain the fractional content of each element.""" 

321 N_all = sum(self.count.values()) 

322 N_elem = sum([self.count.get(e, 0) for e in elements]) 

323 return N_elem / N_all 

324 

325 def __repr__(self): 

326 return f'Species({self.name})' 

327 

328 

329class RedOx: 

330 def __init__(self, species, coeffs, 

331 T=298.15, conc=1e-6, 

332 reference='SHE'): 

333 """RedOx class representing an (electro)chemical reaction. 

334 

335 Initialization 

336 -------------- 

337 

338 species: list[Species] 

339 The reactant and products excluded H2O, protons and electrons. 

340 

341 coeffs: list[float] 

342 The stoichiometric coefficients of the above species. 

343 Positive coefficients are associated with the products, 

344 negative coefficients with the reactants. 

345 

346 T: float 

347 The temperature in Kelvin. Default: 298.15 K. 

348 

349 conc: float 

350 The concentration of ionic species. Default: 1.0e-6 M. 

351 

352 reference: str 

353 The reference electrode. Default: SHE. 

354 available options: SHE, RHE, AgCl, Pt. 

355 

356 

357 Relevant methods 

358 ---------------- 

359 

360 from_species(reactant, products, T, conc, reference) 

361 Initialize the class from the reactant as a Species object 

362 and the product(s) a list of Species objects. 

363 

364 equation() 

365 Print the chemical equation of the reaction. 

366 

367 get_free_energy(U, pH): 

368 Obtain the reaction free energy at a given 

369 applied potential (U) and pH. 

370 

371 """ 

372 self.species = species 

373 self.coeffs = coeffs 

374 self.T = T 

375 self.conc = conc 

376 self.reference = reference 

377 self.count = Counter() 

378 

379 alpha = CONST * T # 0.059 eV @ T=298.15K 

380 const_term = 0 

381 pH_term = 0 

382 U_term = 0 

383 

384 for spec, coef in zip(species, coeffs): 

385 self.count[spec.name] = coef 

386 amounts = spec.balance_electrochemistry() 

387 

388 const_term += coef * ( 

389 spec.energy + alpha * (spec.aq * np.log10(conc)) 

390 ) 

391 pH_term += - coef * alpha * amounts[1] 

392 U_term += - coef * amounts[2] 

393 

394 for name, n in zip(['H2O', 'H+', 'e-'], amounts): 

395 const_term += coef * n * PREDEF_ENERGIES[name] 

396 self.count[name] += coef * n 

397 

398 const_corr, pH_corr = self.get_ref_correction(reference, alpha) 

399 self._vector = [ 

400 float(const_term + const_corr), 

401 float(U_term), 

402 float(pH_term + pH_corr) 

403 ] 

404 

405 @classmethod 

406 def from_species(cls, reactant, products, 

407 T: float = 298.15, conc: float = 1e-6, 

408 reference: str = 'SHE', tol: float = 1e-3): 

409 """Initialize the class from a combination of 

410 reactant and products. The stoichiometric 

411 coefficients are automatically determined. 

412 

413 reactant: Species 

414 products: list[Species] 

415 """ 

416 reac_elem = reactant._main_elements 

417 reac_count = [-reactant._count_array(reac_elem)] 

418 prod_count = [p._count_array(reac_elem) for p in products] 

419 elem_matrix = np.array(reac_count + prod_count).T 

420 coeffs = null_space(elem_matrix).flatten() 

421 

422 if len(coeffs) > 0 and all(coeffs > tol): 

423 coeffs /= coeffs[0] 

424 coeffs[0] = -coeffs[0] 

425 species = (reactant, *products) 

426 phase = cls(species, coeffs, T, conc, reference) 

427 return phase 

428 

429 return None 

430 

431 def get_ref_correction(self, reference: str, 

432 alpha: float) -> Tuple[float, float]: 

433 """Correct the constant and pH contributions to the reaction free energy 

434 based on the reference electrode of choice and the temperature 

435 (alpha=k_B*T*ln(10)) 

436 """ 

437 n_e = self.count['e-'] 

438 gibbs_corr = 0.0 

439 pH_corr = 0.0 

440 if reference in ['RHE', 'Pt']: 

441 pH_corr += n_e * alpha 

442 if reference == 'Pt' and n_e < 0: 

443 gibbs_corr += n_e * 0.5 * PREDEF_ENERGIES['H2O'] 

444 if reference == 'AgCl': 

445 gibbs_corr -= n_e * U_STD_AGCL 

446 if reference == 'SCE': 

447 gibbs_corr -= n_e * U_STD_SCE 

448 

449 return gibbs_corr, pH_corr 

450 

451 def equation(self, max_denom: int = 50) -> str: 

452 """Print the chemical reaction.""" 

453 

454 reactants = [] 

455 products = [] 

456 for s, n in self.count.items(): 

457 if abs(n) <= 1e-6: 

458 continue 

459 nice_coeff = make_coeff_nice(n, max_denom) 

460 substr = f'{nice_coeff}{s}' 

461 if n > 0: 

462 products.append(substr) 

463 else: 

464 reactants.append(substr) 

465 

466 return " ➜ ".join([" + ".join(reactants), " + ".join(products)]) 

467 

468 def get_free_energy(self, U: float, pH: float) -> float: 

469 """Evaluate the reaction free energy 

470 at a given applied potential U and pH""" 

471 return self._vector[0] + self._vector[1] * U + self._vector[2] * pH 

472 

473 

474class Pourbaix: 

475 """Pourbaix class for acqueous stability evaluations. 

476 

477 Allows to determine the most stable phase in a given set 

478 of pH and potential conditions and to evaluate a complete diagram. 

479 

480 Initialization 

481 -------------- 

482 

483 material_name: str 

484 The formula of the target material. It is preferrable 

485 to provide the reduced formula (e.g. RuO2 instad of Ru2O4). 

486 

487 refs_dct: dict 

488 A dictionary containing the formulae of the target material 

489 and its competing phases (solid and/or ionic) as keys, 

490 and their formation energies as values. 

491 

492 T: float 

493 Temperature in Kelvin. Default: 298.15 K. 

494 

495 conc: float 

496 Concentration of the ionic species. Default: 1e-6 mol/L. 

497 

498 reference: str 

499 The reference electrode. Default: SHE. 

500 available options: SHE, RHE, AgCl, Pt. 

501 

502 

503 Relevant methods 

504 ---------------- 

505 

506 get_pourbaix_energy(U, pH) 

507 obtain the energy of the target material 

508 relative to the most stable phase at a given potential U and pH. 

509 If negative, the target material can be regarded as stable. 

510 plot(...) 

511 plot a complete Pourbaix diagram in a given pH and potential window. 

512 

513 

514 Relevant attributes 

515 ------------------- 

516 

517 material: Species 

518 the target material as a Species object 

519 

520 phases: list[RedOx] 

521 the available decomposition reactions of the target material 

522 into its competing phases as a list of RedOx objects. 

523 

524 """ 

525 def __init__(self, 

526 material_name: str, 

527 refs_dct: dict, 

528 T: float = 298.15, 

529 conc: float = 1.0e-6, 

530 reference: str = 'SHE'): 

531 

532 self.material_name = material_name 

533 self.refs = refs_dct 

534 self.T = T 

535 self.conc = conc 

536 self.reference = reference 

537 

538 refs = initialize_refs(refs_dct) 

539 self.material = refs.pop(material_name) 

540 self.phases, phase_matrix = get_phases( 

541 self.material, refs, T, conc, reference 

542 ) 

543 self._const = phase_matrix[:, 0] 

544 self._var = phase_matrix[:, 1:] 

545 

546 def _decompose(self, U, pH): 

547 """Evaluate the reaction energy for decomposing 

548 the target material into each of the available products 

549 at a given pH and applied potential. 

550 """ 

551 return self._const + np.dot(self._var, [U, pH]) 

552 

553 def _get_pourbaix_energy(self, U, pH): 

554 """Evaluate the Pourbaix energy""" 

555 energies = self._decompose(U, pH) 

556 i_min = np.argmin(energies) 

557 return -energies[i_min], i_min 

558 

559 def get_pourbaix_energy(self, U, pH, verbose=False): 

560 """Evaluate the Pourbaix energy, print info about 

561 the most stable phase and the corresponding energy at 

562 a given potential U and pH. 

563 

564 The Pourbaix energy represents the energy of the target material 

565 relative to the most stable competing phase. If negative, 

566 the target material can be considered as stable. 

567 """ 

568 energy, index = self._get_pourbaix_energy(U, pH) 

569 phase = self.phases[index] 

570 if verbose: 

571 if energy <= 0.0: 

572 print(f'{self.material.name} is stable.') 

573 else: 

574 print(f'Stable phase: \n{phase.equation()}') 

575 print(f'Energy: {energy:.3f} eV') 

576 return energy, phase 

577 

578 def get_equations(self, contains: Union[str, None] = None): 

579 """Print the chemical reactions of the available phases. 

580 

581 the argument `contains' allows to filter for phases containing a 

582 particular species 

583 e.g. get_equations(contains='ZnO') 

584 """ 

585 equations = [] 

586 for i, p in enumerate(self.phases): 

587 if contains is not None and contains not in p.count: 

588 continue 

589 equations.append(f'{i}) {p.equation()}') 

590 return equations 

591 

592 def diagram(self, U=None, pH=None): 

593 """Actual evaluation of the complete diagram 

594 

595 Returns 

596 ------- 

597 

598 pour: 

599 The stability domains of the diagram on the pH vs. U grid. 

600 domains are represented by indexes (as integers) 

601 that map to Pourbaix.phases 

602 the target material is stable (index=-1) 

603 

604 meta: 

605 The Pourbaix energy on the pH vs. U grid. 

606 

607 text: 

608 The coordinates and phases information for 

609 text placement on the diagram. 

610 

611 domains: 

612 ... 

613 

614 """ 

615 pour = np.zeros((len(U), len(pH))) 

616 meta = pour.copy() 

617 

618 for i, u in enumerate(U): 

619 for j, p in enumerate(pH): 

620 meta[i, j], pour[i, j] = self._get_pourbaix_energy(u, p) 

621 

622 # Identifying the region where the target material 

623 # is stable and updating the diagram accordingly 

624 where_stable = (meta <= 0) 

625 pour[where_stable] = -1 

626 

627 text = [] 

628 domains = [int(i) for i in np.unique(pour)] 

629 for phase_id in domains: 

630 if phase_id == -1: 

631 where = where_stable 

632 txt = {self.material.name: 1} 

633 else: 

634 where = (pour == phase_id) 

635 phase = self.phases[int(phase_id)] 

636 txt = phase.count 

637 x = np.dot(where.sum(1), U) / where.sum() 

638 y = np.dot(where.sum(0), pH) / where.sum() 

639 text.append((x, y, txt)) 

640 

641 return PourbaixDiagram(self, U, pH, pour, meta, text, domains) 

642 

643 def get_phase_boundaries(self, phmin, phmax, umin, umax, domains, tol=1e-6): 

644 """Plane intersection method for finding 

645 the boundaries between phases seen in the final plot. 

646 

647 Returns a list of tuples, each representing a single boundary, 

648 of the form ([[x1, x2], [y1, y2]], [id1, id2]), namely the 

649 x and y coordinates of the simplices connected by the boundary 

650 and the id's of the phases at each side of the boundary. 

651 """ 

652 from collections import defaultdict 

653 from itertools import combinations 

654 

655 # Planes identifying the diagram frame 

656 planes = [(np.array([0.0, 1.0, 0.0]), umin, 'bottom'), 

657 (np.array([0.0, 1.0, 0.0]), umax, 'top'), 

658 (np.array([1.0, 0.0, 0.0]), phmin, 'left'), 

659 (np.array([1.0, 0.0, 0.0]), phmax, 'right')] 

660 

661 # Planes associated with the stable domains of the diagram. 

662 # Given x=pH, y=U, z=E_pbx=-DeltaG, each plane has expression: 

663 # _vector[2]*x + _vector[1]*y + z = -_vector[0] 

664 # The region where the target material is stable 

665 # (id=-1, if present) is delimited by the xy plane (Epbx=0) 

666 for d in domains: 

667 if d == -1: 

668 plane = np.array([0.0, 0.0, 1.0]) 

669 const = 0.0 

670 else: 

671 pvec = self.phases[d]._vector 

672 plane = np.array([pvec[2], pvec[1], 1]) 

673 const = -pvec[0] 

674 planes.append((plane, const, d)) 

675 

676 # The simplices are found from the intersection points between 

677 # all possible plane triplets. If the z coordinate of the point 

678 # matches the corresponding pourbaix energy, 

679 # then the point is a simplex. 

680 simplices = [] 

681 for (p1, c1, id1), (p2, c2, id2), (p3, c3, id3) in \ 

682 combinations(planes, 3): 

683 A = np.vstack((p1, p2, p3)) 

684 c = np.array([c1, c2, c3]) 

685 ids = (id1, id2, id3) 

686 

687 try: 

688 # triplets containing parallel planes raise a LinAlgError 

689 # and are automatically excluded. 

690 invA = np.linalg.inv(A) 

691 except np.linalg.LinAlgError: 

692 continue 

693 

694 pt = np.dot(invA, c) 

695 Epbx = self._get_pourbaix_energy(pt[1], pt[0])[0] 

696 if pt[2] >= -tol and \ 

697 np.isclose(Epbx, pt[2], rtol=0, atol=tol) and \ 

698 phmin - tol <= pt[0] <= phmax + tol and \ 

699 umin - tol <= pt[1] <= umax + tol: 

700 simplex = np.round(pt[:2], 3) 

701 simplices.append((simplex, ids)) 

702 

703 # The final segments to plot on the diagram are found from 

704 # the pairs of unique simplices that have two neighboring phases 

705 # in common, diagram frame excluded. 

706 duplicate_filter = defaultdict(int) 

707 

708 def are_boundaries_there(comm): 

709 bounds = ['bottom', 'top', 'left', 'right'] 

710 for c in comm: 

711 if c in bounds: 

712 return True 

713 return False 

714 

715 for (s1, id1), (s2, id2) in combinations(simplices, 2): 

716 

717 # common neighboring phases, diagram frame excluded 

718 common = {*id1} & {*id2} 

719 

720 simplices_are_distinct = not np.allclose(s1, s2, rtol=0, atol=tol) 

721 only_two_phases_in_common = (len(common) == 2) 

722 diagram_frame_excluded = not are_boundaries_there(common) 

723 

724 # Filtering out duplicates 

725 if all([simplices_are_distinct, only_two_phases_in_common, 

726 diagram_frame_excluded]): 

727 testarray = sorted([s1, s2], key=lambda s: (s[0], s[1])) 

728 testarray.append(sorted(common)) 

729 testarray = np.array(testarray).flatten() 

730 duplicate_filter[tuple(testarray)] += 1 

731 

732 segments = [] 

733 for segment in duplicate_filter: 

734 coords, phases = np.split(np.array(segment), [4]) 

735 segments.append(( 

736 coords.reshape(2, 2).T, 

737 list(phases.astype(int)) 

738 )) 

739 return segments 

740 

741 

742@dataclass 

743class PourbaixDiagram: 

744 pbx: Pourbaix 

745 U: np.ndarray 

746 pH: np.ndarray 

747 pour: np.ndarray 

748 meta: np.ndarray 

749 text: List[Tuple[float, float, str]] 

750 domains: List[int] 

751 

752 def __post_init__(self): 

753 def _issorted(array): 

754 return all(np.diff(array) > 0) 

755 

756 # We might as well require the input domains to be sorted: 

757 if not _issorted(self.U): 

758 raise ValueError('U must be sorted') 

759 

760 if not _issorted(self.pH): 

761 raise ValueError('pH must be sorted') 

762 

763 @property 

764 def pHrange(self): 

765 return (self.pH[0], self.pH[-1]) 

766 

767 @property 

768 def Urange(self): 

769 return (self.U[0], self.U[-1]) 

770 

771 def _draw_diagram_axes( 

772 self, 

773 cap, 

774 normalize, 

775 include_text, 

776 include_water, 

777 labeltype, cmap, *, 

778 ax): 

779 """Backend for drawing Pourbaix diagrams.""" 

780 

781 meta = self.meta.copy() 

782 if normalize: 

783 meta /= self.pbx.material.natoms 

784 cbarlabel = r'$\Delta G_{pbx}$ (eV/atom)' 

785 else: 

786 cbarlabel = r'$\Delta G_{pbx}$ (eV)' 

787 

788 fig = ax.get_figure() 

789 extent = [*self.pHrange, *self.Urange] 

790 

791 fig.subplots_adjust( 

792 left=0.1, right=0.97, 

793 top=0.97, bottom=0.14 

794 ) 

795 

796 if isinstance(cap, list): 

797 vmin = cap[0] 

798 vmax = cap[1] 

799 else: 

800 vmin = -cap 

801 vmax = cap 

802 

803 colorplot = ax.imshow( 

804 meta, cmap=cmap, 

805 extent=extent, 

806 vmin=vmin, vmax=vmax, 

807 origin='lower', aspect='auto', 

808 interpolation='gaussian' 

809 ) 

810 

811 cbar = fig.colorbar( 

812 colorplot, 

813 ax=ax, 

814 pad=0.02 

815 ) 

816 

817 bounds = self.pbx.get_phase_boundaries( 

818 *self.pHrange, *self.Urange, self.domains 

819 ) 

820 for coords, _ in bounds: 

821 ax.plot(coords[0], coords[1], '-', c='k', lw=1.0) 

822 

823 if labeltype == 'numbers': 

824 add_numbers(ax, self.text) 

825 elif labeltype == 'phases': 

826 add_labels(ax, self.text) 

827 elif labeltype is None: 

828 pass 

829 else: 

830 raise ValueError("The provided label type doesn't exist") 

831 

832 if include_water: 

833 add_redox_lines(ax, self.pH, self.pbx.reference, 'w')