Worked example

A realistic refinement of an amorphous framework material (a-ZIF-like: C, H, N, Zn) against three datasets simultaneously — X-ray \(F_X(Q)\), neutron \(F_N(Q)\), and neutron \(G_N(r)\) — with an ACE potential, molecule moves, and volume moves at ambient pressure. This is a trimmed-down version of a production script; every piece is optional, and the quickstart scaffold underneath is unchanged.

Setup

from mosaic_materials.engine import MCEngine
from mosaic_materials.constraints import (
    EnthalpyConstraint,
    RDFConstraint,
    XrayStructureFactorConstraint,
    NeutronStructureFactorConstraint,
    NeutronTotalRDFConstraint,
)
from mosaic_materials.data import ScatteringData
from mosaic_materials.monte_carlo.adapt import MoveSpec
from mosaic_materials.monte_carlo.callbacks import (
    summary_logger,
    move_magnitude_tuning_logger,
)
from mosaic_materials.monte_carlo.driver import MCDriver
from mosaic_materials.monte_carlo.schedules import ConstantSchedule
from mosaic_materials.moves import (
    TranslateMove,
    TranslateRotateMove,
    VolumeMove,
)

TEMPERATURE = 593.15   # K
DEUTERATION = 5.6 / 6  # fraction of exchangeable H that is D

engine = MCEngine(cores=32, seed=42)
engine.load_lammps_data("azif_4x4x4.data")
engine.add_force_field(
    "pace recursive", "* * zif-ace-24.yace C H N Zn"
)

Constraints

The energy constraint and the RDF come first; the RDF's r_max sets the real-space window every scattering constraint sees, so it is reused for the finite-box correction of the reciprocal-space data.

engine.add_constraint(EnthalpyConstraint(temperature_Kelvin=TEMPERATURE))

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

Each dataset becomes one constraint. sigmas weights the per-point χ² (smaller σ → stronger constraint); the neutron constraints take the deuteration ratio, which mixes the ¹H and ²H coherent scattering lengths (see Total scattering).

xfq = ScatteringData.from_file("azif.xfq", sigmas=0.01)
xfq.apply_finite_box_correction(r_max=r_max, inplace=True)
engine.add_constraint(
    XrayStructureFactorConstraint.from_scattering_data("c_x_fq", xfq)
)

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

ngr = ScatteringData.from_file("azif.ngr", sigmas=0.03)
engine.add_constraint(
    NeutronTotalRDFConstraint.from_scattering_data(
        "c_n_gr", ngr, deuteration_ratio=DEUTERATION
    )
)

X-ray scale

XrayStructureFactorConstraint applies an analytical scale minimisation by default (minimise_differences=True), since X-ray patterns often carry an arbitrary overall scale. The neutron constraints compare absolutely by default.

Moves

Three move types, each with adaptive amplitude tuning: every interval trials the tuned parameters are nudged towards target_acceptance, within [min_values, max_values]. Selectors decide what each move acts on: a random atom, a whole random molecule (the mol-ids in the data file), or nothing (volume moves act on the box).

move_specs = [
    MoveSpec(
        name="translate",
        move_class=TranslateMove,
        selector=lambda system, rng: [system.pick_random_atom(rng)],
        kwargs={"max_disp": 0.1},
        target_acceptance=0.3,
        adjustment_factor=1.05,
        interval=200,
        min_values={"max_disp": 0.01},
        max_values={"max_disp": 0.50},
    ),
    MoveSpec(
        name="translate_rotate",
        move_class=TranslateRotateMove,
        selector=lambda system, rng: system.pick_random_molecule(rng)[1],
        kwargs={"max_disp": 0.1, "max_angle": 2.5, "pivot": "com"},
        target_acceptance=0.3,
        adjustment_factor=1.05,
        interval=200,
        min_values={"max_disp": 0.01, "max_angle": 0.1},
        max_values={"max_disp": 0.50, "max_angle": 3.0},
    ),
    MoveSpec(
        name="volume",
        move_class=VolumeMove,
        selector=lambda system, rng: [],
        kwargs={"max_dlnv": 1e-4, "rigid_molecules": True},
        target_acceptance=0.3,
        adjustment_factor=1.05,
        interval=200,
        min_values={"max_dlnv": 1e-5},
        max_values={"max_dlnv": 1e-3},
    ),
]

rigid_molecules=True scales molecular centres of mass with the box and keeps intramolecular geometry fixed — the acceptance rule then counts molecules, not atoms, in its volume-entropy term (see Theory).

Driver, logging, and the run loop

Volume moves sample the isothermal–isobaric ensemble, so a pressure schedule is supplied (in GPa; ambient ≈ 0).

driver = MCDriver(
    engine,
    move_specs,
    temperature_schedule=ConstantSchedule(TEMPERATURE),
    pressure_schedule=ConstantSchedule(0.0),
    callbacks=[
        summary_logger("summary.csv", interval=100),
        move_magnitude_tuning_logger(move_specs),
    ],
)

for step, move, accepted, chi2, U, V in driver.run(1_000_000):
    if step % 10_000 == 0:
        atoms = engine.to_ase_atoms()
        atoms.info.update(step=step, chi2=chi2, U=U)
        atoms.write("snapshots.xyz", format="extxyz", append=True)

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

summary_logger appends a CSV row every 100 steps (step, temperature, pressure, volume, acceptance counts, per-constraint costs, timing); move_magnitude_tuning_logger logs every amplitude adjustment.

Inspecting the fit

The driver keeps the last accepted SimulationState, which caches every constraint's computed pattern under its label:

import matplotlib.pyplot as plt

ax = driver.state.plot_scattering("c_n_fq", x=nfq.x, y=nfq.y)
ax.set_xlabel(r"$Q$ (Å$^{-1}$)")
ax.set_ylabel(r"$F_{\mathrm{N}}(Q)$")
plt.savefig("fit_nfq.pdf")

Restarting

driver.run(n_steps, start=...) offsets the reported step and the schedules, so a continuation behaves as one long run:

for step, *rest in driver.run(500_000, start=1_000_000):
    ...

Reload the structure from final.data (or the last snapshot) when constructing the engine, and re-use the same constraint and move setup.

Finally, close the engine (or construct it in a with block, as in the quickstart):

engine.close()