Quickstart

A complete, minimal refinement: one structure, one potential, one neutron \(F(Q)\) dataset, single-atom translation moves. Everything else in the package is a variation on this scaffold.

from mosaic_materials.engine import MCEngine
from mosaic_materials.constraints import (
    EnthalpyConstraint,
    RDFConstraint,
    NeutronStructureFactorConstraint,
)
from mosaic_materials.data import ScatteringData
from mosaic_materials.monte_carlo.adapt import MoveSpec
from mosaic_materials.monte_carlo.driver import MCDriver
from mosaic_materials.monte_carlo.schedules import ConstantSchedule
from mosaic_materials.moves import TranslateMove

with MCEngine(cores=4, seed=42) as engine:

    # --- structure and potential ---
    engine.load_lammps_data("structure.data")
    engine.add_force_field(
        "pace recursive", "* * potential.yace C H N Zn"
    )

    # --- constraints: energy, RDF, and one experimental dataset ---
    engine.add_constraint(EnthalpyConstraint(temperature_Kelvin=300.0))

    r_max = min(engine.system.cell_lengths) / 2 - 1.0
    engine.add_constraint(RDFConstraint(r_max=r_max))

    fq = ScatteringData.from_file("sample.nfq", sigmas=0.01)
    fq.apply_finite_box_correction(r_max=r_max, inplace=True)
    engine.add_constraint(
        NeutronStructureFactorConstraint.from_scattering_data("c_n_fq", fq)
    )

    # --- moves ---
    moves = [
        MoveSpec(
            name="translate",
            move_class=TranslateMove,
            selector=lambda system, rng: [system.pick_random_atom(rng)],
            kwargs={"max_disp": 0.1},
            target_acceptance=0.3,
        ),
    ]

    # --- run ---
    driver = MCDriver(engine, moves, ConstantSchedule(300.0))
    for step, move, accepted, chi2, U, V in driver.run(100_000):
        if step % 1_000 == 0:
            print(f"[{step:6d}] chi2={chi2:.3f}  U={U:.3f} eV")

    engine.lmp.command("write_data final.data")

What each block does

Engine. MCEngine starts LAMMPS in-process (cores MPI ranks) and owns the random-number generator. Using it as a context manager guarantees LAMMPS is shut down cleanly. Load a structure from a LAMMPS data file (load_lammps_data) or an ASE Atoms object (load_ase_atoms — a mol-id array, if present, defines molecules for molecule-aware moves).

Constraints. Each constraint scores the current structure. EnthalpyConstraint returns the raw potential energy in eV — the temperature factor \(\beta\) enters exactly once, in the acceptance rule (see Theory). RDFConstraint installs a LAMMPS compute for the partial radial distribution functions \(g_{ij}(r)\); it must be added before any scattering constraint, since those transform its \(g_{ij}(r)\) into \(F(Q)\) or \(G(r)\) (see Total scattering). The χ² against experiment is summed over all scattering constraints.

Data. ScatteringData.from_file reads two- or three-column text (Q F(Q) [sigma]). apply_finite_box_correction convolves the experimental pattern with the simulation box's sinc window so that model and experiment are broadened identically — use the same r_max as the RDF constraint.

Moves. A MoveSpec pairs a move class with a selector that picks which atoms it acts on each step, and tunes the proposal amplitude (max_disp here) towards target_acceptance as the run proceeds.

Driver. MCDriver.run(n) is a generator: it yields (step, move_name, accepted, chi2, U, V) every step, so logging, snapshots and plotting stay in your hands. Nothing happens until you iterate it.

Where next

The worked example extends this scaffold to multiple datasets, molecule moves, volume moves at constant pressure, adaptive tuning bounds, CSV logging, and snapshots.