SimCell¶

SimCell is the main class for building MD simulation boxes. Add layers in the order they should appear along the z axis (you can rotate it later), then call .build() to assemble and pack everything.

from mdinterface import SimCell

Creating a SimCell¶

simbox = SimCell(xysize=[15, 15], verbose=True)

Parameter	Description
`xysize`	Target XY dimensions in Angstroms. Used as a tiling guide for slabs; the actual cell XY is set by the largest fitted slab.
`verbose`	Verbosity: `False`/`0` = quiet, `True`/`1` = INFO, `2` = DEBUG, or a string like `"DEBUG"`.

Adding layers¶

You can add any type of layer to your system with the following methods. Call them in the order you want them along the stacking axis.

`add_slab(species, nlayers)`¶

Adds a periodic solid layer by tiling a unit-cell Specie in XY. The number of repeats is chosen to cover the requested xysize; the actual XY cell is then set by the tiled footprint (an integer multiple of the unit cell). nlayers controls how many unit-cell layers are stacked along Z.

from mdinterface.database import Metal111
gold = Metal111("Au")
simbox.add_slab(gold, nlayers=4)

When multiple slabs are present, each gets distinct atom-type labels (_s0, _s1, …) so their force-field parameters remain independent.

`add_prebuilt(species, nlayers)`¶

Like add_slab, but skips XY tiling — the species positions are used exactly as provided. Intended for structures taken from a previous MD run or geometry optimisation (e.g. a relaxed electrode loaded with hijack and then saved back as a Specie).

simbox.add_prebuilt(relaxed_electrode, nlayers=1)

Any centering or trimming should be done on the species before passing it here.

`add_solvent(solvent, zdim, ...)`¶

Adds a PACKMOL-packed solvent (liquid) layer. Accepts a single Specie, a list of species for mixtures, or None for a solute-only region.

from mdinterface.database import Water, Ion
water = Water(model="ewald")
na    = Ion("Na", ffield="Cheatham")

simbox.add_solvent(
    water,
    zdim=25,
    density=1.0,
    solute=[na],
    nsolute=5,
    solute_pos="right",   # "left", "right", "center", or None (full box)
)

Mixing solvents:

simbox.add_solvent(
    [water, methanol],
    ratio=[3, 1],     # molar ratio
    density=0.95,
    zdim=30,
)

Spatially heterogeneous regions (pockets):

Carve a Sphere, Box, or Cylinder out of the layer and fill it with different content - e.g. an argon gas pocket embedded in bulk water:

from mdinterface.build.regions import Sphere
from mdinterface.database import Argon

argon  = Argon()
bubble = Sphere(center=(7.5, 7.5, 15.0), radius=5.0)

simbox.add_solvent(
    water,
    zdim=30,
    density=1.0,
    regions=[bubble.fill(argon, density=0.8)],
)

A FilledRegion can itself carry regions=[...] for further nested sub-regions, to arbitrary depth - e.g. a neon bubble nested inside the argon pocket:

from mdinterface.build.regions import Sphere
from mdinterface.database import Argon, Neon

argon  = Argon()
neon   = Neon()
pocket = Sphere(center=(7.5, 7.5, 15.0), radius=5.0)
bubble = Sphere(center=(7.5, 7.5, 15.0), radius=1.0)

simbox.add_solvent(
    water,
    zdim=30,
    density=1.0,
    regions=[pocket.fill(argon, density=0.8, regions=[bubble.fill(neon, density=0.1)])],
)

regions isn't limited to one nested pocket - multiple independent, non-nested regions carved out of the same layer are just as common, e.g. two separate gas pockets on opposite sides of the box:

from mdinterface.build.regions import Sphere
from mdinterface.database import Argon, Krypton

argon   = Argon()
krypton = Krypton()

simbox.add_solvent(
    water,
    zdim=30,
    density=1.0,
    regions=[
        Sphere(center=(4.0, 7.5, 15.0), radius=3.0).fill(argon, density=0.8),
        Sphere(center=(11.0, 7.5, 15.0), radius=3.0).fill(krypton, density=0.8),
    ],
)

A region can also confine a solute rather than swap the solvent - e.g. pinning a handful of Na+ ions inside a box-shaped subvolume near one edge instead of letting them roam the whole layer:

from mdinterface.build.regions import Box

corner = Box.from_bounds(0.0, 0.0, 5.0, 6.0, 6.0, 15.0)

simbox.add_solvent(
    water,
    zdim=30,
    density=1.0,
    regions=[corner.fill(solute=[na], nsolute=5)],
)

Region.fill() accepts the same content parameters as add_solvent itself (solvent, solute, nsolute, density, nsolvent, concentration, ratio), minus zdim/xysize since the region defines its own extent. Content parameters are validated at every nesting level: count lists must match the species lists, mixtures require per-species counts or a mixing ratio, and nsolute and concentration are mutually exclusive. conmodel is not supported inside a region.

Both bulk solvent and bulk solute stay outside the top-level regions. Density and concentration use the remaining bulk volume; nested fills likewise subtract their children's volumes. Explicit molecule counts are unchanged. Combining regions with conmodel or solute_pos="center" raises an error because these fixed placements cannot enforce pocket exclusions; use Region.fill(solute=..., nsolute=...) for confined solutes. Sibling regions should not intersect: overlapping bounding boxes produce a warning, and volume accounting assumes disjoint regions.

Every nested region must fit inside its parent's actual shape, not just its bounding box. A region's center can be "random" to choose a placement inside its parent while avoiding sibling bounding boxes. Pass seed to add_solvent to reproduce the region placement; it does not seed PACKMOL's molecular coordinates:

simbox.add_solvent(
    water,
    zdim=30,
    density=1.0,
    regions=[Sphere(center="random", radius=3.0).fill(argon, density=0.8)],
    seed=42,
)

Key parameters:

Parameter	Description
`solvent`	`Specie`, list of `Specie`, or `None` (solute-only)
`zdim`	Z thickness of the layer in Angstroms
`density`	Target density in g/cm³
`nsolvent`	Explicit molecule count (int or list)
`ratio`	Molar ratio for multi-solvent mixtures
`solute`	`Specie` or list of `Specie` to dissolve
`nsolute`	Molecule count per solute species
`solute_pos`	Placement region: `None`, `"center"`, or a `Region` instance. `"left"`/`"right"` strings are deprecated - pass an equivalent `Box` instead (e.g. `Box.from_bounds(...)`)
`regions`	List of `FilledRegion` (from `SomeRegion(...).fill(...)`) for spatially heterogeneous sub-regions within the layer
`seed`	RNG seed for resolving `center="random"` regions, for reproducible placement
`dilate`	Expand the PACKMOL box by this fraction to help convergence
`packmol_tolerance`	PACKMOL distance tolerance in Angstroms (default 2.0)

`add_vacuum(zdim)`¶

Adds an empty vacuum gap.

simbox.add_vacuum(zdim=10)

Building¶

simbox.build(padding=0.5, center=False, stack_axis="z")

Parameter	Description
`padding`	Extra space (Angstroms) added above/below each PACKMOL region
`center`	Default `False`. If `True`, translate the first layer to the box midpoint using its allocated thickness, then wrap individual atoms into the box.
`layered`	Keep per-layer residue numbering instead of merging
`match_cell`	`True` (default): stretch all slabs to the largest XY cell, ensuring a consistent solid/liquid interface. `False`: each slab keeps its natural tiled XY. A `Specie`: lock XY to that species' cell and stretch everything else to match — useful when a pre-relaxed slab or polymer defines the cell.
`hijack`	Replace positions and cell with a prebuilt `ase.Atoms` object
`stack_axis`	Stacking direction: `"z"` (default), `"x"`, or `"y"`

Starting in 2.0.0, center=True places the first layer at the box midpoint. In 1.5.4 it placed that layer across the periodic boundary. For a 100 Å box whose first layer is 20 Å thick, the layer now occupies 40-60 Å instead of 90-100 Å and 0-10 Å. Update coordinate-based selections, restraints, and analysis bins accordingly. center=False is unchanged. Centering does not guarantee that the periodic seam falls in vacuum, and atom-wise wrapping can split other molecules across the boundary. stack_axis moves the centering behavior to the selected axis; hijack overrides the centered coordinates.

Stacking axis¶

By default the system is assembled along Z. Use stack_axis to permute coordinates at the end:

simbox.build(stack_axis="x")   # stacking direction becomes X

This is a coordinate permutation applied after assembly, but the build process itself always runs along Z.

Output¶

Export optional machine-readable structural metadata alongside LAMMPS data:

simbox.write_lammps("data.lammps", metadata="structure.json")

structure.json uses schema mdinterface.structure, version 1. It records the final exported atom, molecule and type IDs, element/species labels, masses, charges, coordinates, species groups, all exported bonded connectivity, cell bounds/tilts, coefficient tokens/annotations, mdinterface version and the SHA-256 checksum of data.lammps. IDs and numerical values follow the written file, including type renumbering and output rounding. Check the schema and checksum before consuming it. Molecule ID zero is a valid exported ID; it is not remapped in JSON.

The optional export requires explicit elements and charges for full style. atomic style records molecule IDs and charges as null because those fields are not exported. write_coeff=False describes only what remains in the resulting file, not omitted topology or coefficients. Metadata paths must differ from the data path and must not already exist. Omitting metadata preserves existing behavior.

Coefficient arrays preserve native tokens because their meaning depends on the force-field style. Units follow the writer's real-unit convention. Full parameter provenance and boundary conditions are null when not known from this export. Do not infer them from an element or filename. Pair style, mixing, constraints and ensemble are deliberately unset: a simulation application must select and validate them. In particular, O-H connectivity is structural information; SHAKE is a protocol decision. No cluster paths or scheduler configuration are stored.

# LAMMPS data file
simbox.write_lammps("data.lammps", atom_style="full", write_coeff=True)

# ASE Atoms object
atoms = simbox.to_ase()

# MDAnalysis Universe
universe = simbox.universe

GROMACS output (experimental)¶

Warning

GROMACS output assumes harmonic bonds and angles, four-term OPLS proper torsions, CVFF impropers, geometric Lennard-Jones mixing, and 0.5 scaling of both LJ and Coulomb 1-4 interactions. It is not a general force-field converter. Constraints, rigid-water settings, electrostatics methods, cutoffs, and simulation protocols must be chosen separately for your model.

write_gromacs() produces these files in the directory of your choice:

{prefix}.gro - atomic coordinates and box vectors, with standard GRO precision (0.001 nm).
{resname}.itp - per-species molecule definitions, bonded terms, and explicit 1-4 pairs derived from shortest bond-graph distances.
{prefix}.top - shared atom types before all molecule includes, OPLS defaults, and consecutive molecule counts matching the GRO coordinate order.

simbox.write_gromacs(prefix="system", outdir="gromacs_input/")

The exporter validates force-field completeness before creating files. Missing or nonfinite parameters, unsupported fifth torsion coefficients, invalid CVFF parameters, and conflicting atom types raise ValueError. Species sharing a residue name must have identical atom order, charges, masses, and topology parameters; use distinct name= values for different models. Assembled residues must match their species definitions, including charges and atom types.

Individual ITP files include atom types by default:

water.write_gromacs_itp("water.itp")

For manual multi-species assembly, keep molecule files separate from shared atom-type definitions:

from mdinterface.io import write_gromacs_top

species = [water, sodium, chloride]
itp_files = []
for specie in species:
    filename = f"{specie.resname}.itp"
    specie.write_gromacs_itp(filename, include_atomtypes=False)
    itp_files.append(filename)
write_gromacs_top(universe, itp_files, species=species)

GROMACS preprocessing tests cover mixtures, chains, and rings. Single-point energy and force comparisons against LAMMPS cover harmonic terms, OPLS torsions, scaled 1-4 interactions, and CVFF impropers, plus water, benzene, ethanol, methane, dimers, mixed molecular systems, and a packed water/NaCl box. Strict packed-system comparisons use double-precision GROMACS to limit rounding of small force components. These checks do not validate arbitrary force fields or production dynamics. Compare engines at identical coordinates because standard GRO rounding can change bonded energies.

Verbosity¶

import mdinterface
mdinterface.set_verbosity("DEBUG")   # package-wide
mdinterface.set_verbosity(1)         # INFO
mdinterface.set_verbosity(0)         # quiet

See the Logging guide for details.

Validate classical exports¶

write_lammps(write_coeff=True) validates parameters for species present in the assembled system before opening the output file. Missing pair or bonded parameters and nonfinite charges or coefficients raise an error. This distinguishes a chemically valid structure from a parameterized classical model. Explicit zero-valued coefficients are retained. Molecular completeness uses stored chemistry or bonded topology; it does not invent bonded interactions for nonbonded solids or infer universally required improper terms.

Use write_lammps("data.lammps", expected_charge=0.0) when the assembled system must be neutral, or provide a nonzero expected total for an intentionally charged system. The absolute tolerance is 1e-5 elementary charges. Omitting expected_charge imposes no neutrality condition. Geometry-only workflows can continue to use to_ase(); exports without coefficient blocks do not require complete force-field coefficients.

Imported LAMMPS bond connectivity is retained when coordinates are replaced from a trajectory, reordered with ASE, or repeated through Specie.repeat(), rather than being inferred again from interatomic distances.