diff --git a/docs/source/user_guide/benchmarks/conformers.rst b/docs/source/user_guide/benchmarks/conformers.rst index 8c5e0c30e..e7a25a8c8 100644 --- a/docs/source/user_guide/benchmarks/conformers.rst +++ b/docs/source/user_guide/benchmarks/conformers.rst @@ -93,3 +93,52 @@ Reference data: * Same as input data * :math:`DLPNO-CCSD(T)` level of theory: a local coupled-cluster method. + + +TorsionNet500CCSDT +=================== + +Summary +------- + +Performance in predicting torsional energy profiles for 500 diverse organic +molecular fragments. Reference data from DLPNO-CCSD(T)/CBS calculations. + +Metrics +------- + +1. RMSE of relative torsional energy profile +2. MAE of relative torsional energy profile + +For each fragment, a torsion scan samples the energy at a series of dihedral +angles. Both the reference and predicted energies are mean-centered per scan, +so only the shape of the profile is compared rather than absolute energy +offsets. RMSE and MAE are calculated between the reference and predicted +profiles for each scan, then averaged across all 500 fragments. + +Computational cost +------------------ + +Medium: tests are likely to take minutes to run on GPU, or less than an hour on +CPU for each model, since each fragment requires a single-point energy +calculation across ~20 conformers in its torsion scan, repeated for all 500 +fragments. + +Data availability +----------------- + +Input structures: + +* B. K. Rai, V. Sresht, Q. Yang, R. Unwalla, M. Tu, A. M. Mathiowetz, and + G. A. Bakken, TorsionNet: A Deep Neural Network to Rapidly Predict + Small-Molecule Torsional Energy Profiles with the Accuracy of Quantum + Mechanics, Journal of Chemical Information and Modeling 62 (2022), 785-800. + PMID: 35119861. + +Reference data: + +* J. L. Weber, R. D. Guha, G. Agarwal, Y. Wei, A. A. Fike, X. Xie, + J. Stevenson, B. Santra, R. A. Friesner, K. Leswing, M. D. Halls, R. Abel, + and L. D. Jacobson, Efficient Long-Range Machine Learning Force Fields for + Liquid and Materials Properties, arXiv:2505.06462 (2025). +* DLPNO-CCSD(T)/CBS level of theory. diff --git a/ml_peg/analysis/conformers/TorsionNet500CCSDT/analyse_TorsionNet500CCSDT.py b/ml_peg/analysis/conformers/TorsionNet500CCSDT/analyse_TorsionNet500CCSDT.py new file mode 100644 index 000000000..9ce413dad --- /dev/null +++ b/ml_peg/analysis/conformers/TorsionNet500CCSDT/analyse_TorsionNet500CCSDT.py @@ -0,0 +1,368 @@ +"""Analyse TorsionNet500CCSDT benchmark.""" + +from __future__ import annotations + +from pathlib import Path + +from ase.io import read, write +import numpy as np +import plotly.graph_objects as go +import pytest + +from ml_peg.analysis.utils.decorators import build_table +from ml_peg.analysis.utils.utils import get_struct_info, load_metrics_config, mae, 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 / "conformers" / "TorsionNet500CCSDT" / "outputs" +OUT_PATH = APP_ROOT / "data" / "conformers" / "TorsionNet500CCSDT" + +METRICS_CONFIG_PATH = Path(__file__).with_name("metrics.yml") +DEFAULT_THRESHOLDS, DEFAULT_TOOLTIPS, DEFAULT_WEIGHTS = load_metrics_config( + METRICS_CONFIG_PATH +) + +# Elemental info for the app's filtering, saved per-fragment (rather than just the +# combined set) in preparation for future partial filtering. +get_struct_info( + calc_path=CALC_PATH, + glob_pattern="*.xyz", + index="0", + include_filenames=True, + write_info=True, + write_structs=False, + out_path=OUT_PATH, +) + + +@pytest.fixture +def fragment_rmse() -> dict[str, dict[str, list]]: + """ + Get per-fragment RMSE of the relative torsional energy profile for each model. + + Returns + ------- + dict[str, dict[str, list]] + Per model, the fragment labels, the corresponding RMSE and MAE for each + torsion scan, and the mean-centered angle/energy profile for each scan + (for plotting torsion curves). + """ + results = {} + + for model_name in MODELS: + model_dir = CALC_PATH / model_name + + labels = [] + rmse_scans = [] + mae_scans = [] + profiles = [] + + if model_dir.exists(): + xyz_files = sorted(model_dir.glob("*.xyz")) + + for xyz_file in xyz_files: + atoms = read(xyz_file, ":") + + angle = [a.info["torsion_angle"] for a in atoms] + ref_energy = np.array([a.info["ref_energy"] for a in atoms]) + model_energy = np.array([a.info["model_energy"] for a in atoms]) + + # Mean-center so only the relative torsional profile is compared. + ref_rel_energy = ref_energy - ref_energy.mean() + model_rel_energy = model_energy - model_energy.mean() + + labels.append(xyz_file.stem) + rmse_scans.append(rmse(ref_rel_energy, model_rel_energy)) + mae_scans.append(mae(ref_rel_energy, model_rel_energy)) + + order = np.argsort(angle) + profiles.append( + { + "angle": np.asarray(angle)[order].tolist(), + "ref_rel_energy": ref_rel_energy[order].tolist(), + "model_rel_energy": model_rel_energy[order].tolist(), + } + ) + + results[model_name] = { + "labels": labels, + "rmse": rmse_scans, + "mae": mae_scans, + "profiles": profiles, + } + + return results + + +@pytest.fixture +def get_rmse(fragment_rmse: dict[str, dict[str, list]]) -> dict[str, float | None]: + """ + Get mean RMSE across TorsionNet500CCSDT torsion scans. + + Parameters + ---------- + fragment_rmse + Per-fragment RMSE and labels for each model. + + Returns + ------- + dict[str, float | None] + Mean RMSE per model, ignoring fragments where the calculation failed. + """ + results = {} + + for model_name in MODELS: + rmse_scans = fragment_rmse[model_name]["rmse"] + results[model_name] = float(np.nanmean(rmse_scans)) if rmse_scans else None + + return results + + +@pytest.fixture +def get_mae(fragment_rmse: dict[str, dict[str, list]]) -> dict[str, float | None]: + """ + Get mean MAE across TorsionNet500CCSDT torsion scans. + + Parameters + ---------- + fragment_rmse + Per-fragment RMSE, MAE, and labels for each model. + + Returns + ------- + dict[str, float | None] + Mean MAE per model, ignoring fragments where the calculation failed. + """ + results = {} + + for model_name in MODELS: + mae_scans = fragment_rmse[model_name]["mae"] + results[model_name] = float(np.nanmean(mae_scans)) if mae_scans else None + + return results + + +def plot_fragment_metric_figure( + model_name: str, metric_name: str, labels: list[str], metric_values: list[float] +) -> go.Figure: + """ + Build a scatter plot of a per-fragment metric for one model. + + Parameters + ---------- + model_name + Name of the model the scatter plot is for. + metric_name + Name of the metric being plotted, e.g. ``"RMSE"`` or ``"MAE"``. + labels + Fragment labels, in scatter point order. + metric_values + Per-fragment metric values, in the same order as ``labels``. + + Returns + ------- + go.Figure + Scatter plot of the metric against fragment. + """ + fig = go.Figure() + fig.add_trace( + go.Scatter( + x=list(range(len(labels))), + y=metric_values, + mode="markers", + customdata=labels, + hovertemplate=( + "Fragment: %{customdata}
" + f"{metric_name}: " + "%{y:.4f} eV
" + ), + showlegend=False, + ) + ) + fig.update_layout( + title={"text": f"Per-fragment {metric_name} - {model_name}"}, + xaxis={"title": {"text": "Fragment"}}, + yaxis={"title": {"text": f"{metric_name} / eV"}}, + ) + return fig + + +@pytest.fixture +def fragment_scatter_figures(fragment_rmse: dict[str, dict[str, list]]) -> None: + """ + Save per-model scatter plots of per-fragment RMSE and MAE for the app. + + Parameters + ---------- + fragment_rmse + Per-fragment RMSE, MAE, and labels for each model. + """ + for model_name in MODELS: + labels = fragment_rmse[model_name]["labels"] + + if not labels: + continue + + out_dir = OUT_PATH / model_name + out_dir.mkdir(parents=True, exist_ok=True) + + rmse_fig = plot_fragment_metric_figure( + model_name, "RMSE", labels, fragment_rmse[model_name]["rmse"] + ) + rmse_fig.write_json(out_dir / "fragment_rmse_scatter.json") + + mae_fig = plot_fragment_metric_figure( + model_name, "MAE", labels, fragment_rmse[model_name]["mae"] + ) + mae_fig.write_json(out_dir / "fragment_mae_scatter.json") + + +def plot_torsion_curve_figure( + model_name: str, label: str, profile: dict[str, list] +) -> go.Figure: + """ + Build a torsion energy profile plot for one fragment. + + Parameters + ---------- + model_name + Name of the model the profile was predicted with. + label + Fragment label the profile belongs to. + profile + Mean-centered ``angle``, ``ref_rel_energy``, and ``model_rel_energy`` for + the scan. + + Returns + ------- + go.Figure + Line plot of relative energy against dihedral angle. + """ + fig = go.Figure() + fig.add_trace( + go.Scatter( + x=profile["angle"], + y=profile["ref_rel_energy"], + mode="lines+markers", + name="CCSD(T)", + ) + ) + fig.add_trace( + go.Scatter( + x=profile["angle"], + y=profile["model_rel_energy"], + mode="lines+markers", + name=model_name, + ) + ) + fig.update_layout( + title={"text": f"{label} - {model_name}"}, + xaxis={"title": {"text": "Dihedral angle / deg"}}, + yaxis={"title": {"text": "Relative energy / eV"}}, + ) + return fig + + +@pytest.fixture +def torsion_curve_figures(fragment_rmse: dict[str, dict[str, list]]) -> None: + """ + Save per-fragment torsion energy profile plots for the app. + + Parameters + ---------- + fragment_rmse + Per-fragment RMSE, labels, and profiles for each model. + """ + for model_name in MODELS: + labels = fragment_rmse[model_name]["labels"] + profiles = fragment_rmse[model_name]["profiles"] + + if not labels: + continue + + out_dir = OUT_PATH / model_name / "torsion_curves" + out_dir.mkdir(parents=True, exist_ok=True) + + for label, profile in zip(labels, profiles, strict=True): + fig = plot_torsion_curve_figure(model_name, label, profile) + fig.write_json(out_dir / f"{label}.json") + + +@pytest.fixture +def torsion_trajectories() -> None: + """ + Save per-fragment torsion scan trajectories for the app's structure viewer. + + Geometries are identical across all models, since only single-point energies + are calculated on the same reference conformers, so this only needs writing + once, from the mock calculator's output, matching the ``info.json`` elemental + info above. + """ + mock_dir = CALC_PATH / "mock" + if not mock_dir.exists(): + return + + out_dir = OUT_PATH / "torsion_trajectories" + out_dir.mkdir(parents=True, exist_ok=True) + + for xyz_file in sorted(mock_dir.glob("*.xyz")): + atoms = read(xyz_file, ":") + angle = [a.info["torsion_angle"] for a in atoms] + order = np.argsort(angle) + write(out_dir / f"{xyz_file.stem}.xyz", [atoms[i] for i in order]) + + +@pytest.fixture +@build_table( + filename=OUT_PATH / "torsionnet500ccsdt_metrics_table.json", + metric_tooltips=DEFAULT_TOOLTIPS, + thresholds=DEFAULT_THRESHOLDS, +) +def metrics( + get_rmse: dict[str, float | None], get_mae: dict[str, float | None] +) -> dict[str, dict]: + """ + Get all TorsionNet500CCSDT metrics. + + Parameters + ---------- + get_rmse + Mean RMSE per model. + get_mae + Mean MAE per model. + + Returns + ------- + dict[str, dict] + Metric names and values for all models. + """ + return { + "RMSE": get_rmse, + "MAE": get_mae, + } + + +def test_torsionnet500ccsdt( + metrics: dict[str, dict], + fragment_scatter_figures: None, + torsion_curve_figures: None, + torsion_trajectories: None, +) -> None: + """ + Run TorsionNet500CCSDT analysis. + + Parameters + ---------- + metrics + All TorsionNet500CCSDT metrics. + fragment_scatter_figures + Per-model fragment-RMSE/MAE scatter figures (side-effect only). + torsion_curve_figures + Per-fragment torsion curve figures (side-effect only). + torsion_trajectories + Per-fragment torsion scan trajectories (side-effect only). + """ + return diff --git a/ml_peg/analysis/conformers/TorsionNet500CCSDT/metrics.yml b/ml_peg/analysis/conformers/TorsionNet500CCSDT/metrics.yml new file mode 100644 index 000000000..c4daf8e92 --- /dev/null +++ b/ml_peg/analysis/conformers/TorsionNet500CCSDT/metrics.yml @@ -0,0 +1,13 @@ +metrics: + RMSE: + good: 0.0 + bad: 0.2 + unit: eV + tooltip: "Mean RMSE of relative torsional energy profiles" + level_of_theory: CCSD(T)/CBS + MAE: + good: 0.0 + bad: 0.2 + unit: eV + tooltip: "Mean MAE of relative torsional energy profiles" + level_of_theory: CCSD(T)/CBS diff --git a/ml_peg/app/conformers/TorsionNet500CCSDT/app_TorsionNet500CCSDT.py b/ml_peg/app/conformers/TorsionNet500CCSDT/app_TorsionNet500CCSDT.py new file mode 100644 index 000000000..cc515ceee --- /dev/null +++ b/ml_peg/app/conformers/TorsionNet500CCSDT/app_TorsionNet500CCSDT.py @@ -0,0 +1,256 @@ +"""Run TorsionNet500CCSDT app.""" + +from __future__ import annotations + +from pathlib import Path + +from dash import Dash, Input, Output, State, callback +from dash.dcc import Store +from dash.html import Div, Iframe + +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, + plot_with_download_controls, +) +from ml_peg.app.utils.load import read_plot +from ml_peg.app.utils.weas import generate_weas_html +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 = "TorsionNet500CCSDT" +DOCS_URL = ( + "https://ddmms.github.io/ml-peg/user_guide/benchmarks/conformers.html" + "#torsionnet500ccsdt" +) +DATA_PATH = APP_ROOT / "data" / "conformers" / "TorsionNet500CCSDT" +INFO_PATH = DATA_PATH / "info.json" +ASSETS_DIR = "/assets/conformers/TorsionNet500CCSDT" + + +# TODO: if another benchmark needs this same lazy-load-on-click pattern, move +# this into ml_peg/app/utils/build_callbacks.py alongside plot_from_scatter/ +# struct_from_scatter instead of duplicating it here. +def _register_curve_callback( + scatter_id: str, + plot_id: str, + curve_dir: Path, + labels: list[str], + struct_plot_id: str, +) -> None: + """ + Attach callbacks that show a torsion curve and structure on point clicks. + + Unlike `plot_from_scatter`, this reads the clicked fragment's curve JSON from + disk on demand instead of requiring every curve to be pre-loaded into a + `plots_list` up front - with up to 500 fragments per model, pre-loading all of + them for every model would be needlessly expensive when only one is ever shown + at a time. + + A second callback shows the 3D structure at the torsion angle of whichever + point on the curve is clicked, using WEAS trajectory mode. The curve's + fragment label is tracked in a `Store` so the structure callback (triggered by + clicks on the curve, whose id is fixed and reused across fragments) knows which + fragment's trajectory file to load. + + Parameters + ---------- + scatter_id + ID of the per-model fragment-RMSE scatter plot. + plot_id + ID of the shared placeholder Div where curves are rendered. + curve_dir + Directory containing this model's per-fragment torsion curve JSON files. + labels + Fragment labels, in the same order as the scatter's points. + struct_plot_id + ID of the shared placeholder Div where structures are rendered. + """ + label_store_id = f"{scatter_id}-curve-label" + + @callback( + Output(plot_id, "children", allow_duplicate=True), + Output(label_store_id, "data", allow_duplicate=True), + Output(struct_plot_id, "children", allow_duplicate=True), + Input(scatter_id, "clickData"), + prevent_initial_call="initial_duplicate", + ) + def show_curve(click_data): + """ + Register callback to show a torsion curve when a scatter point is clicked. + + Also clears any structure shown from a previously displayed curve, so a + stale structure from a different fragment/model doesn't linger on screen + until the new curve is itself clicked. + + Parameters + ---------- + click_data + Clicked data point in the fragment-RMSE scatter plot. + + Returns + ------- + tuple[Div, str | None, Div] + Torsion curve plot on scatter click, the clicked fragment's label, and + an empty Div to clear any previously displayed structure. + """ + if not click_data: + return Div(), None, Div() + + idx = click_data["points"][0]["pointNumber"] + if idx < 0 or idx >= len(labels): + return Div(), None, Div() + + label = labels[idx] + curve_path = curve_dir / f"{label}.json" + graph = read_plot(curve_path, id=f"{scatter_id}-curve") + return plot_with_download_controls(graph), label, Div() + + @callback( + Output(struct_plot_id, "children", allow_duplicate=True), + Input(f"{scatter_id}-curve", "clickData"), + State(label_store_id, "data"), + prevent_initial_call="initial_duplicate", + ) + def show_struct(click_data, label): + """ + Register callback to show a structure when a point on the curve is clicked. + + Parameters + ---------- + click_data + Clicked data point in the torsion curve plot. + label + Fragment label of the currently displayed curve. + + Returns + ------- + Div + Structure at the clicked dihedral angle, in WEAS trajectory mode. + """ + if not click_data or not label: + return Div() + + idx = click_data["points"][0]["pointNumber"] + traj_path = f"{ASSETS_DIR}/torsion_trajectories/{label}.xyz" + + return Div( + Iframe( + srcDoc=generate_weas_html(traj_path, mode="traj", index=idx), + style={ + "height": "550px", + "width": "100%", + "border": "1px solid #ddd", + "borderRadius": "5px", + }, + ) + ) + + +class TorsionNet500CCSDTApp(BaseApp): + """TorsionNet500CCSDT benchmark app layout and callbacks.""" + + def register_callbacks(self) -> None: + """Register callbacks to app.""" + # Build an RMSE and an MAE scatter plot per model (fragment index vs + # that fragment's metric value), and record fragment labels in + # scatter-point order so a later click on a point can be mapped back + # to its curve file. + scatter_cell_to_plot: dict[str, dict[str, Div]] = {} + fragment_labels: dict[str, list[str]] = {} + + for model_name in MODELS: + rmse_scatter_path = DATA_PATH / model_name / "fragment_rmse_scatter.json" + mae_scatter_path = DATA_PATH / model_name / "fragment_mae_scatter.json" + curve_dir = DATA_PATH / model_name / "torsion_curves" + + if ( + not rmse_scatter_path.exists() + or not mae_scatter_path.exists() + or not curve_dir.exists() + ): + continue + + scatter_cell_to_plot[model_name] = { + "RMSE": read_plot( + rmse_scatter_path, id=f"{BENCHMARK_NAME}-{model_name}-rmse-figure" + ), + "MAE": read_plot( + mae_scatter_path, id=f"{BENCHMARK_NAME}-{model_name}-mae-figure" + ), + } + fragment_labels[model_name] = [ + curve_file.stem for curve_file in sorted(curve_dir.glob("*.json")) + ] + + # Clicking a model's RMSE or MAE cell shows that model's corresponding + # fragment scatter. + plot_from_table_cell( + table_id=self.table_id, + plot_id=f"{BENCHMARK_NAME}-scatter-placeholder", + cell_to_plot=scatter_cell_to_plot, + ) + + # Clicking a point on either scatter shows that fragment's curve - the + # curve itself doesn't depend on which metric was clicked through, so + # both scatters for a model share the same labels/curve_dir. + for model_name, labels in fragment_labels.items(): + curve_dir = DATA_PATH / model_name / "torsion_curves" + for metric in ("rmse", "mae"): + _register_curve_callback( + scatter_id=f"{BENCHMARK_NAME}-{model_name}-{metric}-figure", + plot_id=f"{BENCHMARK_NAME}-curve-placeholder", + curve_dir=curve_dir, + labels=labels, + struct_plot_id=f"{BENCHMARK_NAME}-struct-placeholder", + ) + + +def get_app() -> TorsionNet500CCSDTApp: + """ + Get TorsionNet500CCSDT benchmark app layout and callback registration. + + Returns + ------- + TorsionNet500CCSDTApp + Benchmark layout and callback registration. + """ + return TorsionNet500CCSDTApp( + name=BENCHMARK_NAME, + description=( + "Accuracy of predicted torsional energy profiles for 500 molecular " + "fragments, benchmarked against CCSD(T)/CBS reference scans." + ), + docs_url=DOCS_URL, + table_path=DATA_PATH / "torsionnet500ccsdt_metrics_table.json", + extra_components=[ + Div(id=f"{BENCHMARK_NAME}-scatter-placeholder"), + Div(id=f"{BENCHMARK_NAME}-curve-placeholder"), + Div(id=f"{BENCHMARK_NAME}-struct-placeholder"), + # One Store per (model, metric) scatter, remembering which fragment's + # curve is currently displayed, so a click on the curve (whose id is + # fixed and reused across fragments) can be mapped back to the right + # trajectory file. + *[ + Store(id=f"{BENCHMARK_NAME}-{model_name}-{metric}-figure-curve-label") + for model_name in MODELS + for metric in ("rmse", "mae") + ], + ], + info_path=INFO_PATH, + ) + + +if __name__ == "__main__": + # Create Dash app + full_app = Dash(__name__, assets_folder=DATA_PATH.parent.parent) + + # Construct layout and register callbacks + torsionnet500ccsdt_app = get_app() + full_app.layout = torsionnet500ccsdt_app.layout + torsionnet500ccsdt_app.register_callbacks() + + # Run app + full_app.run(port=8072, debug=True) diff --git a/ml_peg/calcs/conformers/TorsionNet500CCSDT/calc_TorsionNet500CCSDT.py b/ml_peg/calcs/conformers/TorsionNet500CCSDT/calc_TorsionNet500CCSDT.py new file mode 100644 index 000000000..f6fed08b6 --- /dev/null +++ b/ml_peg/calcs/conformers/TorsionNet500CCSDT/calc_TorsionNet500CCSDT.py @@ -0,0 +1,73 @@ +"""Run calculations for the TorsionNet500CCSDT benchmark.""" + +from __future__ import annotations + +from pathlib import Path +from typing import Any +from warnings import warn + +from ase import units +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" + +# Unit conversion +HARTREE_TO_EV = units.Hartree + + +@pytest.mark.parametrize("mlip", MODELS.items()) +def test_torsionnet500ccsdt(mlip: tuple[str, Any]) -> None: + """ + Run the TorsionNet500CCSDT benchmark. + + Parameters + ---------- + mlip + Name of model and model to get calculator. + """ + model_name, model = mlip + try: + calc = model.get_calculator(precision="high") + except ModuleNotFoundError as exc: + pytest.skip(f"Skipping {model_name}: {exc}") + + data_path = ( + download_s3_data( + filename="TorsionNet500CCSDT.zip", + key="inputs/conformers/TorsionNet500CCSDT/TorsionNet500CCSDT.zip", + ) + / "TorsionNet500CCSDT" + ) + + xyz_files = sorted(data_path.glob("*.xyz")) + + for o in tqdm(xyz_files): + atoms = read(o, ":") + + for a in atoms: + # Reference energy from the dataset + a.info["ref_energy"] = a.info["E_CCSDT"] * HARTREE_TO_EV + a.info["charge"] = int(a.info.get("charge", 0)) + a.info["spin"] = int(a.info.get("spin", 1)) + + # Model energy + a.calc = calc + + try: + a.info["model_energy"] = a.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy for {o.name}: {exc}", stacklevel=2) + a.info["model_energy"] = np.nan + + write_dir = OUT_PATH / model_name + write_dir.mkdir(parents=True, exist_ok=True) + write(write_dir / o.name, atoms)