Coverage for /builds/debichem-team/python-ase/ase/pourbaix.py: 94.24%

1import functools 

2import re 

3from collections import Counter 

4from dataclasses import dataclass 

5from fractions import Fraction 

6from itertools import chain, combinations, product 

7from typing import List, Tuple, Union 

8 

9import numpy as np 

10from scipy.linalg import null_space 

11 

12from ase.formula import Formula 

13from ase.units import kB 

14 

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

16 

17PREDEF_ENERGIES = { # Default chemical potentials 

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

19 'e-': 0.0, 

20 'H2O': -2.4583 

21} 

22 

23U_STD_AGCL = 0.222 # Standard redox potential of AgCl electrode 

24U_STD_SCE = 0.244 # Standard redox potential of SCE electrode 

25 

26 

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

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

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

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

31 formula_obj = Formula(formula_strip, format=fmt) 

32 return formula_obj, charge, aq 

33 

34 

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

36 """Convert dictionary entries to Species instances""" 

37 refs = {} 

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

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

40 refs[label] = spec 

41 return refs 

42 

43 

44def get_product_combos(reactant, refs): 

45 """Obtain all possible combinations of products. 

46 

47 Obtain - from the available references and based 

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

49 different combinations of products to be inserted 

50 in (electro)chemical reactions. 

51 """ 

52 ref_elem = reactant._main_elements 

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

54 for ref in refs.values(): 

55 contained = ref.contains(ref_elem) 

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

57 allcombos[w].append(ref) 

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

59 

60 

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

62 """Obtain all the possible decomposition pathways 

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

64 """ 

65 phases = [] 

66 phase_matrix = [] 

67 

68 for products in get_product_combos(reactant, refs): 

69 phase = RedOx.from_species( 

70 reactant, products, T, conc, reference, tol 

71 ) 

72 if phase is not None: 

73 phases.append(phase) 

74 phase_matrix.append(phase._vector) 

75 

76 if len(phase_matrix) == 0: 

77 raise ValueError( 

78 "No valid decomposition pathways have been found" + 

79 " given this set of references." 

80 ) 

81 

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

83 

84 

85def get_main_products(count): 

86 """Obtain the reaction products excluded protons, 

87 water and electrons. 

88 """ 

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

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

91 

92 

93def format_label(count) -> str: 

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

95 formatted = [] 

96 for prod in get_main_products(count): 

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

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

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

100 count = label.count(symbol) 

101 if count > 1: 

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

103 if count == 1: 

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

105 formatted.append(label) 

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

107 

108 

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

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

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

112 if frac.numerator == frac.denominator: 

113 return '' 

114 return str(frac) 

115 

116 

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

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

119 import matplotlib.patheffects as pfx 

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

121 txt = ax.text( 

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

123 fontsize=20, 

124 horizontalalignment='center' 

125 ) 

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

127 

128 

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

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

131 import matplotlib.patheffects as pfx 

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

133 label = format_label(species) 

134 annotation = ax.annotate( 

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

136 fontsize=16, horizontalalignment='center' 

137 ) 

138 annotation.set_path_effects( 

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

140 ) 

141 annotation.draggable() 

142 ax.add_artist(annotation) 

143 

144 

145def wrap_text(text) -> str: 

146 import textwrap 

147 

148 textlines = [] 

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

150 label = format_label(species) 

151 textlines.append( 

152 textwrap.fill( 

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

154 width=40, 

155 subsequent_indent=' ' 

156 ) 

157 ) 

158 return '\n'.join(textlines) 

159 

160 

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

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

163 

164 fig.text( 

165 0.75 + offset, 0.5, 

166 wrap_text(text), 

167 fontsize=16, 

168 va='center', 

169 ha='left') 

170 

171 

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

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

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

175 corr = { 

176 'SHE': 0, 

177 'RHE': 0, 

178 'Pt': 0, 

179 'AgCl': -U_STD_AGCL, 

180 'SCE': -U_STD_SCE, 

181 } 

182 kwargs = { 

183 'c': color, 

184 'ls': '--', 

185 'zorder': 2 

186 } 

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

188 slope = -59.2e-3 

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

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

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

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

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

194 else: 

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

196 

197 

198@functools.total_ordering 

199class Species: 

200 """Class representing an individual chemical species, 

201 grouping relevant properties. 

202 

203 Initialization 

204 -------------- 

205 name: str 

206 A label representing the species 

207 

208 formula: Formula 

209 

210 charge: int 

211 the electric charge of the species, if ionic 

212 

213 aq: bool 

214 whether the species is solid (False) 

215 or acqueous (True) 

216 

217 energy: float 

218 the chemical potential of the species 

219 """ 

220 def __init__(self, 

221 name: str, 

222 formula: Formula, 

223 charge: int, 

224 aq: bool, 

225 energy: float): 

226 

227 self.name = name 

228 self.formula = formula 

229 self.energy = energy 

230 self.charge = charge 

231 self.aq = aq 

232 

233 self.count = formula.count() 

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

235 self._main_elements = [ 

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

237 ] 

238 

239 def __eq__(self, other): 

240 if not isinstance(other, Species): 

241 return NotImplemented 

242 return self.name == other.name 

243 

244 def __lt__(self, other): 

245 if not isinstance(other, Species): 

246 return NotImplemented 

247 return self.name < other.name 

248 

249 @classmethod 

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

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

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

253 

254 string: str 

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

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

257 unit formula, and the chemical potential normalized accordingly. 

258 Acqueous species are specified by expliciting 

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

260 Parentheses for grouping functional groups are acceptable. 

261 e.g. 

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

263 Be3(OH)3+++ ➜ wrong 

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

265 Be3O3H3+++(aq) ➜ correct 

266 

267 energy: float 

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

269 

270 reduce: bool 

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

272 Formulae and energies of acqueous species are never reduced. 

273 

274 fmt: str 

275 Formula formatting according to the available options in 

276 ase.formula.Formula 

277 """ 

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

279 if not aq: 

280 if reduce: 

281 formula, n_fu = formula.reduce() 

282 energy /= n_fu 

283 name = str(formula) 

284 else: 

285 name = string 

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

287 

288 def get_chemsys(self): 

289 """Get the possible combinations of elements based 

290 on the stoichiometry. Useful for database queries. 

291 """ 

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

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

294 chemsys = list( 

295 chain.from_iterable( 

296 [combinations(elements, i + 1) 

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

298 ) 

299 ) 

300 return chemsys 

301 

302 def balance_electrochemistry(self): 

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

304 in electrochemical reactions. 

305 """ 

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

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

308 n_e = n_H + self.charge 

309 return n_H2O, n_H, n_e 

310 

311 def _count_array(self, elements): 

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

313 

314 def contains(self, elements): 

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

316 

317 def get_fractional_composition(self, elements): 

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

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

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

321 return N_elem / N_all 

322 

323 def __repr__(self): 

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

325 

326 

327class RedOx: 

328 def __init__(self, species, coeffs, 

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

330 reference='SHE'): 

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

332 

333 Initialization 

334 -------------- 

335 

336 species: list[Species] 

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

338 

339 coeffs: list[float] 

340 The stoichiometric coefficients of the above species. 

341 Positive coefficients are associated with the products, 

342 negative coefficients with the reactants. 

343 

344 T: float 

345 The temperature in Kelvin. Default: 298.15 K. 

346 

347 conc: float 

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

349 

350 reference: str 

351 The reference electrode. Default: SHE. 

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

353 

354 

355 Relevant methods 

356 ---------------- 

357 

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

359 Initialize the class from the reactant as a Species object 

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

361 

362 equation() 

363 Print the chemical equation of the reaction. 

364 

365 get_free_energy(U, pH): 

366 Obtain the reaction free energy at a given 

367 applied potential (U) and pH. 

368 

369 """ 

370 self.species = species 

371 self.coeffs = coeffs 

372 self.T = T 

373 self.conc = conc 

374 self.reference = reference 

375 self.count = Counter() 

376 

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

378 const_term = 0 

379 pH_term = 0 

380 U_term = 0 

381 

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

383 self.count[spec.name] = coef 

384 amounts = spec.balance_electrochemistry() 

385 

386 const_term += coef * ( 

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

388 ) 

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

390 U_term += - coef * amounts[2] 

391 

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

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

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

395 

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

397 self._vector = [ 

398 float(const_term + const_corr), 

399 float(U_term), 

400 float(pH_term + pH_corr) 

401 ] 

402 

403 @classmethod 

404 def from_species(cls, reactant, products, 

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

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

407 """Initialize the class from a combination of 

408 reactant and products. The stoichiometric 

409 coefficients are automatically determined. 

410 

411 reactant: Species 

412 products: list[Species] 

413 """ 

414 reac_elem = reactant._main_elements 

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

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

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

418 coeffs = null_space(elem_matrix).flatten() 

419 

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

421 coeffs /= coeffs[0] 

422 coeffs[0] = -coeffs[0] 

423 species = (reactant, *products) 

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

425 return phase 

426 

427 return None 

428 

429 def get_ref_correction(self, reference: str, 

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

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

432 based on the reference electrode of choice and the temperature 

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

434 """ 

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

436 gibbs_corr = 0.0 

437 pH_corr = 0.0 

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

439 pH_corr += n_e * alpha 

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

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

442 if reference == 'AgCl': 

443 gibbs_corr -= n_e * U_STD_AGCL 

444 if reference == 'SCE': 

445 gibbs_corr -= n_e * U_STD_SCE 

446 

447 return gibbs_corr, pH_corr 

448 

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

450 """Print the chemical reaction.""" 

451 

452 reactants = [] 

453 products = [] 

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

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

456 continue 

457 nice_coeff = make_coeff_nice(n, max_denom) 

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

459 if n > 0: 

460 products.append(substr) 

461 else: 

462 reactants.append(substr) 

463 

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

465 

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

467 """Evaluate the reaction free energy 

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

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

470 

471 

472class Pourbaix: 

473 """Pourbaix class for acqueous stability evaluations. 

474 

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

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

477 

478 Initialization 

479 -------------- 

480 

481 material_name: str 

482 The formula of the target material. It is preferrable 

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

484 

485 refs_dct: dict 

486 A dictionary containing the formulae of the target material 

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

488 and their formation energies as values. 

489 

490 T: float 

491 Temperature in Kelvin. Default: 298.15 K. 

492 

493 conc: float 

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

495 

496 reference: str 

497 The reference electrode. Default: SHE. 

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

499 

500 

501 Relevant methods 

502 ---------------- 

503 

504 get_pourbaix_energy(U, pH) 

505 obtain the energy of the target material 

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

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

508 plot(...) 

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

510 

511 

512 Relevant attributes 

513 ------------------- 

514 

515 material: Species 

516 the target material as a Species object 

517 

518 phases: list[RedOx] 

519 the available decomposition reactions of the target material 

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

521 

522 """ 

523 def __init__(self, 

524 material_name: str, 

525 refs_dct: dict, 

526 T: float = 298.15, 

527 conc: float = 1.0e-6, 

528 reference: str = 'SHE'): 

529 

530 self.material_name = material_name 

531 self.refs = refs_dct 

532 self.T = T 

533 self.conc = conc 

534 self.reference = reference 

535 

536 refs = initialize_refs(refs_dct) 

537 self.material = refs.pop(material_name) 

538 self.phases, phase_matrix = get_phases( 

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

540 ) 

541 self._const = phase_matrix[:, 0] 

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

543 

544 def _decompose(self, U, pH): 

545 """Evaluate the reaction energy for decomposing 

546 the target material into each of the available products 

547 at a given pH and applied potential. 

548 """ 

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

550 

551 def _get_pourbaix_energy(self, U, pH): 

552 """Evaluate the Pourbaix energy""" 

553 energies = self._decompose(U, pH) 

554 i_min = np.argmin(energies) 

555 return -energies[i_min], i_min 

556 

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

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

559 the most stable phase and the corresponding energy at 

560 a given potential U and pH. 

561 

562 The Pourbaix energy represents the energy of the target material 

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

564 the target material can be considered as stable. 

565 """ 

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

567 phase = self.phases[index] 

568 if verbose: 

569 if energy <= 0.0: 

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

571 else: 

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

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

574 return energy, phase 

575 

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

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

578 

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

580 particular species 

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

582 """ 

583 equations = [] 

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

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

586 continue 

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

588 return equations 

589 

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

591 """Actual evaluation of the complete diagram 

592 

593 Returns 

594 ------- 

595 

596 pour: 

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

598 domains are represented by indexes (as integers) 

599 that map to Pourbaix.phases 

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

601 

602 meta: 

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

604 

605 text: 

606 The coordinates and phases information for 

607 text placement on the diagram. 

608 

609 domains: 

610 ... 

611 

612 """ 

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

614 meta = pour.copy() 

615 

616 for i, u in enumerate(U): 

617 for j, p in enumerate(pH): 

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

619 

620 # Identifying the region where the target material 

621 # is stable and updating the diagram accordingly 

622 where_stable = (meta <= 0) 

623 pour[where_stable] = -1 

624 

625 text = [] 

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

627 for phase_id in domains: 

628 if phase_id == -1: 

629 where = where_stable 

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

631 else: 

632 where = (pour == phase_id) 

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

634 txt = phase.count 

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

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

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

638 

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

640 

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

642 """Plane intersection method for finding 

643 the boundaries between phases seen in the final plot. 

644 

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

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

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

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

649 """ 

650 from collections import defaultdict 

651 from itertools import combinations 

652 

653 # Planes identifying the diagram frame 

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

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

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

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

658 

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

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

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

662 # The region where the target material is stable 

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

664 for d in domains: 

665 if d == -1: 

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

667 const = 0.0 

668 else: 

669 pvec = self.phases[d]._vector 

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

671 const = -pvec[0] 

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

673 

674 # The simplices are found from the intersection points between 

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

676 # matches the corresponding pourbaix energy, 

677 # then the point is a simplex. 

678 simplices = [] 

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

680 combinations(planes, 3): 

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

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

683 ids = (id1, id2, id3) 

684 

685 try: 

686 # triplets containing parallel planes raise a LinAlgError 

687 # and are automatically excluded. 

688 invA = np.linalg.inv(A) 

689 except np.linalg.LinAlgError: 

690 continue 

691 

692 pt = np.dot(invA, c) 

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

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

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

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

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

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

699 simplices.append((simplex, ids)) 

700 

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

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

703 # in common, diagram frame excluded. 

704 duplicate_filter = defaultdict(int) 

705 

706 def are_boundaries_there(comm): 

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

708 for c in comm: 

709 if c in bounds: 

710 return True 

711 return False 

712 

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

714 

715 # common neighboring phases, diagram frame excluded 

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

717 

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

719 only_two_phases_in_common = (len(common) == 2) 

720 diagram_frame_excluded = not are_boundaries_there(common) 

721 

722 # Filtering out duplicates 

723 if all([simplices_are_distinct, only_two_phases_in_common, 

724 diagram_frame_excluded]): 

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

726 testarray.append(sorted(common)) 

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

728 duplicate_filter[tuple(testarray)] += 1 

729 

730 segments = [] 

731 for segment in duplicate_filter: 

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

733 segments.append(( 

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

735 list(phases.astype(int)) 

736 )) 

737 return segments 

738 

739 

740@dataclass 

741class PourbaixDiagram: 

742 pbx: Pourbaix 

743 U: np.ndarray 

744 pH: np.ndarray 

745 pour: np.ndarray 

746 meta: np.ndarray 

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

748 domains: List[int] 

749 

750 def __post_init__(self): 

751 def _issorted(array): 

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

753 

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

755 if not _issorted(self.U): 

756 raise ValueError('U must be sorted') 

757 

758 if not _issorted(self.pH): 

759 raise ValueError('pH must be sorted') 

760 

761 @property 

762 def pHrange(self): 

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

764 

765 @property 

766 def Urange(self): 

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

768 

769 def _draw_diagram_axes( 

770 self, 

771 cap, 

772 normalize, 

773 include_text, 

774 include_water, 

775 labeltype, cmap, *, 

776 ax): 

777 """Backend for drawing Pourbaix diagrams.""" 

778 

779 meta = self.meta.copy() 

780 if normalize: 

781 meta /= self.pbx.material.natoms 

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

783 else: 

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

785 

786 fig = ax.get_figure() 

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

788 

789 fig.subplots_adjust( 

790 left=0.1, right=0.97, 

791 top=0.97, bottom=0.14 

792 ) 

793 

794 if isinstance(cap, list): 

795 vmin = cap[0] 

796 vmax = cap[1] 

797 else: 

798 vmin = -cap 

799 vmax = cap 

800 

801 colorplot = ax.imshow( 

802 meta, cmap=cmap, 

803 extent=extent, 

804 vmin=vmin, vmax=vmax, 

805 origin='lower', aspect='auto', 

806 interpolation='gaussian' 

807 ) 

808 

809 cbar = fig.colorbar( 

810 colorplot, 

811 ax=ax, 

812 pad=0.02 

813 ) 

814 

815 bounds = self.pbx.get_phase_boundaries( 

816 *self.pHrange, *self.Urange, self.domains 

817 ) 

818 for coords, _ in bounds: 

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

820 

821 if labeltype == 'numbers': 

822 add_numbers(ax, self.text) 

823 elif labeltype == 'phases': 

824 add_labels(ax, self.text) 

825 elif labeltype is None: 

826 pass 

827 else: 

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

829 

830 if include_water: 

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