Coverage for ase / pourbaix.py: 94.22%

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 

9 

10import numpy as np 

11from scipy.linalg import null_space 

12 

13from ase.formula import Formula 

14from ase.units import kB 

15 

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

17 

18PREDEF_ENERGIES = { # Default chemical potentials 

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

20 'e-': 0.0, 

21 'H2O': -2.4583 

22} 

23 

24U_STD_AGCL = 0.222 # Standard redox potential of AgCl electrode 

25U_STD_SCE = 0.244 # Standard redox potential of SCE electrode 

26 

27 

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

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

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

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

32 formula_obj = Formula(formula_strip, format=fmt) 

33 return formula_obj, charge, aq 

34 

35 

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

37 """Convert dictionary entries to Species instances""" 

38 refs = {} 

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

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

41 refs[label] = spec 

42 return refs 

43 

44 

45def get_product_combos(reactant, refs): 

46 """Obtain all possible combinations of products. 

47 

48 Obtain - from the available references and based 

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

50 different combinations of products to be inserted 

51 in (electro)chemical reactions. 

52 """ 

53 ref_elem = reactant._main_elements 

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

55 for ref in refs.values(): 

56 contained = ref.contains(ref_elem) 

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

58 allcombos[w].append(ref) 

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

60 

61 

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

63 """Obtain all the possible decomposition pathways 

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

65 """ 

66 phases = [] 

67 phase_matrix = [] 

68 

69 for products in get_product_combos(reactant, refs): 

70 phase = RedOx.from_species( 

71 reactant, products, T, conc, reference, tol 

72 ) 

73 if phase is not None: 

74 phases.append(phase) 

75 phase_matrix.append(phase._vector) 

76 

77 if len(phase_matrix) == 0: 

78 raise ValueError( 

79 "No valid decomposition pathways have been found" + 

80 " given this set of references." 

81 ) 

82 

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

84 

85 

86def get_main_products(count): 

87 """Obtain the reaction products excluded protons, 

88 water and electrons. 

89 """ 

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

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

92 

93 

94def format_label(count) -> str: 

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

96 formatted = [] 

97 for prod in get_main_products(count): 

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

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

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

101 count = label.count(symbol) 

102 if count > 1: 

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

104 if count == 1: 

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

106 formatted.append(label) 

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

108 

109 

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

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

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

113 if frac.numerator == frac.denominator: 

114 return '' 

115 return str(frac) 

116 

117 

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

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

120 import matplotlib.patheffects as pfx 

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

122 txt = ax.text( 

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

124 fontsize=20, 

125 horizontalalignment='center' 

126 ) 

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

128 

129 

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

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

132 import matplotlib.patheffects as pfx 

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

134 label = format_label(species) 

135 annotation = ax.annotate( 

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

137 fontsize=16, horizontalalignment='center' 

138 ) 

139 annotation.set_path_effects( 

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

141 ) 

142 annotation.draggable() 

143 ax.add_artist(annotation) 

144 

145 

146def wrap_text(text) -> str: 

147 import textwrap 

148 

149 textlines = [] 

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

151 label = format_label(species) 

152 textlines.append( 

153 textwrap.fill( 

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

155 width=40, 

156 subsequent_indent=' ' 

157 ) 

158 ) 

159 return '\n'.join(textlines) 

160 

161 

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

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

164 

165 fig.text( 

166 0.75 + offset, 0.5, 

167 wrap_text(text), 

168 fontsize=16, 

169 va='center', 

170 ha='left') 

171 

172 

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

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

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

176 corr = { 

177 'SHE': 0, 

178 'RHE': 0, 

179 'Pt': 0, 

180 'AgCl': -U_STD_AGCL, 

181 'SCE': -U_STD_SCE, 

182 } 

183 kwargs = { 

184 'c': color, 

185 'ls': '--', 

186 'zorder': 2 

187 } 

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

189 slope = -59.2e-3 

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

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

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

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

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

195 else: 

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

197 

198 

199@functools.total_ordering 

200class Species: 

201 """Class representing an individual chemical species, 

202 grouping relevant properties. 

203 

204 Initialization 

205 -------------- 

206 name: str 

207 A label representing the species 

208 

209 formula: Formula 

210 

211 charge: int 

212 the electric charge of the species, if ionic 

213 

214 aq: bool 

215 whether the species is solid (False) 

216 or acqueous (True) 

217 

218 energy: float 

219 the chemical potential of the species 

220 """ 

221 def __init__(self, 

222 name: str, 

223 formula: Formula, 

224 charge: int, 

225 aq: bool, 

226 energy: float): 

227 

228 self.name = name 

229 self.formula = formula 

230 self.energy = energy 

231 self.charge = charge 

232 self.aq = aq 

233 

234 self.count = formula.count() 

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

236 self._main_elements = [ 

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

238 ] 

239 

240 def __eq__(self, other): 

241 if not isinstance(other, Species): 

242 return NotImplemented 

243 return self.name == other.name 

244 

245 def __lt__(self, other): 

246 if not isinstance(other, Species): 

247 return NotImplemented 

248 return self.name < other.name 

249 

250 @classmethod 

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

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

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

254 

255 string: str 

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

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

258 unit formula, and the chemical potential normalized accordingly. 

259 Acqueous species are specified by expliciting 

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

261 Parentheses for grouping functional groups are acceptable. 

262 e.g. 

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

264 Be3(OH)3+++ ➜ wrong 

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

266 Be3O3H3+++(aq) ➜ correct 

267 

268 energy: float 

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

270 

271 reduce: bool 

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

273 Formulae and energies of acqueous species are never reduced. 

274 

275 fmt: str 

276 Formula formatting according to the available options in 

277 ase.formula.Formula 

278 """ 

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

280 if not aq: 

281 if reduce: 

282 formula, n_fu = formula.reduce() 

283 energy /= n_fu 

284 name = str(formula) 

285 else: 

286 name = string 

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

288 

289 def get_chemsys(self): 

290 """Get the possible combinations of elements based 

291 on the stoichiometry. Useful for database queries. 

292 """ 

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

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

295 chemsys = list( 

296 chain.from_iterable( 

297 [combinations(elements, i + 1) 

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

299 ) 

300 ) 

301 return chemsys 

302 

303 def balance_electrochemistry(self): 

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

305 in electrochemical reactions. 

306 """ 

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

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

309 n_e = n_H + self.charge 

310 return n_H2O, n_H, n_e 

311 

312 def _count_array(self, elements): 

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

314 

315 def contains(self, elements): 

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

317 

318 def get_fractional_composition(self, elements): 

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

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

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

322 return N_elem / N_all 

323 

324 def __repr__(self): 

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

326 

327 

328class RedOx: 

329 def __init__(self, species, coeffs, 

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

331 reference='SHE'): 

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

333 

334 Initialization 

335 -------------- 

336 

337 species: list[Species] 

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

339 

340 coeffs: list[float] 

341 The stoichiometric coefficients of the above species. 

342 Positive coefficients are associated with the products, 

343 negative coefficients with the reactants. 

344 

345 T: float 

346 The temperature in Kelvin. Default: 298.15 K. 

347 

348 conc: float 

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

350 

351 reference: str 

352 The reference electrode. Default: SHE. 

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

354 

355 

356 Relevant methods 

357 ---------------- 

358 

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

360 Initialize the class from the reactant as a Species object 

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

362 

363 equation() 

364 Print the chemical equation of the reaction. 

365 

366 get_free_energy(U, pH): 

367 Obtain the reaction free energy at a given 

368 applied potential (U) and pH. 

369 

370 """ 

371 self.species = species 

372 self.coeffs = coeffs 

373 self.T = T 

374 self.conc = conc 

375 self.reference = reference 

376 self.count = Counter() 

377 

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

379 const_term = 0 

380 pH_term = 0 

381 U_term = 0 

382 

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

384 self.count[spec.name] = coef 

385 amounts = spec.balance_electrochemistry() 

386 

387 const_term += coef * ( 

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

389 ) 

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

391 U_term += - coef * amounts[2] 

392 

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

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

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

396 

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

398 self._vector = [ 

399 float(const_term + const_corr), 

400 float(U_term), 

401 float(pH_term + pH_corr) 

402 ] 

403 

404 @classmethod 

405 def from_species(cls, reactant, products, 

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

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

408 """Initialize the class from a combination of 

409 reactant and products. The stoichiometric 

410 coefficients are automatically determined. 

411 

412 reactant: Species 

413 products: list[Species] 

414 """ 

415 reac_elem = reactant._main_elements 

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

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

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

419 coeffs = null_space(elem_matrix).flatten() 

420 

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

422 coeffs /= coeffs[0] 

423 coeffs[0] = -coeffs[0] 

424 species = (reactant, *products) 

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

426 return phase 

427 

428 return None 

429 

430 def get_ref_correction(self, reference: str, 

431 alpha: float) -> tuple[float, float]: 

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

433 based on the reference electrode of choice and the temperature 

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

435 """ 

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

437 gibbs_corr = 0.0 

438 pH_corr = 0.0 

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

440 pH_corr += n_e * alpha 

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

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

443 if reference == 'AgCl': 

444 gibbs_corr -= n_e * U_STD_AGCL 

445 if reference == 'SCE': 

446 gibbs_corr -= n_e * U_STD_SCE 

447 

448 return gibbs_corr, pH_corr 

449 

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

451 """Print the chemical reaction.""" 

452 

453 reactants = [] 

454 products = [] 

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

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

457 continue 

458 nice_coeff = make_coeff_nice(n, max_denom) 

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

460 if n > 0: 

461 products.append(substr) 

462 else: 

463 reactants.append(substr) 

464 

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

466 

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

468 """Evaluate the reaction free energy 

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

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

471 

472 

473class Pourbaix: 

474 """Pourbaix class for acqueous stability evaluations. 

475 

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

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

478 

479 Initialization 

480 -------------- 

481 

482 material_name: str 

483 The formula of the target material. It is preferrable 

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

485 

486 refs_dct: dict 

487 A dictionary containing the formulae of the target material 

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

489 and their formation energies as values. 

490 

491 T: float 

492 Temperature in Kelvin. Default: 298.15 K. 

493 

494 conc: float 

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

496 

497 reference: str 

498 The reference electrode. Default: SHE. 

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

500 

501 

502 Relevant methods 

503 ---------------- 

504 

505 get_pourbaix_energy(U, pH) 

506 obtain the energy of the target material 

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

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

509 plot(...) 

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

511 

512 

513 Relevant attributes 

514 ------------------- 

515 

516 material: Species 

517 the target material as a Species object 

518 

519 phases: list[RedOx] 

520 the available decomposition reactions of the target material 

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

522 

523 """ 

524 def __init__(self, 

525 material_name: str, 

526 refs_dct: dict, 

527 T: float = 298.15, 

528 conc: float = 1.0e-6, 

529 reference: str = 'SHE'): 

530 

531 self.material_name = material_name 

532 self.refs = refs_dct 

533 self.T = T 

534 self.conc = conc 

535 self.reference = reference 

536 

537 refs = initialize_refs(refs_dct) 

538 self.material = refs.pop(material_name) 

539 self.phases, phase_matrix = get_phases( 

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

541 ) 

542 self._const = phase_matrix[:, 0] 

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

544 

545 def _decompose(self, U, pH): 

546 """Evaluate the reaction energy for decomposing 

547 the target material into each of the available products 

548 at a given pH and applied potential. 

549 """ 

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

551 

552 def _get_pourbaix_energy(self, U, pH): 

553 """Evaluate the Pourbaix energy""" 

554 energies = self._decompose(U, pH) 

555 i_min = np.argmin(energies) 

556 return -energies[i_min], i_min 

557 

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

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

560 the most stable phase and the corresponding energy at 

561 a given potential U and pH. 

562 

563 The Pourbaix energy represents the energy of the target material 

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

565 the target material can be considered as stable. 

566 """ 

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

568 phase = self.phases[index] 

569 if verbose: 

570 if energy <= 0.0: 

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

572 else: 

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

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

575 return energy, phase 

576 

577 def get_equations(self, contains: str | None = None): 

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

579 

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

581 particular species 

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

583 """ 

584 equations = [] 

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

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

587 continue 

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

589 return equations 

590 

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

592 """Actual evaluation of the complete diagram 

593 

594 Returns 

595 ------- 

596 

597 pour: 

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

599 domains are represented by indexes (as integers) 

600 that map to Pourbaix.phases 

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

602 

603 meta: 

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

605 

606 text: 

607 The coordinates and phases information for 

608 text placement on the diagram. 

609 

610 domains: 

611 ... 

612 

613 """ 

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

615 meta = pour.copy() 

616 

617 for i, u in enumerate(U): 

618 for j, p in enumerate(pH): 

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

620 

621 # Identifying the region where the target material 

622 # is stable and updating the diagram accordingly 

623 where_stable = (meta <= 0) 

624 pour[where_stable] = -1 

625 

626 text = [] 

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

628 for phase_id in domains: 

629 if phase_id == -1: 

630 where = where_stable 

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

632 else: 

633 where = (pour == phase_id) 

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

635 txt = phase.count 

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

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

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

639 

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

641 

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

643 """Plane intersection method for finding 

644 the boundaries between phases seen in the final plot. 

645 

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

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

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

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

650 """ 

651 from collections import defaultdict 

652 from itertools import combinations 

653 

654 # Planes identifying the diagram frame 

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

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

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

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

659 

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

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

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

663 # The region where the target material is stable 

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

665 for d in domains: 

666 if d == -1: 

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

668 const = 0.0 

669 else: 

670 pvec = self.phases[d]._vector 

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

672 const = -pvec[0] 

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

674 

675 # The simplices are found from the intersection points between 

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

677 # matches the corresponding pourbaix energy, 

678 # then the point is a simplex. 

679 simplices = [] 

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

681 combinations(planes, 3): 

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

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

684 ids = (id1, id2, id3) 

685 

686 try: 

687 # triplets containing parallel planes raise a LinAlgError 

688 # and are automatically excluded. 

689 invA = np.linalg.inv(A) 

690 except np.linalg.LinAlgError: 

691 continue 

692 

693 pt = np.dot(invA, c) 

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

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

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

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

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

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

700 simplices.append((simplex, ids)) 

701 

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

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

704 # in common, diagram frame excluded. 

705 duplicate_filter = defaultdict(int) 

706 

707 def are_boundaries_there(comm): 

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

709 for c in comm: 

710 if c in bounds: 

711 return True 

712 return False 

713 

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

715 

716 # common neighboring phases, diagram frame excluded 

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

718 

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

720 only_two_phases_in_common = (len(common) == 2) 

721 diagram_frame_excluded = not are_boundaries_there(common) 

722 

723 # Filtering out duplicates 

724 if all([simplices_are_distinct, only_two_phases_in_common, 

725 diagram_frame_excluded]): 

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

727 testarray.append(sorted(common)) 

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

729 duplicate_filter[tuple(testarray)] += 1 

730 

731 segments = [] 

732 for segment in duplicate_filter: 

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

734 segments.append(( 

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

736 list(phases.astype(int)) 

737 )) 

738 return segments 

739 

740 

741@dataclass 

742class PourbaixDiagram: 

743 pbx: Pourbaix 

744 U: np.ndarray 

745 pH: np.ndarray 

746 pour: np.ndarray 

747 meta: np.ndarray 

748 text: list[tuple[float, float, str]] 

749 domains: list[int] 

750 

751 def __post_init__(self): 

752 def _issorted(array): 

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

754 

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

756 if not _issorted(self.U): 

757 raise ValueError('U must be sorted') 

758 

759 if not _issorted(self.pH): 

760 raise ValueError('pH must be sorted') 

761 

762 @property 

763 def pHrange(self): 

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

765 

766 @property 

767 def Urange(self): 

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

769 

770 def _draw_diagram_axes( 

771 self, 

772 cap, 

773 normalize, 

774 include_text, 

775 include_water, 

776 labeltype, cmap, *, 

777 ax): 

778 """Backend for drawing Pourbaix diagrams.""" 

779 

780 meta = self.meta.copy() 

781 if normalize: 

782 meta /= self.pbx.material.natoms 

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

784 else: 

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

786 

787 fig = ax.get_figure() 

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

789 

790 fig.subplots_adjust( 

791 left=0.1, right=0.97, 

792 top=0.97, bottom=0.14 

793 ) 

794 

795 if isinstance(cap, list): 

796 vmin = cap[0] 

797 vmax = cap[1] 

798 else: 

799 vmin = -cap 

800 vmax = cap 

801 

802 colorplot = ax.imshow( 

803 meta, cmap=cmap, 

804 extent=extent, 

805 vmin=vmin, vmax=vmax, 

806 origin='lower', aspect='auto', 

807 interpolation='gaussian' 

808 ) 

809 

810 cbar = fig.colorbar( 

811 colorplot, 

812 ax=ax, 

813 pad=0.02 

814 ) 

815 

816 bounds = self.pbx.get_phase_boundaries( 

817 *self.pHrange, *self.Urange, self.domains 

818 ) 

819 for coords, _ in bounds: 

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

821 

822 if labeltype == 'numbers': 

823 add_numbers(ax, self.text) 

824 elif labeltype == 'phases': 

825 add_labels(ax, self.text) 

826 elif labeltype is None: 

827 pass 

828 else: 

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

830 

831 if include_water: 

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