diff --git a/.gitignore b/.gitignore index 3f5f8c439..bba1f72af 100644 --- a/.gitignore +++ b/.gitignore @@ -15,6 +15,8 @@ !dvc.yaml *.pdf *.hdf5 +*.mda +*.mdt *.traj *.ipynb *.jpg diff --git a/docs/source/user_guide/benchmarks/index.rst b/docs/source/user_guide/benchmarks/index.rst index ea0dcafb2..7c4e59526 100644 --- a/docs/source/user_guide/benchmarks/index.rst +++ b/docs/source/user_guide/benchmarks/index.rst @@ -19,3 +19,4 @@ Benchmarks conformers molecular_dynamics defect + superacids diff --git a/docs/source/user_guide/benchmarks/superacids.rst b/docs/source/user_guide/benchmarks/superacids.rst new file mode 100644 index 000000000..3268b12f4 --- /dev/null +++ b/docs/source/user_guide/benchmarks/superacids.rst @@ -0,0 +1,98 @@ +========== +Superacids +========== + +HF/SbF5 mixture densities +========================= + +Summary +------- + +Performance in predicting the liquid density of HF/SbF5 mixtures, the archetypal +superacid system, at three compositions: pure HF, a 10 mol % SbF5 mixture, and pure +SbF5. Each system consists of about 200 atoms, simulated at 288.65 K and 1 atm. + +Metrics +------- + +1. MAPE + +For each composition, the density is calculated from the average volume of an NPT +molecular dynamics run of 100 ps. The first half of the simulation is discarded as +equilibration. The mean absolute percentage error over the three compositions is +compared to the reference densities, obtained from experiment. + +Computational cost +------------------ + +High: about 8 GPU hours per model for all three compositions, for a model that +includes a dispersion correction. + +Data availability +----------------- + +Input structures: + +* Generated with Packmol, available from the ML-PEG S3 bucket + +Reference data: + +* Shair and Schurig, Vapor-Liquid Equilibrium of Antimony Pentafluoride-Hydrogen + Fluoride. Ind. Eng. Chem. 43, 1624 (1951). https://doi.org/10.1021/ie50499a042 +* Experimental + + +HF structure factor +=================== + +Summary +------- + +Performance in predicting the total neutron structure factor of liquid HF, which probes +the hydrogen-bonded chain structure of the liquid. A system of 100 HF molecules is +simulated at 296 K and 1.2 bar for 50 ps of NPT molecular dynamics. + +Since the experimental reference is measured on the deuterated liquid, all H are +transmuted to D before the structure factor is computed. S(q) is obtained with MDANSE, +as the Fourier transform of the pair distribution function weighted by coherent +scattering lengths, using the second half of the trajectory sampled every other frame. +The real-space cutoff is kept the same for every model, so that the transform is +truncated identically. + +Metrics +------- + +1. S(q) R-factor + +The relative deviation of the calculated structure factor from experiment, +``sum|S_exp - S_calc| / sum|S_exp|``, evaluated from the first experimental point to +4 1/A. Calculated and experimental structure factors are computed on the same grid of +scattering vectors. + +2. First peak position error + +The absolute error in the position of the maximum of S(q), located over the same range +as the R-factor. The position is read directly off the grid of scattering vectors, +whose spacing is 0.05 1/A. + +The "good" threshold of the R-factor is the statistical noise of the protocol itself, +measured by splitting the production window of one model into two halves. + +Computational cost +------------------ + +High: about 6 GPU hours per model, for a model that includes a dispersion correction. + +Data availability +----------------- + +Input structures: + +* Generated with Packmol, available from the ML-PEG S3 bucket + +Reference data: + +* McLain, Benmore, Siewenie, Urquidi and Turner, On the Structure of Liquid Hydrogen + Fluoride. Angew. Chem. Int. Ed. 43, 1952 (2004). + https://doi.org/10.1002/anie.200353289 +* Experimental diff --git a/ml_peg/analysis/superacids/HF_SbF5_density/analyse_HF_SbF5_density.py b/ml_peg/analysis/superacids/HF_SbF5_density/analyse_HF_SbF5_density.py new file mode 100644 index 000000000..075b065c7 --- /dev/null +++ b/ml_peg/analysis/superacids/HF_SbF5_density/analyse_HF_SbF5_density.py @@ -0,0 +1,312 @@ +"""Analyse HF/SbF5 density benchmark.""" + +from __future__ import annotations + +from pathlib import Path +from warnings import warn + +from ase import units +from ase.io import read +import numpy as np +import pytest + +from ml_peg.analysis.utils.decorators import build_table, plot_parity, plot_scatter +from ml_peg.analysis.utils.utils import ( + build_dispersion_name_map, + get_struct_info, + load_metrics_config, +) +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) +DISPERSION_MODEL_NAMES = build_dispersion_name_map(MODELS) +CALC_PATH = CALCS_ROOT / "superacids" / "HF_SbF5_density" / "outputs" +OUT_PATH = APP_ROOT / "data" / "superacids" / "HF_SbF5_density" + +METRICS_CONFIG_PATH = Path(__file__).with_name("metrics.yml") +DEFAULT_THRESHOLDS, DEFAULT_TOOLTIPS, DEFAULT_WEIGHTS = load_metrics_config( + METRICS_CONFIG_PATH +) + +# Experimental reference densities, from Shair and Schurig, +# Ind. Eng. Chem. 43, 1624 (1951), https://doi.org/10.1021/ie50499a042 +REF_DENSITIES = { + "X_0": 0.989, + "X_10": 1.677, + "X_100": 3.141, +} + +SYSTEMS = sorted(REF_DENSITIES) + +# amu to g conversion factor +AMU_TO_G = 1000 / units.kg +A3_TO_CM3 = 1e-24 + +# Minimum number of production samples for a density to be considered valid +MIN_SAMPLES = 2 + +# MD timestep in fs and number of steps run, to convert steps of volume.dat +# into a time axis. These match the values set in the calculation. +DT_FS = 0.5 +N_NPT_STEPS = 200000 +FS_TO_PS = 1e-3 + +# Composition of each system, as the mol % of SbF5 in the mixture +SYSTEM_COMPOSITIONS = { + "X_0": 0, + "X_10": 10, + "X_100": 100, +} + + +def compute_density(traj_path: Path, volume_path: Path) -> float: + """ + Compute average density from volume.dat and atomic masses. + + Parameters + ---------- + traj_path + Path to trajectory of this system (to get atomic masses). + volume_path + Path to volume.dat file (columns: step, volume_A3). + + Returns + ------- + float + Average density in g/cm³, or NaN if there are too few samples. + """ + # Read total mass from the first frame, so that runs in progress can be analysed + atoms = read(traj_path, index=0) + total_mass_amu = np.sum(atoms.get_masses()) + + # Read volume time series, skip header + data = np.loadtxt(volume_path, comments="#", ndmin=2) + # Take second half as production (discard equilibration) + production = data[len(data) // 2 :, 1] if data.size else np.array([]) + if production.size < MIN_SAMPLES: + return np.nan + avg_volume = np.mean(production) + + return (total_mass_amu * AMU_TO_G) / (avg_volume * A3_TO_CM3) + + +def compute_density_series( + traj_path: Path, volume_path: Path +) -> tuple[list[float], list[float]]: + """ + Compute the instantaneous density over the whole trajectory. + + Note that the average of this series is not exactly the density reported in + the metrics table, which is computed from the average volume. + + Parameters + ---------- + traj_path + Path to trajectory of this system (to get atomic masses). + volume_path + Path to volume.dat file (columns: step, volume_A3). + + Returns + ------- + tuple[list[float], list[float]] + Time in ps, and density in g/cm³. + """ + atoms = read(traj_path, index=0) + total_mass_amu = np.sum(atoms.get_masses()) + + data = np.loadtxt(volume_path, comments="#", ndmin=2) + if not data.size: + return [], [] + + time_ps = data[:, 0] * DT_FS * FS_TO_PS + densities = (total_mass_amu * AMU_TO_G) / (data[:, 1] * A3_TO_CM3) + + return time_ps.tolist(), densities.tolist() + + +def plot_density_series(system: str) -> None: + """ + Plot the density of all models against time for one system. + + Parameters + ---------- + system + System identifier (X_0, X_10, X_100). + """ + composition = SYSTEM_COMPOSITIONS[system] + total_time_ps = N_NPT_STEPS * DT_FS * FS_TO_PS + + @plot_scatter( + filename=OUT_PATH / f"figure_density_time_{system}.json", + title=f"HF/SbF5 mixture with x = {composition}% SbF5", + x_label="Time / ps", + y_label="Density / g/cm³", + show_line=True, + show_markers=False, + hlines={"Target": REF_DENSITIES[system]}, + highlight_range={"Production": [total_time_ps / 2, total_time_ps]}, + ) + def density_series() -> dict[str, list]: + """ + Get the density of all models against time. + + Returns + ------- + dict[str, list] + Times and densities for all models with a trajectory. + """ + results = {} + + for model_name in MODELS: + system_dir = CALC_PATH / model_name / system + traj_path = system_dir / f"{system}.traj" + volume_path = system_dir / "volume.dat" + + # Missing systems are left out of the plot + if not traj_path.exists() or not volume_path.exists(): + continue + + try: + results[model_name] = list( + compute_density_series(traj_path, volume_path) + ) + except Exception as exc: + warn( + f"Error computing density series for {model_name} {system}: {exc}", + stacklevel=2, + ) + + return results + + density_series() + + +@pytest.fixture +def density_series() -> None: + """Plot the density of all models against time, for every system.""" + for system in SYSTEMS: + plot_density_series(system) + + +@pytest.fixture +@plot_parity( + filename=OUT_PATH / "figure_density.json", + title="HF/SbF5 Mixture Densities", + x_label="Predicted density / g/cm³", + y_label="Experimental density / g/cm³", + hoverdata={ + "System": SYSTEMS, + }, +) +def densities() -> dict[str, list]: + """ + Get predicted and reference densities for all systems. + + Returns + ------- + dict[str, list] + Dictionary of reference and predicted densities. + """ + results = {"ref": [REF_DENSITIES[system] for system in SYSTEMS]} | { + mlip: [np.nan] * len(SYSTEMS) for mlip in MODELS + } + + for model_name in MODELS: + for index, system in enumerate(SYSTEMS): + system_dir = CALC_PATH / model_name / system + traj_path = system_dir / f"{system}.traj" + volume_path = system_dir / "volume.dat" + + # Missing systems are left as NaN, to keep systems aligned with `ref` + if not traj_path.exists() or not volume_path.exists(): + continue + + try: + results[model_name][index] = compute_density(traj_path, volume_path) + except Exception as exc: + warn( + f"Error computing density for {model_name} {system}: {exc}", + stacklevel=2, + ) + + return results + + +@pytest.fixture +def density_errors(densities) -> dict[str, float]: + """ + Get mean absolute percentage error for densities. + + Parameters + ---------- + densities + Dictionary of reference and predicted densities. + + Returns + ------- + dict[str, float] + Dictionary of density MAPE for all models. + """ + results = {} + refs = np.array(densities["ref"], dtype=float) + + for model_name in MODELS: + preds = np.array(densities[model_name], dtype=float) + # Models missing any system are scored as None, as in `mae` + if np.isnan(np.sum(preds)): + results[model_name] = None + else: + results[model_name] = float(np.mean(np.abs(preds - refs) / refs) * 100) + + return results + + +@pytest.fixture +@build_table( + filename=OUT_PATH / "hf_sbf5_density_metrics_table.json", + metric_tooltips=DEFAULT_TOOLTIPS, + thresholds=DEFAULT_THRESHOLDS, + mlip_name_map=DISPERSION_MODEL_NAMES, +) +def metrics(density_errors: dict[str, float]) -> dict[str, dict]: + """ + Get all HF/SbF5 density metrics. + + Parameters + ---------- + density_errors + Mean absolute errors for all systems. + + Returns + ------- + dict[str, dict] + Metric names and values for all models. + """ + return { + "MAPE": density_errors, + } + + +def test_hf_sbf5_density(metrics: dict[str, dict], density_series: None) -> None: + """ + Run HF/SbF5 density test. + + Parameters + ---------- + metrics + All HF/SbF5 density metrics. + density_series + Density against time plots for all systems. + """ + # Elemental info for filtering, from the mock calculation + get_struct_info( + calc_path=CALC_PATH, + glob_pattern="*/*.traj", + index=0, + include_dirs=True, + write_structs=False, + out_path=OUT_PATH, + ) diff --git a/ml_peg/analysis/superacids/HF_SbF5_density/metrics.yml b/ml_peg/analysis/superacids/HF_SbF5_density/metrics.yml new file mode 100644 index 000000000..da45319c5 --- /dev/null +++ b/ml_peg/analysis/superacids/HF_SbF5_density/metrics.yml @@ -0,0 +1,7 @@ +metrics: + MAPE: + good: 3.0 + bad: 20.0 + unit: "%" + tooltip: "Mean Absolute Percentage Error in liquid density vs experiment. Experimental reference: Shair and Schurig, Ind. Eng. Chem. 43, 1624 (1951)" + level_of_theory: Experimental diff --git a/ml_peg/analysis/superacids/HF_structure/analyse_HF_structure.py b/ml_peg/analysis/superacids/HF_structure/analyse_HF_structure.py new file mode 100644 index 000000000..4694b4a98 --- /dev/null +++ b/ml_peg/analysis/superacids/HF_structure/analyse_HF_structure.py @@ -0,0 +1,391 @@ +"""Analyse HF neutron structure factor benchmark.""" + +from __future__ import annotations + +import json +from pathlib import Path +from warnings import warn + +from ase.io import iread +import h5py +from MDANSE.Framework.Converters.Converter import Converter +from MDANSE.Framework.Jobs.IJob import IJob +from MDANSE.MolecularDynamics.Trajectory import Trajectory +import numpy as np +import pytest + +from ml_peg.analysis.utils.decorators import build_table, plot_scatter +from ml_peg.analysis.utils.utils import ( + build_dispersion_name_map, + load_metrics_config, + write_struct_info, +) +from ml_peg.app import APP_ROOT +from ml_peg.calcs import CALCS_ROOT +from ml_peg.calcs.utils.utils import download_s3_data +from ml_peg.models import current_models +from ml_peg.models.get_models import get_model_names + +MODELS = get_model_names(current_models) +DISPERSION_MODEL_NAMES = build_dispersion_name_map(MODELS) +CALC_PATH = CALCS_ROOT / "superacids" / "HF_structure" / "outputs" +OUT_PATH = APP_ROOT / "data" / "superacids" / "HF_structure" + +METRICS_CONFIG_PATH = Path(__file__).with_name("metrics.yml") +DEFAULT_THRESHOLDS, DEFAULT_TOOLTIPS, DEFAULT_WEIGHTS = load_metrics_config( + METRICS_CONFIG_PATH +) + +# Time between trajectory frames, from the MD timestep and dump frequency +FRAME_SPACING_FS = 50.0 + +# Second half of the trajectory is production, sampled every STRIDE frames +PRODUCTION_FRACTION = 0.5 +STRIDE = 2 + +# Minimum number of production frames for a structure factor to be computed +MIN_FRAMES = 10 + +# Real space grid for the pair distribution function, in nm. +# The cutoff is kept fixed for every model so that the Fourier transform is +# truncated identically, which would not be the case if it were derived from +# each model's own (NPT) cell. +R_MAX_NM = 0.6 +R_STEP_NM = 0.005 + +# Reciprocal space grid, in 1/nm. Matches the experimental grid once converted +# to 1/A, so calculated and experimental points coincide. +Q_MIN_INV_NM = 7.5 +Q_MAX_INV_NM = 100.0 +Q_STEP_INV_NM = 0.5 + +INV_NM_TO_INV_ANG = 0.1 +NM_TO_ANG = 10.0 + +# Upper bound of the range over which models are scored against experiment, in +# 1/A. The lower bound is the first experimental point. The first peak of S(q) +# is located within the same range. +SCORE_Q_MAX = 4.0 + + +def load_reference_sq() -> tuple[np.ndarray, np.ndarray]: + """ + Load the experimental neutron structure factor. + + Returns + ------- + tuple[numpy.ndarray, numpy.ndarray] + Scattering vector in 1/A, and structure factor. + """ + hf_structure_dir = ( + download_s3_data( + key="inputs/superacids/HF_structure/HF_structure.zip", + filename="HF_structure.zip", + ) + / "HF_structure" + ) + + data = np.loadtxt(hf_structure_dir / "SQ_EXP.dat") + + return data[:, 0], data[:, 1] + + +def compute_sq(traj_path: Path, work_dir: Path) -> tuple[np.ndarray, np.ndarray]: + """ + Compute the total neutron structure factor of a trajectory. + + All H are transmuted to D, since the experimental reference is measured on + the deuterated liquid. S(q) is obtained from the Fourier transform of the + pair distribution function, weighted by coherent scattering lengths. + + Parameters + ---------- + traj_path + Path to the NPT trajectory of this model. + work_dir + Directory to write the converted trajectory and MDANSE output to. + + Returns + ------- + tuple[numpy.ndarray, numpy.ndarray] + Scattering vector in 1/A, and total structure factor. + """ + work_dir.mkdir(parents=True, exist_ok=True) + mdt_path = work_dir / "traj.mdt" + ssf_prefix = work_dir / "ssf" + + # Cell fluctuates in NPT, so the cutoff must fit within every frame analysed + cell_lengths = [ + float(atoms.cell.lengths().min()) for atoms in iread(traj_path, index=":") + ] + n_frames = len(cell_lengths) + first_frame = int(n_frames * PRODUCTION_FRACTION) + + n_production = len(range(first_frame, n_frames, STRIDE)) + if n_production < MIN_FRAMES: + raise ValueError(f"Only {n_production} production frames in {traj_path}") + + min_length = min(cell_lengths[first_frame:]) + if 2 * R_MAX_NM * NM_TO_ANG > min_length: + raise ValueError( + f"Cutoff {R_MAX_NM * NM_TO_ANG} A exceeds half the smallest cell " + f"({min_length:.3f} A) in {traj_path}" + ) + + converter = Converter.create("ASE") + converter.run( + { + "trajectory_file": str(traj_path), + "atom_aliases": "{}", + "time_step": FRAME_SPACING_FS, + "time_unit": "fs", + "n_steps": 0, + "fold": True, + "output_files": (str(work_dir / "traj"), 64, 128, "none", "no logs"), + }, + status=False, + ) + + # Transmute H to D, matching the deuterated experimental sample + trajectory = Trajectory(str(mdt_path)) + transmutation = { + str(index): "H2" + for index, atom_type in enumerate(trajectory.atom_types) + if atom_type == "H" + } + trajectory.close() + + ssf = IJob.create("StaticStructureFactor") + ssf.run( + { + "trajectory": str(mdt_path), + "frames": (first_frame, n_frames, STRIDE), + "r_values": (0.0, R_MAX_NM, R_STEP_NM), + "q_values": (Q_MIN_INV_NM, Q_MAX_INV_NM, Q_STEP_INV_NM), + "atom_selection": "{}", + "atom_transmutation": json.dumps(transmutation), + "grouping_level": "atom", + "weights": "b_coherent", + "output_files": (str(ssf_prefix), ["MDAFormat"], "no logs"), + "running_mode": ("single-core", 1), + }, + status=False, + ) + + with h5py.File(f"{ssf_prefix}.mda", "r") as output: + q = output["ssf/axes/q"][:] * INV_NM_TO_INV_ANG + sq = output["ssf/total"][:] + + return q, sq + + +def compute_r_factor( + q_ref: np.ndarray, + sq_ref: np.ndarray, + q_calc: np.ndarray, + sq_calc: np.ndarray, +) -> float: + """ + Compute the R-factor between calculated and experimental S(q). + + The R-factor is ``sum|S_exp - S_calc| / sum|S_exp|``, evaluated on the + experimental grid up to `SCORE_Q_MAX`. + + Parameters + ---------- + q_ref + Experimental scattering vector in 1/A. + sq_ref + Experimental structure factor. + q_calc + Calculated scattering vector in 1/A. + sq_calc + Calculated structure factor. + + Returns + ------- + float + R-factor, or NaN if the curves do not overlap. + """ + mask = (q_ref >= q_calc.min()) & (q_ref <= min(SCORE_Q_MAX, q_calc.max())) + if not mask.any(): + return np.nan + + # Identity when the calculated and experimental grids coincide + sq_interp = np.interp(q_ref[mask], q_calc, sq_calc) + + return float( + np.sum(np.abs(sq_ref[mask] - sq_interp)) / np.sum(np.abs(sq_ref[mask])) + ) + + +def first_peak_position( + q: np.ndarray, sq: np.ndarray, q_min: float, q_max: float +) -> float: + """ + Get the position of the maximum of S(q) over the scored range. + + The position is read directly off the grid of scattering vectors, which is + fine enough that interpolating between points is not necessary. + + Parameters + ---------- + q + Scattering vector in 1/A. + sq + Structure factor. + q_min + Lower bound of the range searched, in 1/A. + q_max + Upper bound of the range searched, in 1/A. + + Returns + ------- + float + Position of the peak in 1/A, or NaN if the range contains no points. + """ + (indices,) = np.nonzero((q >= q_min) & (q <= q_max) & ~np.isnan(sq)) + if indices.size == 0: + return np.nan + + return float(q[indices[np.argmax(sq[indices])]]) + + +@pytest.fixture +@plot_scatter( + filename=OUT_PATH / "figure_sq.json", + title="HF Neutron Structure Factor", + x_label="q / 1/A", + y_label="S(q)", + show_line=True, + show_markers=False, + highlight_range={"Scored": [Q_MIN_INV_NM * INV_NM_TO_INV_ANG, SCORE_Q_MAX]}, +) +def sq_curves() -> dict[str, list]: + """ + Get experimental and predicted structure factors for all models. + + Returns + ------- + dict[str, list] + Scattering vectors and structure factors for the reference and all + models with a trajectory. + """ + q_ref, sq_ref = load_reference_sq() + results = {"ref": [q_ref.tolist(), sq_ref.tolist()]} + + for model_name in MODELS: + traj_path = CALC_PATH / model_name / "NPT.traj" + + # Missing models are left out of the plot, and scored as None + if not traj_path.exists(): + continue + + try: + q, sq = compute_sq(traj_path, CALC_PATH / model_name / "sq_mdanse") + except Exception as exc: + warn( + f"Error computing structure factor for {model_name}: {exc}", + stacklevel=2, + ) + continue + + results[model_name] = [q.tolist(), sq.tolist()] + + sq_out = OUT_PATH / model_name + sq_out.mkdir(parents=True, exist_ok=True) + np.savetxt( + sq_out / "sq.dat", + np.column_stack([q, sq]), + header="q/1/A S(q)", + fmt="%.6f", + ) + + return results + + +@pytest.fixture +def sq_errors(sq_curves: dict[str, list]) -> dict[str, dict]: + """ + Get structure factor errors for all models. + + Parameters + ---------- + sq_curves + Scattering vectors and structure factors for the reference and models. + + Returns + ------- + dict[str, dict] + R-factors and first peak position errors for all models. + """ + q_ref = np.array(sq_curves["ref"][0], dtype=float) + sq_ref = np.array(sq_curves["ref"][1], dtype=float) + + # Peaks are located over the same range as the R-factor is evaluated on + q_min = float(q_ref.min()) + peak_ref = first_peak_position(q_ref, sq_ref, q_min, SCORE_Q_MAX) + + r_factors = {} + peak_errors = {} + + for model_name in MODELS: + # Models without a structure factor are scored as None, as in `mae` + if model_name not in sq_curves: + r_factors[model_name] = None + peak_errors[model_name] = None + continue + + q = np.array(sq_curves[model_name][0], dtype=float) + sq = np.array(sq_curves[model_name][1], dtype=float) + + r_factors[model_name] = compute_r_factor(q_ref, sq_ref, q, sq) + peak_errors[model_name] = abs( + first_peak_position(q, sq, q_min, SCORE_Q_MAX) - peak_ref + ) + + return { + "S(q) R-factor": r_factors, + "First Peak Position Error": peak_errors, + } + + +@pytest.fixture +@build_table( + filename=OUT_PATH / "hf_structure_metrics_table.json", + metric_tooltips=DEFAULT_TOOLTIPS, + thresholds=DEFAULT_THRESHOLDS, + mlip_name_map=DISPERSION_MODEL_NAMES, + weights=DEFAULT_WEIGHTS, +) +def metrics(sq_errors: dict[str, dict]) -> dict[str, dict]: + """ + Get all HF structure factor metrics. + + Parameters + ---------- + sq_errors + R-factors and first peak position errors for all models. + + Returns + ------- + dict[str, dict] + Metric names and values for all models. + """ + return sq_errors + + +def test_hf_structure(metrics: dict[str, dict]) -> None: + """ + Run HF structure factor test. + + Parameters + ---------- + metrics + All HF structure factor metrics. + """ + write_struct_info( + data_path=CALC_PATH / "mock" / "minimised.xyz", + out_path=OUT_PATH, + index=0, + ) diff --git a/ml_peg/analysis/superacids/HF_structure/metrics.yml b/ml_peg/analysis/superacids/HF_structure/metrics.yml new file mode 100644 index 000000000..80a74f833 --- /dev/null +++ b/ml_peg/analysis/superacids/HF_structure/metrics.yml @@ -0,0 +1,19 @@ +# The `good` threshold of the R-factor is the statistical noise of the protocol +# itself, measured by splitting the production window of one model in two halves +# (R-factor 0.022). The peak position is read directly off the S(q) grid, so its +# `good` threshold is one grid spacing, 0.05 1/A. +metrics: + S(q) R-factor: + good: 0.02 + bad: 0.5 + unit: null + tooltip: "sum|S_exp - S_calc| / sum|S_exp| for the total neutron structure factor of liquid DF, from the first experimental point to 4 1/A. Experimental reference: McLain et al., Angew. Chem. Int. Ed. 43, 1952 (2004)" + level_of_theory: Experimental + weight: 1.0 + First Peak Position Error: + good: 0.05 + bad: 0.8 + unit: "1/A" + tooltip: "Absolute error in the position of the maximum of S(q), located over the same range as the R-factor. Experimental reference: McLain et al., Angew. Chem. Int. Ed. 43, 1952 (2004)" + level_of_theory: Experimental + weight: 1.0 diff --git a/ml_peg/analysis/utils/decorators.py b/ml_peg/analysis/utils/decorators.py index 7df47f73d..b53d03754 100644 --- a/ml_peg/analysis/utils/decorators.py +++ b/ml_peg/analysis/utils/decorators.py @@ -498,6 +498,7 @@ def plot_scatter( horizontal_lines: list[float | dict[str, Any]] | None = None, filename: str = "scatter.json", highlight_range: dict = None, + hlines: dict[str, float] | None = None, ) -> Callable: """ Plot scatter plot of MLIP results. @@ -524,6 +525,9 @@ def plot_scatter( Filename to save plot as JSON. Default is "scatter.json". highlight_range Dictionary of rectangle title and x-axis endpoints. + hlines + Dictionary of label and y-axis value, drawn as dashed horizontal + reference lines. Default is `None`. Returns ------- @@ -599,19 +603,31 @@ def plot_scatter_wrapper(*args, **kwargs) -> dict[str, Any]: ) ) - colors = pc.qualitative.Plotly - - if highlight_range: - for i, (h_text, range) in enumerate(highlight_range.items()): - fig.add_vrect( - x0=range[0], - x1=range[1], - annotation_text=h_text, - annotation_position="top", - fillcolor=colors[i], - opacity=0.25, - line_width=0, - ) + colors = pc.qualitative.Plotly + + # Drawn once, rather than once per model, so that the shaded + # regions do not stack up and darken + if highlight_range: + for i, (h_text, range) in enumerate(highlight_range.items()): + fig.add_vrect( + x0=range[0], + x1=range[1], + annotation_text=h_text, + annotation_position="top", + fillcolor=colors[i], + opacity=0.25, + line_width=0, + ) + + if hlines: + for label, value in hlines.items(): + fig.add_hline( + y=value, + line_dash="dash", + line_color="black", + annotation_text=label, + annotation_position="top right", + ) fig.update_layout( title={"text": title}, diff --git a/ml_peg/app/superacids/HF_SbF5_density/app_HF_SbF5_density.py b/ml_peg/app/superacids/HF_SbF5_density/app_HF_SbF5_density.py new file mode 100644 index 000000000..b62a16f93 --- /dev/null +++ b/ml_peg/app/superacids/HF_SbF5_density/app_HF_SbF5_density.py @@ -0,0 +1,85 @@ +"""Run HF/SbF5 density app.""" + +from __future__ import annotations + +from dash import Dash +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 import current_models +from ml_peg.models.get_models import get_model_names + +# Get all models +MODELS = get_model_names(current_models) +BENCHMARK_NAME = "HF/SbF5 Mixture Densities" +DOCS_URL = ( + "https://ddmms.github.io/ml-peg/user_guide/benchmarks/" + "superacids.html#hf-sbf5-mixture-densities" +) +DATA_PATH = APP_ROOT / "data" / "superacids" / "HF_SbF5_density" + +# Systems simulated, labelled by the mol % of SbF5 in the mixture +SYSTEMS = ("X_0", "X_10", "X_100") + + +class HFSbF5DensityApp(BaseApp): + """HF/SbF5 density benchmark app layout and callbacks.""" + + def register_callbacks(self) -> None: + """Register callbacks to app.""" + scatter = read_plot( + DATA_PATH / "figure_density.json", + id=f"{BENCHMARK_NAME}-figure", + ) + + plot_from_table_column( + table_id=self.table_id, + plot_id=f"{BENCHMARK_NAME}-figure-placeholder", + column_to_plot={"MAPE": scatter}, + ) + + +def get_app() -> HFSbF5DensityApp: + """ + Get HF/SbF5 density benchmark app layout and callback registration. + + Returns + ------- + HFSbF5DensityApp + Benchmark layout and callback registration. + """ + return HFSbF5DensityApp( + name=BENCHMARK_NAME, + description=("Liquid densities of HF/SbF5 mixtures at varying compositions."), + docs_url=DOCS_URL, + table_path=DATA_PATH / "hf_sbf5_density_metrics_table.json", + info_path=DATA_PATH / "info.json", + extra_components=[ + Div(id=f"{BENCHMARK_NAME}-figure-placeholder"), + *( + read_plot( + DATA_PATH / f"figure_density_time_{system}.json", + id=f"{BENCHMARK_NAME}-{system}-figure-density-time", + ) + for system in SYSTEMS + ), + ], + ) + + +if __name__ == "__main__": + # Create Dash app + full_app = Dash(__name__, assets_folder=DATA_PATH.parent.parent) + + # Construct layout and register callbacks + app = get_app() + full_app.layout = app.layout + app.register_callbacks() + + # Run app + full_app.run(port=8056, debug=True) diff --git a/ml_peg/app/superacids/HF_structure/app_HF_structure.py b/ml_peg/app/superacids/HF_structure/app_HF_structure.py new file mode 100644 index 000000000..ceaf93a2f --- /dev/null +++ b/ml_peg/app/superacids/HF_structure/app_HF_structure.py @@ -0,0 +1,81 @@ +"""Run HF structure factor app.""" + +from __future__ import annotations + +from dash import Dash +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 import current_models +from ml_peg.models.get_models import get_model_names + +# Get all models +MODELS = get_model_names(current_models) +BENCHMARK_NAME = "HF Structure Factor" +DOCS_URL = ( + "https://ddmms.github.io/ml-peg/user_guide/benchmarks/" + "superacids.html#hf-structure-factor" +) +DATA_PATH = APP_ROOT / "data" / "superacids" / "HF_structure" + + +class HFStructureApp(BaseApp): + """HF structure factor benchmark app layout and callbacks.""" + + def register_callbacks(self) -> None: + """Register callbacks to app.""" + scatter = read_plot( + DATA_PATH / "figure_sq.json", + id=f"{BENCHMARK_NAME}-figure", + ) + + plot_from_table_column( + table_id=self.table_id, + plot_id=f"{BENCHMARK_NAME}-figure-placeholder", + column_to_plot={ + "S(q) R-factor": scatter, + "First Peak Position Error": scatter, + }, + ) + + +def get_app() -> HFStructureApp: + """ + Get HF structure factor benchmark app layout and callback registration. + + Returns + ------- + HFStructureApp + Benchmark layout and callback registration. + """ + return HFStructureApp( + name=BENCHMARK_NAME, + description=( + "Total neutron structure factor of liquid HF, from NPT molecular " + "dynamics with all H transmuted to D." + ), + docs_url=DOCS_URL, + table_path=DATA_PATH / "hf_structure_metrics_table.json", + info_path=DATA_PATH / "info.json", + extra_components=[ + Div(id=f"{BENCHMARK_NAME}-figure-placeholder"), + ], + ) + + +if __name__ == "__main__": + # Create Dash app + full_app = Dash(__name__, assets_folder=DATA_PATH.parent.parent) + + # Construct layout and register callbacks + app = get_app() + full_app.layout = app.layout + app.register_callbacks() + + # Run app + full_app.run(port=8057, debug=True) diff --git a/ml_peg/app/superacids/superacids.yml b/ml_peg/app/superacids/superacids.yml new file mode 100644 index 000000000..a709b1953 --- /dev/null +++ b/ml_peg/app/superacids/superacids.yml @@ -0,0 +1,3 @@ +title: Superacids +description: Structural and thermophysical properties of superacids and related liquids +weight: 0 diff --git a/ml_peg/calcs/superacids/HF_SbF5_density/calc_HF_SbF5_density.py b/ml_peg/calcs/superacids/HF_SbF5_density/calc_HF_SbF5_density.py new file mode 100644 index 000000000..eebfe42d7 --- /dev/null +++ b/ml_peg/calcs/superacids/HF_SbF5_density/calc_HF_SbF5_density.py @@ -0,0 +1,208 @@ +"""Run calculations for HF/SbF5 density tests.""" + +from __future__ import annotations + +from pathlib import Path +from typing import Any +from warnings import warn + +from ase import Atoms, units +from ase.io import Trajectory, read, write +from ase.md.logger import MDLogger +from ase.md.nose_hoover_chain import IsotropicMTKNPT +from ase.md.velocitydistribution import ( + MaxwellBoltzmannDistribution, + Stationary, + ZeroRotation, +) +from ase.optimize import FIRE +import pytest + +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" + +# Simulation parameters +TEMPERATURE_K = 288.65 +PRESSURE_ATM = 1 +DT_FS = 0.5 # DT in femtoseconds +N_MIN_STEPS = 300 # maximum minimization steps +N_NPT_STEPS = 200000 # NPT production steps +OUT_FREQ = 200 + +# Conversions +ATM_TO_GPA = 1.01325e-4 # 1 atm = 0.000101325 GPa +PRESSURE_AU = PRESSURE_ATM * ATM_TO_GPA * units.GPa +DT = 0.5 * units.fs +TDAMP = 100 * DT_FS * units.fs +PDAMP = 1000 * DT_FS * units.fs + +# Systems +SYSTEMS = ["X_0", "X_10", "X_100"] + + +def read_restart(traj_path: Path) -> tuple[Atoms | None, int]: + """ + Read the last frame of a previous run of this system, if there is one. + + Parameters + ---------- + traj_path + Path to the trajectory of the NPT run. + + Returns + ------- + tuple[Atoms | None, int] + Last frame written, and the step it was written at, or `(None, 0)` if + there is no trajectory to restart from. + """ + if not traj_path.exists(): + return None, 0 + + try: + traj = Trajectory(str(traj_path)) + atoms = traj[-1] + nsteps = (len(traj) - 1) * OUT_FREQ + except Exception as exc: + warn(f"Ignoring unreadable trajectory {traj_path}: {exc}", stacklevel=2) + return None, 0 + + return atoms, nsteps + + +@pytest.mark.very_slow +@pytest.mark.parametrize("mlip", MODELS.items(), ids=lambda x: x[0]) +@pytest.mark.parametrize("system", SYSTEMS) +def test_hf_sbf5_density(mlip: tuple[str, Any], system: str) -> None: + """ + Run HF/SbF5 mixture density test. + + Interrupted runs are resumed from the last frame of the trajectory, so + minimisation and the initial velocity distribution are only applied when + starting from scratch. + + Parameters + ---------- + mlip + Name of model and model to get calculator. + system + System identifier (X_0, X_10, X_100). + """ + model_name, model = mlip + calc = model.get_calculator(precision="low") + + # Add D3 calculator for this test + calc = model.add_d3_calculator(calc) + + write_dir = OUT_PATH / model_name / system + write_dir.mkdir(parents=True, exist_ok=True) + traj_path = write_dir / f"{system}.traj" + + atoms, nsteps = read_restart(traj_path) + restarting = atoms is not None + + if restarting: + print(f"Resuming {system} with model {model_name} from step {nsteps}") + else: + print(f"Simulating {system} with model {model_name}") + + # Download dataset + hf_sbf5_density_dir = ( + download_s3_data( + key="inputs/superacids/HF_SbF5_density/HF_SbF5_density.zip", + filename="HF_SbF5_density.zip", + ) + / "HF_SbF5_density" + ) + + atoms = read(hf_sbf5_density_dir / system / "start.xyz") + + atoms.calc = calc + + if not restarting: + # Minimization + opt = FIRE(atoms, logfile=str(write_dir / "opt.log")) + try: + opt.run(fmax=0.05, steps=N_MIN_STEPS) + except Exception as exc: + warn(f"Error minimising {system}: {exc}", stacklevel=2) + write(write_dir / "minimised.xyz", atoms) + + MaxwellBoltzmannDistribution(atoms, temperature_K=TEMPERATURE_K) + Stationary(atoms) + ZeroRotation(atoms) + + dyn = IsotropicMTKNPT( + atoms=atoms, + timestep=DT, + temperature_K=TEMPERATURE_K, + pressure_au=PRESSURE_AU, + tdamp=TDAMP, + pdamp=PDAMP, + ) + + dyn.nsteps = nsteps + + dyn.attach( + MDLogger( + dyn, + atoms, + str(write_dir / "md.log"), + header=not restarting, + mode="a" if restarting else "w", + ), + interval=OUT_FREQ, + ) + + traj_file = Trajectory(str(traj_path), "a" if restarting else "w", atoms) + vol_file = open(write_dir / "volume.dat", "a" if restarting else "w") + if not restarting: + vol_file.write("# step volume_A3\n") + + last_written = nsteps if restarting else -1 + + def write_frame(_dyn=dyn, _atoms=atoms) -> None: + """ + Append the current frame to the trajectory, and its volume to file. + + Parameters + ---------- + _dyn : IsotropicMTKNPT + The dynamics object. + _atoms : Atoms + The ASE atoms object. + """ + nonlocal last_written + + step = _dyn.nsteps + if step <= last_written: + # Resuming: this frame was already written by the previous run. + return + + traj_file.write() + vol_file.write(f"{step} {_atoms.get_volume():.6f}\n") + vol_file.flush() + last_written = step + + write_frame() # step 0 + dyn.attach(write_frame, interval=OUT_FREQ) + + # Run NPT + if nsteps < N_NPT_STEPS: + try: + dyn.run(N_NPT_STEPS - nsteps) + except Exception as exc: + warn(f"Error running MD for {system}: {exc}", stacklevel=2) + + vol_file.close() + traj_file.close() + + # Save final structure + atoms.info["system"] = system + write(write_dir / f"{system}.xyz", atoms) + + print(f" {system} done") diff --git a/ml_peg/calcs/superacids/HF_structure/calc_HF_structure.py b/ml_peg/calcs/superacids/HF_structure/calc_HF_structure.py new file mode 100644 index 000000000..97bfdabf2 --- /dev/null +++ b/ml_peg/calcs/superacids/HF_structure/calc_HF_structure.py @@ -0,0 +1,193 @@ +"""Run calculations for HF structure factor tests.""" + +from __future__ import annotations + +from pathlib import Path +from typing import Any +from warnings import warn + +from ase import Atoms, units +from ase.io import Trajectory, read, write +from ase.md.logger import MDLogger +from ase.md.nose_hoover_chain import IsotropicMTKNPT +from ase.md.velocitydistribution import ( + MaxwellBoltzmannDistribution, + Stationary, + ZeroRotation, +) +from ase.optimize import FIRE +import pytest + +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" + +# Simulation parameters +TEMPERATURE_K = 296 +PRESSURE_BAR = 1.2 +DT_FS = 0.5 # femtoseconds +N_MIN_STEPS = 300 # maximum minimization steps +N_NPT_STEPS = 100000 # NPT production steps +OUT_FREQ = 100 # trajectory dump frequency + +# Conversions +PRESSURE_AU = PRESSURE_BAR * units.bar +DT = DT_FS * units.fs +TDAMP = 100 * DT_FS * units.fs +PDAMP = 1000 * DT_FS * units.fs + + +def read_restart(traj_path: Path) -> tuple[Atoms | None, int]: + """ + Read the last frame of a previous run, if there is one. + + Parameters + ---------- + traj_path + Path to the restart trajectory of the NPT run. + + Returns + ------- + tuple[Atoms | None, int] + Last frame written, and the step it was written at, or `(None, 0)` if + there is no trajectory to restart from. + """ + if not traj_path.exists(): + return None, 0 + + try: + traj = Trajectory(str(traj_path)) + atoms = traj[-1] + nsteps = (len(traj) - 1) * OUT_FREQ + except Exception as exc: + warn(f"Ignoring unreadable trajectory {traj_path}: {exc}", stacklevel=2) + return None, 0 + + return atoms, nsteps + + +@pytest.mark.very_slow +@pytest.mark.parametrize("mlip", MODELS.items(), ids=lambda x: x[0]) +def test_hf_structure(mlip: tuple[str, Any]) -> None: + """ + Run HF structure factor NPT simulation. + + Parameters + ---------- + mlip + Name of model and model to get calculator. + """ + model_name, model = mlip + calc = model.get_calculator(precision="low") + + calc = model.add_d3_calculator(calc) + + write_dir = OUT_PATH / model_name + write_dir.mkdir(parents=True, exist_ok=True) + + # Restart trajectory, written alongside NPT.xyz for post-processing + restart_path = write_dir / "NPT.traj" + traj_path = write_dir / "NPT.xyz" + + atoms, nsteps = read_restart(restart_path) + restarting = atoms is not None + + if restarting: + print(f"Resuming HF structure with model {model_name} from step {nsteps}") + else: + print(f"Simulating HF structure with model {model_name}") + + # Download dataset + hf_structure_dir = ( + download_s3_data( + key="inputs/superacids/HF_structure/HF_structure.zip", + filename="HF_structure.zip", + ) + / "HF_structure" + ) + + atoms = read(hf_structure_dir / "start.xyz") + traj_path.unlink(missing_ok=True) + + atoms.calc = calc + + if not restarting: + # Minimization + opt = FIRE(atoms, logfile=str(write_dir / "opt.log")) + try: + opt.run(fmax=0.05, steps=N_MIN_STEPS) + except Exception as exc: + warn(f"Error minimising HF structure: {exc}", stacklevel=2) + write(write_dir / "minimised.xyz", atoms) + + MaxwellBoltzmannDistribution(atoms, temperature_K=TEMPERATURE_K) + Stationary(atoms) + ZeroRotation(atoms) + + dyn = IsotropicMTKNPT( + atoms=atoms, + timestep=DT, + temperature_K=TEMPERATURE_K, + pressure_au=PRESSURE_AU, + tdamp=TDAMP, + pdamp=PDAMP, + ) + + dyn.nsteps = nsteps + + dyn.attach( + MDLogger( + dyn, + atoms, + str(write_dir / "md.log"), + header=not restarting, + mode="a" if restarting else "w", + ), + interval=OUT_FREQ, + ) + + restart_file = Trajectory(str(restart_path), "a" if restarting else "w", atoms) + + last_written = nsteps if restarting else -1 + + def write_frame(_dyn=dyn, _atoms=atoms, _path=traj_path) -> None: + """ + Append the current frame to the restart and NPT trajectories. + + Parameters + ---------- + _dyn : IsotropicMTKNPT + The dynamics object. + _atoms : Atoms + The ASE atoms object. + _path : Path + Path to the NPT trajectory file. + """ + nonlocal last_written + + step = _dyn.nsteps + if step <= last_written: + # Resuming: this frame was already written by the previous run. + return + + restart_file.write() + write(_path, _atoms, append=True) + last_written = step + + write_frame() # step 0 + dyn.attach(write_frame, interval=OUT_FREQ) + + # Run NPT + if nsteps < N_NPT_STEPS: + try: + dyn.run(N_NPT_STEPS - nsteps) + except Exception as exc: + warn(f"Error running MD for HF structure: {exc}", stacklevel=2) + + restart_file.close() + + print(f" HF structure done ({model_name})") diff --git a/pyproject.toml b/pyproject.toml index 5d8a1ea1f..6d24d368a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -32,6 +32,7 @@ dependencies = [ "matcalc<0.5,>=0.4.7; python_version >= '3.11'", "matminer<0.10,>=0.9.3", "mdanalysis<3,>=2.9.0", + "mdanse<3,>=2.0.1", "mlipx<0.2,>=0.1.5", "molify<0.3,>=0.2.2", "openpyxl<4,>=3.1.5", diff --git a/tests/test_superacids.py b/tests/test_superacids.py new file mode 100644 index 000000000..3e32e5116 --- /dev/null +++ b/tests/test_superacids.py @@ -0,0 +1,129 @@ +"""Unit tests for superacids benchmark utilities.""" + +from __future__ import annotations + +from pathlib import Path + +from ase import Atoms +from ase.io import Trajectory +import numpy as np +import pytest + +from ml_peg.analysis.superacids.HF_SbF5_density.analyse_HF_SbF5_density import ( + compute_density, +) +from ml_peg.analysis.superacids.HF_structure.analyse_HF_structure import ( + compute_r_factor, +) +from ml_peg.calcs.superacids.HF_SbF5_density.calc_HF_SbF5_density import ( + OUT_FREQ, + read_restart, +) + + +def write_trajectory(traj_path: Path, separations: tuple[float, ...]) -> None: + """ + Write a trajectory of two atoms with a total mass of 100 amu. + + Parameters + ---------- + traj_path + Path to write the trajectory to. + separations + Separation between the two atoms in each frame, in Angstrom. + """ + with Trajectory(str(traj_path), "w") as traj: + for separation in separations: + atoms = Atoms( + "H2", + positions=[[0, 0, 0], [separation, 0, 0]], + cell=[10, 10, 10], + pbc=True, + ) + atoms.set_masses([50.0, 50.0]) + traj.write(atoms) + + +def test_read_restart_without_trajectory(tmp_path: Path) -> None: + """Test nothing is restarted from if the trajectory does not exist.""" + atoms, nsteps = read_restart(tmp_path / "does_not_exist.traj") + + assert atoms is None + assert nsteps == 0 + + +def test_read_restart_returns_last_frame(tmp_path: Path) -> None: + """Test the last frame of the trajectory, and its step, are returned.""" + traj_path = tmp_path / "X_0.traj" + write_trajectory(traj_path, (1.0, 2.0, 3.0)) + + atoms, nsteps = read_restart(traj_path) + + # Frames are written at steps 0, OUT_FREQ and 2 * OUT_FREQ + assert nsteps == 2 * OUT_FREQ + assert atoms.get_distance(0, 1) == pytest.approx(3.0) + + +def test_read_restart_ignores_corrupt_trajectory(tmp_path: Path) -> None: + """Test an unreadable trajectory is warned about, rather than raising.""" + traj_path = tmp_path / "X_0.traj" + traj_path.write_bytes(b"not a trajectory") + + with pytest.warns(UserWarning, match="unreadable trajectory"): + atoms, nsteps = read_restart(traj_path) + + assert atoms is None + assert nsteps == 0 + + +def test_compute_density(tmp_path: Path) -> None: + """Test density is computed from the production half of the volumes.""" + traj_path = tmp_path / "X_0.traj" + write_trajectory(traj_path, (1.0,)) + + # Equilibration at 1000 A^3, production at 200 A^3 + volume_path = tmp_path / "volume.dat" + volume_path.write_text( + "# step volume_A3\n0 1000\n200 1000\n400 200\n600 200\n" + ) + + # 100 amu * 1.66053906660e-24 g/amu / (200 A^3 * 1e-24 cm^3/A^3) + assert compute_density(traj_path, volume_path) == pytest.approx(0.830269533) + + +def test_compute_density_too_few_samples(tmp_path: Path) -> None: + """Test NaN is returned if the production window is too short.""" + traj_path = tmp_path / "X_0.traj" + write_trajectory(traj_path, (1.0,)) + + # Half of two samples is one sample, fewer than MIN_SAMPLES + volume_path = tmp_path / "volume.dat" + volume_path.write_text("# step volume_A3\n0 1000\n200 200\n") + + assert np.isnan(compute_density(traj_path, volume_path)) + + +def test_compute_r_factor_identical_curves() -> None: + """Test the R-factor of a structure factor against itself is zero.""" + q = np.arange(1.0, 4.01, 0.05) + sq = np.exp(-q) + + assert compute_r_factor(q, sq, q, sq) == pytest.approx(0.0) + + +def test_compute_r_factor_constant_offset() -> None: + """Test the R-factor of a structure factor offset by a known amount.""" + q = np.array([1.0, 2.0, 3.0, 4.0]) + sq_ref = np.ones(4) + sq_calc = np.full(4, 1.5) + + # sum|1 - 1.5| = 2, sum|1| = 4 + assert compute_r_factor(q, sq_ref, q, sq_calc) == pytest.approx(0.5) + + +def test_compute_r_factor_without_overlap() -> None: + """Test NaN is returned if the two structure factors do not overlap.""" + q_ref = np.array([1.0, 2.0, 3.0, 4.0]) + q_calc = np.array([10.0, 11.0, 12.0]) + + assert np.isnan(compute_r_factor(q_ref, np.ones(4), q_calc, np.ones(3)))