[AI-Assisted] Add uncertainty-aware thermodynamic benchmark framework

.github/workflows/kij-calibration.yml

@@ -0,0 +1,25 @@
+name: Temporary kij calibration
+
+on:
+  pull_request:
+    paths:
+      - ".github/workflows/kij-calibration.yml"
+      - "src/main/java/neqsim/thermo/util/benchmark/**"
+      - "src/test/java/neqsim/thermo/util/benchmark/**"
+
+permissions:
+  contents: read
+
+jobs:
+  calibrate:
+    runs-on: ubuntu-latest
+    timeout-minutes: 20
+    steps:
+      - uses: actions/checkout@v7
+      - uses: actions/setup-java@v5
+        with:
+          distribution: temurin
+          java-version: "21"
+          cache: maven
+      - name: Run focused calibration
+        run: ./mvnw -Dtest=ThermodynamicBenchmarkTest#calibratesCubicModelsAgainstPublishedData test

docs/REFERENCE_MANUAL_INDEX.md

@@ -159,6 +159,7 @@ NeqSim is distributed under the Apache-2.0 license and can be used via:
 | Mercury Thermodynamics    | [docs/thermo/mercury_thermodynamics.md](thermo/mercury_thermodynamics) | Mercury-focused SRK-TwuCoon-Statoil-EOS usage, TPflash setup, and thesis-linked BIP/correlation guidance |
 | **Søreide-Whitson Model** | [docs/thermo/SoreideWhitsonModel.md](thermo/SoreideWhitsonModel)   | **Gas solubility in brine** - Modified PR EoS with salinity effects, used in NeqSimLive for produced water emissions |
 | GERG-2008                 | [docs/thermo/gerg2008_eoscg.md](thermo/gerg2008_eoscg)             | GERG-2008, GERG-2008-H2, GERG-2008-NH3, and EOS-CG equations of state                                               |
+| Thermodynamic Benchmarks  | [docs/thermo/thermodynamic_benchmarks.md](thermo/thermodynamic_benchmarks) | Experimental-data benchmarks with provenance, uncertainty-aware metrics, and H2-CO2 phase-equilibrium data |
 | GERG-2008-NH3 Notebook    | [docs/examples/GERG2008_NH3_Ammonia_Properties.ipynb](examples/GERG2008_NH3_Ammonia_Properties.ipynb) | Ammonia properties with Gao EOS — density validation, mixture properties, isotherms |
 | Mixing Rules              | [docs/thermo/mixing_rules_guide.md](thermo/mixing_rules_guide)     | Mixing rules and BIPs                                                                                                |
 | Mixing Rule Package       | [docs/thermo/mixingrule/README.md](thermo/mixingrule/)             | Mixing rule implementations                                                                                          |

docs/thermo/thermodynamic_benchmarks.md

@@ -0,0 +1,72 @@
+---
+title: "Thermodynamic Model Benchmarks"
+description: "Auditable comparison of NeqSim property models with published experimental data."
+---
+
+# Thermodynamic model benchmarks
+
+NeqSim provides a reusable benchmark framework in
+`neqsim.thermo.util.benchmark.ThermodynamicBenchmark`. It keeps state,
+composition, property units, experimental uncertainty, citation, DOI, and data
+reuse information together with every comparison.
+
+The framework reports:
+
+- average absolute relative deviation (AARD);
+- signed bias;
+- root-mean-square relative error;
+- maximum absolute relative error;
+- uncertainty-normalized residuals when experimental uncertainty is available;
+- all point-level experimental and predicted values.
+
+## Hydrogen-containing CO2 benchmark
+
+`H2CO2PhaseEquilibriumData.load()` reads 24 bubble- and dew-point values from
+Tables IV and V of:
+
+> Zhang et al. (2026), *Measurement of phase equilibrium characteristics and
+> equation-of-state applicability for hydrogen-containing CO2 systems*,
+> International Journal of Fluid Engineering 3, 013903.
+> [doi:10.1063/5.0288386](https://doi.org/10.1063/5.0288386)
+
+The systems cover:
+
+- 96 mol% CO2 + 4 mol% H2;
+- 96 mol% CO2 + 2 mol% H2 + 2 mol% N2;
+- temperatures from -30 to 20 °C;
+- bubble pressures from 13.8 to 64.3 bara;
+- dew pressures from 65.3 to 112.1 bara.
+
+The article does not report pointwise pressure uncertainties. NeqSim records
+these as unavailable instead of assigning an unsupported uncertainty.
+
+## Java example
+
+```java
+ThermodynamicBenchmark.Dataset dataset = H2CO2PhaseEquilibriumData.load();
+NeqSimPhaseEquilibriumPrediction prediction =
+    new NeqSimPhaseEquilibriumPrediction(
+        NeqSimPhaseEquilibriumPrediction.Model.GERG_2008_H2);
+
+ThermodynamicBenchmark.Report report =
+    ThermodynamicBenchmark.run("GERG-2008-H2", dataset, prediction);
+
+double aardPercent = report.getAverageAbsoluteRelativeDeviationPercent();
+double biasPercent = report.getBiasPercent();
+```
+
+Supported configurations are SRK, PR, standard GERG-2008, and
+GERG-2008-H2. Cubic models use NeqSim database interaction parameters and
+mixing rule 2. The GERG-2008-H2 configuration enables the hydrogen-enhanced
+binary parameters and departure functions.
+
+The experimental pressure is used only as a numerical starting point. It is
+not returned by the prediction adapter and does not alter model parameters.
+
+## Interpretation
+
+A low aggregate deviation is not sufficient evidence of universal model
+validity. Inspect bubble and dew results separately, examine signed residuals
+against temperature and composition, and retain the source validity range.
+Do not tune binary interaction parameters without recording their provenance
+and validating them on data not used for regression.

src/main/java/neqsim/thermo/util/benchmark/BinaryInteractionParameterFitter.java

@@ -0,0 +1,163 @@
+package neqsim.thermo.util.benchmark;
+
+import java.util.ArrayList;
+import java.util.List;
+import neqsim.thermo.util.benchmark.ThermodynamicBenchmark.Dataset;
+import neqsim.thermo.util.benchmark.ThermodynamicBenchmark.Point;
+
+/**
+ * Bounded one-dimensional regression of a constant cubic-EOS binary interaction parameter.
+ *
+ * <p>
+ * The objective is the mean squared relative error. This gives bubble- and dew-point pressures comparable influence
+ * even when their absolute pressure levels differ.
+ * </p>
+ */
+public final class BinaryInteractionParameterFitter {
+  /** Creates model predictions for a trial binary interaction parameter. */
+  public interface PredictionFactory {
+    /**
+     * Creates a prediction function for one trial parameter.
+     *
+     * @param binaryInteractionParameter trial dimensionless binary interaction parameter
+     * @return model prediction function
+     */
+    ThermodynamicBenchmark.Prediction create(double binaryInteractionParameter);
+  }
+
+  /** Immutable scalar-regression result. */
+  public static final class Result {
+    private final double binaryInteractionParameter;
+    private final double rootMeanSquareRelativeErrorPercent;
+    private final int objectiveEvaluations;
+
+    private Result(double binaryInteractionParameter, double rootMeanSquareRelativeErrorPercent,
+        int objectiveEvaluations) {
+      this.binaryInteractionParameter = binaryInteractionParameter;
+      this.rootMeanSquareRelativeErrorPercent = rootMeanSquareRelativeErrorPercent;
+      this.objectiveEvaluations = objectiveEvaluations;
+    }
+
+    /** @return fitted dimensionless binary interaction parameter */
+    public double getBinaryInteractionParameter() {
+      return binaryInteractionParameter;
+    }
+
+    /** @return in-sample root-mean-square relative error in percent */
+    public double getRootMeanSquareRelativeErrorPercent() {
+      return rootMeanSquareRelativeErrorPercent;
+    }
+
+    /** @return number of objective evaluations */
+    public int getObjectiveEvaluations() {
+      return objectiveEvaluations;
+    }
+  }
+
+  private final Dataset dataset;
+  private final PredictionFactory predictionFactory;
+
+  /**
+   * Creates a scalar binary-interaction-parameter regression.
+   *
+   * @param dataset experimental calibration dataset
+   * @param predictionFactory model adapter factory
+   */
+  public BinaryInteractionParameterFitter(Dataset dataset, PredictionFactory predictionFactory) {
+    if (dataset == null || predictionFactory == null) {
+      throw new IllegalArgumentException("Dataset and prediction factory are required");
+    }
+    this.dataset = dataset;
+    this.predictionFactory = predictionFactory;
+  }
+
+  /**
+   * Fits within a closed interval using golden-section minimization.
+   *
+   * @param lowerBound inclusive lower parameter bound
+   * @param upperBound inclusive upper parameter bound
+   * @param tolerance parameter-space convergence tolerance
+   * @param maximumEvaluations maximum objective evaluations
+   * @return fitted parameter and in-sample objective
+   * @throws Exception when a model prediction fails
+   */
+  public Result fit(double lowerBound, double upperBound, double tolerance, int maximumEvaluations) throws Exception {
+    if (!Double.isFinite(lowerBound) || !Double.isFinite(upperBound) || lowerBound >= upperBound
+        || !Double.isFinite(tolerance) || tolerance <= 0.0 || maximumEvaluations < 3) {
+      throw new IllegalArgumentException("Valid bounds, tolerance, and evaluation limit are required");
+    }
+
+    double left = lowerBound;
+    double right = upperBound;
+    double inverseGoldenRatio = (Math.sqrt(5.0) - 1.0) / 2.0;
+    double innerLeft = right - inverseGoldenRatio * (right - left);
+    double innerRight = left + inverseGoldenRatio * (right - left);
+    double objectiveLeft = objective(innerLeft);
+    double objectiveRight = objective(innerRight);
+    int evaluations = 2;
+
+    while (right - left > tolerance && evaluations < maximumEvaluations) {
+      if (objectiveLeft <= objectiveRight) {
+        right = innerRight;
+        innerRight = innerLeft;
+        objectiveRight = objectiveLeft;
+        innerLeft = right - inverseGoldenRatio * (right - left);
+        objectiveLeft = objective(innerLeft);
+      } else {
+        left = innerLeft;
+        innerLeft = innerRight;
+        objectiveLeft = objectiveRight;
+        innerRight = left + inverseGoldenRatio * (right - left);
+        objectiveRight = objective(innerRight);
+      }
+      evaluations++;
+    }
+
+    double fittedParameter;
+    double fittedObjective;
+    if (objectiveLeft <= objectiveRight) {
+      fittedParameter = innerLeft;
+      fittedObjective = objectiveLeft;
+    } else {
+      fittedParameter = innerRight;
+      fittedObjective = objectiveRight;
+    }
+    return new Result(fittedParameter, 100.0 * Math.sqrt(fittedObjective), evaluations);
+  }
+
+  /**
+   * Returns a dataset containing only selected point indexes.
+   *
+   * @param source source dataset
+   * @param indexes zero-based point indexes
+   * @param name descriptive subset name
+   * @return immutable dataset subset retaining source provenance
+   */
+  public static Dataset subset(Dataset source, List<Integer> indexes, String name) {
+    if (source == null || indexes == null || indexes.isEmpty() || name == null || name.trim().isEmpty()) {
+      throw new IllegalArgumentException("Source, indexes, and subset name are required");
+    }
+    List<Point> points = new ArrayList<Point>();
+    for (Integer index : indexes) {
+      if (index == null || index < 0 || index >= source.getPoints().size()) {
+        throw new IllegalArgumentException("Subset index is outside the source dataset");
+      }
+      points.add(source.getPoints().get(index));
+    }
+    return new Dataset(name, source.getCitation(), source.getDoi(), source.getLicense(), points);
+  }
+
+  private double objective(double binaryInteractionParameter) throws Exception {
+    ThermodynamicBenchmark.Prediction prediction = predictionFactory.create(binaryInteractionParameter);
+    double squaredRelativeErrorSum = 0.0;
+    for (Point point : dataset.getPoints()) {
+      double predictedValue = prediction.predict(point);
+      if (!Double.isFinite(predictedValue)) {
+        throw new IllegalStateException("Non-finite prediction during parameter regression");
+      }
+      double relativeError = (predictedValue - point.getExperimentalValue()) / point.getExperimentalValue();
+      squaredRelativeErrorSum += relativeError * relativeError;
+    }
+    return squaredRelativeErrorSum / dataset.getPoints().size();
+  }
+}

src/main/java/neqsim/thermo/util/benchmark/H2CO2PhaseEquilibriumData.java

@@ -0,0 +1,73 @@
+package neqsim.thermo.util.benchmark;
+
+import java.io.BufferedReader;
+import java.io.IOException;
+import java.io.InputStream;
+import java.io.InputStreamReader;
+import java.nio.charset.StandardCharsets;
+import java.util.ArrayList;
+import java.util.LinkedHashMap;
+import java.util.List;
+import java.util.Map;
+import neqsim.thermo.util.benchmark.ThermodynamicBenchmark.Dataset;
+import neqsim.thermo.util.benchmark.ThermodynamicBenchmark.Point;
+import neqsim.thermo.util.benchmark.ThermodynamicBenchmark.Property;
+
+/** Published H2-CO2 and H2-N2-CO2 phase-equilibrium benchmark datasets. */
+public final class H2CO2PhaseEquilibriumData {
+  private static final String RESOURCE = "/data/thermo/benchmark/zhang2026_h2_co2_phase_equilibrium.csv";
+  private static final String DOI = "10.1063/5.0288386";
+  private static final String CITATION = "Zhang et al. (2026), Measurement of phase equilibrium characteristics and "
+      + "equation-of-state applicability for hydrogen-containing CO2 systems, "
+      + "International Journal of Fluid Engineering 3, 013903.";
+  private static final String LICENSE = "Values transcribed from Tables IV-V of the cited open article; cite the original source.";
+
+  private H2CO2PhaseEquilibriumData() {
+  }
+
+  /**
+   * Loads all published binary and ternary bubble- and dew-point data.
+   *
+   * <p>
+   * The article does not report pointwise pressure uncertainty. The framework therefore records uncertainty as
+   * unavailable rather than inventing a value.
+   * </p>
+   *
+   * @return phase-equilibrium dataset containing 24 points
+   * @throws IOException if the packaged resource cannot be read
+   */
+  public static Dataset load() throws IOException {
+    InputStream input = H2CO2PhaseEquilibriumData.class.getResourceAsStream(RESOURCE);
+    if (input == null) {
+      throw new IOException("Missing benchmark resource " + RESOURCE);
+    }
+    List<Point> points = new ArrayList<Point>();
+    try (BufferedReader reader = new BufferedReader(new InputStreamReader(input, StandardCharsets.UTF_8))) {
+      String line = reader.readLine();
+      if (line == null) {
+        throw new IOException("Empty benchmark resource " + RESOURCE);
+      }
+      while ((line = reader.readLine()) != null) {
+        if (line.trim().isEmpty() || line.trim().startsWith("#")) {
+          continue;
+        }
+        String[] values = line.split(",", -1);
+        if (values.length != 8) {
+          throw new IOException("Expected 8 CSV columns but found " + values.length);
+        }
+        double temperatureK = Double.parseDouble(values[1]) + 273.15;
+        double pressureBara = Double.parseDouble(values[3]) * 10.0;
+        Map<String, Double> composition = new LinkedHashMap<String, Double>();
+        composition.put("CO2", Double.parseDouble(values[4]));
+        composition.put("hydrogen", Double.parseDouble(values[5]));
+        double nitrogenFraction = Double.parseDouble(values[6]);
+        if (nitrogenFraction > 0.0) {
+          composition.put("nitrogen", nitrogenFraction);
+        }
+        Property property = "bubble".equals(values[2]) ? Property.BUBBLE_POINT_PRESSURE : Property.DEW_POINT_PRESSURE;
+        points.add(new Point(property, temperatureK, pressureBara, pressureBara, Double.NaN, "bara", composition));
+      }
+    }
+    return new Dataset("Zhang 2026 hydrogen-containing CO2 phase equilibrium", CITATION, DOI, LICENSE, points);
+  }
+}