From 466d85c88d6b1c593f8f5a0fe621d5f5b85772ad Mon Sep 17 00:00:00 2001 From: Prabhath C Date: Tue, 28 Jul 2026 10:00:59 +0200 Subject: [PATCH] Add energy-diff and dissipation columns to gather_results Surfaces free_energy_error, dissipation/ts_dissipation, and the raw forward/backward switching energy-diff (+ lambda) arrays so phase transitions during a ts/tscale sweep can be diagnosed straight from the gathered DataFrame instead of re-parsing each folder by hand. Co-Authored-By: Claude Sonnet 5 --- calphy/postprocessing.py | 72 ++++++++++++++- tests/test_postprocessing.py | 173 +++++++++++++++++++++++++++++++++++ 2 files changed, 243 insertions(+), 2 deletions(-) create mode 100644 tests/test_postprocessing.py diff --git a/calphy/postprocessing.py b/calphy/postprocessing.py index a9f5802..1269ab3 100644 --- a/calphy/postprocessing.py +++ b/calphy/postprocessing.py @@ -64,7 +64,20 @@ def gather_results(mainfolder, reduce_composition=True, extract_phase_prefix=Fal Returns ------- df: pandas DataFrame - DataFrame with results + DataFrame with results. In addition to the columns produced + previously, this also includes: + + - free_energy_error: array (ts/tscale, from temperature_sweep.dat) + or 0.0 (fe/alchemy/composition_scaling); the statistical + standard error of the mean free energy, not a hysteresis check. + - dissipation: mean switching dissipation (fe/alchemy) or NaN (ts/tscale) + - ts_dissipation: max hysteresis over a ts/tscale sweep, or NaN otherwise + - forward_energy_diff / backward_energy_diff: list of arrays, one per + reversible-scaling replica (ts/tscale only, else None); the raw + per-lambda energy differential, useful to see *where* along the + sweep forward/backward diverge (phase-transition diagnostic). + - forward_lambda / backward_lambda: matching lambda arrays for the + above (ts/tscale only, else None) """ try: import pandas as pd @@ -78,6 +91,9 @@ def gather_results(mainfolder, reduce_composition=True, extract_phase_prefix=Fal datadict["temperature"] = [] datadict["pressure"] = [] datadict["free_energy"] = [] + datadict["free_energy_error"] = [] + datadict["dissipation"] = [] + datadict["ts_dissipation"] = [] datadict["reference_phase"] = [] datadict["error_code"] = [] datadict["composition"] = [] @@ -85,6 +101,10 @@ def gather_results(mainfolder, reduce_composition=True, extract_phase_prefix=Fal datadict["ideal_entropy"] = [] datadict["phase_name"] = [] datadict["reference_composition"] = [] + datadict["forward_energy_diff"] = [] + datadict["backward_energy_diff"] = [] + datadict["forward_lambda"] = [] + datadict["backward_lambda"] = [] folders = next(os.walk(mainfolder))[1] for folder in folders: @@ -119,6 +139,13 @@ def gather_results(mainfolder, reduce_composition=True, extract_phase_prefix=Fal # check output file outfile = os.path.join(mainfolder, folder, "report.yaml") datadict["error_code"].append(None) + datadict["free_energy_error"].append(np.nan) + datadict["dissipation"].append(np.nan) + datadict["ts_dissipation"].append(np.nan) + datadict["forward_energy_diff"].append(None) + datadict["backward_energy_diff"].append(None) + datadict["forward_lambda"].append(None) + datadict["backward_lambda"].append(None) # print(inpfile) if not os.path.exists(outfile): @@ -131,12 +158,17 @@ def gather_results(mainfolder, reduce_composition=True, extract_phase_prefix=Fal if mode in ["fe", "alchemy", "composition_scaling"]: datadict["status"].append("True") + datadict["free_energy_error"][-1] = 0.0 # ok, valid calculation, try to parse input file to get info with open(outfile, "r") as fin: out = yaml.safe_load(fin) datadict["free_energy"].append(out["results"]["free_energy"]) + # scalar quality metrics calphy already computes: qdiss (mean switching + # dissipation, fe/alchemy) and max hysteresis over a ts/tscale sweep + datadict["dissipation"][-1] = out["results"].get("dissipation", np.nan) + datadict["ts_dissipation"][-1] = out["results"].get("ts_dissipation", np.nan) # add normal composition el_arr = np.array(out["input"]["element"].split(" ")).astype(str) @@ -170,9 +202,45 @@ def gather_results(mainfolder, reduce_composition=True, extract_phase_prefix=Fal ) if os.path.exists(datafile): datadict["status"].append("True") - t, f = np.loadtxt(datafile, unpack=True, usecols=(0, 1)) + t, f, ferr = np.loadtxt(datafile, unpack=True, usecols=(0, 1, 2)) datadict["temperature"][-1] = t datadict["free_energy"][-1] = f + datadict["free_energy_error"][-1] = ferr + + # raw per-lambda energy differential of each forward/backward + # switching replica; diverging forward vs. backward curves + # signal hysteresis/a structural change during the sweep + # (ts_dissipation above is just the max of this difference) + f_ediffs = [] + b_ediffs = [] + f_lambdas = [] + b_lambdas = [] + i = 1 + while True: + fwdfile = os.path.join( + mainfolder, folder, f"ts.forward_{i}.dat" + ) + bkdfile = os.path.join( + mainfolder, folder, f"ts.backward_{i}.dat" + ) + if not (os.path.exists(fwdfile) and os.path.exists(bkdfile)): + break + fdx, _fp, _fvol, flambda = np.loadtxt( + fwdfile, unpack=True, comments="#" + ) + bdx, _bp, _bvol, blambda = np.loadtxt( + bkdfile, unpack=True, comments="#" + ) + f_ediffs.append(fdx / flambda) + b_ediffs.append(bdx / blambda) + f_lambdas.append(flambda) + b_lambdas.append(blambda) + i += 1 + + datadict["forward_energy_diff"][-1] = f_ediffs if f_ediffs else None + datadict["backward_energy_diff"][-1] = b_ediffs if b_ediffs else None + datadict["forward_lambda"][-1] = f_lambdas if f_lambdas else None + datadict["backward_lambda"][-1] = b_lambdas if b_lambdas else None else: datadict["status"].append("False") errfile = os.path.join(os.getcwd(), mainfolder, folder + ".sub.err") diff --git a/tests/test_postprocessing.py b/tests/test_postprocessing.py new file mode 100644 index 0000000..3676e97 --- /dev/null +++ b/tests/test_postprocessing.py @@ -0,0 +1,173 @@ +"""Tests for calphy.postprocessing.gather_results. + +Builds minimal calculation folders by hand (input_file.yaml + report.yaml, +plus temperature_sweep.dat / ts.forward_i.dat / ts.backward_i.dat for the +ts mode) rather than running real calphy jobs, mirroring the on-disk shape +Calculation.dump() produces -- see calphy/input.py for the phase_name +(default "") and reference_composition (default 0.0) field defaults. +""" +import os + +import numpy as np +import yaml + +from calphy.postprocessing import gather_results + +INPUT_TS = { + "calculations": [ + { + "mode": "ts", + "temperature": [500, 600], + "pressure": 0, + "reference_phase": "solid", + "phase_name": "", + "reference_composition": 0.0, + } + ] +} + +INPUT_FE = { + "calculations": [ + { + "mode": "fe", + "temperature": 500, + "pressure": 0, + "reference_phase": "solid", + "phase_name": "", + "reference_composition": 0.0, + } + ] +} + + +def _write_yaml(path, data): + with open(path, "w") as fh: + yaml.dump(data, fh) + + +def _make_ts_calc(folder, with_replicas=True): + # "replicas" = independent forward/backward switching runs (nsims in + # calphy/integrators.py), each written as a ts.forward_i.dat/ts.backward_i.dat + # pair; with_replicas=False omits them to test the no-file fallback. + os.makedirs(folder, exist_ok=True) + _write_yaml(os.path.join(folder, "input_file.yaml"), INPUT_TS) + _write_yaml( + os.path.join(folder, "report.yaml"), + { + "input": {"element": "Al Cu", "concentration": "0.5 0.5"}, + "results": {"free_energy": -3.5, "ts_dissipation": 0.002}, + }, + ) + t = np.array([500.0, 600.0]) + f = np.array([-3.5, -3.6]) + ferr = np.array([0.001, 0.0012]) + np.savetxt(os.path.join(folder, "temperature_sweep.dat"), np.column_stack((t, f, ferr))) + + if with_replicas: + lam = np.linspace(1, 0.5, 5) + dx = np.arange(5, dtype=float) + 1.0 + p = np.zeros(5) + vol = np.ones(5) + np.savetxt( + os.path.join(folder, "ts.forward_1.dat"), + np.column_stack((dx, p, vol, lam)), + header="x", + ) + np.savetxt( + os.path.join(folder, "ts.backward_1.dat"), + np.column_stack((dx, p, vol, lam)), + header="x", + ) + return t, f, ferr, lam if with_replicas else None + + +def _make_fe_calc(folder, with_report=True): + os.makedirs(folder, exist_ok=True) + _write_yaml(os.path.join(folder, "input_file.yaml"), INPUT_FE) + if with_report: + _write_yaml( + os.path.join(folder, "report.yaml"), + { + "input": {"element": "Al Cu", "concentration": "0.5 0.5"}, + "results": {"free_energy": -3.4, "dissipation": 0.0005}, + }, + ) + + +def test_gather_results_ts_mode_with_replicas(tmp_path): + """ts mode picks up free_energy_error, ts_dissipation, and per-replica + forward/backward energy-diff + lambda arrays.""" + mainfolder = tmp_path / "calcs" + t, f, ferr, lam = _make_ts_calc(str(mainfolder / "run1")) + + df = gather_results(str(mainfolder)) + row = df.loc[df.calculation == "run1"].iloc[0] + + assert row.status == "True" + assert row.calculation_mode == "ts" + np.testing.assert_allclose(row.free_energy_error, ferr) + assert np.isnan(row.dissipation) + assert row.ts_dissipation == 0.002 + + assert row.forward_energy_diff is not None + assert row.backward_energy_diff is not None + assert len(row.forward_energy_diff) == 1 + dx = np.arange(5, dtype=float) + 1.0 + np.testing.assert_allclose(row.forward_energy_diff[0], dx / lam) + np.testing.assert_allclose(row.forward_lambda[0], lam) + np.testing.assert_allclose(row.backward_lambda[0], lam) + + +def test_gather_results_ts_mode_without_replicas(tmp_path): + """Without ts.forward_i.dat/ts.backward_i.dat files, the diff/lambda + columns fall back to None but free_energy_error still parses.""" + mainfolder = tmp_path / "calcs" + t, f, ferr, _ = _make_ts_calc(str(mainfolder / "run1"), with_replicas=False) + + df = gather_results(str(mainfolder)) + row = df.loc[df.calculation == "run1"].iloc[0] + + assert row.status == "True" + np.testing.assert_allclose(row.free_energy_error, ferr) + assert row.forward_energy_diff is None + assert row.backward_energy_diff is None + assert row.forward_lambda is None + assert row.backward_lambda is None + + +def test_gather_results_fe_mode(tmp_path): + """fe mode has no temperature sweep: free_energy_error is fixed at 0.0 + and dissipation is read straight from report.yaml; no diff arrays.""" + mainfolder = tmp_path / "calcs" + _make_fe_calc(str(mainfolder / "run2")) + + df = gather_results(str(mainfolder)) + row = df.loc[df.calculation == "run2"].iloc[0] + + assert row.status == "True" + assert row.calculation_mode == "fe" + assert row.free_energy == -3.4 + assert row.free_energy_error == 0.0 + assert row.dissipation == 0.0005 + assert np.isnan(row.ts_dissipation) + assert row.forward_energy_diff is None + assert row.backward_energy_diff is None + + +def test_gather_results_missing_report(tmp_path): + """A calculation folder with no report.yaml (e.g. still running or + failed) reports status False and leaves every new column at its + NaN/None default -- no regression vs. the pre-existing columns.""" + mainfolder = tmp_path / "calcs" + _make_fe_calc(str(mainfolder / "run3"), with_report=False) + + df = gather_results(str(mainfolder)) + row = df.loc[df.calculation == "run3"].iloc[0] + + assert row.status == "False" + assert np.isnan(row.free_energy) + assert np.isnan(row.free_energy_error) + assert np.isnan(row.dissipation) + assert np.isnan(row.ts_dissipation) + assert row.forward_energy_diff is None + assert row.backward_energy_diff is None