Skip to content

Add HPHT_CH4_H2O free energy profile of proton hopping benchmark#436

Open
ThomasThevenet wants to merge 1 commit intoddmms:mainfrom
ThomasThevenet:feature/hpht_ch4_h2o-benchmark
Open

Add HPHT_CH4_H2O free energy profile of proton hopping benchmark#436
ThomasThevenet wants to merge 1 commit intoddmms:mainfrom
ThomasThevenet:feature/hpht_ch4_h2o-benchmark

Conversation

@ThomasThevenet
Copy link
Copy Markdown

Pre-review checklist for PR author

PR author must check the checkboxes below when creating the PR.

Summary

"HPHT_CH4_H2O" is a new benchmark of the "Molecular Reactions" category which tests the free energy profile of proton hopping reaction in CH4/H2O liquid from molecular dynamics simulations at 3000 K and 22-69 GPa. The associated analysis performs a molecular recognition algorithm to collect reaction coordinate values and compute the free energy profile associated to the reaction : H3O+ + CH4 = H2O + CH5+. Reference profiles were computed at the DFT (GGA) level (PBE+D3).
Currently three metrics are available :

  1. Global MAE of predicted profiles with respect to reference ones.
  2. MAE on free energy of reaction (F[product] - F[reactant])
  3. MAE on foward free energy barrier (F[transition state] - F[reactant])

The application displays parity plot for MAE of both free energy of reaction and barrier. Interaction of the user with parity plots dynamically displays the predicted and reference free energy profiles of the selected structure.

Linked issue

None.

Progress

  • Calculations
  • Analysis
  • Application
  • Documentation

to do:

  • Prepare reference data

  • Upload reference data to S3 bucket

  • Write documentation

  • Talk about the full pipeline test

  • Testing

Analysis and application components have been tested locally from precomputed molecular dynamics trajectories. Full MD-based pipeline is not manageable with local ressources and has not been tested.

New decorators/callbacks

No new decorators were introduced. A new Dash callback leveraging clickData was implemented to enable interactive visualization of free energy profiles. Selecting a point in the parity plot dynamically displays the corresponding structure's profiles for all MLIP models.

@ElliottKasoar ElliottKasoar added the new benchmark Proposals and suggestions for new benchmarks label Apr 2, 2026
Copy link
Copy Markdown
Collaborator

@ElliottKasoar ElliottKasoar left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @ThomasThevenet, thanks for this, and apologies for the delay in reviewing. It's looking great so far!

If you can please share the data needed to run this, or remind me where you have if you already have, I can test it more thoroughly, but I've left some initial comments/thoughts/questions.

Please also try to add the documentation when you can too.

Let us know if anything is unclear!

write_dir = OUT_PATH / model_name
write_dir.mkdir(parents=True, exist_ok=True)

calc = model.get_calculator()
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
calc = model.get_calculator()
calc = model.get_calculator(precision="low")

You may need to rebase for this to work, but we've added a change where the default precision is "high" (usually float64) since we generally want this for static tests, but previously if you didn't specify a precision explicitly, it would default to the dtypes we set in models.yml, which tended to be equivalent to "low", and I think make sense since we're doing MD here

Comment on lines +8 to +22
level_of_theory: PBE + D3

DF_MAE:
good: 0
bad: 50
unit: kJ/mol
tooltip: "Mean absolute error of reaction free energies"
level_of_theory: PBE + D3

DF#_MAE:
good: 0
bad: 50
unit: kJ/mol
tooltip: "Mean absolute error of free energy barriers"
level_of_theory: PBE + D3
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
level_of_theory: PBE + D3
DF_MAE:
good: 0
bad: 50
unit: kJ/mol
tooltip: "Mean absolute error of reaction free energies"
level_of_theory: PBE + D3
DF#_MAE:
good: 0
bad: 50
unit: kJ/mol
tooltip: "Mean absolute error of free energy barriers"
level_of_theory: PBE + D3
level_of_theory: PBE+D3
DF_MAE:
good: 0
bad: 50
unit: kJ/mol
tooltip: "Mean absolute error of reaction free energies"
level_of_theory: PBE+D3
DF#_MAE:
good: 0
bad: 50
unit: kJ/mol
tooltip: "Mean absolute error of free energy barriers"
level_of_theory: PBE+D3

These need to be without spaces to be matched up correctly

Comment on lines +5 to +19
bad: 50
unit: kJ/mol
tooltip: "Mean absolute error of free energy profiles"
level_of_theory: PBE + D3

DF_MAE:
good: 0
bad: 50
unit: kJ/mol
tooltip: "Mean absolute error of reaction free energies"
level_of_theory: PBE + D3

DF#_MAE:
good: 0
bad: 50
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Genuine question: How much have you considered these thresholds, and are you happy with the idea behind what they represent?

If you're happy with them, great, but if not, we're very happy to discuss their intent further!

Comment on lines +3 to +17
FEP_MAE:
good: 0
bad: 50
unit: kJ/mol
tooltip: "Mean absolute error of free energy profiles"
level_of_theory: PBE + D3

DF_MAE:
good: 0
bad: 50
unit: kJ/mol
tooltip: "Mean absolute error of reaction free energies"
level_of_theory: PBE + D3

DF#_MAE:
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would consider making these slightly more 'human readable'. E.g. you don't need underscores (if it doesn't work, especially with a #, you may need to wrap it in ".

I'd potentially even writing out "Free Energy Profile MAE" etc. since otherwise it's very unlikely people will have any dense of what they are without the tooltips (which of course are there to give more detail, but I think for example DF and DF# are quite confusing on their own).

import pytest
from ase import units
from ase.io import write, read
#from ase.io.extxyz import read
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
#from ase.io.extxyz import read


starting_frames_dir= (
download_s3_data(
key="inputs/molecular_reactions/HPHT_CH4_H2O/HPHT_CH4_H2O.zip",
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Apologies if I missed this, but have you shared this with us? If not, could you so that we can upload it?

Comment on lines +67 to +93
opt = LBFGS(atoms, logfile= write_dir / f"{structure_name}.opt")
opt.run(fmax=0.2)
#VELOCITIES INITIALISATION
rng = np.random.default_rng(seed=13)
MaxwellBoltzmannDistribution(atoms, temperature_K=3000, rng=rng)
Stationary(atoms)
ZeroRotation(atoms)

#OUTPUT SETUP


traj_path = write_dir / f"{structure_name}.extxyz"
log_path = write_dir / f"{structure_name}.log"

def write_frame():
write(traj_path, atoms, format="extxyz", append=True)

#MOLECULAR DYNAMICS
dyn = NoseHooverChainNVT(
atoms,
timestep=0.5 * units.fs,
temperature_K=3000,
tdamp=100 * units.fs,
tchain=3,
tloop=2,
logfile=str(log_path),
)
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's not a requirement, but can I suggest looking at swapping this out for janus_core's molecular dynamics (see examples here: https://stfc.github.io/janus-core/tutorials/python/md.html)?

It will essentially do all of this in a single command (geometry optimisation, setting the temperature distribution, removing rotations, logging).

If you're interested and run into any problems with this, let me know!


MODELS = get_model_names(current_models)
BENCHMARK_NAME = "HPHT_CH4_H2O"
DOCS_URL = "https://ddmms.github.io/ml-peg/user_guide/benchmarks/molecular.html#proton"
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The last part proton will need to match the title you create in the molecular.rst file documenting this test (lower case, with spaces replaced by -, so you may need to change this.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

new benchmark Proposals and suggestions for new benchmarks

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants