feat:add parametric melt-quench workflow and docs

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",

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

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

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.

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