Add energy-diff and dissipation columns to gather_results

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,13 +91,20 @@ 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"] = []
     datadict["calculation"] = []
     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")

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