From bfc8706d9ab134c3a096f2e81172ff1c8614a90a Mon Sep 17 00:00:00 2001 From: bpota Date: Fri, 20 Feb 2026 16:10:27 +0000 Subject: [PATCH 1/4] Adding thermal conductivity benchmark --- .../user_guide/benchmarks/bulk_crystal.rst | 78 ++ .../analyse_thermal_conductivity.py | 648 ++++++++++++++ .../calc_thermal_conductivity.py | 356 ++++++++ .../thermal_conductivity/collect_kappas.py | 74 ++ .../data/collect_ref_kappas.py | 73 ++ .../data/reference_generation.py | 360 ++++++++ .../thermal_conductivity.py | 823 ++++++++++++++++++ 7 files changed, 2412 insertions(+) create mode 100644 ml_peg/analysis/bulk_crystal/thermal_conductivity/analyse_thermal_conductivity.py create mode 100644 ml_peg/calcs/bulk_crystal/thermal_conductivity/calc_thermal_conductivity.py create mode 100644 ml_peg/calcs/bulk_crystal/thermal_conductivity/collect_kappas.py create mode 100644 ml_peg/calcs/bulk_crystal/thermal_conductivity/data/collect_ref_kappas.py create mode 100644 ml_peg/calcs/bulk_crystal/thermal_conductivity/data/reference_generation.py create mode 100644 ml_peg/calcs/bulk_crystal/thermal_conductivity/thermal_conductivity.py diff --git a/docs/source/user_guide/benchmarks/bulk_crystal.rst b/docs/source/user_guide/benchmarks/bulk_crystal.rst index 26484416c..80d751589 100644 --- a/docs/source/user_guide/benchmarks/bulk_crystal.rst +++ b/docs/source/user_guide/benchmarks/bulk_crystal.rst @@ -482,6 +482,84 @@ Reference data: `arXiv:2307.10072 `_ +Thermal conductivity +==================== + +Summary +------- + +Performance in evaluating lattice thermal conductivity for 103 binary compounds in +rocksalt, zincblende, and wurtzite structures. + +Metrics +------- + +1. κSRE (PBE) + +Mean symmetric relative error (SRE) between predicted and PBE reference lattice +thermal conductivity values. + +DFT-optimised structures are relaxed with fixed symmetries to 0.0001 eV/Å forces. +Interatomic force constants are calculated using finite differences with a 0.03 Å +displacement. The Wigner transport equation is solved within the single-mode +relaxation-time approximation (RTA) to obtain lattice thermal conductivity. The +symmetric relative error in the thermal conductivity at 300 K for a single material +is calculated as: + +.. math:: + + \operatorname{SRE}(\kappa_i) = + 2\frac{|\kappa_i^{\mathrm{pred}} - \kappa_i^{\mathrm{ref}}|} + {\kappa_i^{\mathrm{pred}} + \kappa_i^{\mathrm{ref}}} + +where :math:`\kappa^{\mathrm{pred}}` and :math:`\kappa^{\mathrm{ref}}` are the +predicted and reference thermal conductivities for material :math:`i`. The mean SRE +is calculated across all materials in the dataset. + +2. κSRME (PBE) + +Mean symmetric relative mean error (SRME) between predicted and PBE reference lattice +phonon-mode contributions to thermal conductivity. Thermal conductivity is calculated +as above. The symmetric relative mean error in the phonon-mode contributions to +thermal conductivity at 300 K for a single material is calculated as: + +.. math:: + + \operatorname{SRME}(\kappa_i) = + 2\frac{\sum_j |\kappa_i^{\mathrm{pred}}(q)_s - + \kappa_i^{\mathrm{ref}}(q)_s|} + {\kappa_i^{\mathrm{pred}} + \kappa_i^{\mathrm{ref}}} + +where :math:`\kappa_i^{\mathrm{pred}}(q)_s` and +:math:`\kappa_i^{\mathrm{ref}}(q)_s` are the predicted and reference thermal +conductivity contributions from phonon mode :math:`s` at wavevector :math:`q` for +material :math:`i`. The mean SRME is calculated across all materials in the dataset. + +Both metrics are based on Póta, B., Ahlawat, P., Csányi, G., & Simoncelli, M. (2024). +Thermal conductivity predictions with foundation atomistic models. arXiv preprint +arXiv:2408.00755. + +Computational cost +------------------ + +High: tests are likely to take hours to run on both CPU and GPU. Thermal conductivity +evaluation requires a CPU. + +Data availability +----------------- + +Input structures: + +* DFT-optimised structures from the reference dataset. + +Reference data: + +* DFT (PBE) data +* A. Togo, L. Chaput, and I. Tanaka, Phys. Rev. B, 91, 094306 (2015), and + A. Seko et al., Phys. Rev. Lett., 115, 205901 (2015). +* https://github.com/atztogo/phonondb/blob/main/mdr/phono3py_103compounds_fd_PBE/README.md + + Phonons ======= diff --git a/ml_peg/analysis/bulk_crystal/thermal_conductivity/analyse_thermal_conductivity.py b/ml_peg/analysis/bulk_crystal/thermal_conductivity/analyse_thermal_conductivity.py new file mode 100644 index 000000000..2960852d6 --- /dev/null +++ b/ml_peg/analysis/bulk_crystal/thermal_conductivity/analyse_thermal_conductivity.py @@ -0,0 +1,648 @@ +"""Analyse thermal conductivity benchmark.""" + +from __future__ import annotations + +import os +from pathlib import Path +import traceback + +from ase.io import read +import h5py +import numpy as np +import pandas as pd +import pytest + +from ml_peg.analysis.utils.decorators import build_table, plot_parity +from ml_peg.analysis.utils.utils import load_metrics_config +from ml_peg.app import APP_ROOT +from ml_peg.calcs import CALCS_ROOT +from ml_peg.calcs.bulk_crystal.thermal_conductivity import thermal_conductivity as tc +from ml_peg.models.get_models import get_model_names +from ml_peg.models.models import current_models + +MODELS = get_model_names(current_models) +CALC_PATH = CALCS_ROOT / "bulk_crystal" / "thermal_conductivity" / "outputs" +REF_PATH = CALCS_ROOT / "bulk_crystal" / "thermal_conductivity" / "data" +OUT_PATH = APP_ROOT / "data" / "bulk_crystal" / "thermal_conductivity" + +METRICS_CONFIG_PATH = Path(__file__).with_name("metrics.yaml") +DEFAULT_THRESHOLDS, DEFAULT_TOOLTIPS, DEFAULT_WEIGHTS = load_metrics_config( + METRICS_CONFIG_PATH +) + +STRUCTURE_FILE = REF_PATH / "phononDB-PBE-structures.extxyz" + + +def get_system_names() -> list[str]: + """ + Get list of thermal conductivity system names. + + Returns + ------- + list[str] + List of system names from structure files. + """ + system_names = [] + system_ids = [] + + atoms_list = read(STRUCTURE_FILE, index=":") + + for atoms in atoms_list: + system_names.append(f"{atoms.info['name']}") + system_ids.append(f"{atoms.info[tc.TCKeys.mat_id]}") + return [x for _, x in sorted(zip(system_ids, system_names, strict=False))] + + +def get_system_ids() -> list[str]: + """ + Get list of thermal conductivity system IDs. + + Returns + ------- + list[str] + List of system IDs from structure files. + """ + system_ids = [] + atoms_list = read(STRUCTURE_FILE, index=":") + for atoms in atoms_list: + system_ids.append(f"{atoms.info[tc.TCKeys.mat_id]}") + return sorted(system_ids) + + +def calc_kappa_metrics_from_dfs( + df_pred: pd.DataFrame, df_true: pd.DataFrame +) -> pd.DataFrame: + """ + Compute per-material thermal-conductivity prediction metrics from two dataframes. + + This function takes raw ML predictions and DFT reference results and computes + benchmark metrics. It handles array-type columns (e.g., stress tensors and + mode-resolved properties), calculates averaged quantities, and computes SRD, + SRE, and SRME. + + Parameters + ---------- + df_pred : pandas.DataFrame + DataFrame containing ML model predictions with columns for thermal + conductivity tensors, mode-resolved properties, and other structural + information. + df_true : pandas.DataFrame + DataFrame containing DFT reference calculations with the same structure + as `df_pred`. + + Returns + ------- + pandas.DataFrame + A copy of `df_pred` with additional columns for benchmark metrics, e.g.: + + - SRD: Symmetric Relative Difference between ML and DFT conductivities + - SRE: Absolute value of SRD + - SRME: Mode-resolved error + - DFT_kappa_tot_avg: Reference DFT conductivity values + """ + # Remove precomputed columns + df_pred[tc.TCKeys.kappa_tot_avg] = df_pred[tc.TCKeys.kappa_tot_rta].map( + tc.calculate_kappa_avg + ) + + df_pred[tc.TCKeys.srd] = ( + 2 + * (df_pred[tc.TCKeys.kappa_tot_avg] - df_true[tc.TCKeys.kappa_tot_avg]) + / (df_pred[tc.TCKeys.kappa_tot_avg] + df_true[tc.TCKeys.kappa_tot_avg]) + ) + + # turn temperature list to the first temperature (300K) TODO: allow multiple + # temperatures to be tested + df_pred[tc.TCKeys.srd] = df_pred[tc.TCKeys.srd].map( + lambda x: x if isinstance(x, float) else x[0] + ) + + # We substitute NaN values with 0 predicted conductivity, yielding -2 for SRD + df_pred[tc.TCKeys.srd] = df_pred[tc.TCKeys.srd].fillna(-2) + + df_pred[tc.TCKeys.sre] = df_pred[tc.TCKeys.srd].abs() + + df_pred[tc.TCKeys.srme] = calc_kappa_srme_dataframes(df_pred, df_true) + + df_pred[tc.TCKeys.true_kappa_tot_avg] = df_true[tc.TCKeys.kappa_tot_avg] + + return df_pred + + +def calc_kappa_srme_dataframes( + df_pred: pd.DataFrame, df_true: pd.DataFrame +) -> list[float]: + """ + Calculate the Symmetric Relative Mean Error (SRME) for each material. + + SRME is a comprehensive metric that evaluates both the overall accuracy of thermal + conductivity predictions and the accuracy of individual phonon mode contributions. + It is symmetric (like SRD) to treat over- and under-predictions equally, and + accounts for the mean error across all phonon modes weighted by their contributions. + + The function handles various edge cases: + - Returns 2.0 for materials with imaginary frequencies (unphysical predictions) + - Returns 2.0 for materials where symmetry is broken during relaxation + - Returns 2.0 for failed calculations or missing data + + Parameters + ---------- + df_pred : pd.DataFrame + ML predictions including mode-resolved properties. + df_true : pd.DataFrame + DFT reference data including mode-resolved properties. + + Returns + ------- + list[float]: + SRME values for each material. Values are between 0 and 2, where: + - 0 indicates perfect agreement in both total κ and mode-resolved properties + - 2 indicates complete failure (imaginary frequencies, broken symmetry, etc.) + - Values in between indicate partial agreement, with lower being better + """ + srme_list: list[float] = [] + for idx, row_pred in df_pred.iterrows(): + row_true = df_true.loc[idx] + + # NOTE code below just until before return used to be wrapped in try/except in + # which case SRME=2 was set for the failing material + if row_pred.get(tc.TCKeys.has_imag_ph_modes) is True: + srme_list.append(2) + continue + if relaxed_space_group_number := row_pred.get(tc.TCKeys.final_spg_num): + if initial_space_group_number := row_pred.get(tc.TCKeys.init_spg_num): + if relaxed_space_group_number != initial_space_group_number: + srme_list.append(2) + continue + elif relaxed_space_group_number != row_true.get(tc.TCKeys.spg_num): + srme_list.append(2) + continue + result = calc_kappa_srme(row_pred, row_true) + srme_list.append(float(result[0])) # append the first temperature's SRME + + return srme_list + + +def calc_kappa_srme(kappas_pred: pd.Series, kappas_true: pd.Series) -> np.ndarray: + """ + Calculate the Symmetric Relative Mean Error (SRME) for a single material. + + SRME = 2 * (sum|κ_pred,i - κ_true,i| * w_i) / (κ_pred,tot + κ_true,tot) + + where: + - κ_pred,i and κ_true,i are mode-resolved conductivities for mode i + - w_i are the mode weights + - κ_pred,tot and κ_true,tot are total conductivities + + The calculation involves: + 1. Computing mode-resolved average conductivities if not pre-computed + 2. Calculating the weighted mean absolute error across all modes + 3. Normalizing by the sum of total conductivities to make it symmetric and relative + + Parameters + ---------- + kappas_pred : Series + Series containing ML predictions including: + - kappa_tot_avg: Average total conductivity + - mode_kappa_tot: Mode-resolved conductivities + - mode_weights: Mode weights for averaging. + kappas_true : Series + Series containing DFT reference data with same structure. + + Returns + ------- + np.ndarray + SRME values per temperature, each between 0 and 2, where: + - 0 indicates perfect agreement in both total κ and mode-resolved properties + - 2 indicates complete disagreement or invalid results + On error conditions (missing data, NaN values), returns np.array([2.0]). + """ + if np.any(np.isnan(kappas_true[tc.TCKeys.kappa_tot_avg])): + raise ValueError("found NaNs in kappa_tot_avg reference values") + if ( # return highest possible SRME=2 if any of these conditions are met: + # only have NaN averaged kappa preds + np.all(np.isnan(kappas_pred[tc.TCKeys.kappa_tot_avg])) + # some mode-resolved kappa preds are NaN + or np.any(np.isnan(kappas_pred[tc.TCKeys.kappa_tot_rta])) + # some mode weights are NaN + or np.any(np.isnan(kappas_pred[tc.TCKeys.mode_weights])) + ): + return np.array([2.0]) + + mode_kappa_tot_avgs = {} # store results for pred and true + # Try different data sources in order of preference for both pred and true data + for label, kappas in {"preds": kappas_pred, "true": kappas_true}.items(): + keys = set(kappas.keys()) + if tc.TCKeys.mode_kappa_tot_avg in kappas: + kappas = kappas[tc.TCKeys.mode_kappa_tot_avg] + elif tc.TCKeys.mode_kappa_tot_rta in kappas: + kappas = tc.calculate_kappa_avg(kappas[tc.TCKeys.mode_kappa_tot_rta]) + elif { + tc.TCKeys.kappa_p_rta, + tc.TCKeys.kappa_c, + tc.TCKeys.heat_capacity, + } <= keys: + kappas = tc.calculate_kappa_avg( + tc.calc_mode_kappa_tot( + kappas[tc.TCKeys.kappa_p_rta], + kappas[tc.TCKeys.kappa_c], + kappas[tc.TCKeys.heat_capacity], + ) + ) + else: + raise ValueError( + f"Neither mode_kappa_tot_avg, mode_kappa_tot nor individual kappa\n" + f"components found in {label}, got\n{keys}" + ) + mode_kappa_tot_avgs[label] = np.asarray(kappas) + + # calculating microscopic error for all temperatures + microscopic_error = ( + np.abs(mode_kappa_tot_avgs["preds"] - mode_kappa_tot_avgs["true"]).sum( + axis=tuple(range(1, np.asarray(mode_kappa_tot_avgs["preds"]).ndim)) + ) + / np.asarray(kappas_pred[tc.TCKeys.mode_weights]).sum() + ) + + denominator = ( + kappas_pred[tc.TCKeys.kappa_tot_avg] + kappas_true[tc.TCKeys.kappa_tot_avg] + ) + return 2 * microscopic_error / denominator + + +@pytest.fixture +def kappa_stats() -> dict[str, pd.DataFrame]: + """ + Get thermal conductivity statistics for all models. + + Returns + ------- + dict[str, pd.DataFrame] + Dictionary of DataFrames containing thermal conductivity statistics for + each model. + """ + results: dict[str, pd.DataFrame] = {} + # ref_df = pd.read_json(REF_PATH / "FT_stats.json", orient="index") + + ref_df = pd.DataFrame.from_dict( + tc.hdf5_to_dict(h5py.File(REF_PATH / "PBE" / "kappas.hdf5", "r")), + orient="index", + ) + + fast_ref_df = pd.DataFrame.from_dict( + tc.hdf5_to_dict(h5py.File(REF_PATH / "PBE" / "fast_kappas.hdf5", "r")), + orient="index", + ) + + results["ref"] = ref_df + + for model_name in MODELS: + model_dir = CALC_PATH / model_name + pred_df = None + fast_pred_df = None + + if not model_dir.exists(): + continue + + if (model_dir / "kappas.hdf5").exists(): + current_pred_dict = tc.hdf5_to_dict( + h5py.File(model_dir / "kappas.hdf5", "r") + ) + pred_df = pd.DataFrame.from_dict(current_pred_dict, orient="index") + elif (model_dir / "kappas.json.gz").exists(): + pred_df = pd.read_json(model_dir / "kappas.json.gz", orient="index") + else: + subdirs = sorted(os.listdir(model_dir)) + + if not any(not (Path(d) / "kappa.hdf5").exists() for d in subdirs): + try: + current_pred_dict = tc.load_hdf5_subdir_dicts( + model_dir, "kappa.hdf5" + ) + except Exception: + print(f"Error loading kappas for {model_name}...") + traceback.print_exc() + + if (model_dir / "fast_kappas.hdf5").exists(): + current_pred_dict = tc.hdf5_to_dict( + h5py.File(model_dir / "fast_kappas.hdf5", "r") + ) + fast_pred_df = pd.DataFrame.from_dict(current_pred_dict, orient="index") + elif (model_dir / "fast_kappas.json.gz").exists(): + fast_pred_df = pd.read_json( + model_dir / "fast_kappas.json.gz", orient="index" + ) + else: + subdirs = sorted(os.listdir(model_dir)) + + if not any(not (Path(d) / "fast_kappa.hdf5").exists() for d in subdirs): + try: + current_pred_dict = tc.load_hdf5_subdir_dicts( + model_dir, "fast_kappa.hdf5" + ) + except Exception: + print(f"Error loading fast_kappas for {model_name}...") + traceback.print_exc() + + if pred_df is None and fast_pred_df is None: + continue + + if pred_df is not None and fast_pred_df is None: + pred_df = calc_kappa_metrics_from_dfs(pred_df, ref_df) + elif pred_df is None: + pred_df = calc_kappa_metrics_from_dfs(fast_pred_df, fast_ref_df) + pred_df["fast_sre"] = pred_df[tc.TCKeys.sre] + pred_df["fast_srme"] = pred_df[tc.TCKeys.srme] + pred_df.drop(columns=[tc.TCKeys.sre, tc.TCKeys.srme], inplace=True) + elif fast_pred_df is not None: + pred_df = calc_kappa_metrics_from_dfs(pred_df, ref_df) + fast_pred_df = calc_kappa_metrics_from_dfs(fast_pred_df, fast_ref_df) + pred_df["fast_sre"] = fast_pred_df[tc.TCKeys.sre] + pred_df["fast_srme"] = fast_pred_df[tc.TCKeys.srme] + else: + print( + f"Unexpected case for model {model_name} with pred_df is not None: " + f"{pred_df is not None} and fast_pred_df is not None: " + f"{fast_pred_df is not None}" + ) + continue + + results[model_name] = pred_df + + return results + + +@pytest.fixture +@plot_parity( + filename=OUT_PATH / "figure_conductivity.json", + title="Thermal Conductivity Parity Plot", + x_label="Predicted κ (W/mK)", + y_label="Reference κ (W/mK)", + hoverdata={ + "System": get_system_names(), + "System ID": get_system_ids(), + }, +) +def conductivity(kappa_stats: dict[str, pd.DataFrame]) -> dict[str, list]: + """ + Get thermal conductivity parity plot data. + + Parameters + ---------- + kappa_stats : dict[str, pd.DataFrame] + Dictionary of DataFrames containing thermal conductivity statistics for + each model. + + Returns + ------- + dict[str, list] + Dictionary mapping model names to lists of thermal conductivity values. + """ + conductivity_data = {} + for model_name in MODELS: + if model_name not in kappa_stats: + continue + df = kappa_stats[model_name] + if tc.TCKeys.kappa_tot_avg not in df: + continue + conductivity_data[model_name] = df[tc.TCKeys.kappa_tot_avg].tolist() + conductivity_data["ref"] = kappa_stats["ref"][tc.TCKeys.kappa_tot_avg].tolist() + return conductivity_data + + +@pytest.fixture +def mean_sre(kappa_stats: dict[str, pd.DataFrame]) -> dict[str, float]: + """ + Calculate mean Symmetric Relative Error (SRE) for each model. + + Parameters + ---------- + kappa_stats : dict[str, pd.DataFrame] + Dictionary of DataFrames containing thermal conductivity statistics for + each model. + + Returns + ------- + dict[str, float] + Dictionary mapping model names to their mean SRE values. + """ + sre_values = {} + for model_name in MODELS: + if model_name not in kappa_stats: + continue + df = kappa_stats[model_name] + if tc.TCKeys.sre not in df: + continue + sre_values[model_name] = df[tc.TCKeys.sre].mean() + return sre_values + + +@pytest.fixture +def mean_srme(kappa_stats: dict[str, pd.DataFrame]) -> dict[str, float]: + """ + Calculate mean Symmetric Relative Mean Error (SRME) for each model. + + Parameters + ---------- + kappa_stats : dict[str, pd.DataFrame] + Dictionary of DataFrames containing thermal conductivity statistics for + each model. + + Returns + ------- + dict[str, float] + Dictionary mapping model names to their mean SRME values. + """ + srme_values = {} + for model_name in MODELS: + if model_name not in kappa_stats: + continue + df = kappa_stats[model_name] + if tc.TCKeys.srme not in df: + continue + srme_values[model_name] = df[tc.TCKeys.srme].mean() + return srme_values + + +@pytest.fixture +def instability(kappa_stats: dict[str, pd.DataFrame]) -> dict[str, float]: + """ + Calculate instability metric for each model. + + Parameters + ---------- + kappa_stats : dict[str, pd.DataFrame] + Dictionary of DataFrames containing thermal conductivity statistics for + each model. + + Returns + ------- + dict[str, float] + Dictionary mapping model names to their instability values + (fraction of unstable materials). + """ + instability_values = {} + ref_length = len(kappa_stats["ref"]) + for model_name in MODELS: + if model_name not in kappa_stats: + continue + df = kappa_stats[model_name] + if tc.TCKeys.has_imag_ph_modes not in df: + continue + has_imag_ph_modes = df[tc.TCKeys.has_imag_ph_modes].values + has_imag_ph_modes_filtered = has_imag_ph_modes[ + np.logical_not(np.isnan(has_imag_ph_modes)) + ] + instability_values[model_name] = has_imag_ph_modes_filtered.sum() / ref_length + + return instability_values + + +@pytest.fixture +def failure(kappa_stats: dict[str, pd.DataFrame]) -> dict[str, float]: + """ + Calculate failure metric for each model. + + Parameters + ---------- + kappa_stats : dict[str, pd.DataFrame] + Dictionary of DataFrames containing thermal conductivity statistics for + each model. + + Returns + ------- + dict[str, float] + Dictionary mapping model names to their failure values + (fraction of failed calculations). + """ + failure_values = {} + ref_length = len(kappa_stats["ref"]) + for model_name in MODELS: + if model_name not in kappa_stats: + continue + df = kappa_stats[model_name] + if tc.TCKeys.has_imag_ph_modes not in df: + continue + has_nan_frequency = df[tc.TCKeys.has_imag_ph_modes].isna().values + failure_values[model_name] = has_nan_frequency.sum() / ref_length + + return failure_values + + +@pytest.fixture +def mean_fast_sre(kappa_stats: dict[str, pd.DataFrame]) -> dict[str, float]: + """ + Calculate mean Symmetric Relative Error (SRE) for each model with "fast" settings. + + In "fast" settings the qpoint mesh is reduced in the phonon interaction strength c + alculation to allow quicker benchmarking. + + Parameters + ---------- + kappa_stats : dict[str, pd.DataFrame] + Dictionary of DataFrames containing thermal conductivity statistics for + each model. + + Returns + ------- + dict[str, float] + Dictionary mapping model names to their mean SRE values. + """ + sre_values = {} + for model_name in MODELS: + if model_name not in kappa_stats: + continue + df = kappa_stats[model_name] + if "fast_sre" not in df: + continue + sre_values[model_name] = df["fast_sre"].mean() + return sre_values + + +@pytest.fixture +def mean_fast_srme(kappa_stats: dict[str, pd.DataFrame]) -> dict[str, float]: + """ + Calculate mean Symmetric Relative Mean Error (SRME) with "fast" settings. + + Calculation processes each model. In "fast" settings the qpoint mesh is reduced + in the phonon interaction strength calculation. + + Parameters + ---------- + kappa_stats : dict[str, pd.DataFrame] + Dictionary of DataFrames containing thermal conductivity statistics for + each model. + + Returns + ------- + dict[str, float] + Dictionary mapping model names to their mean SRME values. + """ + srme_values = {} + for model_name in MODELS: + if model_name not in kappa_stats: + continue + df = kappa_stats[model_name] + if "fast_srme" not in df: + continue + srme_values[model_name] = df["fast_srme"].mean() + return srme_values + + +@pytest.fixture +@build_table( + filename=OUT_PATH / "kappas.json", + metric_tooltips=DEFAULT_TOOLTIPS, + thresholds=DEFAULT_THRESHOLDS, +) +def metrics( + mean_srme: dict[str, float], + mean_sre: dict[str, float], + instability: dict[str, float], + failure: dict[str, float], + mean_fast_sre: dict[str, float], + mean_fast_srme: dict[str, float], +) -> dict[str, dict]: + """ + Get all thermal conductivity metrics. + + Parameters + ---------- + mean_srme : dict[str, float] + Mean Symmetric Relative Mean Error for each model. + mean_sre : dict[str, float] + Mean Symmetric Relative Error for each model. + instability : dict[str, float] + Fraction of unstable materials (has imaginary frequencies) for each model. + failure : dict[str, float] + Fraction of failed calculations (NaN frequencies) for each model. + mean_fast_sre : dict[str, float] + Mean Symmetric Relative Error for each model (fast version). + mean_fast_srme : dict[str, float] + Mean Symmetric Relative Mean Error for each model (fast version). + + Returns + ------- + dict[str, dict] + Metric names and values for all models. + """ + return { + "κSRE": mean_sre, + "κSRME": mean_srme, + "Instability": instability, + "Failure": failure, + "Fast κSRE": mean_fast_sre, + "Fast κSRME": mean_fast_srme, + } + + +def test_thermal_conducticity(metrics: dict[str, dict]): + """ + Run thermal conductivity benchmark tests. + + Parameters + ---------- + metrics : dict[str, dict] + Metric names and values for all models. + """ + return diff --git a/ml_peg/calcs/bulk_crystal/thermal_conductivity/calc_thermal_conductivity.py b/ml_peg/calcs/bulk_crystal/thermal_conductivity/calc_thermal_conductivity.py new file mode 100644 index 000000000..80d15243a --- /dev/null +++ b/ml_peg/calcs/bulk_crystal/thermal_conductivity/calc_thermal_conductivity.py @@ -0,0 +1,356 @@ +""" +Run calculations for thermal conductivity tests. + +Code is adapted from https://github.com/MPA2suite/k_SRME/blob/6ff4c867/k_srme/conductivity.py +by Balázs Póta, Paramvir Ahlawat, Gábor Csányi, Michele Simoncelli +and https://github.com/janosh/matbench-discovery/blob/main/matbench_discovery/phonons/thermal_conductivity.py +by Janosh Riebesell. +See https://arxiv.org/abs/2408.00755 for details. +""" + +from __future__ import annotations + +from copy import deepcopy +from pathlib import Path +import sys +import traceback +from typing import Any, Literal +import warnings + +import ase +from ase.constraints import FixSymmetry +from ase.filters import ExpCellFilter, Filter, FrechetCellFilter +from ase.io import read, write +from ase.optimize.optimize import Optimizer +import h5py +import pandas as pd +import pytest +from tqdm import tqdm + +from ml_peg.calcs.bulk_crystal.thermal_conductivity import thermal_conductivity as tc +from ml_peg.models.get_models import load_models +from ml_peg.models.models import current_models + +MODELS = load_models(current_models) + +DATA_PATH = Path(__file__).parent / "data" +OUT_PATH = Path(__file__).parent / "outputs" + +if len(sys.argv) >= 3: + try: + PARALLEL_TASK_ID = int(sys.argv[1]) + PARALLEL_TASK_NUM = int(sys.argv[2]) + except ValueError: + print( + "Could not parse parallel task arguments, running in serial.", + file=sys.stderr, + ) + PARALLEL_TASK_ID = 0 + PARALLEL_TASK_NUM = 1 +else: + PARALLEL_TASK_ID = 0 + PARALLEL_TASK_NUM = 1 + + +ase_optimizer: Literal[ + "GPMin", + "GOQN", + "BFGSLineSearch", + "QuasiNewton", + "SciPyFminBFGS", + "BFGS", + "LBFGSLineSearch", + "SciPyFminCG", + "FIRE2", + "FIRE", + "LBFGS", +] = "LBFGS" +ase_filter: Literal["frechet", "exp"] = "frechet" +max_steps = 300 +fmax = 1e-4 # Run until the forces are smaller than this in eV/A +enforce_relax_symm = True # Enforce symmetry during relaxation + +# Symmetry parameters +symprec = 1e-5 # symmetry precision for enforcing relaxation and conductivity calcs + +# Conductivity to be calculated if symmetry group changed during relaxation +conductivity_broken_symm = False +save_forces = True # Save force sets to file +temperatures: list[float] = [300] +displacement_distance = 0.03 +ignore_imaginary_freqs = True + +default_dtype = "float64" + +STRUCTURE_PATH = DATA_PATH / "phononDB-PBE-structures.extxyz" + +FAST_ONLY = True + + +@pytest.mark.parametrize("mlip", MODELS.items()) +def test_thermal_conductivity(mlip: tuple[str, Any]) -> None: + """ + Run thermal conductivity test. + + Parameters + ---------- + mlip + Name of model use and model to get calculator. + """ + model_name, model = mlip + model.default_dtype = default_dtype + calculator = model.get_calculator() + + # Download dataset + # thermal_conductivity_dir = ( + # download_s3_data( + # key="inputs/bulk_crystal/thermal_conductivity/thermal_conductivity.zip", + # filename="thermal_conductivity.zip", + # ) + # / "thermal_conductivity" + # ) + + atoms_list = read(STRUCTURE_PATH, format="extxyz", index=":") + atoms_list = atoms_list[PARALLEL_TASK_ID - 0 :: PARALLEL_TASK_NUM] + + kappa_dicts = {} + fast_kappa_dicts = {} + + for i, atoms_input in enumerate(tqdm(atoms_list, desc="Atoms loop")): + atoms = atoms_input.copy() + + atoms.calc = calculator + + structure_id = atoms.info.get(tc.TCKeys.mat_id, f"structure_{i}") + + out_dir = OUT_PATH / model_name / structure_id + out_dir.mkdir(parents=True, exist_ok=True) + + # Relax structure before calculating thermal conductivity + relax_path = out_dir / "relaxed.extxyz" + + formula = atoms.info.get("name", atoms.get_chemical_formula()) + + mat_id = atoms.info[tc.TCKeys.mat_id] + init_info = deepcopy(atoms.info) + info_dict: dict[str, Any] = { + str(tc.TCKeys.mat_id): mat_id, + str(tc.TCKeys.formula): formula, + } + err_dict: dict[str, list[str]] = {"errors": [], "error_traceback": []} + + # Select filter class + if ase_filter in {"frechet", "exp"}: + filter_cls: type[Filter] = { + "frechet": FrechetCellFilter, + "exp": ExpCellFilter, + }[ase_filter] + else: + # Default to FrechetCellFilter if not specified (for MACE compatibility) + filter_cls = FrechetCellFilter + + # Select optimizer class + optimizer_dict = { + "GPMin": ase.optimize.GPMin, + "GOQN": ase.optimize.GoodOldQuasiNewton, + "BFGSLineSearch": ase.optimize.BFGSLineSearch, + "QuasiNewton": ase.optimize.BFGSLineSearch, + "SciPyFminBFGS": ase.optimize.sciopt.SciPyFminBFGS, + "BFGS": ase.optimize.BFGS, + "LBFGSLineSearch": ase.optimize.LBFGSLineSearch, + "SciPyFminCG": ase.optimize.sciopt.SciPyFminCG, + "FIRE2": ase.optimize.FIRE2, + "FIRE": ase.optimize.FIRE, + "LBFGS": ase.optimize.LBFGS, + } + optim_cls: type[Optimizer] = optimizer_dict[ase_optimizer] + + # Initialize variables that might be needed in error handling + relax_dict: dict[str, Any] = { + "max_stress": None, + "reached_max_steps": False, + "broken_symmetry": False, + } + # initial space group for symmetry breaking detection + init_spg_num = tc.get_spacegroup_number_from_atoms(atoms, symprec=symprec) + fast_results_dict = None + + # Relaxation + try: + results_dict = {} + atoms.calc = calculator + if max_steps > 0: + if enforce_relax_symm: + atoms.set_constraint(FixSymmetry(atoms)) + filtered_atoms = filter_cls(atoms, mask=[True] * 3 + [False] * 3) + else: + filtered_atoms = filter_cls(atoms) + + optimizer = optim_cls(filtered_atoms, logfile=out_dir / "relax.log") + optimizer.run(fmax=fmax, steps=max_steps) + + step_count = getattr( + optimizer, "nsteps", None + ) # Get optimizer step count + if ( + step_count is None + ): # fallback to extract from state_dict if available + state = getattr(optimizer, "state_dict", dict)() + step_count = state.get("step", 0) + + reached_max_steps = step_count >= max_steps + if reached_max_steps: + print(f"Material {mat_id=} reached {max_steps=} during relaxation") + + # maximum residual stress component in for xx,yy,zz and xy,yz,xz + # components separately result is a array of 2 elements + max_stress = atoms.get_stress().reshape((2, 3), order="C").max(axis=1) + + atoms.calc = None + atoms.constraints = None + atoms.info = init_info | atoms.info + + # Check if symmetry was broken during relaxation + relaxed_spg_num = tc.get_spacegroup_number_from_atoms( + atoms, symprec=symprec + ) + broken_symmetry = init_spg_num != relaxed_spg_num + + relax_dict = { + "max_stress": max_stress, + "reached_max_steps": reached_max_steps, + "broken_symmetry": broken_symmetry, + tc.TCKeys.final_spg_num: relaxed_spg_num, + tc.TCKeys.init_spg_num: init_spg_num, + } + + except (ValueError, RuntimeError, OSError, KeyError) as exc: + warnings.warn( + f"Failed to relax {formula=}, {mat_id=}: {exc!r}", stacklevel=2 + ) + traceback.print_exc() + err_dict["errors"] += [f"RelaxError: {exc!r}"] + err_dict["error_traceback"] += [traceback.format_exc()] + results_dict = info_dict | relax_dict | err_dict + + write(relax_path, atoms, format="extxyz") + + try: + ph3 = tc.init_phono3py( + atoms, + fc2_supercell=atoms.info["fc2_supercell"], + fc3_supercell=atoms.info["fc3_supercell"], + q_point_mesh=atoms.info["q_point_mesh"], + displacement_distance=displacement_distance, + symprec=symprec, + ) + + ph3, fc2_set, freqs = tc.get_fc2_and_freqs( + ph3, calculator=calculator, pbar_kwargs={"disable": True} + ) + + has_imag_ph_modes = tc.check_imaginary_freqs(freqs) + freqs_dict = { + tc.TCKeys.has_imag_ph_modes: has_imag_ph_modes, + tc.TCKeys.ph_freqs: freqs, + } + + # Determine if conductivity calculation should proceed + if ignore_imaginary_freqs: + # NequIP/Allegro mode: ignore imaginary frequencies + ltc_condition = True + else: + # MACE mode: check both imaginary freqs and broken symmetry + broken_symmetry = relax_dict.get("broken_symmetry", False) + ltc_condition = not has_imag_ph_modes and ( + not broken_symmetry or conductivity_broken_symm + ) + + if ltc_condition: + # fc3_set = tc.calculate_fc3_set( + # ph3, calculator=calculator, pbar_kwargs={"leave": False} + # ) + ph3.produce_fc3(symmetrize_fc3r=True) + else: + reason = [] + if has_imag_ph_modes: + reason.append("imaginary frequencies") + if relax_dict.get("broken_symmetry") and not conductivity_broken_symm: + reason.append("broken symmetry") + warnings.warn( + f"{' and '.join(reason).capitalize()} detected for {mat_id}, " + f"skipping FC3 and LTC calculation!", + stacklevel=2, + ) + + if not ltc_condition: + results_dict = info_dict | relax_dict | freqs_dict | err_dict + + except (ValueError, RuntimeError, OSError, KeyError) as exc: + warnings.warn( + f"Failed to calculate force sets {mat_id}: {exc!r}", stacklevel=2 + ) + traceback.print_exc() + err_dict["errors"] += [f"ForceConstantError: {exc!r}"] + err_dict["error_traceback"] += [traceback.format_exc()] + results_dict = info_dict | relax_dict | err_dict + + # Calculation of conductivity + try: + with tc.tqdm_gridpoints(desc="Conducitivity calc"): + ph3.mesh_numbers = atoms.info["fast_q_point_mesh"] + ph3, fast_kappa_dict, _cond = tc.calculate_conductivity( + ph3, temperatures=temperatures, log_level=2 + ) + fast_results_dict = ( + info_dict | relax_dict | freqs_dict | fast_kappa_dict | err_dict + ) + + if not FAST_ONLY: + with tc.tqdm_gridpoints(desc="Conducitivity calc"): + ph3.mesh_numbers = atoms.info["q_point_mesh"] + ph3, kappa_dict, _cond = tc.calculate_conductivity( + ph3, temperatures=temperatures, log_level=2 + ) + results_dict = ( + info_dict | relax_dict | freqs_dict | kappa_dict | err_dict + ) + + except (ValueError, RuntimeError, OSError, KeyError) as exc: + warnings.warn( + f"Failed to calculate conductivity {mat_id}: {exc!r}", stacklevel=2 + ) + traceback.print_exc() + err_dict["errors"] += [f"ConductivityError: {exc!r}"] + err_dict["error_traceback"] += [traceback.format_exc()] + results_dict = info_dict | relax_dict | freqs_dict | err_dict + + if fast_results_dict is None: + fast_results_dict = relax_dict + results_dict[tc.TCKeys.mat_id] = structure_id + fast_results_dict[tc.TCKeys.mat_id] = structure_id + + if not FAST_ONLY: + results_path = out_dir / "kappa.hdf5" + with h5py.File(results_path, "w") as f: + tc.dict_to_hdf5(results_dict, f) + kappa_dicts[structure_id] = results_dict + + fast_results_path = out_dir / "fast_kappa.hdf5" + with h5py.File(fast_results_path, "w") as f: + tc.dict_to_hdf5(fast_results_dict, f) + fast_kappa_dicts[structure_id] = fast_results_dict + + if PARALLEL_TASK_NUM == 1: + if not FAST_ONLY: + df = pd.DataFrame(kappa_dicts).T + df.index.name = tc.TCKeys.mat_id + df.reset_index().to_json(OUT_PATH / model_name / "kappas.json.gz") + with h5py.File(OUT_PATH / model_name / "kappas.hdf5", "w") as f: + tc.dict_to_hdf5(kappa_dicts, f) + + fast_df = pd.DataFrame(fast_kappa_dicts).T + fast_df.index.name = tc.TCKeys.mat_id + fast_df.reset_index().to_json(OUT_PATH / model_name / "fast_kappas.json.gz") + with h5py.File(OUT_PATH / model_name / "fast_kappas.hdf5", "w") as f: + tc.dict_to_hdf5(fast_kappa_dicts, f) diff --git a/ml_peg/calcs/bulk_crystal/thermal_conductivity/collect_kappas.py b/ml_peg/calcs/bulk_crystal/thermal_conductivity/collect_kappas.py new file mode 100644 index 000000000..fdae9472c --- /dev/null +++ b/ml_peg/calcs/bulk_crystal/thermal_conductivity/collect_kappas.py @@ -0,0 +1,74 @@ +"""Collects thermal conductivity data from HDF5 files.""" + +from __future__ import annotations + +from pathlib import Path +import sys + +import h5py +import pandas as pd + +from ml_peg.calcs.bulk_crystal.thermal_conductivity import thermal_conductivity as tc + +OUT_PATH = Path(__file__).parent / "outputs" + + +if len(sys.argv) == 1: + raise ValueError( + "Please provide the path to the directory containing the kappa.hdf5 files " + "as a command-line argument." + ) + +models = sys.argv[1:] + + +def collect_kappas(parent_dir: Path, filename_no_ext: str, output_name_no_ext: str): + """ + Collect kappa data from HDF5 files in the directory and save as JSON and HDF5. + + Parameters + ---------- + parent_dir : Path + The directory containing the HDF5 files to collect. + filename_no_ext : str + The base name of the HDF5 files to look for (without extension). + output_name_no_ext : str + The base name to use for the output JSON and HDF5 files (without extension). + """ + print(f"Loading {filename_no_ext} from {parent_dir} subdirectories...") + dicts = tc.load_hdf5_subdir_dicts(parent_dir, filename_no_ext + ".hdf5") + print("Loading finished.") + + df = pd.DataFrame(dicts).T + df.index.name = tc.TCKeys.mat_id + df.reset_index().to_json(OUT_PATH / f"{output_name_no_ext}.json.gz") + with h5py.File(parent_dir / f"{output_name_no_ext}.hdf5", "w") as f: + tc.dict_to_hdf5(dicts, f) + + +for model in models: + model_dir = OUT_PATH / model + if not model_dir.exists(): + print(f"Model directory {model_dir} does not exist. Skipping.") + continue + + subdirs = [d for d in model_dir.iterdir() if d.is_dir()] + + if not any(not (d / "fast_kappa.hdf5").exists() for d in subdirs): + try: + collect_kappas(model_dir, "fast_kappa", "fast_kappa") + except Exception as exc: + print(f"Error collecting fast kappas for {model}: {exc}") + else: + print( + "Fast kappa files not found in all subdirectories." + "Skipping fast kappa collection." + ) + + if not any(not (d / "kappa.hdf5").exists() for d in subdirs): + try: + collect_kappas(model_dir, "kappa", "kappa") + except Exception as exc: + print(f"Error collecting kappas for {model}: {exc}") + else: + print("Kappa files not found in all subdirectories. Skipping kappa collection.") diff --git a/ml_peg/calcs/bulk_crystal/thermal_conductivity/data/collect_ref_kappas.py b/ml_peg/calcs/bulk_crystal/thermal_conductivity/data/collect_ref_kappas.py new file mode 100644 index 000000000..3b927a49d --- /dev/null +++ b/ml_peg/calcs/bulk_crystal/thermal_conductivity/data/collect_ref_kappas.py @@ -0,0 +1,73 @@ +""" +Collects reference thermal conductivity data. + +It loads the data from JSON and HDF5 files, processes it into pandas DataFrames, +and saves it back to disk in both formats for easy access in future analyses. +""" + +from __future__ import annotations + +import importlib +from pathlib import Path + +import h5py +import pandas as pd + +DATA_PATH = Path(__file__).parent + +tc = importlib.import_module(DATA_PATH.parent / "thermal_conductivity.py") + +PBE_DATA_PATH = DATA_PATH / "PBE" + +FAST_JSON_PATH = PBE_DATA_PATH / "fast_kappas.json.gz" +FAST_HDF5_PATH = PBE_DATA_PATH / "fast_kappas.hdf5" + +ORIGINAL_JSON_PATH = PBE_DATA_PATH / "kappas.json.gz" +ORIGINAL_HDF5_PATH = PBE_DATA_PATH / "kappas.hdf5" + + +if FAST_JSON_PATH.exists(): + print(f"Loading fast kappa dicts from {FAST_JSON_PATH}...") + df_fast = pd.read_json(FAST_JSON_PATH) + fast_dicts = df_fast.set_index(tc.TCKeys.mat_id).to_dict(orient="index") + print("Loading finished.") + +elif FAST_HDF5_PATH.exists(): + print(f"Loading fast kappa dicts from {FAST_HDF5_PATH}...") + with h5py.File(FAST_HDF5_PATH, "r") as f: + fast_dicts = tc.hdf5_to_dict(f) + print("Loading finished.") +else: + print(f"Loading fast kappa dicts from {PBE_DATA_PATH} subdirectories...") + fast_dicts = tc.load_hdf5_subdir_dicts(PBE_DATA_PATH, "fast_kappa.hdf5") + print("Loading finished.") + +if ORIGINAL_JSON_PATH.exists(): + print(f"Loading original kappa dicts from {ORIGINAL_JSON_PATH}...") + df_original = pd.read_json(ORIGINAL_JSON_PATH) + original_dicts = df_original.set_index(tc.TCKeys.mat_id).to_dict(orient="index") + print("Loading finished.") + +elif ORIGINAL_HDF5_PATH.exists(): + print(f"Loading original kappa dicts from {ORIGINAL_HDF5_PATH}...") + with h5py.File(ORIGINAL_HDF5_PATH, "r") as f: + original_dicts = tc.hdf5_to_dict(f) + print("Loading finished.") +else: + print(f"Loading original kappa dicts from {PBE_DATA_PATH} subdirectories...") + original_dicts = tc.load_hdf5_subdir_dicts(PBE_DATA_PATH, "kappa.hdf5") + print("Loading finished.") + + +df_fast = pd.DataFrame(fast_dicts).T +df_original = pd.DataFrame(original_dicts).T + +df_fast.index.name = tc.TCKeys.mat_id +df_original.index.name = tc.TCKeys.mat_id +df_fast.reset_index().to_json(FAST_JSON_PATH) +df_original.reset_index().to_json(ORIGINAL_JSON_PATH) + +with h5py.File(FAST_HDF5_PATH, "w") as f: + tc.dict_to_hdf5(fast_dicts, f) +with h5py.File(ORIGINAL_HDF5_PATH, "w") as f: + tc.dict_to_hdf5(original_dicts, f) diff --git a/ml_peg/calcs/bulk_crystal/thermal_conductivity/data/reference_generation.py b/ml_peg/calcs/bulk_crystal/thermal_conductivity/data/reference_generation.py new file mode 100644 index 000000000..f0f6d5e25 --- /dev/null +++ b/ml_peg/calcs/bulk_crystal/thermal_conductivity/data/reference_generation.py @@ -0,0 +1,360 @@ +"""Reference data generation for thermal conductivity calculations using phono3py.""" + +from __future__ import annotations + +import importlib +from io import BytesIO +import os +from pathlib import Path +import re +import sys + +from ase import Atoms +import h5py +from phono3py.cui.load import load +import requests +from tqdm import tqdm + +DATA_PATH = Path(__file__).parent + +sys.path.append(str(DATA_PATH.parent)) +tc = importlib.import_module("thermal_conductivity") + +PH3_PARAMS_PATH = DATA_PATH / "PBE" + +MAT_ID_FILE = DATA_PATH / "mat_id-phonondb.txt" +STRUCTURE_FILE = DATA_PATH / "phononDB-PBE-structures.extxyz" +SCRAPED_YAMLS_FILE = PH3_PARAMS_PATH / "scraped_yamls.txt" + +SCRAPE = False + +FAST_ONLY = True + +SKIP_EXISTING = True + +TEMPERATURES = [300] + + +def get_original_qmesh(mat_id_txt: str, name: str, spg_no: int) -> list[int]: + """ + Determine the original q-point mesh for a given material. + + Uses its name and space group number. + + Parameters + ---------- + mat_id_txt : str + The content of the mat_id-phonondb.txt file as a string. + name : str + The name of the material (e.g., "Si", "Ge", etc.). + spg_no : int + The space group number of the material. + + Returns + ------- + list[int] + The original q-point mesh corresponding to the given space group number. + For example: + -For space group numbers 216 and 225, the original qpoint mesh is [19, 19, 19]. + -For space group number 186, the original qpoint mesh is [19, 19, 15]. + -If the space group number is not recognized, a ValueError is raised. + """ + if spg_no == 216 or spg_no == 225: + return [19, 19, 19] + if spg_no == 186: + return [19, 19, 15] + + raise ValueError(f"Unsupported space group number: {spg_no} for {name}") + + +def get_fast_qmesh(mat_id_txt: str, name: str, spg_no: int) -> list[int]: + """ + Determine the "fast" q-point mesh for a given material. + + Uses its name and space group number. The "fast" q-point mesh is a reduced mesh + used for quicker calculations, and it is determined based on the space group number + of the material. This function checks the space group number and returns the + corresponding fast q-point mesh. + + Parameters + ---------- + mat_id_txt : str + The content of the mat_id-phonondb.txt file as a string. + name : str + The name of the material (e.g., "Si", "Ge", etc.). + spg_no : int + The space group number of the material. + + Returns + ------- + list[int] + The fast q-point mesh corresponding to the given space group number. + For example: + - For space group numbers 216 and 225, the fast q-point mesh is [9, 9, 9]. + - For space group number 186, the fast q-point mesh is [9, 9, 7]. + - If the space group number is not recognized, a ValueError is raised. + """ + if spg_no == 216 or spg_no == 225: + return [9, 9, 9] + if spg_no == 186: + return [9, 9, 7] + + raise ValueError(f"Unsupported space group number: {spg_no} for {name}") + + +def get_mat_id(mat_id_txt: str, name: str, spg_no: int) -> str: + """ + Extract the material ID from the mat_id_txt. + + Uses the material name and space group number. + + Parameters + ---------- + mat_id_txt : str + The content of the mat_id-phonondb.txt file as a string. + name : str + The name of the material (e.g., "Si", "Ge", etc.). + spg_no : int + The space group number of the material. + + Returns + ------- + str + The material ID corresponding to the given name and space group number, or "N/A" + if not found. + """ + for line in mat_id_txt: + if " ".join([name, str(spg_no)]) in line: + return line.split()[-1] # Assuming the MP ID is the first element + return "N/A" + + +def create_yaml(mat_id: str, yaml_str: str) -> Path: + """ + Save the phono3py parameters as a YAML file for a given material ID. + + Parameters + ---------- + mat_id : str + Material ID to use for naming the YAML file and directory. + yaml_str : str + The YAML content as a string to be saved to the file. + + Returns + ------- + Path + The path to the saved YAML file. + """ + os.makedirs(PH3_PARAMS_PATH / mat_id, exist_ok=True) + with open(PH3_PARAMS_PATH / mat_id / "phono3py_params.yaml", "w") as f: + f.write(yaml_str) + return PH3_PARAMS_PATH / mat_id / "phono3py_params.yaml" + + +def scrape_phono3py_data( + phono3py_yaml_path: str = "data/DFT", +) -> tuple[list[str], list[str]]: + """ + Scrape phono3py parameters from the phonondb README and save them as YAML files. + + The README contains links to compressed YAML files with phonon properties + and thermal conductivity parameters for a set of materials. This function + extracts those links, downloads the YAML files, decompresses them, and saves + them locally. It also extracts material names and space group numbers for + later use in matching with material IDs. + + Parameters + ---------- + phono3py_yaml_path : str, optional + Local directory to save the scraped YAML files, by default "data/DFT". + + Returns + ------- + tuple[list[str], list[str]] + A tuple containing a list of YAML strings and a corresponding list of material + names. + """ + raw_url = "https://raw.githubusercontent.com/atztogo/phonondb/refs/heads/main/mdr/phono3py_103compounds_fd_PBE/README.md" + params_link_regex = r"\[\s*phono3py_params\.yaml\.xz\s*\]\((https?://[^)\s]+)\)" + name_regex = r"(?m)^\|\s*[A-Za-z0-9]+-([A-Z][a-z]?(?:[A-Z][a-z]?|\d+)*)\s*\|" + + phonondb_readme_path = PH3_PARAMS_PATH / "phonon3db_readme.md" + if os.path.exists(phonondb_readme_path): + print( + "Found existing phonon3db README at " + f"{phonondb_readme_path}. Loading from file instead of downloading." + ) + with open(phonondb_readme_path) as f: + phonondb_readme = f.read() + else: + r = requests.get(raw_url, timeout=30) + r.raise_for_status() + phonondb_readme = r.text + with open(phonondb_readme_path, "w") as f: + f.write(phonondb_readme) + + params_links = re.findall(params_link_regex, phonondb_readme) + name_list = re.findall(name_regex, phonondb_readme) + + if os.path.exists(SCRAPED_YAMLS_FILE): + print( + "Found existing scraped YAMLs at " + f"{SCRAPED_YAMLS_FILE}. Loading from file instead of downloading." + ) + with open(SCRAPED_YAMLS_FILE) as f: + yaml_list = [y for y in f.read().split("END\n") if y.strip()] + + return yaml_list, name_list + + try: + import lzma + except ImportError: + raise ImportError("lzma module is required to decompress .xz files.") from None + + yaml_list = [] + + for link in tqdm(params_links, desc="Download"): + xz_bytes = requests.get(link, timeout=30).content + + yaml_bytes = lzma.decompress(xz_bytes) + + yaml_str = yaml_bytes.decode("utf-8") + + yaml_list.append(yaml_str) + + return yaml_list, name_list + + +yaml_list, name_list = scrape_phono3py_data() + + +with open(SCRAPED_YAMLS_FILE, "w") as f: + for yaml_str in yaml_list: + f.write(yaml_str + "\nEND\n") + + +with open(MAT_ID_FILE) as f: + mat_id_txt = f.readlines() + + +with open(STRUCTURE_FILE, "w") as f: + f.write("") # Clear the file before appending new structures + +N = len(name_list) + +if N != len(yaml_list): + raise ValueError( + "Number of names " + f"({len(name_list)}) does not match number of YAMLs ({len(yaml_list)})!" + ) +pbar = tqdm(range(N), desc="Processing phono3py") +for i in pbar: + yaml_str = yaml_list[i] + name = name_list[i] + pbar.set_postfix_str(name) + + ph3 = load(BytesIO(yaml_str.encode("utf-8")), produce_fc=True, symmetrize_fc=True) + + symm_no = ph3.symmetry.dataset.number + mat_id = get_mat_id(mat_id_txt, name, symm_no) + + original_qmesh = get_original_qmesh(mat_id_txt, name, symm_no) + fast_qmesh = get_fast_qmesh(mat_id_txt, name, symm_no) + + info_dict = { + "name": name, + "primitive_matrix": ph3.primitive_matrix, + "fc2_supercell": ph3.supercell_matrix, + "fc3_supercell": ph3.supercell_matrix, + "symm.no": ph3.symmetry.dataset.number, + tc.TCKeys.mat_id: mat_id, + "q_point_mesh": original_qmesh, + "fast_q_point_mesh": fast_qmesh, + } + + yaml_file = create_yaml(mat_id, yaml_str) + if SKIP_EXISTING: + if FAST_ONLY and os.path.exists(yaml_file.parent / "fast_kappa.hdf5"): + print( + "Skipping {mat_id}: fast_kappa.hdf5 already exists at " + f"{yaml_file.parent / 'fast_kappa.hdf5'}" + ) + continue + elif os.path.exists(yaml_file.parent / "kappa.hdf5") and os.path.exists( + yaml_file.parent / "fast_kappa.hdf5" + ): + print( + "Skipping {mat_id}: kappa.hdf5 and fast_kappa.hdf5 already exist at " + f"{yaml_file.parent / 'kappa.hdf5'}" + ) + continue + pbar.set_postfix_str(f"{name}-{mat_id}") + + atoms = Atoms( + positions=ph3.unitcell.positions, + cell=ph3.unitcell.cell, + symbols=ph3.unitcell.symbols, + pbc=True, + ) + atoms.info.update(info_dict) + + atoms.write(STRUCTURE_FILE, format="extxyz", append=True) + + ph3.mesh_numbers = fast_qmesh + + with tc.tqdm_gridpoints(desc="Conducitivity calc"): + ph3, fast_kappa_dict, _fast_cond = tc.calculate_conductivity( + ph3, temperatures=TEMPERATURES, log_level=2 + ) + + fast_kappa_dict = { + tc.TCKeys.kappa_tot_avg: tc.calculate_kappa_avg( + fast_kappa_dict[tc.TCKeys.kappa_tot_rta] + ), + tc.TCKeys.mode_kappa_tot_avg: tc.calculate_kappa_avg( + fast_kappa_dict[tc.TCKeys.mode_kappa_tot_rta] + ), + tc.TCKeys.q_points: fast_kappa_dict[tc.TCKeys.q_points], + tc.TCKeys.temperatures: TEMPERATURES, + tc.TCKeys.ph_freqs: _fast_cond.frequencies, + tc.TCKeys.heat_capacity: _fast_cond.mode_heat_capacities, + tc.TCKeys.spg_num: ph3.symmetry.dataset.number, + tc.TCKeys.name: name, + "q_mesh": fast_qmesh, + "q_point_mesh": fast_qmesh, + "weights": fast_kappa_dict[tc.TCKeys.mode_weights], + tc.TCKeys.mode_weights: fast_kappa_dict[tc.TCKeys.mode_weights], + } + + with h5py.File(PH3_PARAMS_PATH / mat_id / "fast_kappa.hdf5", "w") as f: + tc.dict_to_hdf5(fast_kappa_dict, f) + + if not FAST_ONLY: + ph3.mesh_numbers = original_qmesh + + with tc.tqdm_gridpoints(desc="Conducitivity calc"): + ph3, original_kappa_dict, _cond = tc.calculate_conductivity( + ph3, temperatures=TEMPERATURES, log_level=2 + ) + + original_kappa_dict = { + tc.TCKeys.kappa_tot_avg: tc.calculate_kappa_avg( + original_kappa_dict[tc.TCKeys.kappa_tot_rta] + ), + tc.TCKeys.mode_kappa_tot_avg: tc.calculate_kappa_avg( + original_kappa_dict[tc.TCKeys.mode_kappa_tot_rta] + ), + tc.TCKeys.q_points: original_kappa_dict[tc.TCKeys.q_points], + tc.TCKeys.temperatures: TEMPERATURES, + tc.TCKeys.ph_freqs: _cond.frequencies, + tc.TCKeys.heat_capacity: _cond.mode_heat_capacities, + tc.TCKeys.spg_num: ph3.symmetry.dataset.number, + tc.TCKeys.name: name, + "q_mesh": original_qmesh, + "q_point_mesh": original_qmesh, + "weights": original_kappa_dict[tc.TCKeys.mode_weights], + tc.TCKeys.mode_weights: original_kappa_dict[tc.TCKeys.mode_weights], + } + + with h5py.File(PH3_PARAMS_PATH / mat_id / "kappa.hdf5", "w") as f: + tc.dict_to_hdf5(original_kappa_dict, f) diff --git a/ml_peg/calcs/bulk_crystal/thermal_conductivity/thermal_conductivity.py b/ml_peg/calcs/bulk_crystal/thermal_conductivity/thermal_conductivity.py new file mode 100644 index 000000000..9b1008d45 --- /dev/null +++ b/ml_peg/calcs/bulk_crystal/thermal_conductivity/thermal_conductivity.py @@ -0,0 +1,823 @@ +""" +Module for calculating thermal conductivity using Phono3py. + +Code is adapted from https://github.com/MPA2suite/k_SRME/blob/6ff4c867/k_srme/conductivity.py +by Balázs Póta, Paramvir Ahlawat, Gábor Csányi, Michele Simoncelli +and https://github.com/janosh/matbench-discovery/blob/main/matbench_discovery/phonons/thermal_conductivity.py +by Janosh Riebesell. +See https://arxiv.org/abs/2408.00755 for details. +""" + +from __future__ import annotations + +from collections.abc import Sequence +from contextlib import contextmanager +from copy import deepcopy +from enum import Enum +from pathlib import Path +import re +import sys +import traceback +from typing import Any +import warnings + +from ase import Atoms +from ase.calculators.calculator import Calculator +from ase.utils import atoms_to_spglib_cell +import h5py +import numpy as np +import pandas as pd +from phono3py.api_phono3py import Phono3py +from phono3py.conductivity.wigner_rta import ConductivityWignerRTA +from phonopy.structure.atoms import PhonopyAtoms +import spglib +from tqdm import tqdm + + +# Backport-ish "StrEnum" behavior that works in Python 3.11 +class StrEnumCompat(str, Enum): + """ + Enum members behave like their underlying string values. + + - hash(TCKeys.mat_id) == hash("material_id") + - TCKeys.mat_id == "material_id" + """ + + def __str__(self) -> str: + """ + Return the underlying string value. + + Returns + ------- + str + The underlying string value. + """ + return str(self.value) + + def __hash__(self) -> int: + """ + Return the hash of the underlying string value. + + Returns + ------- + int + The hash of the underlying string value. + """ + return hash(self.value) + + def __eq__(self, other) -> bool: + """ + Compare against another enum or string value. + + Parameters + ---------- + other : Any + The value to compare against. + + Returns + ------- + bool + True if the values are equal, otherwise False. + """ + if isinstance(other, Enum): + return self.value == other.value + return self.value == other + + def __format__(self, format_spec: str) -> str: + """ + Format the underlying string value. + + Parameters + ---------- + format_spec : str + The format specification. + + Returns + ------- + str + The formatted string value. + """ + return format(self.value, format_spec) + + +# Use the real StrEnum on 3.11+, otherwise the compat one +if sys.version_info >= (3, 11): + from enum import StrEnum as _BaseStrEnum +else: + _BaseStrEnum = StrEnumCompat + + +class TCKeys(_BaseStrEnum): + """Keys for thermal conductivity dictionary.""" + + kappa_tot_rta = "kappa_tot_rta" + kappa_tot_avg = "kappa_tot_avg" + kappa_p_rta = "kappa_p_rta" + kappa_c = "kappa_c" + mode_weights = "mode_weights" + q_points = "q_points" + ph_freqs = "ph_freqs" + mode_kappa_tot_rta = "mode_kappa_tot_rta" + mode_kappa_tot_avg = "mode_kappa_tot_avg" + has_imag_ph_modes = "has_imag_ph_modes" + temperatures = "temperatures" + mat_id = "material_id" + formula = "formula" + heat_capacity = "heat_capacity" + spg_num = "spg_num" + init_spg_num = "initial_spg_num" + relaxed_space_group_number = "relaxed_space_group_number" + final_spg_num = "final_spg_num" + srd = "srd" + sre = "sre" + srme = "srme" + true_kappa_tot_avg = "true_kappa_tot_avg" + name = "name" + stability = "stability" + + +def calculate_fc2_set( + ph3: Phono3py, calculator: Calculator, pbar_kwargs: dict[str, Any] | None = None +) -> np.ndarray: + """ + Calculate 2nd order force constants. + + Requires initializing Phono3py with an FC2 supercell matrix. + + Parameters + ---------- + ph3 : Phono3py + Phono3py object for which to calculate force constants. + calculator : Calculator + ASE calculator to compute forces. + pbar_kwargs : dict[str, Any] | None + Arguments passed to tqdm progress bar. + Defaults to None. + + Returns + ------- + np.ndarray + Array of forces for each displacement. + """ + # print(f"Computing FC2 force set in {ph3.unitcell.formula}.") + + forces: list[np.ndarray] = [] + n_atoms = len(ph3.phonon_supercell) + + displacements = ph3.phonon_supercells_with_displacements + for supercell in tqdm( + displacements, + desc=f"FC2 calculation: {ph3.unitcell.formula}", + **pbar_kwargs or {}, + ): + if supercell is not None: + atoms = Atoms( + supercell.symbols, + cell=supercell.cell, + positions=supercell.positions, + pbc=True, + ) + atoms.calc = calculator + force = atoms.get_forces() + else: + force = np.zeros((n_atoms, 3)) + forces += [force] + + force_set = np.array(forces) + ph3.phonon_forces = force_set + return force_set + + +def calculate_fc3_set( + ph3: Phono3py, + calculator: Calculator, + pbar_kwargs: dict[str, Any] | None = None, +) -> np.ndarray: + """ + Calculate 3rd order force constants. + + Parameters + ---------- + ph3 : Phono3py + Phono3py object for which to calculate force constants. + calculator : Calculator + ASE calculator to compute forces. + pbar_kwargs : dict[str, Any] | None + Passed to tqdm progress bar. + Defaults to None. + + Returns + ------- + np.ndarray + Array of forces for each displacement. + """ + forces: list[np.ndarray] = [] + n_atoms = len(ph3.supercell) + + desc = f"FC3 calculation: {ph3.unitcell.formula}" + task_idx = (pbar_kwargs or {}).get("position") + if task_idx: + desc = f"{task_idx}. {desc}" + displacements = ph3.supercells_with_displacements + for supercell in tqdm(displacements, desc=desc, **pbar_kwargs or {}): + if supercell is None: + forces += [np.zeros((n_atoms, 3))] + else: + atoms = Atoms( + supercell.symbols, + cell=supercell.cell, + positions=supercell.positions, + pbc=True, + ) + atoms.calc = calculator + forces += [atoms.get_forces()] + + force_set = np.array(forces) + ph3.forces = force_set + return force_set + + +def init_phono3py( + atoms: Atoms, + *, + fc2_supercell: np.ndarray, + fc3_supercell: np.ndarray, + q_point_mesh: tuple[int, int, int] = (20, 20, 20), + displacement_distance: float = 0.03, + symprec: float = 1e-5, + **kwargs: Any, +) -> Phono3py: + """ + Initialize Phono3py object from ASE Atoms. + + Parameters + ---------- + atoms : Atoms + ASE Atoms object to initialize from. + fc2_supercell : np.ndarray + Supercell matrix for 2nd order force constants. + fc3_supercell : np.ndarray + Supercell matrix for 3rd order force constants. + q_point_mesh : tuple[int, int, int] + Mesh size for q-point sampling. Defaults + to (20, 20, 20). + displacement_distance : float + Displacement distance for force calculations. + Defaults to 0.03. + symprec : float + Symmetry precision for finding space group. Defaults to 1e-5. + **kwargs : Any + Passed to Phono3py constructor. + + Returns + ------- + Phono3py + Initialized Phono3py object. + + Raises + ------ + ValueError + If required metadata is missing from atoms.info. + """ + unit_cell = PhonopyAtoms(atoms.symbols, cell=atoms.cell, positions=atoms.positions) + ph3 = Phono3py( + unitcell=unit_cell, + supercell_matrix=fc3_supercell, + phonon_supercell_matrix=fc2_supercell, + primitive_matrix="auto", + symprec=symprec, + **kwargs, + ) + ph3.mesh_numbers = q_point_mesh + + ph3.generate_displacements(distance=displacement_distance) + + return ph3 + + +def get_fc2_and_freqs( + ph3: Phono3py, calculator: Calculator, pbar_kwargs: dict[str, Any] | None = None +) -> tuple[Phono3py, np.ndarray, np.ndarray]: + """ + Calculate 2nd order force constants and phonon frequencies. + + Parameters + ---------- + ph3 : Phono3py + Phono3py object for which to calculate force constants. + calculator : Calculator + ASE calculator to compute forces. + pbar_kwargs : dict[str, Any] | None + Arguments passed to tqdm progress bar. + Defaults to None. + + Returns + ------- + tuple[Phono3py, np.ndarray, np.ndarray] + Tuple of (Phono3py object, force + constants array, frequencies array [shape: (n_bz_grid, n_bands)]). + + Raises + ------ + ValueError + If mesh_numbers not set. + """ + if ph3.mesh_numbers is None: + raise ValueError( + "mesh_numbers was not found in phono3py object and was not provided as " + "an argument when calculating phonons from phono3py object." + ) + + pbar_kwargs = {"leave": False} | (pbar_kwargs or {}) + fc2_set = calculate_fc2_set(ph3, calculator, pbar_kwargs=pbar_kwargs) + + ph3.produce_fc2(symmetrize_fc2=True) + ph3.init_phph_interaction(symmetrize_fc3q=False) + ph3.run_phonon_solver() + + freqs, _eigvecs, _grid = ph3.get_phonon_data() + + return ph3, fc2_set, freqs + + +def load_force_sets( + ph3: Phono3py, fc2_set: np.ndarray, fc3_set: np.ndarray +) -> Phono3py: + """ + Load pre-computed force sets into Phono3py object. + + Parameters + ---------- + ph3 : Phono3py + Phono3py object to load force sets into. + fc2_set : np.ndarray + 2nd order force constants array. + fc3_set : np.ndarray + 3rd order force constants array. + + Returns + ------- + Phono3py + Phono3py object with loaded force sets. + """ + ph3.phonon_forces = fc2_set + ph3.forces = fc3_set + ph3.produce_fc2(symmetrize_fc2=True) + ph3.produce_fc3(symmetrize_fc3r=True) + + return ph3 + + +def calculate_conductivity( + ph3: Phono3py, + temperatures: Sequence[float], + boundary_mfp: float = 1e6, + mode_kappa_thresh: float = 1e-6, + **kwargs: Any, +) -> tuple[Phono3py, dict[str, np.ndarray], ConductivityWignerRTA]: + """ + Calculate thermal conductivity. + + Parameters + ---------- + ph3 : Phono3py + Phono3py object for which to calculate conductivity. + temperatures : list[float] + Temperatures to compute conductivity at in Kelvin. + boundary_mfp : float + Mean free path in micrometer to calculate simple boundary + scattering contribution to thermal conductivity. Defaults to 1e6. + mode_kappa_thresh : float + Threshold for mode kappa consistency check. Defaults + to 1e-6. + **kwargs : Any + Passed to Phono3py.run_thermal_conductivity(). + + Returns + ------- + tuple[Phono3py, dict[str, np.ndarray], ConductivityWignerRTA] + (Phono3py object, + conductivity dict, conductivity object). + """ + ph3.init_phph_interaction(symmetrize_fc3q=False) + + ph3.run_thermal_conductivity( + **kwargs, + temperatures=temperatures, + is_isotope=True, + # use type="wigner" to include both wave-like coherence (kappa_c) and + # particle-like (kappa_p) conductivity contributions + conductivity_type="wigner", + boundary_mfp=boundary_mfp, + ) + + kappa = ph3.thermal_conductivity + + kappa_dict = { + TCKeys.kappa_tot_rta: deepcopy(kappa.kappa_TOT_RTA[0]), + TCKeys.kappa_p_rta: deepcopy(kappa.kappa_P_RTA[0]), + TCKeys.kappa_c: deepcopy(kappa.kappa_C[0]), + TCKeys.mode_weights: deepcopy(kappa.grid_weights), + TCKeys.q_points: deepcopy(kappa.qpoints), + TCKeys.ph_freqs: deepcopy(kappa.frequencies), + } + mode_kappa_total = kappa_dict[TCKeys.mode_kappa_tot_rta] = calc_mode_kappa_tot( + deepcopy(kappa.mode_kappa_P_RTA[0]), + deepcopy(kappa.mode_kappa_C[0]), + deepcopy(kappa.mode_heat_capacities), + ) + + sum_mode_kappa_tot = mode_kappa_total.sum( + axis=tuple(range(1, mode_kappa_total.ndim - 1)) + ) / np.sum(kappa_dict[TCKeys.mode_weights]) + + kappa_tot_rta = kappa_dict[TCKeys.kappa_tot_rta] + if np.any(np.abs(sum_mode_kappa_tot - kappa_tot_rta) > mode_kappa_thresh): + warnings.warn( + f"Total mode kappa does not sum to total kappa. {sum_mode_kappa_tot=}, " + f"{kappa_tot_rta=}", + stacklevel=2, + ) + + return ph3, kappa_dict, kappa + + +def calc_mode_kappa_tot( + mode_kappa_p_rta: np.ndarray, + mode_kappa_coherence: np.ndarray, + heat_capacity: np.ndarray, +) -> np.ndarray: + """ + Calculate total mode kappa from particle-like RTA and coherence terms. + + Parameters + ---------- + mode_kappa_p_rta : np.ndarray + Mode kappa from particle-like RTA with shape + (T, q-points, bands, xyz). + mode_kappa_coherence : np.ndarray + Mode kappa from wave-like coherence with + shape (T, q-points, bands, bands, xyz). + heat_capacity : np.ndarray + Mode heat capacities with shape + (T, q-points, bands). + + Returns + ------- + np.ndarray + Total (particle-like + wave-like) thermal conductivity per phonon + mode with shape (T, q-points, bands). + """ + # Temporarily silence divide warnings since we handle NaN values below + with np.errstate(divide="ignore", invalid="ignore"): + mode_kappa_c_per_mode = 2 * ( # None equiv to np.newaxis + (mode_kappa_coherence * heat_capacity[:, :, :, None, None]) + / (heat_capacity[:, :, :, None, None] + heat_capacity[:, :, None, :, None]) + ).sum(axis=2) + + mode_kappa_c_per_mode[np.isnan(mode_kappa_c_per_mode)] = 0 + + return mode_kappa_c_per_mode + mode_kappa_p_rta + + +def check_imaginary_freqs(frequencies: np.ndarray, threshold: float = -0.01) -> bool: + """ + Check if frequencies are imaginary. + + Parameters + ---------- + frequencies : np.ndarray + Frequencies to check. + threshold : float + Threshold for imaginary frequencies. Defaults to -0.01. + + Returns + ------- + bool + True if imaginary frequencies are found. + """ + # Return True if all frequencies are NaN, indicating invalid or missing data + if np.all(pd.isna(frequencies)): + return True + + # Check for imaginary frequencies in non-acoustic modes at gamma point (q=0) + # Indices 3+ correspond to optical modes which should never be negative + if np.any(frequencies[0, 3:] < 0): + return True + + # Check acoustic modes at gamma point against threshold. First 3 modes at q=0 + # are acoustic and may be slightly negative due to numerical noise + if np.any(frequencies[0, :3] < threshold): + return True + + # Check for imaginary frequencies at any q-point except gamma + # All frequencies should be positive away from gamma point + return bool(np.any(frequencies[1:] < 0)) + + +def get_spacegroup_number_from_atoms(atoms: Atoms, symprec: float = 1e-5) -> int: + """ + Get space group number from ASE Atoms object. + + Parameters + ---------- + atoms : Atoms + ASE Atoms object to analyze. + symprec : float + Symmetry precision. Defaults to 1e-5. + + Returns + ------- + int + Space group number. + """ + dataset = spglib.get_symmetry_dataset(atoms_to_spglib_cell(atoms), symprec=symprec) + return dataset.number + + +def calculate_kappa_avg(kappa: np.ndarray) -> np.ndarray: + """ + Calculate directionally averaged trace of the conductivity tensor. + + Takes a thermal conductivity tensor and returns its trace (average of diagonal + components). This represents the average thermal conductivity in the 3 spatial + directions, which is a useful scalar metric for comparing materials. + + Parameters + ---------- + kappa : np.ndarray + Thermal conductivity tensor, typically of shape (..., 3, 3) where + the last two dimensions represent the 3x3 conductivity tensor. + Earlier dimensions may include temperatures or other parameters. + + Returns + ------- + np.ndarray + Average conductivity value(s). Returns np.nan if the input contains + any NaN values or if the calculation fails. For multiple temperatures, + returns an array of averages. + """ + if np.any(pd.isna(kappa)): + return np.array([np.nan]) + try: + return np.asarray(kappa)[..., :3].mean(axis=-1) + except Exception: + warnings.warn( + f"Failed to calculate kappa_avg: {traceback.format_exc()}", stacklevel=2 + ) + return np.array([np.nan]) + + +_GRID_RE = re.compile(r"Grid point\s+\d+\s*\(\s*(\d+)\s*/\s*(\d+)\s*\)") + + +class _TqdmIntercept: + """ + Internal class to intercept stdout and update a tqdm progress bar. + + Based on output lines. Not intended for external use. + + Parameters + ---------- + pbar : tqdm.tqdm + The tqdm progress bar instance to update based on intercepted output. + """ + + def __init__(self, pbar): + """ + Initialize the interceptor with a tqdm progress bar. + + Parameters + ---------- + pbar : tqdm.tqdm + The tqdm progress bar instance to update based on intercepted output. + """ + self.pbar = pbar + self._buf = "" + + def write(self, s: str) -> int: + """ + Intercept writes to stdout, parse for progress, and update the tqdm bar. + + This method looks for lines in the output matching the pattern + "Grid point X/Y" and updates the tqdm progress bar accordingly. + It handles partial lines and ensures that the bar is updated smoothly + without breaking the rendering. + + Parameters + ---------- + s : str + String to write, typically a line of output from the conductivity + calculation. + + Returns + ------- + int + The number of characters written (length of the input string). + """ + if not s: + return 0 + self._buf += s + + while True: + nl = self._buf.find("\n") + if nl < 0: + break + line = self._buf[:nl] + self._buf = self._buf[nl + 1 :] + + m = _GRID_RE.search(line) + if not m: + continue + + cur = int(m.group(1)) + total = int(m.group(2)) + + # Learn total from the output + if self.pbar.total is None or self.pbar.total != total: + self.pbar.total = total + + # Start bar at 0 when cur==1 (or when the first seen cur is 2, show 1, etc.) + if cur >= total: + new_n = total + else: + new_n = max(0, min(total, cur - 1)) + if new_n != self.pbar.n: + self.pbar.n = new_n + self.pbar.refresh() + + return len(s) + + def flush(self): + """ + Flush the buffer and update the progress bar to completion if needed. + + This should be called when the progress is complete to ensure the bar finished. + """ + if self._buf: + self.write("\n") + if self.pbar.n < self.pbar.total: + self.pbar.n = self.pbar.total + self.pbar.refresh() + + # def isatty(self): + # return False + + @property + def encoding(self) -> str: + """ + Return the encoding of the underlying real stdout. + + Or 'utf-8' if it cannot be determined. + + Returns + ------- + str + Encoding of the underlying real stdout, or 'utf-8' if it cannot be + determined. + """ + return getattr(sys.__stdout__, "encoding", "utf-8") + + +@contextmanager +def tqdm_gridpoints(desc="Grid points", intercept_stderr=False, leave=False): + """ + Context manager that creates a ``tqdm`` progress bar for grid-point iteration. + + Context manager that creates a ``tqdm`` progress bar for grid-point iteration + and temporarily redirects standard output (and optionally standard error) so + printed messages are routed through the bar without breaking its rendering. + + Parameters + ---------- + desc : str, optional + Description shown next to the progress bar. Default is ``"Grid points"``. + intercept_stderr : bool, optional + If ``True``, temporarily redirect ``sys.stderr`` to the same interceptor as + ``sys.stdout`` while inside the context. Default is ``False``. + leave : bool, optional + Whether to leave the progress bar displayed after completion, passed to + ``tqdm``. Default is ``False``. + + Yields + ------ + tqdm.tqdm + The active progress bar instance, allowing manual updates and metadata + changes within the context. + """ + old_out, old_err = sys.stdout, sys.stderr + real_out = sys.__stdout__ if sys.__stdout__ is not None else old_out + + with tqdm( + total=None, desc=desc, file=real_out, dynamic_ncols=True, leave=leave + ) as pbar: + interceptor = _TqdmIntercept(pbar) + try: + sys.stdout = interceptor + if intercept_stderr: + sys.stderr = interceptor + yield pbar + finally: + interceptor.flush() + sys.stdout = old_out + sys.stderr = old_err + + +def dict_to_hdf5( + data: dict[str, Any], h5file: h5py.File, group_name: str | None = None +) -> None: + """ + Save a dictionary to an HDF5 file. + + Parameters + ---------- + data : dict[str, Any] + Dictionary to save. + h5file : h5py.File + Open HDF5 file object. + group_name : str + Optional name of the group to save the data under. Defaults to None. + """ + if group_name is None: + group = h5file + else: + group = h5file.require_group(group_name) + + for key, value in data.items(): + if isinstance(value, dict): + dict_to_hdf5(value, group, group_name=f"{key}") + else: + if value is None: + value = np.array( + [np.nan] + ) # h5py does not support None, use NaN instead + try: + group.create_dataset(key, data=value) + except TypeError as e: + print( + f"Error: Could not save key '{key}' to HDF5 file. Unsupported " + f"type: {type(value)}, with value: {value}" + ) + # traceback.print_exc() + raise e + + +def hdf5_to_dict(h5file: h5py.File, group_name: str | None = None) -> dict[str, Any]: + """ + Load a dictionary from an HDF5 file. + + Parameters + ---------- + h5file : h5py.File + Open HDF5 file object. + group_name : str + Optional name of the group to load the data from. Defaults to None. + + Returns + ------- + dict[str, Any] + Loaded dictionary. + """ + if group_name is None: + group = h5file + else: + group = h5file[group_name] + + data: dict[str, Any] = {} + for key, item in group.items(): + if isinstance(item, h5py.Group): + data[key] = hdf5_to_dict(h5file, group_name=f"{group.name}/{key}") + else: + data[key] = item[()] + + return data + + +def load_hdf5_subdir_dicts( + path: Path | str, filename: str +) -> dict[str, dict[str, Any]]: + """ + Load dictionaries from HDF5 files in subdirectories. + + Parameters + ---------- + path : str | Path + Path to the directory containing subdirectories with HDF5 files. + filename (str): Name of the HDF5 file to load from each subdirectory. + filename : str + Name of the HDF5 file to load from each subdirectory. + + Returns + ------- + dict[str, dict[str, Any]] + Dictionary with subdirectory names as keys and loaded dictionaries as values. + """ + subdirs = [d for d in Path(path).iterdir() if d.is_dir()] + dicts = {} + for d in subdirs: + try: + with h5py.File(d / filename, "r") as f: + dicts[d.name] = hdf5_to_dict(f) + except FileNotFoundError: + raise FileNotFoundError(f"File not found: {d / filename}") from None + return dicts From 88c58f349c25700388a2f53b229d6fd625cda6f4 Mon Sep 17 00:00:00 2001 From: bpota Date: Thu, 9 Jul 2026 19:51:42 +0100 Subject: [PATCH 2/4] Add thermal conductivity app component and fix analysis MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Add app_thermal_conductivity.py Dash sub-app - Replace unicode κ with k in docstrings/labels for encoding safety - Sort DataFrames by index before metric computation for determinism - Fix output filename: figure_conductivity -> figure_thermal_conductivity - Add OC157 dvc.yaml, update uv.lock, add CLAUDE.md Co-Authored-By: Claude Sonnet 4.6 Claude-Session: https://claude.ai/code/session_01S7VGFQPPbiFQsxbHa1yKbu --- .../analyse_thermal_conductivity.py | 32 ++++---- .../app_thermal_conductivity.py | 76 +++++++++++++++++++ .../calc_thermal_conductivity.py | 6 +- .../thermal_conductivity.py | 6 +- 4 files changed, 99 insertions(+), 21 deletions(-) create mode 100644 ml_peg/app/bulk_crystal/thermal_conductivity/app_thermal_conductivity.py diff --git a/ml_peg/analysis/bulk_crystal/thermal_conductivity/analyse_thermal_conductivity.py b/ml_peg/analysis/bulk_crystal/thermal_conductivity/analyse_thermal_conductivity.py index 2960852d6..8632c175d 100644 --- a/ml_peg/analysis/bulk_crystal/thermal_conductivity/analyse_thermal_conductivity.py +++ b/ml_peg/analysis/bulk_crystal/thermal_conductivity/analyse_thermal_conductivity.py @@ -156,7 +156,7 @@ def calc_kappa_srme_dataframes( ------- list[float]: SRME values for each material. Values are between 0 and 2, where: - - 0 indicates perfect agreement in both total κ and mode-resolved properties + - 0 indicates perfect agreement in both total k and mode-resolved properties - 2 indicates complete failure (imaginary frequencies, broken symmetry, etc.) - Values in between indicate partial agreement, with lower being better """ @@ -187,12 +187,12 @@ def calc_kappa_srme(kappas_pred: pd.Series, kappas_true: pd.Series) -> np.ndarra """ Calculate the Symmetric Relative Mean Error (SRME) for a single material. - SRME = 2 * (sum|κ_pred,i - κ_true,i| * w_i) / (κ_pred,tot + κ_true,tot) + SRME = 2 * (sum|k_pred,i - k_true,i| * w_i) / (k_pred,tot + k_true,tot) where: - - κ_pred,i and κ_true,i are mode-resolved conductivities for mode i + - k_pred,i and k_true,i are mode-resolved conductivities for mode i - w_i are the mode weights - - κ_pred,tot and κ_true,tot are total conductivities + - k_pred,tot and k_true,tot are total conductivities The calculation involves: 1. Computing mode-resolved average conductivities if not pre-computed @@ -213,7 +213,7 @@ def calc_kappa_srme(kappas_pred: pd.Series, kappas_true: pd.Series) -> np.ndarra ------- np.ndarray SRME values per temperature, each between 0 and 2, where: - - 0 indicates perfect agreement in both total κ and mode-resolved properties + - 0 indicates perfect agreement in both total k and mode-resolved properties - 2 indicates complete disagreement or invalid results On error conditions (missing data, NaN values), returns np.array([2.0]). """ @@ -288,11 +288,13 @@ def kappa_stats() -> dict[str, pd.DataFrame]: tc.hdf5_to_dict(h5py.File(REF_PATH / "PBE" / "kappas.hdf5", "r")), orient="index", ) + ref_df.sort_index(inplace=True) fast_ref_df = pd.DataFrame.from_dict( tc.hdf5_to_dict(h5py.File(REF_PATH / "PBE" / "fast_kappas.hdf5", "r")), orient="index", ) + fast_ref_df.sort_index(inplace=True) results["ref"] = ref_df @@ -348,13 +350,17 @@ def kappa_stats() -> dict[str, pd.DataFrame]: continue if pred_df is not None and fast_pred_df is None: + pred_df.sort_index(inplace=True) pred_df = calc_kappa_metrics_from_dfs(pred_df, ref_df) elif pred_df is None: + fast_pred_df.sort_index(inplace=True) pred_df = calc_kappa_metrics_from_dfs(fast_pred_df, fast_ref_df) pred_df["fast_sre"] = pred_df[tc.TCKeys.sre] pred_df["fast_srme"] = pred_df[tc.TCKeys.srme] pred_df.drop(columns=[tc.TCKeys.sre, tc.TCKeys.srme], inplace=True) elif fast_pred_df is not None: + pred_df.sort_index(inplace=True) + fast_pred_df.sort_index(inplace=True) pred_df = calc_kappa_metrics_from_dfs(pred_df, ref_df) fast_pred_df = calc_kappa_metrics_from_dfs(fast_pred_df, fast_ref_df) pred_df["fast_sre"] = fast_pred_df[tc.TCKeys.sre] @@ -374,10 +380,10 @@ def kappa_stats() -> dict[str, pd.DataFrame]: @pytest.fixture @plot_parity( - filename=OUT_PATH / "figure_conductivity.json", + filename=OUT_PATH / "figure_thermal_conductivity.json", title="Thermal Conductivity Parity Plot", - x_label="Predicted κ (W/mK)", - y_label="Reference κ (W/mK)", + x_label="Predicted k (W/mK)", + y_label="Reference k (W/mK)", hoverdata={ "System": get_system_names(), "System ID": get_system_ids(), @@ -591,7 +597,7 @@ def mean_fast_srme(kappa_stats: dict[str, pd.DataFrame]) -> dict[str, float]: @pytest.fixture @build_table( - filename=OUT_PATH / "kappas.json", + filename=OUT_PATH / "thermal_conductivity.json", metric_tooltips=DEFAULT_TOOLTIPS, thresholds=DEFAULT_THRESHOLDS, ) @@ -627,12 +633,12 @@ def metrics( Metric names and values for all models. """ return { - "κSRE": mean_sre, - "κSRME": mean_srme, + "kSRE": mean_sre, + "kSRME": mean_srme, "Instability": instability, "Failure": failure, - "Fast κSRE": mean_fast_sre, - "Fast κSRME": mean_fast_srme, + "Fast kSRE": mean_fast_sre, + "Fast kSRME": mean_fast_srme, } diff --git a/ml_peg/app/bulk_crystal/thermal_conductivity/app_thermal_conductivity.py b/ml_peg/app/bulk_crystal/thermal_conductivity/app_thermal_conductivity.py new file mode 100644 index 000000000..1e1bafc8f --- /dev/null +++ b/ml_peg/app/bulk_crystal/thermal_conductivity/app_thermal_conductivity.py @@ -0,0 +1,76 @@ +"""Run Thermal thermal_Conductivity app.""" + +from __future__ import annotations + +from dash.html import Div + +from ml_peg.app import APP_ROOT +from ml_peg.app.base_app import BaseApp +from ml_peg.app.utils.build_callbacks import ( + plot_from_table_column, +) +from ml_peg.app.utils.load import read_plot +from ml_peg.models.get_models import get_model_names +from ml_peg.models.models import current_models + +# Get all models +MODELS = get_model_names(current_models) +BENCHMARK_NAME = "Thermal Conductivity" +DOCS_URL = ( + "https://ddmms.github.io/ml-peg/user_guide/" + "benchmarks/bulk_crystal.html#thermal_conductivity" +) +DATA_PATH = APP_ROOT / "data" / "bulk_crystal" / "thermal_conductivity" + + +class ThermalConductivityApp(BaseApp): + """Thermal Conductivity benchmark app layout and callbacks.""" + + def register_callbacks(self) -> None: + """Register callbacks to app.""" + scatter = read_plot( + DATA_PATH / "figure_thermal_conductivity.json", + id=f"{BENCHMARK_NAME}-figure", + ) + + # Assets dir will be parent directory - individual files for each system + # structs_dir = DATA_PATH / MODELS[0] + # structs = [ + # f"assets/bulk_crystal/thermal_conductivity/" + # f"{MODELS[0]}/{struct_file.stem}.xyz" + # for struct_file in sorted(structs_dir.glob("*.xyz")) + # ] + + plot_from_table_column( + table_id=self.table_id, + plot_id=f"{BENCHMARK_NAME}-figure-placeholder", + column_to_plot={"MAE": scatter}, + ) + + # struct_from_scatter( + # scatter_id=f"{BENCHMARK_NAME}-figure", + # struct_id=f"{BENCHMARK_NAME}-struct-placeholder", + # structs=structs, + # mode="struct", + # ) + + +def get_app() -> ThermalConductivityApp: + """ + Get Thermal Conductivity benchmark app layout and callback registration. + + Returns + ------- + ThermalConductivityApp + Benchmark layout and callback registration. + """ + return ThermalConductivityApp( + name=BENCHMARK_NAME, + description="Thermal conductivity for 103 binary crystals.", + docs_url=DOCS_URL, + table_path=DATA_PATH / "thermal_conductivity.json", + extra_components=[ + Div(id=f"{BENCHMARK_NAME}-figure-placeholder"), + # Div(id=f"{BENCHMARK_NAME}-struct-placeholder"), + ], + ) diff --git a/ml_peg/calcs/bulk_crystal/thermal_conductivity/calc_thermal_conductivity.py b/ml_peg/calcs/bulk_crystal/thermal_conductivity/calc_thermal_conductivity.py index 80d15243a..239a3051e 100644 --- a/ml_peg/calcs/bulk_crystal/thermal_conductivity/calc_thermal_conductivity.py +++ b/ml_peg/calcs/bulk_crystal/thermal_conductivity/calc_thermal_conductivity.py @@ -267,9 +267,9 @@ def test_thermal_conductivity(mlip: tuple[str, Any]) -> None: ) if ltc_condition: - # fc3_set = tc.calculate_fc3_set( - # ph3, calculator=calculator, pbar_kwargs={"leave": False} - # ) + tc.calculate_fc3_set( + ph3, calculator=calculator, pbar_kwargs={"leave": False} + ) ph3.produce_fc3(symmetrize_fc3r=True) else: reason = [] diff --git a/ml_peg/calcs/bulk_crystal/thermal_conductivity/thermal_conductivity.py b/ml_peg/calcs/bulk_crystal/thermal_conductivity/thermal_conductivity.py index 9b1008d45..1e124d36e 100644 --- a/ml_peg/calcs/bulk_crystal/thermal_conductivity/thermal_conductivity.py +++ b/ml_peg/calcs/bulk_crystal/thermal_conductivity/thermal_conductivity.py @@ -654,11 +654,7 @@ def flush(self): This should be called when the progress is complete to ensure the bar finished. """ - if self._buf: - self.write("\n") - if self.pbar.n < self.pbar.total: - self.pbar.n = self.pbar.total - self.pbar.refresh() + return # def isatty(self): # return False From 9fa6a9e38732d851e06308f169888ed92b33b7be Mon Sep 17 00:00:00 2001 From: bpota Date: Sun, 12 Jul 2026 22:27:44 +0100 Subject: [PATCH 3/4] Fix thermal conductivity analysis and app --- .../analyse_thermal_conductivity.py | 100 ++-- .../thermal_conductivity/metrics.yml | 44 ++ .../app_thermal_conductivity.py | 11 +- .../calc_thermal_conductivity.py | 430 +++++++++--------- .../thermal_conductivity.py | 2 +- 5 files changed, 352 insertions(+), 235 deletions(-) create mode 100644 ml_peg/analysis/bulk_crystal/thermal_conductivity/metrics.yml diff --git a/ml_peg/analysis/bulk_crystal/thermal_conductivity/analyse_thermal_conductivity.py b/ml_peg/analysis/bulk_crystal/thermal_conductivity/analyse_thermal_conductivity.py index 8632c175d..e6ad9e2b9 100644 --- a/ml_peg/analysis/bulk_crystal/thermal_conductivity/analyse_thermal_conductivity.py +++ b/ml_peg/analysis/bulk_crystal/thermal_conductivity/analyse_thermal_conductivity.py @@ -2,9 +2,9 @@ from __future__ import annotations -import os from pathlib import Path import traceback +import warnings from ase.io import read import h5py @@ -25,7 +25,7 @@ REF_PATH = CALCS_ROOT / "bulk_crystal" / "thermal_conductivity" / "data" OUT_PATH = APP_ROOT / "data" / "bulk_crystal" / "thermal_conductivity" -METRICS_CONFIG_PATH = Path(__file__).with_name("metrics.yaml") +METRICS_CONFIG_PATH = Path(__file__).with_name("metrics.yml") DEFAULT_THRESHOLDS, DEFAULT_TOOLTIPS, DEFAULT_WEIGHTS = load_metrics_config( METRICS_CONFIG_PATH ) @@ -166,7 +166,8 @@ def calc_kappa_srme_dataframes( # NOTE code below just until before return used to be wrapped in try/except in # which case SRME=2 was set for the failing material - if row_pred.get(tc.TCKeys.has_imag_ph_modes) is True: + has_imag_ph_modes = row_pred.get(tc.TCKeys.has_imag_ph_modes, False) + if pd.notna(has_imag_ph_modes) and bool(has_imag_ph_modes): srme_list.append(2) continue if relaxed_space_group_number := row_pred.get(tc.TCKeys.final_spg_num): @@ -270,6 +271,42 @@ def calc_kappa_srme(kappas_pred: pd.Series, kappas_true: pd.Series) -> np.ndarra return 2 * microscopic_error / denominator +def _add_missing_error_rows( + df: pd.DataFrame, + ref_df: pd.DataFrame, + model_name: str, + result_type: str = "", +) -> pd.DataFrame: + """ + Warn and add default-error rows for missing result files. + + Parameters + ---------- + df : pd.DataFrame + Prediction results. + ref_df : pd.DataFrame + Reference results defining the required index. + model_name : str + Model name used in the warning. + result_type : str + Optional result-type prefix used in the warning. + + Returns + ------- + pd.DataFrame + Prediction results reindexed to the reference when rows are missing. + """ + missing = ref_df.index.difference(df.index) + if len(missing): + warnings.warn( + f"Missing {len(missing)} {result_type}thermal-conductivity result files " + f"for {model_name}; using default error rows.", + stacklevel=2, + ) + return df.reindex(ref_df.index) + return df + + @pytest.fixture def kappa_stats() -> dict[str, pd.DataFrame]: """ @@ -312,18 +349,17 @@ def kappa_stats() -> dict[str, pd.DataFrame]: ) pred_df = pd.DataFrame.from_dict(current_pred_dict, orient="index") elif (model_dir / "kappas.json.gz").exists(): - pred_df = pd.read_json(model_dir / "kappas.json.gz", orient="index") + pred_df = pd.read_json(model_dir / "kappas.json.gz").set_index( + tc.TCKeys.mat_id + ) else: - subdirs = sorted(os.listdir(model_dir)) - - if not any(not (Path(d) / "kappa.hdf5").exists() for d in subdirs): - try: - current_pred_dict = tc.load_hdf5_subdir_dicts( - model_dir, "kappa.hdf5" - ) - except Exception: - print(f"Error loading kappas for {model_name}...") - traceback.print_exc() + try: + current_pred_dict = tc.load_hdf5_subdir_dicts(model_dir, "kappa.hdf5") + if current_pred_dict: + pred_df = pd.DataFrame.from_dict(current_pred_dict, orient="index") + except Exception: + print(f"Error loading kappas for {model_name}...") + traceback.print_exc() if (model_dir / "fast_kappas.hdf5").exists(): current_pred_dict = tc.hdf5_to_dict( @@ -331,20 +367,28 @@ def kappa_stats() -> dict[str, pd.DataFrame]: ) fast_pred_df = pd.DataFrame.from_dict(current_pred_dict, orient="index") elif (model_dir / "fast_kappas.json.gz").exists(): - fast_pred_df = pd.read_json( - model_dir / "fast_kappas.json.gz", orient="index" + fast_pred_df = pd.read_json(model_dir / "fast_kappas.json.gz").set_index( + tc.TCKeys.mat_id ) else: - subdirs = sorted(os.listdir(model_dir)) - - if not any(not (Path(d) / "fast_kappa.hdf5").exists() for d in subdirs): - try: - current_pred_dict = tc.load_hdf5_subdir_dicts( - model_dir, "fast_kappa.hdf5" + try: + current_pred_dict = tc.load_hdf5_subdir_dicts( + model_dir, "fast_kappa.hdf5" + ) + if current_pred_dict: + fast_pred_df = pd.DataFrame.from_dict( + current_pred_dict, orient="index" ) - except Exception: - print(f"Error loading fast_kappas for {model_name}...") - traceback.print_exc() + except Exception: + print(f"Error loading fast_kappas for {model_name}...") + traceback.print_exc() + + if pred_df is not None: + pred_df = _add_missing_error_rows(pred_df, ref_df, model_name) + if fast_pred_df is not None: + fast_pred_df = _add_missing_error_rows( + fast_pred_df, fast_ref_df, model_name, "fast " + ) if pred_df is None and fast_pred_df is None: continue @@ -411,6 +455,7 @@ def conductivity(kappa_stats: dict[str, pd.DataFrame]) -> dict[str, list]: df = kappa_stats[model_name] if tc.TCKeys.kappa_tot_avg not in df: continue + df = df.reindex(kappa_stats["ref"].index) conductivity_data[model_name] = df[tc.TCKeys.kappa_tot_avg].tolist() conductivity_data["ref"] = kappa_stats["ref"][tc.TCKeys.kappa_tot_avg].tolist() return conductivity_data @@ -496,9 +541,7 @@ def instability(kappa_stats: dict[str, pd.DataFrame]) -> dict[str, float]: if tc.TCKeys.has_imag_ph_modes not in df: continue has_imag_ph_modes = df[tc.TCKeys.has_imag_ph_modes].values - has_imag_ph_modes_filtered = has_imag_ph_modes[ - np.logical_not(np.isnan(has_imag_ph_modes)) - ] + has_imag_ph_modes_filtered = has_imag_ph_modes[pd.notna(has_imag_ph_modes)] instability_values[model_name] = has_imag_ph_modes_filtered.sum() / ref_length return instability_values @@ -600,6 +643,7 @@ def mean_fast_srme(kappa_stats: dict[str, pd.DataFrame]) -> dict[str, float]: filename=OUT_PATH / "thermal_conductivity.json", metric_tooltips=DEFAULT_TOOLTIPS, thresholds=DEFAULT_THRESHOLDS, + weights=DEFAULT_WEIGHTS, ) def metrics( mean_srme: dict[str, float], diff --git a/ml_peg/analysis/bulk_crystal/thermal_conductivity/metrics.yml b/ml_peg/analysis/bulk_crystal/thermal_conductivity/metrics.yml new file mode 100644 index 000000000..98aec599d --- /dev/null +++ b/ml_peg/analysis/bulk_crystal/thermal_conductivity/metrics.yml @@ -0,0 +1,44 @@ +metrics: + kSRE: + good: 0.0 + bad: 1.0 + unit: null + tooltip: "mean of Symmetric Relative Error in thermal conductivity" + level_of_theory: PBE + + kSRME: + good: 0.0 + bad: 1.0 + unit: null + tooltip: "mean of Symmetric Relative Mean Error in mode contributions to thermal conductivity" + level_of_theory: PBE + + Instability: + good: 0.0 + bad: 1.0 + unit: null + tooltip: "Fraction of unstable materials (has imaginary frequencies)" + level_of_theory: PBE + + Failure: + good: 0.0 + bad: 1.0 + unit: null + tooltip: "Fraction of materials for which the calculation failed (NaN frequencies)" + level_of_theory: PBE + + Fast kSRE: + good: 0.0 + bad: 1.0 + unit: null + weight: 0.0 + tooltip: "mean of Symmetric Relative Error in thermal conductivity, with reduced q point mesh" + level_of_theory: PBE + + Fast kSRME: + good: 0.0 + bad: 1.0 + unit: null + weight: 0.0 + tooltip: "mean of Symmetric Relative Mean Error in mode contributions to thermal conductivity, with reduced q point mesh" + level_of_theory: PBE diff --git a/ml_peg/app/bulk_crystal/thermal_conductivity/app_thermal_conductivity.py b/ml_peg/app/bulk_crystal/thermal_conductivity/app_thermal_conductivity.py index 1e1bafc8f..6dfbd64c9 100644 --- a/ml_peg/app/bulk_crystal/thermal_conductivity/app_thermal_conductivity.py +++ b/ml_peg/app/bulk_crystal/thermal_conductivity/app_thermal_conductivity.py @@ -32,6 +32,10 @@ def register_callbacks(self) -> None: DATA_PATH / "figure_thermal_conductivity.json", id=f"{BENCHMARK_NAME}-figure", ) + fast_scatter = read_plot( + DATA_PATH / "figure_fast_thermal_conductivity.json", + id=f"{BENCHMARK_NAME}-fast-figure", + ) # Assets dir will be parent directory - individual files for each system # structs_dir = DATA_PATH / MODELS[0] @@ -44,7 +48,12 @@ def register_callbacks(self) -> None: plot_from_table_column( table_id=self.table_id, plot_id=f"{BENCHMARK_NAME}-figure-placeholder", - column_to_plot={"MAE": scatter}, + column_to_plot={ + "kSRE": scatter, + "kSRME": scatter, + "Fast kSRE": fast_scatter, + "Fast kSRME": fast_scatter, + }, ) # struct_from_scatter( diff --git a/ml_peg/calcs/bulk_crystal/thermal_conductivity/calc_thermal_conductivity.py b/ml_peg/calcs/bulk_crystal/thermal_conductivity/calc_thermal_conductivity.py index 239a3051e..18069c0bb 100644 --- a/ml_peg/calcs/bulk_crystal/thermal_conductivity/calc_thermal_conductivity.py +++ b/ml_peg/calcs/bulk_crystal/thermal_conductivity/calc_thermal_conductivity.py @@ -117,216 +117,15 @@ def test_thermal_conductivity(mlip: tuple[str, Any]) -> None: fast_kappa_dicts = {} for i, atoms_input in enumerate(tqdm(atoms_list, desc="Atoms loop")): - atoms = atoms_input.copy() - - atoms.calc = calculator - - structure_id = atoms.info.get(tc.TCKeys.mat_id, f"structure_{i}") + structure_id = atoms_input.info.get(tc.TCKeys.mat_id, f"structure_{i}") out_dir = OUT_PATH / model_name / structure_id out_dir.mkdir(parents=True, exist_ok=True) - # Relax structure before calculating thermal conductivity - relax_path = out_dir / "relaxed.extxyz" - - formula = atoms.info.get("name", atoms.get_chemical_formula()) - - mat_id = atoms.info[tc.TCKeys.mat_id] - init_info = deepcopy(atoms.info) - info_dict: dict[str, Any] = { - str(tc.TCKeys.mat_id): mat_id, - str(tc.TCKeys.formula): formula, - } - err_dict: dict[str, list[str]] = {"errors": [], "error_traceback": []} - - # Select filter class - if ase_filter in {"frechet", "exp"}: - filter_cls: type[Filter] = { - "frechet": FrechetCellFilter, - "exp": ExpCellFilter, - }[ase_filter] - else: - # Default to FrechetCellFilter if not specified (for MACE compatibility) - filter_cls = FrechetCellFilter - - # Select optimizer class - optimizer_dict = { - "GPMin": ase.optimize.GPMin, - "GOQN": ase.optimize.GoodOldQuasiNewton, - "BFGSLineSearch": ase.optimize.BFGSLineSearch, - "QuasiNewton": ase.optimize.BFGSLineSearch, - "SciPyFminBFGS": ase.optimize.sciopt.SciPyFminBFGS, - "BFGS": ase.optimize.BFGS, - "LBFGSLineSearch": ase.optimize.LBFGSLineSearch, - "SciPyFminCG": ase.optimize.sciopt.SciPyFminCG, - "FIRE2": ase.optimize.FIRE2, - "FIRE": ase.optimize.FIRE, - "LBFGS": ase.optimize.LBFGS, - } - optim_cls: type[Optimizer] = optimizer_dict[ase_optimizer] - - # Initialize variables that might be needed in error handling - relax_dict: dict[str, Any] = { - "max_stress": None, - "reached_max_steps": False, - "broken_symmetry": False, - } - # initial space group for symmetry breaking detection - init_spg_num = tc.get_spacegroup_number_from_atoms(atoms, symprec=symprec) - fast_results_dict = None - - # Relaxation - try: - results_dict = {} - atoms.calc = calculator - if max_steps > 0: - if enforce_relax_symm: - atoms.set_constraint(FixSymmetry(atoms)) - filtered_atoms = filter_cls(atoms, mask=[True] * 3 + [False] * 3) - else: - filtered_atoms = filter_cls(atoms) - - optimizer = optim_cls(filtered_atoms, logfile=out_dir / "relax.log") - optimizer.run(fmax=fmax, steps=max_steps) - - step_count = getattr( - optimizer, "nsteps", None - ) # Get optimizer step count - if ( - step_count is None - ): # fallback to extract from state_dict if available - state = getattr(optimizer, "state_dict", dict)() - step_count = state.get("step", 0) - - reached_max_steps = step_count >= max_steps - if reached_max_steps: - print(f"Material {mat_id=} reached {max_steps=} during relaxation") - - # maximum residual stress component in for xx,yy,zz and xy,yz,xz - # components separately result is a array of 2 elements - max_stress = atoms.get_stress().reshape((2, 3), order="C").max(axis=1) - - atoms.calc = None - atoms.constraints = None - atoms.info = init_info | atoms.info - - # Check if symmetry was broken during relaxation - relaxed_spg_num = tc.get_spacegroup_number_from_atoms( - atoms, symprec=symprec - ) - broken_symmetry = init_spg_num != relaxed_spg_num - - relax_dict = { - "max_stress": max_stress, - "reached_max_steps": reached_max_steps, - "broken_symmetry": broken_symmetry, - tc.TCKeys.final_spg_num: relaxed_spg_num, - tc.TCKeys.init_spg_num: init_spg_num, - } - - except (ValueError, RuntimeError, OSError, KeyError) as exc: - warnings.warn( - f"Failed to relax {formula=}, {mat_id=}: {exc!r}", stacklevel=2 - ) - traceback.print_exc() - err_dict["errors"] += [f"RelaxError: {exc!r}"] - err_dict["error_traceback"] += [traceback.format_exc()] - results_dict = info_dict | relax_dict | err_dict - - write(relax_path, atoms, format="extxyz") - - try: - ph3 = tc.init_phono3py( - atoms, - fc2_supercell=atoms.info["fc2_supercell"], - fc3_supercell=atoms.info["fc3_supercell"], - q_point_mesh=atoms.info["q_point_mesh"], - displacement_distance=displacement_distance, - symprec=symprec, - ) - - ph3, fc2_set, freqs = tc.get_fc2_and_freqs( - ph3, calculator=calculator, pbar_kwargs={"disable": True} - ) - - has_imag_ph_modes = tc.check_imaginary_freqs(freqs) - freqs_dict = { - tc.TCKeys.has_imag_ph_modes: has_imag_ph_modes, - tc.TCKeys.ph_freqs: freqs, - } - - # Determine if conductivity calculation should proceed - if ignore_imaginary_freqs: - # NequIP/Allegro mode: ignore imaginary frequencies - ltc_condition = True - else: - # MACE mode: check both imaginary freqs and broken symmetry - broken_symmetry = relax_dict.get("broken_symmetry", False) - ltc_condition = not has_imag_ph_modes and ( - not broken_symmetry or conductivity_broken_symm - ) - - if ltc_condition: - tc.calculate_fc3_set( - ph3, calculator=calculator, pbar_kwargs={"leave": False} - ) - ph3.produce_fc3(symmetrize_fc3r=True) - else: - reason = [] - if has_imag_ph_modes: - reason.append("imaginary frequencies") - if relax_dict.get("broken_symmetry") and not conductivity_broken_symm: - reason.append("broken symmetry") - warnings.warn( - f"{' and '.join(reason).capitalize()} detected for {mat_id}, " - f"skipping FC3 and LTC calculation!", - stacklevel=2, - ) - - if not ltc_condition: - results_dict = info_dict | relax_dict | freqs_dict | err_dict - - except (ValueError, RuntimeError, OSError, KeyError) as exc: - warnings.warn( - f"Failed to calculate force sets {mat_id}: {exc!r}", stacklevel=2 - ) - traceback.print_exc() - err_dict["errors"] += [f"ForceConstantError: {exc!r}"] - err_dict["error_traceback"] += [traceback.format_exc()] - results_dict = info_dict | relax_dict | err_dict - - # Calculation of conductivity - try: - with tc.tqdm_gridpoints(desc="Conducitivity calc"): - ph3.mesh_numbers = atoms.info["fast_q_point_mesh"] - ph3, fast_kappa_dict, _cond = tc.calculate_conductivity( - ph3, temperatures=temperatures, log_level=2 - ) - fast_results_dict = ( - info_dict | relax_dict | freqs_dict | fast_kappa_dict | err_dict - ) - - if not FAST_ONLY: - with tc.tqdm_gridpoints(desc="Conducitivity calc"): - ph3.mesh_numbers = atoms.info["q_point_mesh"] - ph3, kappa_dict, _cond = tc.calculate_conductivity( - ph3, temperatures=temperatures, log_level=2 - ) - results_dict = ( - info_dict | relax_dict | freqs_dict | kappa_dict | err_dict - ) - - except (ValueError, RuntimeError, OSError, KeyError) as exc: - warnings.warn( - f"Failed to calculate conductivity {mat_id}: {exc!r}", stacklevel=2 - ) - traceback.print_exc() - err_dict["errors"] += [f"ConductivityError: {exc!r}"] - err_dict["error_traceback"] += [traceback.format_exc()] - results_dict = info_dict | relax_dict | freqs_dict | err_dict + results_dict, fast_results_dict = calc_thermal_conductivity_per_structure( + atoms_input, calculator, out_dir + ) - if fast_results_dict is None: - fast_results_dict = relax_dict results_dict[tc.TCKeys.mat_id] = structure_id fast_results_dict[tc.TCKeys.mat_id] = structure_id @@ -354,3 +153,224 @@ def test_thermal_conductivity(mlip: tuple[str, Any]) -> None: fast_df.reset_index().to_json(OUT_PATH / model_name / "fast_kappas.json.gz") with h5py.File(OUT_PATH / model_name / "fast_kappas.hdf5", "w") as f: tc.dict_to_hdf5(fast_kappa_dicts, f) + + +def calc_thermal_conductivity_per_structure( + atoms_input: ase.Atoms, calculator: ase.Calculator, out_dir: Path +) -> tuple[dict[str, Any], dict[str, Any]]: + """ + Calculate thermal conductivity results for a single structure. + + Parameters + ---------- + atoms_input : ase.Atoms + Input atomic structure to evaluate. + calculator : ase.Calculator + Calculator used for the structure relaxation and property evaluation. + out_dir : Path + Directory where intermediate and final results are written. + + Returns + ------- + tuple[dict[str, Any], dict[str, Any]] + A tuple containing the main thermal conductivity result dictionary and + the fast calculation result dictionary. + """ + atoms = atoms_input.copy() + + atoms.calc = calculator + + # Relax structure before calculating thermal conductivity + relax_path = out_dir / "relaxed.extxyz" + + formula = atoms.info.get("name", atoms.get_chemical_formula()) + + mat_id = atoms.info[tc.TCKeys.mat_id] + init_info = deepcopy(atoms.info) + info_dict: dict[str, Any] = { + str(tc.TCKeys.mat_id): mat_id, + str(tc.TCKeys.formula): formula, + } + err_dict: dict[str, list[str]] = {"errors": [], "error_traceback": []} + + # Select filter class + if ase_filter in {"frechet", "exp"}: + filter_cls: type[Filter] = { + "frechet": FrechetCellFilter, + "exp": ExpCellFilter, + }[ase_filter] + else: + # Default to FrechetCellFilter if not specified (for MACE compatibility) + filter_cls = FrechetCellFilter + + # Select optimizer class + optimizer_dict = { + "GPMin": ase.optimize.GPMin, + "GOQN": ase.optimize.GoodOldQuasiNewton, + "BFGSLineSearch": ase.optimize.BFGSLineSearch, + "QuasiNewton": ase.optimize.BFGSLineSearch, + "SciPyFminBFGS": ase.optimize.sciopt.SciPyFminBFGS, + "BFGS": ase.optimize.BFGS, + "LBFGSLineSearch": ase.optimize.LBFGSLineSearch, + "SciPyFminCG": ase.optimize.sciopt.SciPyFminCG, + "FIRE2": ase.optimize.FIRE2, + "FIRE": ase.optimize.FIRE, + "LBFGS": ase.optimize.LBFGS, + } + optim_cls: type[Optimizer] = optimizer_dict[ase_optimizer] + + # Initialize variables that might be needed in error handling + relax_dict: dict[str, Any] = { + "max_stress": None, + "reached_max_steps": False, + "broken_symmetry": False, + } + # initial space group for symmetry breaking detection + init_spg_num = tc.get_spacegroup_number_from_atoms(atoms, symprec=symprec) + + # Relaxation + try: + results_dict = {} + atoms.calc = calculator + if max_steps > 0: + if enforce_relax_symm: + atoms.set_constraint(FixSymmetry(atoms, symprec=symprec)) + filtered_atoms = filter_cls(atoms, mask=[True] * 3 + [False] * 3) + else: + filtered_atoms = filter_cls(atoms) + + optimizer = optim_cls(filtered_atoms, logfile=out_dir / "relax.log") + optimizer.run(fmax=fmax, steps=max_steps) + + step_count = getattr(optimizer, "nsteps", None) # Get optimizer step count + if step_count is None: # fallback to extract from state_dict if available + state = getattr(optimizer, "state_dict", dict)() + step_count = state.get("step", 0) + + reached_max_steps = step_count >= max_steps + if reached_max_steps: + print(f"Material {mat_id=} reached {max_steps=} during relaxation") + + # maximum residual stress component in for xx,yy,zz and xy,yz,xz + # components separately result is a array of 2 elements + max_stress = abs(atoms.get_stress()).reshape((2, 3), order="C").max(axis=1) + + atoms.calc = None + atoms.constraints = None + atoms.info = init_info | atoms.info + + # Check if symmetry was broken during relaxation + relaxed_spg_num = tc.get_spacegroup_number_from_atoms( + atoms, symprec=symprec + ) + broken_symmetry = init_spg_num != relaxed_spg_num + + relax_dict = { + "max_stress": max_stress, + "reached_max_steps": reached_max_steps, + "broken_symmetry": broken_symmetry, + tc.TCKeys.final_spg_num: relaxed_spg_num, + tc.TCKeys.init_spg_num: init_spg_num, + } + + except (ValueError, RuntimeError, OSError, KeyError) as exc: + warnings.warn(f"Failed to relax {formula=}, {mat_id=}: {exc!r}", stacklevel=2) + traceback.print_exc() + ph3 = None + err_dict["errors"] += [f"RelaxError: {exc!r}"] + err_dict["error_traceback"] += [traceback.format_exc()] + results_dict = info_dict | relax_dict | err_dict + return results_dict, results_dict + + write(relax_path, atoms, format="extxyz") + + try: + ph3 = tc.init_phono3py( + atoms, + fc2_supercell=atoms.info["fc2_supercell"], + fc3_supercell=atoms.info["fc3_supercell"], + q_point_mesh=atoms.info["q_point_mesh"], + displacement_distance=displacement_distance, + symprec=symprec, + ) + + ph3, fc2_set, freqs = tc.get_fc2_and_freqs( + ph3, calculator=calculator, pbar_kwargs={"disable": True} + ) + + has_imag_ph_modes = tc.check_imaginary_freqs(freqs) + freqs_dict = { + tc.TCKeys.has_imag_ph_modes: has_imag_ph_modes, + tc.TCKeys.ph_freqs: freqs, + } + + # Determine if conductivity calculation should proceed + if ignore_imaginary_freqs: + # NequIP/Allegro mode: ignore imaginary frequencies + ltc_condition = True + else: + # MACE mode: check both imaginary freqs and broken symmetry + broken_symmetry = relax_dict.get("broken_symmetry", False) + ltc_condition = not has_imag_ph_modes and ( + not broken_symmetry or conductivity_broken_symm + ) + + if ltc_condition: + tc.calculate_fc3_set( + ph3, calculator=calculator, pbar_kwargs={"leave": False} + ) + ph3.produce_fc3(symmetrize_fc3r=True) + else: + reason = [] + if has_imag_ph_modes: + reason.append("imaginary frequencies") + if relax_dict.get("broken_symmetry") and not conductivity_broken_symm: + reason.append("broken symmetry") + warnings.warn( + f"{' and '.join(reason).capitalize()} detected for {mat_id}, " + f"skipping FC3 and LTC calculation!", + stacklevel=2, + ) + + if not ltc_condition: + results_dict = info_dict | relax_dict | freqs_dict | err_dict + return results_dict, results_dict + + except (ValueError, RuntimeError, OSError, KeyError) as exc: + warnings.warn(f"Failed to calculate force sets {mat_id}: {exc!r}", stacklevel=2) + traceback.print_exc() + err_dict["errors"] += [f"ForceConstantError: {exc!r}"] + err_dict["error_traceback"] += [traceback.format_exc()] + results_dict = info_dict | relax_dict | err_dict + return results_dict, results_dict + + # Calculation of conductivity + try: + with tc.tqdm_gridpoints(desc="Conducitivity calc"): + ph3.mesh_numbers = atoms.info["fast_q_point_mesh"] + ph3, fast_kappa_dict, _cond = tc.calculate_conductivity( + ph3, temperatures=temperatures, log_level=2 + ) + fast_results_dict = ( + info_dict | relax_dict | freqs_dict | fast_kappa_dict | err_dict + ) + + if not FAST_ONLY: + with tc.tqdm_gridpoints(desc="Conducitivity calc"): + ph3.mesh_numbers = atoms.info["q_point_mesh"] + ph3, kappa_dict, _cond = tc.calculate_conductivity( + ph3, temperatures=temperatures, log_level=2 + ) + results_dict = info_dict | relax_dict | freqs_dict | kappa_dict | err_dict + + except (ValueError, RuntimeError, OSError, KeyError) as exc: + warnings.warn( + f"Failed to calculate conductivity {mat_id}: {exc!r}", stacklevel=2 + ) + traceback.print_exc() + err_dict["errors"] += [f"ConductivityError: {exc!r}"] + err_dict["error_traceback"] += [traceback.format_exc()] + results_dict = info_dict | relax_dict | freqs_dict | err_dict + return results_dict, results_dict + + return results_dict, fast_results_dict diff --git a/ml_peg/calcs/bulk_crystal/thermal_conductivity/thermal_conductivity.py b/ml_peg/calcs/bulk_crystal/thermal_conductivity/thermal_conductivity.py index 1e124d36e..928d679c2 100644 --- a/ml_peg/calcs/bulk_crystal/thermal_conductivity/thermal_conductivity.py +++ b/ml_peg/calcs/bulk_crystal/thermal_conductivity/thermal_conductivity.py @@ -815,5 +815,5 @@ def load_hdf5_subdir_dicts( with h5py.File(d / filename, "r") as f: dicts[d.name] = hdf5_to_dict(f) except FileNotFoundError: - raise FileNotFoundError(f"File not found: {d / filename}") from None + continue return dicts From 1e903eaff4bed13233006746ccb6f038cf0fcb80 Mon Sep 17 00:00:00 2001 From: bpota Date: Sun, 12 Jul 2026 22:45:41 +0100 Subject: [PATCH 4/4] Update thermal conductivity for current APIs --- .../thermal_conductivity/analyse_thermal_conductivity.py | 8 ++++---- .../thermal_conductivity/app_thermal_conductivity.py | 6 ++++-- .../thermal_conductivity/calc_thermal_conductivity.py | 2 +- 3 files changed, 9 insertions(+), 7 deletions(-) diff --git a/ml_peg/analysis/bulk_crystal/thermal_conductivity/analyse_thermal_conductivity.py b/ml_peg/analysis/bulk_crystal/thermal_conductivity/analyse_thermal_conductivity.py index e6ad9e2b9..48ecd3aa8 100644 --- a/ml_peg/analysis/bulk_crystal/thermal_conductivity/analyse_thermal_conductivity.py +++ b/ml_peg/analysis/bulk_crystal/thermal_conductivity/analyse_thermal_conductivity.py @@ -13,12 +13,12 @@ import pytest from ml_peg.analysis.utils.decorators import build_table, plot_parity -from ml_peg.analysis.utils.utils import load_metrics_config +from ml_peg.analysis.utils.utils import load_metrics_config, write_struct_info from ml_peg.app import APP_ROOT from ml_peg.calcs import CALCS_ROOT from ml_peg.calcs.bulk_crystal.thermal_conductivity import thermal_conductivity as tc +from ml_peg.models import current_models from ml_peg.models.get_models import get_model_names -from ml_peg.models.models import current_models MODELS = get_model_names(current_models) CALC_PATH = CALCS_ROOT / "bulk_crystal" / "thermal_conductivity" / "outputs" @@ -686,7 +686,7 @@ def metrics( } -def test_thermal_conducticity(metrics: dict[str, dict]): +def test_thermal_conductivity(metrics: dict[str, dict]) -> None: """ Run thermal conductivity benchmark tests. @@ -695,4 +695,4 @@ def test_thermal_conducticity(metrics: dict[str, dict]): metrics : dict[str, dict] Metric names and values for all models. """ - return + write_struct_info(data_path=STRUCTURE_FILE, out_path=OUT_PATH) diff --git a/ml_peg/app/bulk_crystal/thermal_conductivity/app_thermal_conductivity.py b/ml_peg/app/bulk_crystal/thermal_conductivity/app_thermal_conductivity.py index 6dfbd64c9..b68827300 100644 --- a/ml_peg/app/bulk_crystal/thermal_conductivity/app_thermal_conductivity.py +++ b/ml_peg/app/bulk_crystal/thermal_conductivity/app_thermal_conductivity.py @@ -10,17 +10,18 @@ plot_from_table_column, ) from ml_peg.app.utils.load import read_plot +from ml_peg.models import current_models from ml_peg.models.get_models import get_model_names -from ml_peg.models.models import current_models # Get all models MODELS = get_model_names(current_models) BENCHMARK_NAME = "Thermal Conductivity" DOCS_URL = ( "https://ddmms.github.io/ml-peg/user_guide/" - "benchmarks/bulk_crystal.html#thermal_conductivity" + "benchmarks/bulk_crystal.html#thermal-conductivity" ) DATA_PATH = APP_ROOT / "data" / "bulk_crystal" / "thermal_conductivity" +INFO_PATH = DATA_PATH / "info.json" class ThermalConductivityApp(BaseApp): @@ -82,4 +83,5 @@ def get_app() -> ThermalConductivityApp: Div(id=f"{BENCHMARK_NAME}-figure-placeholder"), # Div(id=f"{BENCHMARK_NAME}-struct-placeholder"), ], + info_path=INFO_PATH, ) diff --git a/ml_peg/calcs/bulk_crystal/thermal_conductivity/calc_thermal_conductivity.py b/ml_peg/calcs/bulk_crystal/thermal_conductivity/calc_thermal_conductivity.py index 18069c0bb..7acb2227c 100644 --- a/ml_peg/calcs/bulk_crystal/thermal_conductivity/calc_thermal_conductivity.py +++ b/ml_peg/calcs/bulk_crystal/thermal_conductivity/calc_thermal_conductivity.py @@ -28,8 +28,8 @@ from tqdm import tqdm from ml_peg.calcs.bulk_crystal.thermal_conductivity import thermal_conductivity as tc +from ml_peg.models import current_models from ml_peg.models.get_models import load_models -from ml_peg.models.models import current_models MODELS = load_models(current_models)