Added benchmark CRBH20 reaction energy barriers dataset#343
Added benchmark CRBH20 reaction energy barriers dataset#343Quantumplations wants to merge 11 commits into
Conversation
…c_rxn_barriers.py script
…rrect display of model MAEs in browser
|
Hey sorry for the long delay in getting on top of this! I've rebased the PR, changed the naming to CRBH20, did the pre-commit formatting, and also I removed some raw data from the repo. I'm also attaching the data zip file for the S3 bucket. Please review and let me know if there's anything else that needs to be done! |
| if hasattr(model, "add_d3_calculator"): | ||
| calc = model.add_d3_calculator(calc) |
There was a problem hiding this comment.
| if hasattr(model, "add_d3_calculator"): | |
| calc = model.add_d3_calculator(calc) | |
| calc = model.add_d3_calculator(calc) |
There was a problem hiding this comment.
this should be handled
|
|
||
| # ml_peg imports for model management | ||
| from ml_peg.models.get_models import load_models | ||
| from ml_peg.models.models import current_models |
There was a problem hiding this comment.
| from ml_peg.models.models import current_models | |
| from ml_peg.models import current_models |
There was a problem hiding this comment.
due to changes we made
| from ase.io import read, write | ||
| import pytest | ||
|
|
||
| # ml_peg imports for model management |
There was a problem hiding this comment.
| # ml_peg imports for model management | |
| from ml_peg.calcs.utils.utils import download_s3_data |
| # Apply D3 dispersion if the model supports/requires it (standard ml-peg pattern) | ||
| if hasattr(model, "add_d3_calculator"): | ||
| calc = model.add_d3_calculator(calc) | ||
|
|
There was a problem hiding this comment.
| data_dir = ( | |
| download_s3_data( | |
| key="inputs/molecular/CRBH20/CRBH20.zip", | |
| filename="CRBH20.zip", | |
| ) | |
| / "CRBH20" | |
| ) |
| # 2. Loop through Reaction Folders (1 to 20) | ||
| for i in range(1, 21): | ||
| rxn_id = str(i) | ||
| rxn_path = DATA_PATH / rxn_id |
There was a problem hiding this comment.
| rxn_path = DATA_PATH / rxn_id | |
| rxn_path = data_dir / rxn_id |
| # """Run calculations of reaction energy barriers for 20 systems""" | ||
|
|
||
| # from __future__ import annotations | ||
|
|
||
| # from copy import copy | ||
| # from pathlib import Path | ||
| # from typing import Any | ||
|
|
||
| # from ase import units | ||
| # from ase.io import read, write | ||
| # import numpy as np | ||
| # import pytest | ||
|
|
||
| # from ml_peg.calcs.utils.utils import download_s3_data | ||
| # 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" | ||
|
|
||
| # # Unit conversion | ||
| # EV_TO_KJ_PER_MOL = units.mol / units.kJ | ||
|
|
||
| # @pytest.mark.parametrize("mlip", MODELS.items()) | ||
| # def test_rxn_barrier_calculation(mlip: tuple[str, Any]) -> None: | ||
| # """ | ||
| # Run calculations of the reaction energy barriers for the 20 systems in CRBH20. | ||
|
|
||
| # Parameters | ||
| # ---------- | ||
| # mlip | ||
| # Name of model use and model to get calculator. | ||
| # """ | ||
| # model_name, model = mlip | ||
| # calc = model.get_calculator() | ||
|
|
||
| # # Add D3 calculator for this test (for models where applicable) | ||
| # calc = model.add_d3_calculator(calc) | ||
|
|
||
| # # Download X23 dataset | ||
| # lattice_energy_dir = ( | ||
| # download_s3_data( | ||
| # key="inputs/molecular_crystal/X23/X23.zip", | ||
| # filename="lattice_energy.zip", | ||
| # ) | ||
| # / "lattice_energy" | ||
| # ) | ||
|
|
||
| # with open(lattice_energy_dir / "list") as f: | ||
| # systems = f.read().splitlines() | ||
|
|
||
| # for system in systems: | ||
| # molecule_path = lattice_energy_dir / system / "POSCAR_molecule" | ||
| # solid_path = lattice_energy_dir / system / "POSCAR_solid" | ||
| # ref_path = lattice_energy_dir / system / "lattice_energy_DMC" | ||
| # num_molecules_path = lattice_energy_dir / system / "nmol" | ||
|
|
||
| # molecule = read(molecule_path, index=0, format="vasp") | ||
| # molecule.calc = calc | ||
| # molecule.get_potential_energy() | ||
|
|
||
| # solid = read(solid_path, index=0, format="vasp") | ||
| # solid.calc = copy(calc) | ||
| # solid.get_potential_energy() | ||
|
|
||
| # ref = np.loadtxt(ref_path)[0] | ||
| # num_molecules = np.loadtxt(num_molecules_path) | ||
|
|
||
| # solid.info["ref"] = ref | ||
| # solid.info["num_molecules"] = num_molecules | ||
| # solid.info["system"] = system | ||
| # molecule.info["ref"] = ref | ||
| # molecule.info["num_molecules"] = num_molecules | ||
| # molecule.info["system"] = system | ||
|
|
||
| # # Write output structures | ||
| # write_dir = OUT_PATH / model_name | ||
| # write_dir.mkdir(parents=True, exist_ok=True) | ||
| # write(write_dir / f"{system}.xyz", [solid, molecule]) | ||
|
|
||
|
|
||
| # import os | ||
| # from mace.calculators import mace_mp | ||
| # from ase.io import read | ||
|
|
||
| # # 1. Initialize the model ONCE (much faster than reloading it 40 times) | ||
| # # Define the path to your downloaded model | ||
| # # Make sure the filename matches exactly what you downloaded | ||
| # model_path = os.path.abspath("models/mace-mp-0b3-medium.model") | ||
|
|
||
| # print(f"Loading MACE model from: {model_path}") | ||
|
|
||
| # # Load the model using the path | ||
| # macemp = mace_mp( | ||
| # model=model_path, # <--- Point to the local file here | ||
| # dispersion=True, # Keep dispersion on | ||
| # default_dtype="float64", | ||
| # device="cuda" # Use "cpu" if you don't have a GPU | ||
| # ) | ||
|
|
||
| # # Conversion: 1 eV = 23.0605 kcal/mol | ||
| # EV_TO_KCAL = 23.0605 | ||
|
|
||
| # print(f"{'Rxn ID':<8} | {'Barrier (eV)':<12} | {'Barrier (kcal/mol)':<18}") | ||
| # print("="*50) | ||
|
|
||
| # # 2. Loop through all reaction folders (assuming 1 to 20) | ||
| # for i in range(1, 21): | ||
| # rxn_id = str(i) | ||
| # energies = {} | ||
|
|
||
| # # Check both Reactant and Transition State | ||
| # for state in ['react', 'ts']: | ||
| # # Construct the path: e.g., "1/react/POSCAR" | ||
| # path = os.path.join(rxn_id, state, 'POSCAR') | ||
|
|
||
| # if os.path.exists(path): | ||
| # try: | ||
| # # Read the geometry | ||
| # atoms = read(path) | ||
|
|
||
| # # Attach the calculator we loaded earlier | ||
| # atoms.calc = macemp | ||
|
|
||
| # # Calculate Potential Energy | ||
| # e_pot = atoms.get_potential_energy() | ||
| # energies[state] = e_pot | ||
|
|
||
| # # OPTIONAL: Write to file to save progress | ||
| # # This overwrites ('w') to prevent duplicate lines if you re-run | ||
| # output_file = os.path.join(rxn_id, state, 'energy-mace') | ||
| # with open(output_file, "w") as f: | ||
| # print(e_pot, file=f) | ||
|
|
||
| # except Exception as e: | ||
| # print(f"Error calculating {path}: {e}") | ||
| # else: | ||
| # # Silently skip if folder doesn't exist | ||
| # # (useful if you only have some folders) | ||
| # pass | ||
|
|
||
| # # 3. Calculate and Print Barrier | ||
| # if 'react' in energies and 'ts' in energies: | ||
| # barrier_ev = energies['ts'] - energies['react'] | ||
| # barrier_kcal = barrier_ev * EV_TO_KCAL | ||
| # print(f"{rxn_id:<8} | {barrier_ev:.4f} | {barrier_kcal:.4f}") | ||
| # else: | ||
| # print(f"{rxn_id:<8} | {'MISSING DATA':<30}") |
I've uplaoded and written a code suggestion for the download. If you accept hte code suggestion, remember to fetch and pull the changes to your local machine before doing any edits, otherwise they will be overwritten. Thanks! |
|
|
||
| # ml_peg imports | ||
| from ml_peg.analysis.utils.decorators import build_table, plot_parity | ||
| from ml_peg.analysis.utils.utils import build_d3_name_map, load_metrics_config, mae |
There was a problem hiding this comment.
| from ml_peg.analysis.utils.utils import build_d3_name_map, load_metrics_config, mae | |
| from ml_peg.analysis.utils.utils import build_dispersion_name_map, load_metrics_config, mae |
| from ml_peg.app import APP_ROOT | ||
| from ml_peg.calcs import CALCS_ROOT | ||
| from ml_peg.models.get_models import get_model_names | ||
| from ml_peg.models.models import current_models |
There was a problem hiding this comment.
| from ml_peg.models.models import current_models | |
| from ml_peg.models import current_models |
|
|
||
| # --- Configuration --- | ||
| MODELS = get_model_names(current_models) | ||
| D3_MODEL_NAMES = build_d3_name_map(MODELS) |
There was a problem hiding this comment.
| D3_MODEL_NAMES = build_d3_name_map(MODELS) | |
| D3_MODEL_NAMES = build_dispersion_name_map(MODELS) |
| ) | ||
| 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 |
There was a problem hiding this comment.
| from ml_peg.models.models import current_models | |
| from ml_peg.models import current_models |
| OUT_PATH = Path(__file__).parent / "outputs" | ||
|
|
||
| # Unit conversion: 1 eV = 23.0605 kcal/mol | ||
| EV_TO_KCAL = 23.0605 |
There was a problem hiding this comment.
| EV_TO_KCAL = 23.0605 | |
| EV_TO_KCAL = units.mol / units.kcal |
Assuming you also add the units import above. In general it's better to import defined values than defining this ourselves everywhere
| from pathlib import Path | ||
| from typing import Any | ||
|
|
||
| from ase.io import read, write |
There was a problem hiding this comment.
| from ase.io import read, write | |
| from ase import units | |
| from ase.io import read, write |
See below
|
|
||
| # 1. Initialize the Calculator | ||
| # The ml-peg wrapper handles device selection (cuda/cpu) and loading | ||
| calc = model.get_calculator() |
There was a problem hiding this comment.
| calc = model.get_calculator() | |
| calc = model.get_calculator(precision="high") |
This is the default, but it's good to specify that we want 'high' (float64, usually) precision when we get the calculator
| print("=" * 50) | ||
|
|
||
| # 2. Loop through Reaction Folders (1 to 20) | ||
| for i in range(1, 21): |
There was a problem hiding this comment.
It would be nice to wrap this in use tqdm here to show the progress as we go through these iterations. This also potentially reduces the need for the print statements - we don't really need to print the barrier to the terminal, but is useful to see progress
| atoms.calc = calc | ||
|
|
||
| # Compute Energy | ||
| e_pot = atoms.get_potential_energy() |
There was a problem hiding this comment.
This isn't something we strictly enforce yet, but to get ahead of things, it would be useful to wrap this in a try...except to catch any errors e.g. if a model doesn't support an element in the dataset. If it fails, we'd probably want the energy to be set to NaN
| atoms.calc = calc | ||
|
|
||
| # Compute Energy | ||
| e_pot = atoms.get_potential_energy() |
There was a problem hiding this comment.
| atoms.calc = calc | |
| # Compute Energy | |
| e_pot = atoms.get_potential_energy() | |
| atoms.calc = calc | |
| atoms.info.setdefault("charge", 0) | |
| atoms.info.setdefault("spin", 1) | |
| # Compute Energy | |
| e_pot = atoms.get_potential_energy() |
Unless these are read in from the input files, generally we need to set these, else certain models fail
| # If the file doesn't exist, we provide defaults, but usually it should exist. | ||
| try: | ||
| DEFAULT_THRESHOLDS, DEFAULT_TOOLTIPS, DEFAULT_WEIGHTS = load_metrics_config( | ||
| METRICS_CONFIG_PATH | ||
| ) | ||
| except FileNotFoundError: | ||
| # Fallback defaults if metrics.yml is missing | ||
| DEFAULT_THRESHOLDS = {"MAE": [1.0, 5.0]} # Green < 1.0, Red > 5.0 | ||
| DEFAULT_TOOLTIPS = {"MAE": "Mean Absolute Error"} | ||
| DEFAULT_WEIGHTS = {} |
There was a problem hiding this comment.
| # If the file doesn't exist, we provide defaults, but usually it should exist. | |
| try: | |
| DEFAULT_THRESHOLDS, DEFAULT_TOOLTIPS, DEFAULT_WEIGHTS = load_metrics_config( | |
| METRICS_CONFIG_PATH | |
| ) | |
| except FileNotFoundError: | |
| # Fallback defaults if metrics.yml is missing | |
| DEFAULT_THRESHOLDS = {"MAE": [1.0, 5.0]} # Green < 1.0, Red > 5.0 | |
| DEFAULT_TOOLTIPS = {"MAE": "Mean Absolute Error"} | |
| DEFAULT_WEIGHTS = {} | |
| DEFAULT_THRESHOLDS, DEFAULT_TOOLTIPS, DEFAULT_WEIGHTS = load_metrics_config( | |
| METRICS_CONFIG_PATH | |
| ) |
If the metrics.yml doesn't exist, I think we want it to raise an error, so we can fix it. This makes it harder to catch
| REF_BARRIERS_KCAL = { | ||
| 1: 1.7194 * 23.0605, | ||
| 2: 1.9241 * 23.0605, | ||
| 3: 1.7499 * 23.0605, | ||
| 4: 1.8238 * 23.0605, | ||
| 5: 1.7237 * 23.0605, | ||
| 6: 1.5653 * 23.0605, | ||
| 7: 1.0911 * 23.0605, | ||
| 8: 1.8983 * 23.0605, | ||
| 9: 1.5477 * 23.0605, | ||
| 10: 1.7115 * 23.0605, | ||
| 11: 1.7379 * 23.0605, | ||
| 12: 2.0361 * 23.0605, | ||
| 13: 1.8739 * 23.0605, | ||
| 14: 1.9760 * 23.0605, | ||
| 15: 1.8865 * 23.0605, | ||
| 16: 1.5741 * 23.0605, | ||
| 17: 1.2587 * 23.0605, | ||
| 18: 1.7497 * 23.0605, | ||
| 19: 1.6989 * 23.0605, | ||
| 20: 1.7654 * 23.0605, | ||
| } |
There was a problem hiding this comment.
As in the calc, I would set a conversion factor using units.mol / units.kcal (probably as a single constant, not written out each time).
It's also ok to stick to eV if you want, it doesn't have to be kcal/mol
|
|
||
|
|
||
| def numeric_sort_key(filepath: Path): | ||
| """Sort helper to ensure files 1, 2, ... 10 come in numerical order, not alpha.""" |
There was a problem hiding this comment.
Alpha? As in alphanumerically? I think I'd either use the full word, or leave is sorting numerically, I think it's clear enough
| # --- DEBUGGING START --- | ||
| print(f"\nDEBUG: CALC_PATH is set to: {CALC_PATH.resolve()}") | ||
| print(f"DEBUG: OUT_PATH is set to: {OUT_PATH.resolve()}") | ||
| # --- DEBUGGING END --- |
There was a problem hiding this comment.
I don't think we really need these debugging prints?
| # --- DEBUGGING START --- | |
| print(f"\nDEBUG: CALC_PATH is set to: {CALC_PATH.resolve()}") | |
| print(f"DEBUG: OUT_PATH is set to: {OUT_PATH.resolve()}") | |
| # --- DEBUGGING END --- |
| if not model_dir.exists(): | ||
| print(f"[MISSING] -> Skipped {model_dir}") | ||
| continue | ||
| print("[FOUND]") | ||
| # --- DEBUGGING END --- | ||
| if not model_dir.exists(): | ||
| continue |
There was a problem hiding this comment.
| if not model_dir.exists(): | |
| print(f"[MISSING] -> Skipped {model_dir}") | |
| continue | |
| print("[FOUND]") | |
| # --- DEBUGGING END --- | |
| if not model_dir.exists(): | |
| continue | |
| if not model_dir.exists(): | |
| continue |
I don't think we need debugging in this form. We will make it clearer which models were successful in other ways in future.
You also have duplicated if not model_dir.exists(): checks, so the second is never reached?
| # Handle missing data (e.g., if calc failed) | ||
| # For parity plots, lists must be equal length. | ||
| # We append None or NaN, though ml-peg might prefer dropping the point. | ||
| # Here we assume completeness or append 0.0 with a warning. | ||
| model_barriers.append(None) |
There was a problem hiding this comment.
| # Handle missing data (e.g., if calc failed) | |
| # For parity plots, lists must be equal length. | |
| # We append None or NaN, though ml-peg might prefer dropping the point. | |
| # Here we assume completeness or append 0.0 with a warning. | |
| model_barriers.append(None) | |
| model_barriers.append(np.nan) |
Generally we're going to NaNs for missing/failed results
|
|
||
| # Extract ML Barrier (calculated in the previous script) | ||
| # stored as "barrier_kcal" | ||
| barrier_ml = reactant.info.get("barrier_kcal", 0.0) |
There was a problem hiding this comment.
| barrier_ml = reactant.info.get("barrier_kcal", 0.0) | |
| barrier_ml = reactant.info.get("barrier_kcal", np.nan) |
It this key doesn't exist, we should be consistent with the above failure output, so either NaN or None if we didn't change above
|
|
||
| # Store reference energies (only once, during the first model loop) | ||
| if not ref_stored: | ||
| ref_val = REF_BARRIERS_KCAL.get(rxn_id, 0.0) |
There was a problem hiding this comment.
| ref_val = REF_BARRIERS_KCAL.get(rxn_id, 0.0) | |
| ref_val = REF_BARRIERS_KCAL.get(rxn_id, np.nan) |
As above, although this shouldn't really happen. 0 will give confusing results
| results[model_name] = model_barriers | ||
|
|
||
| # Mark reference as stored so we don't duplicate it | ||
| if any(x is not None for x in model_barriers): |
There was a problem hiding this comment.
Generally it's easier to do this within if not ref_stored:, so we're not doing every checks for every model.
To be honest, if the reference is missing so gets set to None/NaN, I think something has probably gone wrong that continuing isn't necessarily useful anyway, so you could probably simplify REF_BARRIERS_KCAL.get(, as that's entirely in your control, and never set a None/NaN for the reference
| # Filter out None values in case of failed calculations | ||
| y_true = [] | ||
| y_pred = [] | ||
| for r, p in zip(reaction_barriers["ref"], reaction_barriers[model_name], strict=False): | ||
| if r is not None and p is not None: | ||
| y_true.append(r) | ||
| y_pred.append(p) |
There was a problem hiding this comment.
We don't want to filter out None/NaNs really. This is also an advantage to NaN over None as generally it's better handled - the mean of something with a NaN in it is NaN.
If all but one calculation failed, we don't want the score to be that error, the test has failed. In future, with element filtering, we'll use partial results more effectively
| good: 1.0 # Any error below 1.0 will be Green (Chemical Accuracy) | ||
| bad: 3.0 # Any error above 3.0 will be Red |
There was a problem hiding this comment.
| good: 1.0 # Any error below 1.0 will be Green (Chemical Accuracy) | |
| bad: 3.0 # Any error above 3.0 will be Red | |
| good: 1.0 | |
| bad: 3.0 |
You don't really need to explain what these do, and the colours are also not necessarily correct.
Genuine question: How much have you considered these thresholds/do you understand what they should represent?
| @@ -0,0 +1,7 @@ | |||
| metrics: | |||
| MAE: | |||
| tooltip: "Mean Absolute Error (kcal/mol)" | |||
There was a problem hiding this comment.
| tooltip: "Mean Absolute Error (kcal/mol)" | |
| tooltip: "Mean Absolute Error" |
This should be taken automatically fro the unit
| # Placeholder docs link | ||
| DOCS_URL = "https://github.com/ddmms/ml-peg" |
There was a problem hiding this comment.
Can you add documentation for this, and then update this please?
| # Dash Asset Paths: The 'assets' folder is mounted at APP_ROOT/data | ||
| # So we construct the URL path starting from there. | ||
| structs = [ | ||
| f"assets/reaction_barriers/CRBH20/{valid_model}/{f.name}" for f in files |
There was a problem hiding this comment.
| f"assets/reaction_barriers/CRBH20/{valid_model}/{f.name}" for f in files | |
| f"/assets/reaction_barriers/CRBH20/{valid_model}/{f.name}" for f in files |
Minor fix due to backend changes
| # Initialize Dash | ||
| # IMPORTANT: We set assets_folder to 'ml_peg/app/data' so it finds your JSONs/XYZs |
There was a problem hiding this comment.
| # Initialize Dash | |
| # IMPORTANT: We set assets_folder to 'ml_peg/app/data' so it finds your JSONs/XYZs |
We don't really need these sorts of comments
| # Placeholder docs link | ||
| DOCS_URL = "https://github.com/ddmms/ml-peg" | ||
|
|
||
| # Point to where we saved the JSONs (ml_peg/app/data/reaction_barriers/CRBH20) |
There was a problem hiding this comment.
| # Point to where we saved the JSONs (ml_peg/app/data/reaction_barriers/CRBH20) |
This doesn't really add anything
There was a problem hiding this comment.
I believe this has been uploaded, so we shouldn't need this as part of the pull request anymore
Pre-review checklist for PR author
PR author must check the checkboxes below when creating the PR.
Summary
Linked issue
Resolves 1 issue.
Resolves #291
Progress
Testing
New decorators/callbacks