diff --git a/docs/source/benchmarks/molecular_liquids/radial_distribution.rst b/docs/source/benchmarks/molecular_liquids/radial_distribution.rst index bd8f3a5d..7598fe4f 100644 --- a/docs/source/benchmarks/molecular_liquids/radial_distribution.rst +++ b/docs/source/benchmarks/molecular_liquids/radial_distribution.rst @@ -16,12 +16,15 @@ and emergent properties of liquid systems. Description ----------- -The benchmark performs an **MD** simulation using the **MLIP** model in the **NVT** ensemble at -**300 K** for **500,000 steps**, leveraging the `jax-md `_ engine -from the `mlip `_ library. The starting configuration is already +The benchmark performs an **MD** simulation using the **MLIP** model in the **NPT** ensemble for **500,000 steps**, +leveraging the `jax-md `_ engine from the +`mlip `_ library. Water is run at **295.15 K** and **1 atm**, +while all other solvents are run at **293.15 K** and **1 atm**. The starting configuration is already equilibrated. For every specific atom pair (e.g., **oxygen-oxygen** in water) the radial distribution function (**RDF** or **g(r)**) is calculated from the simulation. +# TODO: Add details of density computation + details RE experimental density to sections below + .. figure:: img/rdf.png :figwidth: 35% :align: center diff --git a/src/mlipaudit/benchmarks/solvent_radial_distribution/solvent_radial_distribution.py b/src/mlipaudit/benchmarks/solvent_radial_distribution/solvent_radial_distribution.py index b8b72c32..175c5f26 100644 --- a/src/mlipaudit/benchmarks/solvent_radial_distribution/solvent_radial_distribution.py +++ b/src/mlipaudit/benchmarks/solvent_radial_distribution/solvent_radial_distribution.py @@ -20,6 +20,7 @@ from ase import Atoms, units from ase.io import read as ase_read from mlip.simulation import SimulationState +from mlip.simulation.enums import MDIntegrator from pydantic import BaseModel, ConfigDict, NonNegativeFloat from mlipaudit.benchmark import ( @@ -35,6 +36,7 @@ create_mdtraj_trajectory_from_simulation_state, run_simulation, ) +from mlipaudit.utils.molecular_liquids import compute_densities from mlipaudit.utils.stability import is_simulation_stable logger = logging.getLogger("mlipaudit") @@ -43,29 +45,26 @@ "num_steps": 500_000, "snapshot_interval": 500, "num_episodes": 1000, - "temperature_kelvin": 295.15, + "temperature_kelvin": 293.15, + "pressure_bar": 1.01325, } SIMULATION_CONFIG_DEV = { "num_steps": 5, "snapshot_interval": 1, "num_episodes": 1, - "temperature_kelvin": 295.15, + "temperature_kelvin": 293.15, + "pressure_bar": 1.01325, } SIMULATION_CONFIG_FAST = { "num_steps": 250_000, "snapshot_interval": 250, "num_episodes": 1000, - "temperature_kelvin": 295.15, + "temperature_kelvin": 293.15, + "pressure_bar": 1.01325, } NUM_DEV_SYSTEMS = 1 -BOX_CONFIG = { # In Angstrom - "CCl4": 28.575, - "methanol": 29.592, - "acetonitrile": 27.816, -} - SYSTEM_ATOM_OF_INTEREST = { "CCl4": "C", "methanol": "O", @@ -74,6 +73,11 @@ MIN_RADII, MAX_RADII = 0.0, 20.0 # In Angstrom +MOLECULE_CONFIG = { + "CCl4": {"molecule_weight": 153.823, "atoms_per_molecule": 5}, + "methanol": {"molecule_weight": 32.042, "atoms_per_molecule": 6}, + "acetonitrile": {"molecule_weight": 41.053, "atoms_per_molecule": 6}, +} REFERENCE_MAXIMA = { "CCl4": {"type": "C-C", "distance": 5.9, "range": (0.0, 20.0)}, "acetonitrile": {"type": "N-N", "distance": 4.0, "range": (3.5, 4.5)}, @@ -85,6 +89,12 @@ "methanol": (0.0, 20.0), } +REFERENCE_DENSITIES = { + "CCl4": 1.594, + "acetonitrile": 0.786, + "methanol": 0.791, +} + class SolventRadialDistributionModelOutput(ModelOutput): """Model output containing the final simulation states for @@ -108,19 +118,23 @@ class SolventRadialDistributionStructureResult(BaseModel): Attributes: structure_name: The structure name. + densities: List of densities in g/cm3. + average_density: Average density over the final 4 fifths of the frames. + density_deviation: Deviation of the average density from the reference. radii: The radii values in Angstrom. - rdf: The radial distribution function values at the - radii. - first_solvent_peak: The first solvent peak, i.e. - the radius at which the rdf is the maximum. - peak_deviation: The deviation of the - first solvent peak from the reference. + rdf: The radial distribution function values at the radii. + first_solvent_peak: The first solvent peak, i.e. the radius at which the + rdf is the maximum. + peak_deviation: The deviation of the first solvent peak from the reference. failed: Whether the simulation was successful. If unsuccessful, the other attributes will be not be set. score: The score for the molecule. """ structure_name: str + densities: list[float] | None = None + average_density: float | None = None + density_deviation: NonNegativeFloat | None = None radii: list[float] | None = None rdf: list[float] | None = None first_solvent_peak: float | None = None @@ -136,7 +150,8 @@ class SolventRadialDistributionResult(BenchmarkResult): Attributes: structure_names: The names of the structures. structures: List of per structure results. - avg_peak_deviation: The average deviation across all structures. + avg_density_deviation: The average density deviation across all structures. + avg_peak_deviation: The average peak deviation across all structures. failed: Whether all the simulations failed and no analysis could be performed. Defaults to False. score: The final score for the benchmark between @@ -145,6 +160,7 @@ class SolventRadialDistributionResult(BenchmarkResult): structure_names: list[str] structures: list[SolventRadialDistributionStructureResult] + avg_density_deviation: NonNegativeFloat | None = None avg_peak_deviation: NonNegativeFloat | None = None @@ -179,10 +195,11 @@ class SolventRadialDistributionBenchmark(Benchmark): required_elements = {"N", "H", "O", "C", "Cl"} def run_model(self) -> None: - """Run an MD simulation for each structure. + """Run an MD simulation for each structure using the NPT ensemble. + The MD simulation is performed using the JAX MD engine and starts from - the reference structure. NOTE: This benchmark runs a simulation in the - NVT ensemble, which is not recommended for a water RDF calculation. + the reference structure. The NPT integrator uses Langevin dynamics with + a Monte Carlo barostat. """ if self.run_mode == RunMode.DEV: md_kwargs = SIMULATION_CONFIG_DEV @@ -195,13 +212,13 @@ def run_model(self) -> None: for system_name in self._system_names: logger.info("Running MD for %s radial distribution function.", system_name) - md_kwargs["box"] = BOX_CONFIG[system_name] simulation_state = run_simulation( atoms=self._load_system(system_name), force_field=self.force_field, + md_integrator=MDIntegrator.NPT_MC_LANGEVIN, + molecule_indices=self._load_molecule_indices(system_name), **md_kwargs, ) - simulation_states.append(simulation_state) self.model_output = SolventRadialDistributionModelOutput( @@ -237,15 +254,21 @@ def analyze(self) -> SolventRadialDistributionResult: continue num_succeeded += 1 - - box_length = BOX_CONFIG[system_name] + mol_config = MOLECULE_CONFIG[system_name] + densities = compute_densities( + simulation_state, + mol_config["molecule_weight"], + int(mol_config["atoms_per_molecule"]), + ) + n_frames_equilibration = len(densities) // 5 + average_density = np.mean(densities[n_frames_equilibration:]) + density_deviation = abs(average_density - REFERENCE_DENSITIES[system_name]) traj = create_mdtraj_trajectory_from_simulation_state( simulation_state=simulation_state, topology_path=self.data_input_dir / self.name / self._get_pdb_file_name(system_name), - cell_lengths=(box_length, box_length, box_length), ) pair_indices = traj.top.select( f"symbol == {SYSTEM_ATOM_OF_INTEREST[system_name]}" @@ -283,12 +306,16 @@ def analyze(self) -> SolventRadialDistributionResult: peak_deviation = abs( first_solvent_peak - REFERENCE_MAXIMA[system_name]["distance"] ) + score = math.exp( -ALPHA * peak_deviation / REFERENCE_MAXIMA[system_name]["distance"] ) structure_result = SolventRadialDistributionStructureResult( structure_name=system_name, + densities=densities, + average_density=average_density, + density_deviation=density_deviation, radii=radii.tolist(), rdf=rdf, first_solvent_peak=first_solvent_peak, @@ -309,6 +336,11 @@ def analyze(self) -> SolventRadialDistributionResult: return SolventRadialDistributionResult( structure_names=self.model_output.structure_names, structures=structure_results, + avg_density_deviation=statistics.mean( + structure.density_deviation + for structure in structure_results + if structure.density_deviation is not None + ), avg_peak_deviation=statistics.mean( structure.peak_deviation for structure in structure_results @@ -322,10 +354,10 @@ def analyze(self) -> SolventRadialDistributionResult: @property def _system_names(self) -> list[str]: if self.run_mode == RunMode.STANDARD: - return list(BOX_CONFIG.keys()) + return list(SYSTEM_ATOM_OF_INTEREST.keys()) # reduced number of cases for DEV and FAST run mode - return list(BOX_CONFIG.keys())[:NUM_DEV_SYSTEMS] + return list(SYSTEM_ATOM_OF_INTEREST.keys())[:NUM_DEV_SYSTEMS] def _load_system(self, system_name) -> Atoms: atoms = ase_read( @@ -335,6 +367,17 @@ def _load_system(self, system_name) -> Atoms: atoms.info["spin"] = DEFAULT_SPIN return atoms + def _load_molecule_indices(self, system_name) -> np.ndarray: + molecule_indices_filename = self._get_molecule_indices_file_name(system_name) + molecule_indices = np.load( + self.data_input_dir / self.name / molecule_indices_filename + ) + return molecule_indices + @staticmethod def _get_pdb_file_name(system_name: str) -> str: return f"{system_name}_eq.pdb" + + @staticmethod + def _get_molecule_indices_file_name(system_name: str) -> str: + return f"{system_name}_molecule_indices.npy" diff --git a/src/mlipaudit/benchmarks/water_radial_distribution/water_radial_distribution.py b/src/mlipaudit/benchmarks/water_radial_distribution/water_radial_distribution.py index 23d1a0ac..a35dd237 100644 --- a/src/mlipaudit/benchmarks/water_radial_distribution/water_radial_distribution.py +++ b/src/mlipaudit/benchmarks/water_radial_distribution/water_radial_distribution.py @@ -21,6 +21,7 @@ from ase import Atoms, units from ase.io import read as ase_read from mlip.simulation import SimulationState +from mlip.simulation.enums import MDIntegrator from pydantic import ConfigDict, NonNegativeFloat from sklearn.metrics import mean_absolute_error, root_mean_squared_error @@ -34,6 +35,7 @@ from mlipaudit.run_mode import RunMode from mlipaudit.scoring import ALPHA, compute_metric_score from mlipaudit.utils import run_simulation +from mlipaudit.utils.molecular_liquids import compute_densities from mlipaudit.utils.stability import is_simulation_stable from mlipaudit.utils.trajectory_helpers import ( create_mdtraj_trajectory_from_simulation_state, @@ -46,7 +48,7 @@ "snapshot_interval": 500, "num_episodes": 1000, "temperature_kelvin": 295.15, - "box": 24.772, + "pressure_bar": 1.01325, } SIMULATION_CONFIG_FAST = { @@ -54,7 +56,7 @@ "snapshot_interval": 250, "num_episodes": 1000, "temperature_kelvin": 295.15, - "box": 24.772, + "pressure_bar": 1.01325, } SIMULATION_CONFIG_DEV = { @@ -62,20 +64,25 @@ "snapshot_interval": 1, "num_episodes": 1, "temperature_kelvin": 295.15, - "box": 24.772, + "pressure_bar": 1.01325, } WATERBOX_N500 = "water_box_n500_eq.pdb" +MOLECULE_INDICES_PATH = "water_box_n500_molecule_indices.npy" REFERENCE_DATA = "experimental_reference.npz" +MOLECULE_WEIGHT = 18.01528 # g/mol +ATOMS_PER_MOLECULE = 3 +REFERENCE_PEAK_DISTANCE = 2.80 # A RMSE_SCORE_THRESHOLD = 0.1 SOLVENT_PEAK_RANGE = (2.8, 3.0) RADII_RANGE = (2.5, 10.0) +REFERENCE_DENSITY = 0.997773 # g/cm3 + class WaterRadialDistributionModelOutput(ModelOutput): - """Model output containing the final simulation state of - the water box. + """Model output containing the final simulation state of the water box. Attributes: simulation_state: The final simulation state of the water @@ -93,9 +100,11 @@ class WaterRadialDistributionResult(BenchmarkResult): """Result object for the water radial distribution benchmark. Attributes: + densities: List of densities in g/cm3. + average_density: Average density over the final 4 fifths of the frames. + density_deviation: Deviation of the average density from the reference. radii: The radii values in Angstrom. - rdf: The radial distribution function values at the - radii. + rdf: The radial distribution function values at the radii. mae: The MAE of the radial distribution function values. rmse: The RMSE of the radial distribution function values. first_solvent_peak: The first solvent peak, i.e. @@ -106,10 +115,12 @@ class WaterRadialDistributionResult(BenchmarkResult): radial distribution function error metrics. failed: Whether all the simulations failed and no analysis could be performed. Defaults to False. - score: The final score for the benchmark between - 0 and 1. + score: The final score for the benchmark between 0 and 1. """ + densities: list[float] | None = None + average_density: float | None = None + density_deviation: NonNegativeFloat | None = None radii: list[float] | None = None rdf: list[float] | None = None mae: float | None = None @@ -150,17 +161,19 @@ class WaterRadialDistributionBenchmark(Benchmark): required_elements = {"H", "O"} def run_model(self) -> None: - """Run an MD simulation for each structure. + """Run an MD simulation for the water box system using the NPT ensemble. The MD simulation is performed using the JAX MD engine and starts from - the reference structure. NOTE: This benchmark runs a simulation in the - NVT ensemble, which is not recommended for a water RDF calculation. + the reference structure. The NPT integrator uses Langevin dynamics with + a Monte Carlo barostat. """ - logger.info("Running MD for for water radial distribution function.") + logger.info("Running MD for water radial distribution function.") simulation_state = run_simulation( atoms=self._water_box_n500, force_field=self.force_field, + md_integrator=MDIntegrator.NPT_MC_LANGEVIN, + molecule_indices=self._molecule_indices, **self._md_kwargs, ) @@ -180,17 +193,21 @@ def analyze(self) -> WaterRadialDistributionResult: if self.model_output is None: raise RuntimeError("Must call run_model() first.") - if self.model_output.failed or not is_simulation_stable( - self.model_output.simulation_state - ): + simulation_state = self.model_output.simulation_state + + if self.model_output.failed or not is_simulation_stable(simulation_state): return WaterRadialDistributionResult(failed=True, score=0.0) - box_length = self._md_kwargs["box"] + densities = compute_densities( + simulation_state, MOLECULE_WEIGHT, ATOMS_PER_MOLECULE + ) + n_frames_equilibration = len(densities) // 5 + average_density = np.mean(densities[n_frames_equilibration:]) + density_deviation = abs(average_density - REFERENCE_DENSITY) traj = create_mdtraj_trajectory_from_simulation_state( - self.model_output.simulation_state, + simulation_state, self.data_input_dir / self.name / WATERBOX_N500, - cell_lengths=(box_length, box_length, box_length), ) oxygen_indices = traj.top.select("symbol == O") @@ -233,23 +250,23 @@ def analyze(self) -> WaterRadialDistributionResult: first_solvent_peak = radii[np.argmax(g_r)].item() - peak_deviation = max( - 0, - SOLVENT_PEAK_RANGE[0] - first_solvent_peak, - first_solvent_peak - SOLVENT_PEAK_RANGE[1], - ) + peak_deviation = abs(first_solvent_peak - REFERENCE_PEAK_DISTANCE) + peak_deviation_score = math.exp( - -ALPHA - * peak_deviation - / ((SOLVENT_PEAK_RANGE[0] + SOLVENT_PEAK_RANGE[1]) / 2) + -ALPHA * peak_deviation / REFERENCE_PEAK_DISTANCE ) rmse_score = compute_metric_score( np.array([rmse]), RMSE_SCORE_THRESHOLD, ALPHA ).item() + score = (peak_deviation_score + rmse_score) / 2 + # TODO: Remove `range_of_interest=SOLVENT_PEAK_RANGE`? return WaterRadialDistributionResult( + densities=densities, + average_density=average_density, + density_deviation=density_deviation, radii=radii.tolist(), rdf=rdf, mae=mae, @@ -276,9 +293,17 @@ def _water_box_n500(self) -> Atoms: atoms.info["spin"] = DEFAULT_SPIN return atoms + @functools.cached_property + def _molecule_indices(self) -> np.ndarray: + molecule_indices = np.load( + self.data_input_dir / self.name / MOLECULE_INDICES_PATH + ) + return molecule_indices + @functools.cached_property def _reference_data(self): """The experimental reference data for the water RDF benchmark. + Contains keys 'r_OO' and 'g_OO', the radii and RDF values. The radii are in Angstrom. """ diff --git a/src/mlipaudit/ui/solvent_radial_distribution.py b/src/mlipaudit/ui/solvent_radial_distribution.py index 4a52ff4d..e1bf07bc 100644 --- a/src/mlipaudit/ui/solvent_radial_distribution.py +++ b/src/mlipaudit/ui/solvent_radial_distribution.py @@ -57,6 +57,7 @@ def _process_data_into_dataframe( "Model name": model_name, "Score": result.score, "Average peak deviation (Å)": result.avg_peak_deviation, + "Average density deviation (g/cm3)": result.avg_density_deviation, } for structure_res in result.structures: if structure_res.failed: @@ -65,6 +66,10 @@ def _process_data_into_dataframe( model_data_converted[ f"{structure_res.structure_name} peak deviation (Å)" ] = structure_res.peak_deviation + + model_data_converted[ + f"{structure_res.structure_name} density deviation (g/cm3)" + ] = structure_res.density_deviation converted_data_scores.append(model_data_converted) df = pd.DataFrame(converted_data_scores) return df diff --git a/src/mlipaudit/ui/water_radial_distribution.py b/src/mlipaudit/ui/water_radial_distribution.py index d9a4b5bb..1e2aeef1 100644 --- a/src/mlipaudit/ui/water_radial_distribution.py +++ b/src/mlipaudit/ui/water_radial_distribution.py @@ -25,9 +25,6 @@ WaterRadialDistributionBenchmark, WaterRadialDistributionResult, ) -from mlipaudit.benchmarks.water_radial_distribution.water_radial_distribution import ( - SOLVENT_PEAK_RANGE, -) from mlipaudit.ui.page_wrapper import UIPageWrapper from mlipaudit.ui.utils import ( display_failed_models, @@ -47,6 +44,7 @@ ] +# TODO: Update this file with density scores @st.cache_resource def _load_tip3p() -> NpzFile: return np.load(WATER_RADIAL_DISTRIBUTION_DATA_DIR / "tip3p_500ps.npz") @@ -74,9 +72,9 @@ def _process_data_into_dataframe( "RMSE (Å)": result.rmse, "MAE (Å)": result.mae, "First solvent peak (Å)": result.first_solvent_peak, - "Solvent peak acceptable minimum (Å)": SOLVENT_PEAK_RANGE[0], - "Solvent peak acceptable maximum (Å)": SOLVENT_PEAK_RANGE[1], "Peak deviation (Å)": result.peak_deviation, + "Equilibrium density (g/cm3)": result.average_density, + "Density deviation (g/cm3)": result.density_deviation, } converted_data_scores.append(model_data_converted) model_names.append(model_name) diff --git a/src/mlipaudit/utils/molecular_liquids.py b/src/mlipaudit/utils/molecular_liquids.py new file mode 100644 index 00000000..2b3620d4 --- /dev/null +++ b/src/mlipaudit/utils/molecular_liquids.py @@ -0,0 +1,45 @@ +# Copyright 2025 InstaDeep Ltd +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +import numpy as np +from ase import units +from mlip.simulation import SimulationState + +ANGSTROM3_TO_CM3 = 1e-24 + + +def compute_densities( + simulation_state: SimulationState, molecule_weight: float, atoms_per_molecule: int +) -> np.ndarray: + """Compute the density (g/cm3) for each frame of the simulation. + + Used by `WaterRadialDistributionFunction` and `SolventRadialDistributionFunction`. + TODO: Should be moved into benchmark file when a density benchmark is created. + + Args: + simulation_state: The final simulation state. + molecule_weight: Molecular weight of each solvent molecule. + atoms_per_molecule: Number of atoms in each solvent molecule. + + Returns: + densities: Computed density (g/cm3) for each frame of the simulation. + """ + n_molecules = simulation_state.positions.shape[1] / atoms_per_molecule + volumes = np.abs(np.linalg.det(simulation_state.cell)) + + density_numerator = molecule_weight * n_molecules + density_denominator = units.mol * volumes * ANGSTROM3_TO_CM3 + + densities = density_numerator / density_denominator + return densities diff --git a/src/mlipaudit/utils/trajectory_helpers.py b/src/mlipaudit/utils/trajectory_helpers.py index c43a9019..a6be14e1 100644 --- a/src/mlipaudit/utils/trajectory_helpers.py +++ b/src/mlipaudit/utils/trajectory_helpers.py @@ -19,6 +19,7 @@ import mdtraj as md import numpy as np from ase import Atoms, units +from ase.cell import Cell from ase.io import write as ase_write from mlip.simulation import SimulationState @@ -35,11 +36,16 @@ def create_mdtraj_trajectory_from_simulation_state( to save the trajectory as an xyz file. All input values should be in Angstrom units. Note that the resulting trajectory uses nm as units. + Note that if `simulation_state` carries `cell` information (e.g. in the case of an + NPT simulation), then the `cell_lengths` and `cell_angles` inputs are ignored. + Args: simulation_state: The state containing the trajectory. topology_path: The path towards the topology file. Typically, a pdb file. - cell_lengths: The lengths of the unit cell in Angstrom. Default is `None`. - cell_angles: The angles of the unit cell in degrees. Default is `(90, 90, 90)`. + cell_lengths: The lengths of the unit cell in Angstrom. Ignored if + `simulation_state` contains `cell` information. Default is `None`. + cell_angles: The angles of the unit cell in degrees. Ignored if + `simulation_state` contains `cell` information. Default is `(90, 90, 90)`. Returns: The converted trajectory. @@ -49,7 +55,12 @@ def create_mdtraj_trajectory_from_simulation_state( _tmp_path = Path(tmpdir) ase_write(_tmp_path / "traj.xyz", ase_traj) traj = md.load(_tmp_path / "traj.xyz", top=topology_path) - if cell_lengths is not None: + if simulation_state.cell is not None: + # cell shape: (n_frames, 3, 3) — extract cellpar per frame + cellpars = np.array([Cell(c).cellpar() for c in simulation_state.cell]) + traj.unitcell_lengths = cellpars[:, :3] * (units.Angstrom / units.nm) + traj.unitcell_angles = cellpars[:, 3:] + elif cell_lengths is not None: # converting length units to nm for mdtraj cell_lengths_converted = [ cell_length * (units.Angstrom / units.nm) @@ -76,10 +87,14 @@ def create_ase_trajectory_from_simulation_state( An ASE trajectory as a list of `ase.Atoms`. """ num_frames = simulation_state.positions.shape[0] + cell = simulation_state.cell + if cell is None: + cell = [None] * num_frames trajectory = [ Atoms( numbers=simulation_state.atomic_numbers, positions=simulation_state.positions[frame], + cell=cell[frame], ) for frame in range(num_frames) ] diff --git a/tests/conftest.py b/tests/conftest.py index c739978f..50a61bd0 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -128,6 +128,7 @@ def _factory(simulation_state: SimulationState | None = None): positions=np.random.rand(10, 2, 3), forces=np.random.rand(10, 2, 3), temperature=np.random.rand(10), + cell=np.abs(np.random.randn(10, 1, 1)) * np.eye(3), ) mock_engine.configure_mock(state=state) return mock_engine diff --git a/tests/data/solvent_radial_distribution/CCl4_molecule_indices.npy b/tests/data/solvent_radial_distribution/CCl4_molecule_indices.npy new file mode 100644 index 00000000..07f37e11 Binary files /dev/null and b/tests/data/solvent_radial_distribution/CCl4_molecule_indices.npy differ diff --git a/tests/data/solvent_radial_distribution/acetonitrile_molecule_indices.npy b/tests/data/solvent_radial_distribution/acetonitrile_molecule_indices.npy new file mode 100644 index 00000000..29fad1c0 Binary files /dev/null and b/tests/data/solvent_radial_distribution/acetonitrile_molecule_indices.npy differ diff --git a/tests/data/solvent_radial_distribution/methanol_molecule_indices.npy b/tests/data/solvent_radial_distribution/methanol_molecule_indices.npy new file mode 100644 index 00000000..5abd8bf0 Binary files /dev/null and b/tests/data/solvent_radial_distribution/methanol_molecule_indices.npy differ diff --git a/tests/data/water_radial_distribution/water_box_n500_molecule_indices.npy b/tests/data/water_radial_distribution/water_box_n500_molecule_indices.npy new file mode 100644 index 00000000..b03b57ae Binary files /dev/null and b/tests/data/water_radial_distribution/water_box_n500_molecule_indices.npy differ diff --git a/tests/solvent_radial_distribution/test_solvent_radial_distribution.py b/tests/solvent_radial_distribution/test_solvent_radial_distribution.py index 3f5f6b28..e96970a6 100644 --- a/tests/solvent_radial_distribution/test_solvent_radial_distribution.py +++ b/tests/solvent_radial_distribution/test_solvent_radial_distribution.py @@ -73,15 +73,15 @@ def test_full_run_with_mocked_engine( atoms = ase_read(INPUT_DATA_DIR / benchmark.name / "CCl4_eq.pdb") num_frames = 2 - positions = np.tile( - np.array(atoms.positions), - reps=(num_frames, 1, 1), - ) + positions = np.tile(np.array(atoms.positions), reps=(num_frames, 1, 1)) + cells = np.tile(np.array(atoms.get_cell()), reps=(num_frames, 1, 1)) benchmark.model_output = SolventRadialDistributionModelOutput( structure_names=["CCl4"], simulation_states=[ - SimulationState(positions=positions, temperature=np.ones(10)) + SimulationState( + positions=positions, temperature=np.ones(num_frames), cell=cells + ) ], ) @@ -93,9 +93,12 @@ def test_full_run_with_mocked_engine( # Some leeway around the maximum of 5.9 assert 5.4 < result.structures[0].first_solvent_peak < 6.4 - assert result.structures[0].peak_deviation < 0.5 + # Target density = 1.594, initial density = 1.368 + assert 1.3 < result.structures[0].average_density < 1.6 + assert result.structures[0].density_deviation < 0.3 + def test_analyze_raises_error_if_run_first(solvent_radial_distribution_benchmark): """Verifies the RuntimeError using the new fixture.""" diff --git a/tests/water_radial_distribution/test_water_radial_distribution.py b/tests/water_radial_distribution/test_water_radial_distribution.py index c60045e1..37eda854 100644 --- a/tests/water_radial_distribution/test_water_radial_distribution.py +++ b/tests/water_radial_distribution/test_water_radial_distribution.py @@ -77,10 +77,11 @@ def test_full_run_with_mocked_engine( np.load(INPUT_DATA_DIR / benchmark.name / "positions.npy"), reps=(num_frames, 1, 1), ) + cells = np.tile(24.772 * np.eye(3), reps=(num_frames, 1, 1)) benchmark.model_output = WaterRadialDistributionModelOutput( simulation_state=SimulationState( - positions=positions, temperature=np.ones(10) + positions=positions, temperature=np.ones(num_frames), cell=cells ) ) @@ -96,6 +97,10 @@ def test_full_run_with_mocked_engine( max_radii = np.array(result.radii)[np.argmax(np.array(result.rdf))] assert 2.5 < max_radii < 3.0 + # Target density = 0.997773, initial density = 0.9859266 + assert 0.9 < result.average_density < 1.1 + assert result.density_deviation < 0.1 + def test_analyze_raises_error_if_run_first(water_radial_distribution_benchmark): """Verifies the RuntimeError using the new fixture."""