Build a periodic solvent box¶

mlmd.solvate places one solute, an optional single counterion, and an integer number of identical solvent molecules in an orthorhombic periodic cell. It chooses the solvent count from total mass divided by cell volume, then uses Packmol to avoid close intermolecular contacts. The result is a starting geometry for relaxation and equilibration, not an equilibrated liquid.

This tutorial starts with a reproducible packing plan, runs Packmol locally, and passes the resulting state to the simulation workflow. The browser examples inspect configuration and counting only. ASE and the Packmol executable run in your local Python environment.

1. Install the packing dependencies¶

From a checkout of matmalib, install the small solvation extra:

python -m pip install -e '.[solvation]'

This installs ASE and Packmol wheels. The package requires Packmol 21.2.1 or newer; the underlying periodic input format requires at least 20.15.0. If you installed Packmol separately, set packmol="/absolute/path/to/packmol". The default also checks for an executable alongside the active Python interpreter.

Solvation does not evaluate a potential, so it does not download AIMNet weights or require a GPU. Later AIMNet simulations need the MLMD dependencies described in the simulation tutorial.

2. Prepare and inspect a packing plan¶

Save this as prepare_box.py and run python prepare_box.py from a new working directory. It generates methane and water templates with ASE, then prepares a 15 Å box at a requested total density of 1.0 g/mL:

from ase.build import molecule
from matmalib import mlmd

simulation = mlmd()
solute = molecule("CH4")
solvent = molecule("H2O")

config = mlmd.SolvationConfig(
    cell=15.0,
    density_g_ml=1.0,
    tolerance=2.0,
    seed=12345,
    volume_samples=20000,
    work_dir="packing_plan",
    dry_run=True,
)
report = simulation.solvate(solute, solvent, config)

print("Status:", report.status)
for entry in report.molecules:
    print(entry["role"], entry["formula"], entry["count"])
print("Atoms:", report.atom_count)
print("Actual density:", report.achieved_density_g_ml, "g/mL")
print("Warnings:", report.warnings)
print("Packmol input:", report.files["packmol.inp"])

Expect one methane and 112 water molecules, or 341 atoms, using ASE's ordinary atomic masses. The achieved density is about 1.0006 g/mL. Rounding the solvent count explains the small difference from the target. A dry run returns Data, also available as simulation.data, with status "prepared"; there is no packed geometry yet. It requires work_dir, and you must omit output and output_format.

The directory contains normalized molecule templates, packmol.inp, and summary.json. Inspect the molecule counts, cell, tolerances, and warnings before packing. Dry runs work even when Packmol itself is unavailable, but still require ASE.

Inputs can instead be XYZ, extended XYZ, or PDB paths, matmalib State/Glycan objects, or ASE Atoms. Provide one molecule per component. A PDB must contain one model with no alternate conformations and full or blank occupancies; the reader does not infer missing atoms or add hydrogens. Template cell records are discarded because the new cell comes from SolvationConfig.cell.

How the solvent count is chosen¶

The implementation calculates the target mass in atomic mass units, subtracts the solute and counterion masses, then rounds to the nearest positive integer number of solvent molecules:

V_mL = Lx * Ly * Lz * 1e-24
target_mass_u = density_g_ml * V_mL * N_A
N_solvent = floor((target_mass_u - fixed_mass_u) / solvent_mass_u + 0.5)
actual_density = (fixed_mass_u + N_solvent * solvent_mass_u) / (N_A * V_mL)

The volume is the whole cell, not a solvent-accessible volume obtained by subtracting the solute's van der Waals volume. Custom or isotope masses on input states/ASE atoms are retained and affect the count.

Change the box length or density in this browser exercise to see the integer jumps:

import math
from matmalib import mlmd

config = mlmd.SolvationConfig(cell=15.0, density_g_ml=1.0)
NA = 6.02214076e23
methane_mass, water_mass = 16.043, 18.015
volume_ml = config.cell ** 3 * 1e-24
count = math.floor(
    (config.density_g_ml * volume_ml * NA - methane_mass) / water_mass + 0.5
)
if count < 1:
    print("This example needs a larger box or greater target density.")
else:
    density = (methane_mass + count * water_mass) / (NA * volume_ml)
    print("Water molecules:", count)
    print("Total atoms:", 5 + 3 * count)
    print("Achieved density:", round(density, 6), "g/mL")

This cell illustrates the counting formula with stated masses; it does not run a packing calculation. The local API obtains masses from your actual inputs.

3. Run Packmol and keep the final state¶

Append this to prepare_box.py after inspecting the plan:

from dataclasses import replace

packed = simulation.solvate(
    solute,
    solvent,
    replace(config, dry_run=False, work_dir="packing_run"),
    output="solvated.xyz",
    output_format="state",
)
packing = simulation.data
print("Status:", packing.status)
print("Periodic axes:", packed.pbc)
print("Cell:", packed.cell)
print("Atoms:", len(packed.atoms))

The new work directory preserves the plan for comparison. output_format="state" writes extended XYZ and returns a State. Without output, the method also returns a State; with output="solvated.xyz" alone, it returns the output Path. The final state includes masses, charge, multiplicity, cell, full three-dimensional PBC, and solvation_report.

Packmol holds the solute's geometric center at the cell center and, when present, holds the counterion at a nearby nonoverlapping location. It translates and rotates solvent molecules while preserving each template's internal geometry. These are packing constraints only: they do not freeze the solute or counterion in subsequent dynamics. MD constraints belong in MinimizationConfig, EquilibrationConfig, or MDConfig.

After Packmol finishes, matmalib verifies atom order, count, elements, finite coordinates, bounds, fixed placements, and intermolecular distances including periodic images. A close contact triggers one retry with face padding of tolerance / 2; the initial files are retained with an initial_ prefix. A second failure raises and leaves a report/log for inspection. Packmol's own success message alone does not bypass these checks.

The generated packmol.inp can also be run manually from its directory using packmol < packmol.inp, as described in the Packmol user guide. That command produces the raw packed.pdb; use simulation.solvate for matmalib's final validation, metadata, and cell.xyz output.

4. Add a counterion and control electronic state¶

The current builder supports exactly one optional counterion molecule, which may contain one or several atoms. It does not infer charge balance, add multiple ions, choose a salt concentration, or support multiple solvent species.

from ase import Atoms
from matmalib import mlmd

# Replace these paths and charge with your prepared molecular inputs.
simulation = mlmd()
ion = Atoms("Cl", positions=[[0.0, 0.0, 0.0]])
ion.info["charge"] = -1.0

plan = simulation.solvate(
    "cation.xyz", "water.xyz",
    mlmd.SolvationConfig(
        cell=(25.0, 25.0, 30.0),
        density_g_ml=1.0,
        counterion=ion,
        counterion_gap=0.5,
        charge=0.0,
        multiplicity=1,
        dry_run=True,
        work_dir="ionic_plan",
    ),
)
print(plan.fixed_centers_angstrom)
print(plan.molecules)

counterion_gap is the minimum requested surface-to-surface van der Waals gap between the solute and counterion, in Å. The placement must also satisfy the generic atomic tolerance, fit inside the cell, and avoid periodic overlaps. Neither input chemistry nor charge balance is inferred from a filename.

For completed packing, charge=None sums component charge metadata: one solute, one counterion if present, and N_solvent solvent charges. An explicit charge overrides that sum. Set the combined system's multiplicity directly; molecular multiplicities are not additive. Ordinary XYZ files frequently lack this metadata, so supply explicit values when needed.

5. Choose controls and interpret diagnostics¶

Control	Default	Effect
`cell`	`30.0` Å	One cubic length or `(Lx, Ly, Lz)`; no triclinic packing. Lengths are rounded to 0.001 Å for PDB consistency.
`density_g_ml`	`1.0`	Total target density; must leave enough target mass for at least one solvent molecule.
`tolerance`	`2.0` Å	Minimum intermolecular atom-center separation used by packing and collision checks; not a van der Waals radius.
`counterion_gap`	`0.5` Å	Additional solute/counterion surface gap.
`nloop` / `timeout`	`1000` / `600` s	Packmol optimization loops and subprocess time limit.
`seed`	`12345`	Positive integer used for packing and volume sampling, at most 2147483647.
`volume_samples`	`200000`	Samples per template for the diagnostic sphere-union volume estimate; at least 2.
`vdw_radii`	`{}`	Element-to-radius overrides in Å; otherwise use ASE's van der Waals radii.
`work_dir`	`None`	Preserve inputs, logs, cells, and JSON here; otherwise use temporary storage. Required for a dry run.
`overwrite`	`False`	Protect output and previously generated work-directory files.

Each molecule's volume is estimated as the union of its atomic van der Waals spheres, including overlaps; single-atom volumes are analytic. vdw_volume_standard_error_angstrom3 measures Monte Carlo sampling error, not uncertainty in the radii or a liquid's thermodynamic volume. Increasing volume_samples refines this diagnostic; it does not change the mass-based solvent count.

estimated_vdw_packing_fraction is the sum of those molecular volumes divided by cell volume. A value above one adds a warning that packing may be difficult. It is not a measured occupied fraction of the final periodic configuration. The chosen radii also affect whether the centered solute and counterion fit, so missing radii must be supplied, for example vdw_radii={"Cl": 1.8}.

PDB intermediates limit the box to 99,999 atoms and cell lengths to the range that rounds to 0.001–9999.999 Å. Reusing a work directory requires overwrite=True; only known generated files are replaced, and unrelated files are preserved. Avoid using input templates at the generated solute.pdb, solvent.pdb, or counterion.pdb paths in that directory.

Outputs to retain¶

File or object	Purpose
`summary.json`, `simulation.data.summary`	Counts, masses, target/achieved density, cell, volume estimates, warnings, status, and failure details.
`solute.pdb`, `solvent.pdb`, optional `counterion.pdb`	Centered templates used for packing.
`packmol.inp`, `packmol.log`, `packed.pdb`	Reproducible input, subprocess output, and raw packing result.
`cell.xyz`	Final extended XYZ with periodic cell, masses, and electronic-state metadata.
`cell.pdb`	Final geometry/cell for inspection; prefer extended XYZ to preserve simulation metadata.
`State.solvation_report`	Packing report attached to the returned state.

Save packing = simulation.data before the next operation: the interface replaces .data each time. If packing raises, inspect .data.status, .data.error, and the preserved work-directory log. Increase the cell or reduce target density when the fixed molecule cannot fit or the requested packing is too dense; increasing iteration limits does not make an impossible geometry feasible.

6. Relax and equilibrate the cell¶

With the MLMD environment installed, continue from the completed packed state:

from matmalib import mlmd

simulation = mlmd()  # AIMNet defaults; resolves model weights when first used
relaxed = simulation.minimize("solvated.xyz", mlmd.MinimizationConfig())
ready = simulation.equilibrate(
    relaxed,
    mlmd.EquilibrationConfig(
        steps=2000,
        timestep_fs=0.5,
        ensemble="NPT",
        pressure_bar=1.0,
        temperature_K=300.0,
        log_file="box_equilibration.log",
    ),
    output="equilibrated.xyz",
    output_format="state",
)
print(simulation.data.summary)

This is a workflow example, not a convergence criterion. Evaluate whether the potential describes your solute, solvent, ions, and thermodynamic conditions; monitor density, temperature, forces, and structural relaxation. NPT needs a stress-capable calculator and does not support fixed atoms, SHAKE, or persistent center-of-mass constraints. NVT Langevin dynamics supports those constraints when appropriate. Continue with production MD or metadynamics after assessing equilibration.