Skills › Scientific Agent Skills › pycalphad › equilibrate.py
equilibrate.py
184 lines · skills/pycalphad/scripts/equilibrate.py at 5f2bc03
1#!/usr/bin/env python3
2"""Compute TDB-based equilibria at one bulk mole composition over temperatures."""
3from __future__ import annotations
5import argparse
6import csv
7import hashlib
8import json
9from pathlib import Path
11import numpy as np
12import pycalphad as pc
13from pycalphad import variables as v
14from pycalphad.core.utils import filter_phases, unpack_species
17def finite(value, name, positive=False):
18 if isinstance(value, bool) or not isinstance(value, (int, float)) or not np.isfinite(value):
19 raise ValueError(f"{name} must be a finite number")
20 if positive and value <= 0:
21 raise ValueError(f"{name} must be positive")
22 return float(value)
25def validate(db, settings):
26 keys = {"components", "phases", "composition_basis", "dependent_component", "independent_mole_fractions",
27 "temperatures_k", "pressure_pa", "database_temperature_range_k", "pdens",
28 "mass_balance_tolerance", "phase_fraction_tolerance", "gibbs_energy_tolerance_j_per_mol"}
29 if set(settings) != keys:
30 raise ValueError(f"Missing settings: {sorted(keys-set(settings))}; unknown: {sorted(set(settings)-keys)}")
31 for key in ("components", "phases"):
32 values = settings[key]
33 if not isinstance(values, list) or not values or any(not isinstance(s, str) for s in values) or len(set(values)) != len(values):
34 raise ValueError(f"{key} must be a nonempty list of unique names")
35 comps = settings["components"]
36 phases = settings["phases"]
37 if set(comps) - db.elements:
38 raise ValueError(f"Components absent from database elements: {sorted(set(comps) - db.elements)}")
39 if set(phases) - set(db.phases):
40 raise ValueError(f"Unknown phases: {sorted(set(phases) - set(db.phases))}")
41 active = filter_phases(db, unpack_species(db, comps), phases)
42 if set(active) != set(phases):
43 raise ValueError(f"Incompatible or order-disorder-filtered phases: {sorted(set(phases) - set(active))}")
44 if settings["composition_basis"] != "mole_fraction":
45 raise ValueError("Only elemental mole_fraction conditions are accepted; convert mass fractions explicitly")
46 material = set(comps) - {"VA"}
47 dependent = settings["dependent_component"]
48 if dependent not in material:
49 raise ValueError("dependent_component must be a non-vacancy selected element")
50 fractions = settings["independent_mole_fractions"]
51 if not isinstance(fractions, dict) or set(fractions) != material - {dependent}:
52 raise ValueError("Specify exactly one mole fraction for each nondependent, non-vacancy component")
53 fractions = {name: finite(value, name) for name, value in fractions.items()}
54 if any(value < 0 or value > 1 for value in fractions.values()) or sum(fractions.values()) > 1:
55 raise ValueError("Mole fractions must be in [0,1] and sum to at most one")
56 bulk = {**fractions, dependent: 1 - sum(fractions.values())}
57 temperatures = settings["temperatures_k"]
58 if not isinstance(temperatures, list) or not temperatures or len(temperatures) > 1000:
59 raise ValueError("temperatures_k must contain between 1 and 1000 explicit temperatures")
60 for temperature in temperatures:
61 finite(temperature, "temperature_k", True)
62 if len(set(temperatures)) != len(temperatures):
63 raise ValueError("Duplicate temperatures are not accepted")
64 bounds = settings["database_temperature_range_k"]
65 if bounds is not None:
66 if not isinstance(bounds, list) or len(bounds) != 2:
67 raise ValueError("database_temperature_range_k must be [low, high] or null if unknown")
68 low, high = [finite(x, "database temperature bound", True) for x in bounds]
69 if low >= high or min(temperatures) < low or max(temperatures) > high:
70 raise ValueError("Requested temperatures lie outside the declared database range")
71 for name in ("pressure_pa", "mass_balance_tolerance", "phase_fraction_tolerance", "gibbs_energy_tolerance_j_per_mol"):
72 finite(settings[name], name, True)
73 pdens = settings["pdens"]
74 if isinstance(pdens, bool) or not isinstance(pdens, int) or not 10 <= pdens <= 10000:
75 raise ValueError("pdens must be an integer from 10 through 10000")
76 return bulk
79def solve_point(db, settings, bulk, temperature, pdens):
80 conditions = {v.T: temperature, v.P: settings["pressure_pa"], v.N: 1}
81 conditions.update({v.X(name): value for name, value in settings["independent_mole_fractions"].items()})
82 result = pc.equilibrium(db, settings["components"], settings["phases"], conditions, calc_opts={"pdens": pdens})
83 # pycalphad clips independent composition conditions near zero and one.
84 # Preserve the actual imposed coordinates without replacing the requested bulk.
85 solver_bulk = {name: float(result.coords[str(v.X(name))].values.item())
86 for name in settings["independent_mole_fractions"]}
87 solver_bulk[settings["dependent_component"]] = 1 - sum(solver_bulk.values())
88 condition_adjustment = max(abs(solver_bulk[name] - bulk[name]) for name in bulk)
89 energy = float(result.GM.values.item())
90 names = result.Phase.values.ravel()
91 amounts = result.NP.values.ravel()
92 components = [str(c) for c in result.component.values]
93 compositions = result.X.values.reshape(-1, len(components))
94 if not np.isfinite(energy) or set(components) != set(bulk):
95 raise ValueError("Equilibrium failed: nonfinite energy or unexpected output components")
96 vertices = []
97 reconstructed = {name: 0.0 for name in bulk}
98 phase_sums = {name: 0.0 for name in settings["phases"]}
99 for index, (name, amount, composition) in enumerate(zip(names, amounts, compositions)):
100 if name == "":
101 continue # pycalphad pads unused vertices with blank names and NaN amounts.
102 if name not in phase_sums or not np.isfinite(amount) or not np.isfinite(composition).all():
103 raise ValueError("Equilibrium contains an invalid or fictitious stable vertex")
104 if amount < -1e-10 or amount > 1 + 1e-10 or np.min(composition) < -1e-10 or abs(composition.sum() - 1) > 1e-7:
105 raise ValueError("Equilibrium has invalid phase fraction or composition")
106 fractions = dict(zip(components, map(float, composition)))
107 vertices.append({"vertex": index, "phase": str(name), "mole_phase_fraction": float(amount), "mole_fractions": fractions})
108 phase_sums[str(name)] += float(amount)
109 for component, fraction in fractions.items():
110 reconstructed[component] += float(amount) * fraction
111 if not vertices:
112 raise ValueError("Equilibrium produced no stable phases")
113 amount_error = abs(sum(phase_sums.values()) - 1)
114 mass_error = max(abs(reconstructed[c] - bulk[c]) for c in bulk)
115 return {"temperature_k": temperature, "gibbs_energy_j_per_mol": energy, "vertices": vertices,
116 "phase_totals": phase_sums, "reconstructed_bulk_mole_fractions": reconstructed,
117 "solver_bulk_mole_fractions": solver_bulk,
118 "composition_condition_adjustment_absolute_error": condition_adjustment,
119 "phase_fraction_sum_error": amount_error, "mass_balance_absolute_error": mass_error,
120 "balance_passed": amount_error <= settings["mass_balance_tolerance"] and mass_error <= settings["mass_balance_tolerance"]}
123def calculate(db, settings):
124 bulk = validate(db, settings)
125 points = []
126 for temperature in settings["temperatures_k"]:
127 baseline = solve_point(db, settings, bulk, temperature, settings["pdens"])
128 refined = solve_point(db, settings, bulk, temperature, 2 * settings["pdens"])
129 fraction_change = max(abs(baseline["phase_totals"][p] - refined["phase_totals"][p]) for p in settings["phases"])
130 energy_change = abs(baseline["gibbs_energy_j_per_mol"] - refined["gibbs_energy_j_per_mol"])
131 points.append({"baseline": baseline, "refined": refined,
132 "maximum_phase_fraction_change": fraction_change,
133 "gibbs_energy_change_j_per_mol": energy_change,
134 "checks_passed": baseline["balance_passed"] and refined["balance_passed"] and
135 fraction_change <= settings["phase_fraction_tolerance"] and energy_change <= settings["gibbs_energy_tolerance_j_per_mol"]})
136 return {"bulk_mole_fractions": bulk, "points": points,
137 "all_checks_passed": all(point["checks_passed"] for point in points),
138 "excluded_database_phases": sorted(set(db.phases) - set(settings["phases"]))}
141def run(database_path, settings_path, output_dir):
142 database_path, settings_path, output_dir = map(Path, (database_path, settings_path, output_dir))
143 if output_dir.exists():
144 raise ValueError("Output directory already exists; choose a new directory")
145 settings = json.loads(settings_path.read_text())
146 db = pc.Database(str(database_path))
147 report = calculate(db, settings)
148 report.update(versions={"pycalphad": pc.__version__, "numpy": np.__version__}, settings=settings,
149 database_sha256=hashlib.sha256(database_path.read_bytes()).hexdigest(),
150 settings_sha256=hashlib.sha256(settings_path.read_bytes()).hexdigest(),
151 phase_fraction_basis="molar, N=1; vacancies excluded from bulk composition",
152 database_experimental_validity="not_evaluated_by_helper",
153 refinement="pdens doubled; phase totals and GM compared; not a proof of the global minimum")
154 output_dir.mkdir(parents=True)
155 (output_dir / "report.json").write_text(json.dumps(report, indent=2, allow_nan=False) + "\n")
156 with (output_dir / "phase-equilibria.csv").open("w", newline="") as handle:
157 writer = csv.writer(handle)
158 components = sorted(report["bulk_mole_fractions"])
159 writer.writerow(["run", "temperature_k", "vertex", "phase", "mole_phase_fraction", "gibbs_energy_j_per_mol", *[f"X_{c}" for c in components]])
160 for point in report["points"]:
161 for variant in ("baseline", "refined"):
162 result = point[variant]
163 for vertex in result["vertices"]:
164 writer.writerow([variant, result["temperature_k"], vertex["vertex"], vertex["phase"],
165 vertex["mole_phase_fraction"], result["gibbs_energy_j_per_mol"],
166 *[vertex["mole_fractions"][c] for c in components]])
167 return report
170def main():
171 parser = argparse.ArgumentParser(description=__doc__)
172 parser.add_argument("database", type=Path, help="Thermodynamic TDB database")
173 parser.add_argument("settings", type=Path, help="Conditions JSON")
174 parser.add_argument("output", type=Path, help="New output directory")
175 args = parser.parse_args()
176 try:
177 report = run(args.database, args.settings, args.output)
178 except (ValueError, OSError, TypeError, KeyError) as error:
179 parser.error(str(error))
180 print(json.dumps({"output": str(args.output), "all_checks_passed": report["all_checks_passed"]}))
183if __name__ == "__main__":
184 main()