diff --git a/amorphouspy/src/amorphouspy/__init__.py b/amorphouspy/src/amorphouspy/__init__.py index 682a604..51fd35a 100644 --- a/amorphouspy/src/amorphouspy/__init__.py +++ b/amorphouspy/src/amorphouspy/__init__.py @@ -42,6 +42,7 @@ from amorphouspy.workflows.elastic_mod import elastic_simulation from amorphouspy.workflows.md import md_simulation from amorphouspy.workflows.meltquench import melt_quench_simulation +from amorphouspy.workflows.parametric import parametric_melt_quench from amorphouspy.workflows.structural_analysis import analyze_structure, find_rdf_minimum, plot_analysis_results_plotly from amorphouspy.workflows.viscosity import fit_vft, get_viscosity, viscosity_ensemble, viscosity_simulation @@ -80,6 +81,7 @@ "load_lammps_dump", "md_simulation", "melt_quench_simulation", + "parametric_melt_quench", "parse_formula", "plan_system", "plot_analysis_results_plotly", diff --git a/amorphouspy/src/amorphouspy/workflows/parametric.py b/amorphouspy/src/amorphouspy/workflows/parametric.py new file mode 100644 index 0000000..7d38c54 --- /dev/null +++ b/amorphouspy/src/amorphouspy/workflows/parametric.py @@ -0,0 +1,123 @@ +"""Parametric study workflow for system size and cooling rate sweeps. + +Author: Achraf Atila (achraf.atila@bam.de) +""" + +from pathlib import Path + +from amorphouspy.potentials.potential import generate_potential +from amorphouspy.structure import get_ase_structure, get_structure_dict +from amorphouspy.workflows.meltquench import melt_quench_simulation + + +def parametric_melt_quench( + composition: dict[str, float], + potential_type: str, + study: dict[int, list[float]], + temperature_high: float | None = None, + temperature_low: float = 300.0, + timestep: float = 1.0, + heating_rate: float | None = 1e13, + n_print: int = 1000, + equilibration_steps: int | None = None, + *, + server_kwargs: dict | None = None, + seed: int = 12345, + tmp_working_directory: str | Path | None = None, +) -> list[dict]: + """Run a parametric melt-quench study over system sizes and cooling rates. + + Generates one initial structure per unique system size, then runs a melt-quench + simulation for each cooling rate assigned to that size. Each size can have its + own subset of cooling rates, so the run matrix does not have to be a full grid. + + Args: + composition: Oxide glass composition mapping oxide formula to mol%. + Example: ``{"SiO2": 70, "Na2O": 15, "CaO": 15}``. + potential_type: Potential name (``"pmmcs"``, ``"bjp"``, or ``"shik"``). + study: Mapping of target atom count to the cooling rates (K/s) to run for + that size. Example:: + + { + 1000: [1e11, 1e12, 1e13], + 10000: [1e11], + } + + temperature_high: Melt temperature in K. Defaults to the protocol-specific + value when ``None`` (4000 K for SHIK, 5000 K for PMMCS/BJP). + temperature_low: Final glass temperature in K (default 300 K). + timestep: MD timestep in femtoseconds (default 1.0 fs). + heating_rate: Heating rate in K/s (default 1e13 K/s). Pass ``None`` to + defer to the protocol default inside ``melt_quench_simulation``. + n_print: Output frequency in simulation steps (default 1000). + equilibration_steps: Override for all fixed equilibration stages inside the + protocol. ``None`` uses protocol-specific defaults. + server_kwargs: LAMMPS server kwargs (e.g. ``{"cores": 4}``). + seed: Random seed for velocity initialisation (default 12345). + tmp_working_directory: Directory for LAMMPS temporary files. + + Returns: + Flat list of result dicts, one per (n_atoms, cooling_rate) pair, in the + order they appear in ``study``. Each dict contains: + + - ``n_atoms`` — actual atom count after integer formula-unit rounding + - ``target_n_atoms`` — the requested target (may differ slightly) + - ``cooling_rate`` — cooling rate used (K/s) + - ``heating_rate`` — heating rate used (K/s) + - ``structure`` — final quenched ASE ``Atoms`` object + - ``result`` — per-stage thermodynamic history from ``melt_quench_simulation`` + + Example: + >>> results = parametric_melt_quench( + ... composition={"SiO2": 70, "Na2O": 15, "CaO": 15}, + ... potential_type="pmmcs", + ... study={ + ... 1000: [1e11, 1e12, 1e13], + ... 10000: [1e11], + ... }, + ... ) + >>> for r in results: + ... print(r["target_n_atoms"], r["cooling_rate"], len(r["structure"])) + + """ + heating_rate = heating_rate if heating_rate is not None else 1e13 + + # Build one (structure, potential) per unique system size. + structure_cache: dict[int, tuple] = {} + for n_atoms in study: + atoms_dict = get_structure_dict(composition=composition, target_atoms=n_atoms) + structure = get_ase_structure(atoms_dict=atoms_dict) + potential = generate_potential(atoms_dict=atoms_dict, potential_type=potential_type) + structure_cache[n_atoms] = (structure, potential) + + results = [] + for n_atoms, cooling_rates in study.items(): + structure, potential = structure_cache[n_atoms] + for cooling_rate in cooling_rates: + mq = melt_quench_simulation( + structure=structure, + potential=potential, + temperature_high=temperature_high, + temperature_low=temperature_low, + timestep=timestep, + heating_rate=heating_rate, + cooling_rate=cooling_rate, + n_print=n_print, + equilibration_steps=equilibration_steps, + server_kwargs=server_kwargs, + seed=seed, + tmp_working_directory=tmp_working_directory, + ) + + results.append( + { + "n_atoms": len(structure), + "target_n_atoms": n_atoms, + "cooling_rate": cooling_rate, + "heating_rate": heating_rate, + "structure": mq["structure"], + "result": mq["result"], + } + ) + + return results diff --git a/amorphouspy/src/tests/test_parametric.py b/amorphouspy/src/tests/test_parametric.py new file mode 100644 index 0000000..5341d95 --- /dev/null +++ b/amorphouspy/src/tests/test_parametric.py @@ -0,0 +1,80 @@ +"""Unit tests for parametric_melt_quench workflow.""" + +from unittest.mock import MagicMock, patch + +import pytest +from amorphouspy.workflows.parametric import parametric_melt_quench +from ase import Atoms + +COMPOSITION = {"SiO2": 100} +MOCK_STRUCTURE = Atoms("SiO2", positions=[[0, 0, 0], [1, 0, 0], [0, 1, 0]]) + + +@pytest.fixture(autouse=True) +def _patch_deps(): + """Patch all external dependencies of parametric_melt_quench for every test. + + Replaces get_structure_dict, get_ase_structure, generate_potential, and + melt_quench_simulation with lightweight fakes so no LAMMPS installation is needed. + Yields the get_ase_structure mock for tests that need to inspect call counts. + """ + with ( + patch("amorphouspy.workflows.parametric.get_structure_dict", return_value={"Si": 1, "O": 2}), + patch("amorphouspy.workflows.parametric.get_ase_structure", return_value=MOCK_STRUCTURE) as mock_structure, + patch("amorphouspy.workflows.parametric.generate_potential", return_value=MagicMock()), + patch( + "amorphouspy.workflows.parametric.melt_quench_simulation", + return_value={"structure": MOCK_STRUCTURE, "result": {}}, + ), + ): + yield mock_structure + + +@pytest.fixture +def mock_structure(_patch_deps): + """Expose the get_ase_structure mock from _patch_deps for call-count assertions.""" + return _patch_deps + + +def test_structure_generated_once_per_size(mock_structure): + """Structure generation is cached per unique size, not repeated per cooling rate. + + With 2 distinct sizes and 3 total (size, rate) pairs, get_ase_structure must be + called exactly twice — once per size entry in the study dict. + """ + study = {500: [1e11, 1e12], 1000: [1e13]} # 2 sizes, 3 total runs + + parametric_melt_quench(COMPOSITION, "pmmcs", study) + + assert mock_structure.call_count == 2 # NOT 3 + + +def test_results_length_matches_total_runs(): + """Returned list has one entry per (n_atoms, cooling_rate) pair, not per size.""" + study = {500: [1e11, 1e12], 1000: [1e13]} # 3 (n_atoms, cooling_rate) pairs + + results = parametric_melt_quench(COMPOSITION, "pmmcs", study) + + assert len(results) == 3 + + +def test_result_dict_has_correct_keys_and_metadata(): + """Each result dict contains all expected keys with values from the call arguments.""" + study = {500: [1e12]} + + results = parametric_melt_quench(COMPOSITION, "pmmcs", study, heating_rate=1e14) + + r = results[0] + assert set(r.keys()) == {"n_atoms", "target_n_atoms", "cooling_rate", "heating_rate", "structure", "result"} + assert r["target_n_atoms"] == 500 + assert r["cooling_rate"] == 1e12 + assert r["heating_rate"] == 1e14 + + +def test_heating_rate_none_defaults_to_1e13(): + """Passing heating_rate=None substitutes the protocol default of 1e13 K/s.""" + study = {500: [1e12]} + + results = parametric_melt_quench(COMPOSITION, "pmmcs", study, heating_rate=None) + + assert results[0]["heating_rate"] == 1e13 diff --git a/docs/theory/workflows/parametric.md b/docs/theory/workflows/parametric.md new file mode 100644 index 0000000..3dba77e --- /dev/null +++ b/docs/theory/workflows/parametric.md @@ -0,0 +1,124 @@ +# Parametric Study + +The parametric study workflow sweeps over **system size** (`n_atoms`) and **cooling rate** in a single call to identify converged simulation settings before committing to expensive production runs. + +--- + +## Why a Parametric Study? + +Two parameters dominate the cost and quality of a melt-quench glass simulation: + +| Parameter | Effect on quality | Effect on cost | +|---|---|---| +| **System size** (`n_atoms`) | Larger boxes suppress finite-size artefacts in structural properties (RDF, coordination, ring statistics) | Scales as O(N log N) per step; larger boxes also need more steps to equilibrate | +| **Cooling rate** | Slower cooling → more relaxed glass, closer to experimental fictive temperature | Steps scale linearly with 1/rate | + +The `study` dict lets you assign different cooling rates to each system size — run a full sweep on small systems to establish trends, then limit large expensive systems to only the rates that matter. + +--- + +## Process Overview + +```mermaid +graph TD + A[study dict: n_atoms → cooling_rates] --> B{For each unique n_atoms} + B --> C[Generate random structure] + C --> D[Generate potential] + D --> E{For each assigned cooling_rate} + E --> F[Run melt_quench_simulation] + F --> G[Collect structure + thermo] + G --> H[Flat results list] +``` + +Structures are cached per system size, so each unique system size is generated only once regardless of how many cooling rates are assigned to it. + +--- + +## Basic Usage + +### `parametric_melt_quench(composition, potential_type, study, ...)` + +```python +from amorphouspy import parametric_melt_quench + +results = parametric_melt_quench( + composition={"SiO2": 70, "Na2O": 15, "CaO": 15}, + potential_type="pmmcs", + study={ + 300: [1e11, 1e12, 1e13], # small system — full rate sweep + 1000: [1e11, 1e12, 1e13], # medium system — full rate sweep + 10000: [1e13], # large system — production rate only + 100000: [1e13], # large system — production rate only + }, + temperature_high=5000.0, + temperature_low=300.0, + timestep=1.0, + # equilibration_steps=5000, # reduced for fast screening; remove for production + seed=42, +) +``` + +**Parameters:** + +| Parameter | Type | Default | Description | +|---|---|---|---| +| `composition` | `dict[str, float]` | — | Oxide mol% composition, e.g. `{"SiO2": 70, "Na2O": 15, "CaO": 15}` | +| `potential_type` | `str` | — | Potential name: `"pmmcs"`, `"bjp"`, or `"shik"` | +| `study` | `dict[int, list[float]]` | — | Mapping of target atom count to the cooling rates (K/s) to run for that size | +| `temperature_high` | `float | None` | `None` | Melt temperature in K. `None` uses the protocol default (5000 K for PMMCS/BJP, 4000 K for SHIK) | +| `temperature_low` | `float` | `300.0` | Final glass temperature in K | +| `timestep` | `float` | `1.0` | MD timestep in femtoseconds | +| `heating_rate` | `float | None` | `1e13` | Heating rate in K/s. `None` defers to the protocol default inside `melt_quench_simulation` | +| `n_print` | `int` | `1000` | Output frequency in simulation steps | +| `equilibration_steps` | `int | None` | `None` | Override for all fixed equilibration stages. `None` uses protocol defaults | +| `server_kwargs` | `dict | None` | `None` | LAMMPS server kwargs, e.g. `{"cores": 4}` | +| `seed` | `int` | `12345` | Random seed for velocity initialisation | +| `tmp_working_directory` | `str | Path | None` | `None` | Directory for LAMMPS temporary files | + +**Returns:** Flat list of dicts in the order they appear in `study` (outer loop over sizes, inner loop over each size's cooling rates): + +| Key | Type | Description | +|---|---|---| +| `"n_atoms"` | `int` | Actual atom count after integer formula-unit rounding | +| `"target_n_atoms"` | `int` | Requested target (may differ slightly) | +| `"cooling_rate"` | `float` | Cooling rate used (K/s) | +| `"heating_rate"` | `float` | Heating rate used (K/s) | +| `"structure"` | `Atoms` | Final quenched ASE `Atoms` object | +| `"result"` | `list[dict | None]` | Per-stage thermodynamic history from `melt_quench_simulation` | + +--- + +## Convergence Analysis + +### System-size convergence + +Plot a structural property (e.g. density, first-peak position in the RDF) against `n_atoms` at the slowest cooling rate. A plateau indicates that finite-size effects are small. Depending on the property of interest, different sizes may be needed — short-range order converges with ~1000 atoms, while the thermal expansion coefficient requires much larger systems. + +### Cooling-rate convergence + +Plot the same property against cooling rate for a selected system. Density typically increases as the cooling rate decreases (more relaxed glass). Extrapolating to zero rate (e.g. Vogel-Fulcher-Tammann fit) gives an estimate of the equilibrium glass density. + +### Example workflow + +```python +# 1. Screen — fast rates on small systems to identify the trend +study = {200: [1e13, 1e14, 1e15], 500: [1e13, 1e14, 1e15]} +results = parametric_melt_quench(..., study=study, equilibration_steps=100000) + +# 2. Refine — moderate rates on medium systems, protocol equilibration defaults +study = {500: [1e11, 1e12, 1e13], 1000: [1e11, 1e12, 1e13], 2000: [1e11, 1e12]} +results = parametric_melt_quench(..., study=study) + +# 3. Production — single run at the converged (n_atoms, cooling_rate) setting +results = parametric_melt_quench(..., study={2000: [1e11]}) +# Pass structure to analyze_structure(), compute_rdf(), etc. +``` + +--- + +## Tips + +- **Heating rate**: defaults to `1e13 K/s`, matching the standard melt-quench convention. Override with `heating_rate=...` only when you need a specific ratio. +- **Structure caching**: the same random initial structure is reused across all cooling rates for a given `n_atoms`. This keeps comparisons clean but means all cooling rates for a given size start from the same initial configuration. +- **Reduced equilibration for screening**: pass `equilibration_steps=5000` to cut run times dramatically during the sweep. Remove it (or set `None`) for production. +- **Independent samples**: the parametric grid uses a single random seed per system size. For statistical averaging, run the grid multiple times with different `seed` values and average results. diff --git a/mkdocs.yml b/mkdocs.yml index 3232d03..8aa2a19 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -43,16 +43,13 @@ nav: - how_to_guides/index.md - how_to_guides/installation.md - how_to_guides/running_api.md - - MD Tutorial: notebooks/MDTutorial.ipynb - Meltquench Tutorial: notebooks/MeltquenchTutorial.ipynb - # - Meltquench PMMCS: notebooks/Meltquench_PMMCS.ipynb - # - Meltquench BJP: notebooks/Meltquench_BJP.ipynb - # - Meltquench SHIK: notebooks/Meltquench_SHIK.ipynb - Structure Tutorial: notebooks/StructureTutorial.ipynb - - Elastic Properties: notebooks/ElasticProperties.ipynb - Viscosity: notebooks/ViscosityTutorial.ipynb - CTE: notebooks/CTE.ipynb - # - Meltquench + Analysis: notebooks/meltquench_Structure_analysis.ipynb + - Elastic Properties: notebooks/ElasticProperties.ipynb + - MD Tutorial: notebooks/MDTutorial.ipynb + #- Parametric Study: notebooks/ParametricStudyTutorial.ipynb - Theory: - theory/index.md - Structure generation: theory/structure.md @@ -68,6 +65,7 @@ nav: - Elastic Moduli: theory/workflows/elastic.md - Viscosity: theory/workflows/viscosity.md - CTE: theory/workflows/cte.md + - Parametric Study: theory/workflows/parametric.md - Structural Analysis: - theory/analysis/index.md - RDF & Coordination: theory/analysis/rdf.md