diff --git a/docs/source/user_guide/benchmarks/electrolytes.rst b/docs/source/user_guide/benchmarks/electrolytes.rst new file mode 100644 index 000000000..6d16498c9 --- /dev/null +++ b/docs/source/user_guide/benchmarks/electrolytes.rst @@ -0,0 +1,83 @@ +==================== +Electrolytes +==================== + +LIB Electrolyte Inter-Intra Properties +====================================== + +Summary +------- + +These tests examine the model's ability to accurately predict inter- and intra- molecular properties of a common LIB battery electrolyte and solvent. +Inter-molecular forces heavily influence the prediction of density and diffusivity for liquids. Evaluate models on a mix of 200 LIB full electrolyte +and neat solvent configs across a range of densities. The following predicted properties will be tested against PBE DFT: + +Intra-forces +Inter-forces +Inter-energy +Intra-virial +Inter-virial + + +Metrics +------- + +1. RMSE (PBE) + +Root mean square errors for each predicted property compared to PBE data. + +All properties listed above are calculated for each structure. The intra decomposition for a frame is achieved by isolating each molecule and evaluating it with PBE. +The intra properties are then calculated by subtracting the intra property from the total property prediction. The D3 correction is applied both on the models and the PBE functional. + +Computational cost +------------------ + +Small: tests are likely to take seconds to 10 minutes of GPU time per model. + + +Data availability +----------------- + +Input structures: + +* Built from LIB ful electrolyte (LiPF6 EC:EMC) and neat solvent (EC:EMC) configs. + +Reference data: + +* DFT data + + * PBE-D3(BJ) + + +LIB Electrolyte Volume Scans +============================ + +Summary +------- + +Evaluate model energy predictions across battery solvent and battery electrolyte Volume Scans. + +Metrics +------- + +(1) Energy RMSE + +Root mean square error (RMSE) between predicted and reference energy values for each volume scan config. + +Volume scans consist of an initial config, where the molecules are frozen and the volume isotropically expanded or contracted. +The resulting set of configurations represent a scan across different electrolyte densities with all intra properties remaining unchanged. +The relative energy difference between densities is fully dependent on inter-molecular interactions, which heavily influence the density and diffusivity of an electrolyte. +The D3 correction is applied both on the models and the PBE functional. + +Computational cost +------------------ + +Small: tests are likely to take seconds to 10 minutes of GPU time per model. + +Data availability +----------------- + +Input structures: + +* Constructed using the aseMolec package https://github.com/imagdau/aseMolec.git +* PBE diff --git a/docs/source/user_guide/benchmarks/index.rst b/docs/source/user_guide/benchmarks/index.rst index f968dcc35..97e468b48 100644 --- a/docs/source/user_guide/benchmarks/index.rst +++ b/docs/source/user_guide/benchmarks/index.rst @@ -8,6 +8,7 @@ Benchmarks bulk_crystal conformers defect + electrolytes f_block molecular molecular_crystal diff --git a/ml_peg/analysis/electrolytes/LIB_electrolyte_inter_intra/analyse_LIB_electrolyte_inter_intra.py b/ml_peg/analysis/electrolytes/LIB_electrolyte_inter_intra/analyse_LIB_electrolyte_inter_intra.py new file mode 100644 index 000000000..8dea777c1 --- /dev/null +++ b/ml_peg/analysis/electrolytes/LIB_electrolyte_inter_intra/analyse_LIB_electrolyte_inter_intra.py @@ -0,0 +1,222 @@ +"""Analyse LIB electrolyte inter intra benchmark.""" + +from __future__ import annotations + +from pathlib import Path + +from ase.io import read +import numpy as np +import pytest + +from ml_peg.analysis.utils.decorators import ( + build_table, + plot_density_scatter, + plot_parity, +) +from ml_peg.analysis.utils.utils import get_struct_info, load_metrics_config, rmse +from ml_peg.app import APP_ROOT +from ml_peg.calcs import CALCS_ROOT +from ml_peg.models import current_models +from ml_peg.models.get_models import get_model_names + +MODELS = get_model_names(current_models) + +CALC_PATH = CALCS_ROOT / "electrolytes" / "LIB_electrolyte_inter_intra" / "outputs" +REF_PATH = CALC_PATH / "ref" +OUT_PATH = APP_ROOT / "data" / "electrolytes" / "LIB_electrolyte_inter_intra" +METRICS_CONFIG_PATH = Path(__file__).with_name("metrics.yml") +DEFAULT_THRESHOLDS, DEFAULT_TOOLTIPS, _ = load_metrics_config(METRICS_CONFIG_PATH) + +SYSTEM_INFO = get_struct_info( + calc_path=CALC_PATH, + glob_pattern="*.xyz", + index=":", + include_filenames=True, + write_structs=False, + out_path=OUT_PATH, + info_keys=["sys_formula"], +) + +property_metadata = { + "Intra-Forces": ["arrays", "forces_intram"], + "Inter-Forces": ["arrays", "forces_interm"], + "Inter-Energy": ["info", "energy_interm"], + "Intra-Virial": ["info", "virial_intram"], + "Inter-Virial": ["info", "virial_interm"], +} + + +def get_property_results(prop_key: str) -> dict[str, float]: + """ + Get inter-intra results for a specific property. + + Parameters + ---------- + prop_key + String of property name. + + Returns + ------- + dict[str, list] + Dictionary of reference and predicted inter-intra property. + """ + results = {"ref": []} | {mlip: [] for mlip in MODELS} + + stored, property = property_metadata[prop_key] + + for model in results.keys(): + if model == "ref": + configs = read(REF_PATH / "intrainter_PBED3.xyz", ":") + + else: + configs = read(CALC_PATH / model / f"intrainter_{model}_D3.xyz", ":") + + for frame in configs: + frame_data = getattr(frame, stored) + property_data = frame_data[property] + results[model].append(property_data.tolist()) + + if "Forces" in prop_key: + results[model] = np.concatenate(results[model]).flatten() + results[model] = results[model].tolist() + + if "Virial" in prop_key: + results[model] = np.array(results[model]).flatten() + results[model] = results[model].tolist() + + return results + + +def plot_parity_results(prop_key: str, results: dict[str, float]) -> None: + """ + Plot inter-intra property parity plots. + + Parameters + ---------- + prop_key + Name of inter-intra property to be plotted. + results + Results from all models for a single property. + """ + + @plot_parity( + filename=OUT_PATH / f"{prop_key.lower()}_parity.json", + title=prop_key, + x_label=f"Predicted {prop_key} / {DEFAULT_THRESHOLDS[prop_key]['unit']}", + y_label=f"DFT {prop_key} / {DEFAULT_THRESHOLDS[prop_key]['unit']}", + plot_combined=False, + ) + def plot_parity_result() -> dict[str, list[float]]: + """ + Plot the inter-intra propery parity plots. + + Returns + ------- + dict[str, tuple[list[float]]] + Dictionary of reference and predicted inter-intra property. + """ + return results + + plot_parity_result() + + +def plot_density_parity_results(prop_key: str, results: dict[str, float]) -> None: + """ + Plot inter-intra property parity density plots. + + Parameters + ---------- + prop_key + Name of inter-intra property to be plotted. + results + Results from all models for a single property. + """ + + @plot_density_scatter( + filename=OUT_PATH / f"{prop_key.lower()}_density_parity.json", + title=prop_key, + x_label=f"Predicted {prop_key} / {DEFAULT_THRESHOLDS[prop_key]['unit']}", + y_label=f"DFT {prop_key} / {DEFAULT_THRESHOLDS[prop_key]['unit']}", + ) + def plot_density_parity_result() -> dict[str, list[float]]: + """ + Plot the inter-intra propery density parity plots. + + Returns + ------- + dict[str, tuple[list[float]]] + Dictionary of reference and predicted inter-intra property. + """ + results_formatted: dict[str, dict] = {} + ref_vals = results["ref"] + for model, model_pred in results.items(): + if model != "ref": + results_formatted[model] = { + "ref": ref_vals, + "pred": model_pred, + } + return results_formatted + + plot_density_parity_result() + + +@pytest.fixture +def get_property_rmses() -> dict[str, dict]: + """ + Get model prediction RMSEs for all inter-intra properties. + + Returns + ------- + dict[str, dict] + Dictionary of inter-intra properties and the respective RMSE per model. + """ + property_rmse = {prop_key: {} for prop_key in property_metadata.keys()} + + for prop_key in property_metadata.keys(): + results = get_property_results(prop_key) + if "Forces" in prop_key: + plot_density_parity_results(prop_key, results) + else: + plot_parity_results(prop_key, results) + for model in MODELS: + model_rmse = rmse(results["ref"], results[model]) + property_rmse[prop_key][model] = model_rmse + + return property_rmse + + +@pytest.fixture +@build_table( + filename=OUT_PATH / "inter_intra_metrics_table.json", + metric_tooltips=DEFAULT_TOOLTIPS, + thresholds=DEFAULT_THRESHOLDS, +) +def rmse_metrics(get_property_rmses: dict[str, dict]) -> dict[str, dict]: + """ + Get all inter intra RMSE metrics. + + Parameters + ---------- + get_property_rmses + Dictionary for every property containing each model's RMSE. + + Returns + ------- + dict[str, dict] + Dictionary for every property containing each model's RMSE. + """ + return get_property_rmses + + +def test_rmse_metrics( + rmse_metrics: dict[str, dict], +) -> None: + """ + Run inter-intra property test. + + Parameters + ---------- + rmse_metrics + All inter-intra metrics. + """ + return diff --git a/ml_peg/analysis/electrolytes/LIB_electrolyte_inter_intra/metrics.yml b/ml_peg/analysis/electrolytes/LIB_electrolyte_inter_intra/metrics.yml new file mode 100644 index 000000000..92ec8fca4 --- /dev/null +++ b/ml_peg/analysis/electrolytes/LIB_electrolyte_inter_intra/metrics.yml @@ -0,0 +1,36 @@ +metrics: + Intra-Forces: + good: 0.1 + bad: 0.5 + unit: eV/Å + tooltip: "Weighted Root Mean Square Deviation of Intra-Forces" + level_of_theory: PBE + weight: 1 + Inter-Forces: + good: 0.1 + bad: 0.5 + unit: eV/Å + tooltip: "Weighted Root Mean Square Deviation of Inter-Forces" + level_of_theory: PBE + weight: 1 + Inter-Energy: + good: 5 + bad: 60.0 + unit: meV/atom + tooltip: "Weighted Root Mean Square Deviation of Inter-Energy" + level_of_theory: PBE + weight: 1 + Intra-Virial: + good: 10.0 + bad: 50.0 + unit: meV + tooltip: "Weighted Root Mean Square Deviation of Intra-Virial" + level_of_theory: PBE + weight: 1 + Inter-Virial: + good: 1.0 + bad: 50.0 + unit: meV + tooltip: "Weighted Root Mean Square Deviation of Inter-Virial" + level_of_theory: PBE + weight: 1 diff --git a/ml_peg/analysis/electrolytes/LIB_electrolyte_volume_scans/analyse_LIB_electrolyte_volume_scans.py b/ml_peg/analysis/electrolytes/LIB_electrolyte_volume_scans/analyse_LIB_electrolyte_volume_scans.py new file mode 100644 index 000000000..c271bcaba --- /dev/null +++ b/ml_peg/analysis/electrolytes/LIB_electrolyte_volume_scans/analyse_LIB_electrolyte_volume_scans.py @@ -0,0 +1,180 @@ +"""Analyse LIB electrolyte Volume Scans benchmark.""" + +from __future__ import annotations + +from pathlib import Path + +from ase.io import read, write +from aseMolec import extAtoms +import pytest + +from ml_peg.analysis.utils.decorators import build_table, plot_scatter +from ml_peg.analysis.utils.utils import get_struct_info, load_metrics_config, rmse +from ml_peg.app import APP_ROOT +from ml_peg.calcs import CALCS_ROOT +from ml_peg.models import current_models +from ml_peg.models.get_models import get_model_names + +MODELS = get_model_names(current_models) + +CALC_PATH = CALCS_ROOT / "electrolytes" / "LIB_electrolyte_volume_scans" / "outputs" +REF_PATH = CALC_PATH / "ref" +OUT_PATH = APP_ROOT / "data" / "electrolytes" / "LIB_electrolyte_volume_scans" + +METRICS_CONFIG_PATH = Path(__file__).with_name("metrics.yml") +DEFAULT_THRESHOLDS, DEFAULT_TOOLTIPS, _ = load_metrics_config(METRICS_CONFIG_PATH) + +SYSTEM_INFO = get_struct_info( + calc_path=CALC_PATH, + glob_pattern="*.xyz", + index=":", + include_filenames=True, + write_structs=False, + out_path=OUT_PATH, + info_keys=["sys_formula"], +) + +conf_types = ["Solvent", "Electrolyte"] + + +def get_volscan_results( + conf_type: str, +) -> tuple[dict[str, list[float]], dict[str, list[float]]]: + """ + Get relative energies per atom for a type of Volume Scan. + + Parameters + ---------- + conf_type + Name of Volume Scan type to be plotted. + + Returns + ------- + results + Relative energy of each model per Volume Scan config density. + """ + results = {"ref": []} | {mlip: [] for mlip in MODELS} + densities = {"ref": []} | {mlip: [] for mlip in MODELS} + + for model in results: + if model == "ref": + configs = read(REF_PATH / f"{conf_type.lower()}_VS_PBED3.xyz", ":") + + else: + model_path = CALC_PATH / model / f"{conf_type.lower()}_VS_{model}_D3.xyz" + if not model_path.exists(): + continue + + configs = read(model_path, ":") + structs_dir = OUT_PATH / model + structs_dir.mkdir(parents=True, exist_ok=True) + write(structs_dir / f"{model}-{conf_type.lower()}-volscan.extxyz", configs) + + energies = [ + frame.calc.__dict__["results"]["energy"] * 1000 / len(frame) + for frame in configs + ] + relative_energies = energies - min(energies) + densities = [extAtoms.get_density_gcm3(frame) for frame in configs] + + results[model].append(densities) + results[model].append(relative_energies.tolist()) + + return results + + +def plot_volscans(conf_type: str, model: str, results: dict[str, float]) -> None: + """ + Plot Volume Scan scatter plots. + + Parameters + ---------- + conf_type + Name of Volume Scan type to be plotted. + model + Name of MLIP. + results + Results from all models for a single Volume Scan. + """ + + @plot_scatter( + filename=OUT_PATH / f"{conf_type.lower()}_{model}_volscan_scatter.json", + title=f"{conf_type} Volume Scan", + x_label="Density / g cm-3", + y_label="Energy wrt Minimum Energy / meV/atom", + show_line=True, + ) + def plot_result() -> dict[str, list[float]]: + """ + Plot the Volume Scan plots. + + Returns + ------- + model_results + Dictionary of reference and a specific MLIP's Volume Scan energies. + """ + return {"ref": results["ref"], model: results[model]} + + plot_result() + + +@pytest.fixture +def get_volscan_rmses() -> dict[str, dict]: + """ + Get model prediction RMSEs for all volume scan energies. + + Returns + ------- + volscan_rmse + Dictionary of energy RMSE per model for each volume scan. + """ + volscan_rmse = {conf_type: {} for conf_type in conf_types} + + for conf_type in conf_types: + results = get_volscan_results(conf_type) + for model in MODELS: + try: + model_rmse = rmse(results["ref"][1], results[model][1]) + volscan_rmse[conf_type][model] = model_rmse + plot_volscans(conf_type, model, results) + except IndexError: + volscan_rmse[conf_type][model] = None + + return volscan_rmse + + +@pytest.fixture +@build_table( + filename=OUT_PATH / "vol_scan_rmses_table.json", + metric_tooltips=DEFAULT_TOOLTIPS, + thresholds=DEFAULT_THRESHOLDS, +) +def vs_rmse_metrics(get_volscan_rmses: dict[str, dict]) -> dict[str, dict]: + """ + Get all Volume Scan RMSE metrics. + + Parameters + ---------- + get_volscan_rmses + Dictionary of each model's RMSE for all Volume Scans. + + Returns + ------- + dict[str, dict] + Dictionary of each model's RMSE for all Volume Scans. + """ + return get_volscan_rmses + + +def test_vs_rmse_metrics( + vs_rmse_metrics: dict[str, dict], +) -> None: + """ + Run Volume Scans test. + + Parameters + ---------- + vs_rmse_metrics + All Volume Scan metrics. + """ + return diff --git a/ml_peg/analysis/electrolytes/LIB_electrolyte_volume_scans/metrics.yml b/ml_peg/analysis/electrolytes/LIB_electrolyte_volume_scans/metrics.yml new file mode 100644 index 000000000..b84b4ed6a --- /dev/null +++ b/ml_peg/analysis/electrolytes/LIB_electrolyte_volume_scans/metrics.yml @@ -0,0 +1,15 @@ +metrics: + Solvent: + good: 10.0 + bad: 1000.0 + unit: meV/atom + tooltip: "Weighted Root Mean Square Deviation of Solvent Volume Scan energies" + level_of_theory: PBE + weight: 1 + Electrolyte: + good: 10.0 + bad: 1000.0 + unit: meV/atom + tooltip: "Weighted Root Mean Square Deviation of Electrolyte Volume Scan energies" + level_of_theory: PBE + weight: 1 diff --git a/ml_peg/analysis/utils/decorators.py b/ml_peg/analysis/utils/decorators.py index 5305d1c49..2c1c79c69 100644 --- a/ml_peg/analysis/utils/decorators.py +++ b/ml_peg/analysis/utils/decorators.py @@ -39,9 +39,10 @@ def plot_parity( filename: str = "parity.json", symbol_by: list | None = None, symbol_labels: dict[str, str] | None = None, + plot_combined: bool = True, ) -> Callable: """ - Plot parity plot of MLIP results against reference data. + Plot parity plots of MLIP results against reference data. Parameters ---------- @@ -62,6 +63,8 @@ def plot_parity( symbol_labels Optional mapping from ``symbol_by`` values to shorter display names used in the legend. Values absent from this dict are shown as-is. + plot_combined + Option to plot data from all models in a single parity plot. Returns ------- @@ -123,11 +126,12 @@ def plot_parity_wrapper(*args, **kwargs) -> dict[str, Any]: marker_kwargs = { "marker": {"symbol": [group_symbol[g] for g in symbol_by]} } + traces = [] for mlip, value in results.items(): if mlip == "ref": continue - fig.add_trace( + traces.append( go.Scatter( x=value, y=ref, @@ -158,42 +162,84 @@ def plot_parity_wrapper(*args, **kwargs) -> dict[str, Any]: full_fig = fig.full_figure_for_development() x_range = full_fig.layout.xaxis.range y_range = full_fig.layout.yaxis.range + if not plot_combined: + for trace in traces: + fig = go.Figure() + fig.add_trace(trace) + full_fig = fig.full_figure_for_development() + x_range = full_fig.layout.xaxis.range + y_range = full_fig.layout.yaxis.range - lims = [ - np.min([x_range, y_range]), # min of both axes - np.max([x_range, y_range]), # max of both axes - ] + lims = [ + np.min([x_range, y_range]), # min of both axes + np.max([x_range, y_range]), # max of both axes + ] - fig.add_trace( - go.Scatter( - x=lims, - y=lims, - mode="lines", - showlegend=False, + fig.add_trace( + go.Scatter( + x=lims, + y=lims, + mode="lines", + showlegend=False, + ) + ) + + fig.update_traces() + fig.update_layout( + title={"text": title}, + xaxis={"title": {"text": x_label}}, + yaxis={"title": {"text": y_label}}, + ) + Path(filename).parent.mkdir(parents=True, exist_ok=True) + out = Path(filename).with_stem( + f"{Path(filename).stem}_{trace.name}" + ) + fig.write_json(out) + else: + fig = go.Figure() + + for trace in traces: + fig.add_trace(trace) + + full_fig = fig.full_figure_for_development() + x_range = full_fig.layout.xaxis.range + y_range = full_fig.layout.yaxis.range + + lims = [ + np.min([x_range, y_range]), # min of both axes + np.max([x_range, y_range]), # max of both axes + ] + + fig.add_trace( + go.Scatter( + x=lims, + y=lims, + mode="lines", + showlegend=False, + ) ) - ) - fig.update_layout( - title={"text": title}, - xaxis={"title": {"text": x_label}}, - yaxis={"title": {"text": y_label}}, - ) - if symbol_by: fig.update_layout( - legend2={ - "orientation": "h", - "yanchor": "bottom", - "y": 1.02, - "xanchor": "left", - "x": 0, - } + title={"text": title}, + xaxis={"title": {"text": x_label}}, + yaxis={"title": {"text": y_label}}, ) + if symbol_by: + fig.update_layout( + legend2={ + "orientation": "h", + "yanchor": "bottom", + "y": 1.02, + "xanchor": "left", + "x": 0, + } + ) - fig.update_traces() + fig.update_traces() - # Write to file - Path(filename).parent.mkdir(parents=True, exist_ok=True) - fig.write_json(filename) + # Write to file + Path(filename).parent.mkdir(parents=True, exist_ok=True) + fig.write_json(filename) return results diff --git a/ml_peg/app/electrolytes/LIB_electrolyte_inter_intra/app_LIB_electrolyte_inter_intra.py b/ml_peg/app/electrolytes/LIB_electrolyte_inter_intra/app_LIB_electrolyte_inter_intra.py new file mode 100644 index 000000000..747699214 --- /dev/null +++ b/ml_peg/app/electrolytes/LIB_electrolyte_inter_intra/app_LIB_electrolyte_inter_intra.py @@ -0,0 +1,86 @@ +"""Run LIB electrolyte inter intra benchmark 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_cell +from ml_peg.app.utils.load import read_density_plot_for_model, read_plot +from ml_peg.models import current_models +from ml_peg.models.get_models import get_model_names + +MODELS = get_model_names(current_models) + +BENCHMARK_NAME = "LIB electrolyte Inter-Intra Properties" +DOCS_URL = "https://ddmms.github.io/ml-peg/user_guide/benchmarks/electrolytes.html#lib-electrolyte-inter-intra-properties" +DATA_PATH = APP_ROOT / "data" / "electrolytes" / "LIB_electrolyte_inter_intra" +INFO_PATH = ( + APP_ROOT / "data" / "electrolytes" / "LIB_electrolyte_inter_intra" / "info.json" +) + + +class LIBelectrolyteInterIntraApp(BaseApp): + """LIB electrolyte inter intra benchmark app layout and callbacks.""" + + def register_callbacks(self) -> None: + """Register callbacks to app.""" + plots = { + model: { + "Intra-Forces": read_density_plot_for_model( + DATA_PATH / "intra-forces_density_parity.json", + model=model, + id=f"{BENCHMARK_NAME}-{model}-figure", + ), + "Inter-Forces": read_density_plot_for_model( + DATA_PATH / "inter-forces_density_parity.json", + model=model, + id=f"{BENCHMARK_NAME}-{model}-figure", + ), + "Inter-Energy": read_plot( + DATA_PATH / f"inter-energy_parity_{model}.json", + id=f"{BENCHMARK_NAME}-{model}-figure", + ), + "Intra-Virial": read_plot( + DATA_PATH / f"intra-virial_parity_{model}.json", + id=f"{BENCHMARK_NAME}-{model}-figure", + ), + "Inter-Virial": read_plot( + DATA_PATH / f"inter-virial_parity_{model}.json", + id=f"{BENCHMARK_NAME}-{model}-figure", + ), + } + for model in MODELS + } + + plot_from_table_cell( + table_id=self.table_id, + plot_id=f"{BENCHMARK_NAME}-figure-placeholder", + cell_to_plot=plots, + ) + + +def get_app() -> LIBelectrolyteInterIntraApp: + """ + Get LIB electrolyte inter intra benchmark app layout and callback registration. + + Returns + ------- + LIBelectrolyteInterIntraApp + Benchmark layout and callback registration. + """ + return LIBelectrolyteInterIntraApp( + name=BENCHMARK_NAME, + description=( + "Evaluate model inter/intra property prediction " + "for different densities of LIB electrolyte" + " and neat solvent configs" + ), + # docs_url=DOCS_URL, + table_path=DATA_PATH / "inter_intra_metrics_table.json", + extra_components=[ + Div(id=f"{BENCHMARK_NAME}-figure-placeholder"), + ], + info_path=INFO_PATH, + ) diff --git a/ml_peg/app/electrolytes/LIB_electrolyte_volume_scans/app_LIB_electrolyte_volume_scans.py b/ml_peg/app/electrolytes/LIB_electrolyte_volume_scans/app_LIB_electrolyte_volume_scans.py new file mode 100644 index 000000000..150fda95e --- /dev/null +++ b/ml_peg/app/electrolytes/LIB_electrolyte_volume_scans/app_LIB_electrolyte_volume_scans.py @@ -0,0 +1,95 @@ +"""Run LIB electrolyte Volume Scans Benchmark 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_cell, + struct_from_multi_scatters, +) +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 + +MODELS = get_model_names(current_models) + +BENCHMARK_NAME = "LIB Electrolyte Volume-Scans" +DOCS_URL = "https://ddmms.github.io/ml-peg/user_guide/benchmarks/electrolytes.html#lib-electrolyte-volume-scans" +DATA_PATH = APP_ROOT / "data" / "electrolytes" / "LIB_electrolyte_volume_scans" +INFO_PATH = DATA_PATH / "info.json" + + +class LIBelectrolyteVolumeScansApp(BaseApp): + """LIB Electrolyte Volume Scans benchmark app layout and callbacks.""" + + def register_callbacks(self) -> None: + """Register callbacks to app.""" + scatter_plots = { + model: { + "Solvent": read_plot( + DATA_PATH / f"solvent_{model}_volscan_scatter.json", + id=f"{BENCHMARK_NAME}-{model}-figure-solventVS", + ), + "Electrolyte": read_plot( + DATA_PATH / f"electrolyte_{model}_volscan_scatter.json", + id=f"{BENCHMARK_NAME}-{model}-figure-electrolyteVS", + ), + } + for model in MODELS + } + + assets_dir = "/assets/electrolytes/LIB_electrolyte_volume_scans" + structs = { + model: { + "Solvent": f"{assets_dir}/{model}/{model}-solvent-volscan.extxyz", + "Electrolyte": f"{assets_dir}/{model}/" + f"{model}-electrolyte-volscan.extxyz", + } + for model in MODELS + } + + plot_from_table_cell( + table_id=self.table_id, + plot_id=f"{BENCHMARK_NAME}-figure-placeholder", + cell_to_plot=scatter_plots, + ) + + for model in MODELS: + for volscan in ("solvent", "electrolyte"): + struct_from_multi_scatters( + scatter_id=f"{BENCHMARK_NAME}-{model}-figure-{volscan}VS", + struct_id=f"{BENCHMARK_NAME}-struct-placeholder", + structs=[ + structs[model][volscan.capitalize()], + structs[model][volscan.capitalize()], + ], + mode="traj", + ) + + +def get_app() -> LIBelectrolyteVolumeScansApp: + """ + Get Volume Scan benchmark app layout and callback registration. + + Returns + ------- + LIBelectrolyteVolumeScansApp + Benchmark layout and callback registration. + """ + return LIBelectrolyteVolumeScansApp( + name=BENCHMARK_NAME, + description=( + "Evaluate model energy predictions on " + "battery solvent and electrolyte Volume Scans" + ), + docs_url=DOCS_URL, + table_path=DATA_PATH / "vol_scan_rmses_table.json", + extra_components=[ + Div(id=f"{BENCHMARK_NAME}-figure-placeholder"), + Div(id=f"{BENCHMARK_NAME}-struct-placeholder"), + ], + info_path=INFO_PATH, + ) diff --git a/ml_peg/calcs/electrolytes/LIB_electrolyte_inter_intra/calc_LIB_electrolyte_inter_intra.py b/ml_peg/calcs/electrolytes/LIB_electrolyte_inter_intra/calc_LIB_electrolyte_inter_intra.py new file mode 100644 index 000000000..539249fd7 --- /dev/null +++ b/ml_peg/calcs/electrolytes/LIB_electrolyte_inter_intra/calc_LIB_electrolyte_inter_intra.py @@ -0,0 +1,91 @@ +"""Run calculations for LIB electrolyte inter intra benchmark.""" + +from __future__ import annotations + +from copy import copy +from pathlib import Path +import shutil +from typing import Any + +from ase.io import read, write +import pytest + +pytest.importorskip("aseMolec", reason="Please install `asemolec` extra") +from aseMolec import anaAtoms +import numpy as np +import pytest +from tqdm import tqdm + +from ml_peg.calcs.utils.utils import download_s3_data +from ml_peg.models import current_models +from ml_peg.models.get_models import load_models + +MODELS = load_models(current_models) + +OUT_PATH = Path(__file__).parent / "outputs" + + +@pytest.mark.parametrize("mlip", MODELS.items()) +def test_intra_inter(mlip: tuple[str, Any]) -> None: + """ + Run calculations required for intra/inter molecule property comparison. + + Parameters + ---------- + mlip + Name of model use and model to get calculator. + """ + model_name, model = mlip + calc = model.get_calculator(precision="high") + # Add D3 calculator for this test + calc = model.add_d3_calculator(calc) + + out_dir = OUT_PATH / model_name + ref_dir = OUT_PATH / "ref" + out_dir.mkdir(parents=True, exist_ok=True) + ref_dir.mkdir(parents=True, exist_ok=True) + + data_path = ( + download_s3_data( + key="inputs/electrolytes/LIB_electrolyte_inter_intra/LIB_electrolyte_inter_intra.zip", + filename="LIB_electrolyte_inter_intra.zip", + ) + / "LIB_electrolyte_inter_intra" + ) + + structure_paths = data_path.glob("*.xyz") + + for struct_path in tqdm(structure_paths, total=6): + # Copy structure files to output directory + shutil.copy(struct_path, ref_dir / struct_path.name) + + file_prefix = out_dir / f"{struct_path.stem[:-6]}_{model_name}_D3.xyz" + configs = read(struct_path, ":") + for struct in configs: + struct.calc = copy(calc) + struct.info["spin"] = 0 + struct.info["charge"] = 1 + try: + struct.info["energy"] = struct.get_potential_energy() + struct.arrays["forces"] = struct.get_forces() + struct.info["virial"] = ( + -struct.get_stress(voigt=False) * struct.get_volume() + ) + except Exception as e: + print(f"Error calculating energy for {struct_path.name}: {e}") + struct.info["energy"] = float("nan") + struct.arrays["forces"] = np.full((len(struct), 3), np.nan) + struct.info["virial"] = np.full((3, 3), np.nan) + + struct.calc = None + write(file_prefix, configs) + + eval_file_prefix = out_dir + test = read(eval_file_prefix / f"output_{model_name}_D3.xyz", ":") + single_molecule_test = [] + for molsym in ["EMC", "EC", "PF6", "Li"]: + single_molecule_test += read( + eval_file_prefix / f"output{molsym}_{model_name}_D3.xyz", ":" + ) + anaAtoms.collect_molec_results_dict(test, single_molecule_test) + write(eval_file_prefix / f"intrainter_{model_name}_D3.xyz", test) diff --git a/ml_peg/calcs/electrolytes/LIB_electrolyte_volume_scans/calc_LIB_electrolyte_volume_scans.py b/ml_peg/calcs/electrolytes/LIB_electrolyte_volume_scans/calc_LIB_electrolyte_volume_scans.py new file mode 100644 index 000000000..f8415d6c2 --- /dev/null +++ b/ml_peg/calcs/electrolytes/LIB_electrolyte_volume_scans/calc_LIB_electrolyte_volume_scans.py @@ -0,0 +1,77 @@ +"""Run calculations for LIB electrolyte Volume Scans.""" + +from __future__ import annotations + +from copy import copy +from pathlib import Path +import shutil +from typing import Any + +from ase.io import read, write +import numpy as np +import pytest +from tqdm import tqdm + +from ml_peg.calcs.utils.utils import download_s3_data +from ml_peg.models import current_models +from ml_peg.models.get_models import load_models + +MODELS = load_models(current_models) +OUT_PATH = Path(__file__).parent / "outputs" + + +@pytest.mark.parametrize("mlip", MODELS.items()) +def test_volume_scans(mlip: tuple[str, Any]) -> None: + """ + Run calculations required for Volume Scan tests. + + Parameters + ---------- + mlip + Name of models and models used for Volume Scan calculations. + """ + model_name, model = mlip + calc = model.get_calculator(precision="high") + # Add D3 calculator for this test + calc = model.add_d3_calculator(calc) + + out_dir = OUT_PATH / model_name + ref_dir = OUT_PATH / "ref" + out_dir.mkdir(parents=True, exist_ok=True) + ref_dir.mkdir(parents=True, exist_ok=True) + + data_path = ( + download_s3_data( + key="inputs/electrolytes/LIB_electrolyte_volume_scans/LIB_electrolyte_volume_scans.zip", + filename="LIB_electrolyte_volume_scans.zip", + ) + / "LIB_electrolyte_volume_scans" + ) + + structure_paths = data_path.glob("*.xyz") + + for struct_path in tqdm(structure_paths, total=2): + # Copy structure files to output directory + shutil.copy(struct_path, ref_dir / struct_path.name) + + file_prefix = out_dir / f"{struct_path.stem[:-6]}_{model_name}_D3.xyz" + configs = read(struct_path, ":") + for struct in configs: + struct.calc = copy(calc) + struct.info["spin"] = 0 + struct.info["charge"] = 1 + try: + struct.calc = copy(calc) + struct.info["energy"] = struct.get_potential_energy() + struct.arrays["forces"] = struct.get_forces() + struct.info["virial"] = ( + -struct.get_stress(voigt=False) * struct.get_volume() + ) + except Exception as e: + print(f"Error calculating energy for {struct_path.name}: {e}") + struct.info["energy"] = float("nan") + struct.arrays["forces"] = np.full((len(struct), 3), np.nan) + struct.info["virial"] = np.full((3, 3), np.nan) + + struct.calc = None + write(file_prefix, configs) diff --git a/ml_peg/calcs/physicality/oxidation_states/calc_oxidation_states.py b/ml_peg/calcs/physicality/oxidation_states/calc_oxidation_states.py index fdfa21992..f45184fbf 100644 --- a/ml_peg/calcs/physicality/oxidation_states/calc_oxidation_states.py +++ b/ml_peg/calcs/physicality/oxidation_states/calc_oxidation_states.py @@ -36,7 +36,7 @@ def test_iron_oxidation_state_md(mlip: tuple[str, Any]) -> None: Name of model used and model. """ model_name, model = mlip - calc = model.get_calculator(precision="low") + calc = model.get_calculator(precision="high") # Add D3 calculator for this test calc = model.add_d3_calculator(calc)