SimCell

mdinterface.build.builder.SimCell

Fluent builder for assembling layered simulation boxes.

Layers are stacked in the order they are added. Slabs (solid interfaces), solvent regions, and vacuum gaps can be freely combined.

Example

simbox = SimCell(xysize=[15, 15]) simbox.add_slab(gold, nlayers=3) simbox.add_solvent(water, solute=[na, cl], nsolute=[5, 5], zdim=25, density=1.0) simbox.add_slab(gold, nlayers=3) simbox.add_vacuum(zdim=5) simbox.build() simbox.write_lammps("data.lammps")

Source code in mdinterface/build/builder.py

class SimCell:
    """
    Fluent builder for assembling layered simulation boxes.

    Layers are stacked in the order they are added. Slabs (solid interfaces),
    solvent regions, and vacuum gaps can be freely combined.

    Example
    -------
    >>> simbox = SimCell(xysize=[15, 15])
    >>> simbox.add_slab(gold, nlayers=3)
    >>> simbox.add_solvent(water, solute=[na, cl], nsolute=[5, 5], zdim=25, density=1.0)
    >>> simbox.add_slab(gold, nlayers=3)
    >>> simbox.add_vacuum(zdim=5)
    >>> simbox.build()
    >>> simbox.write_lammps("data.lammps")
    """

    def __init__(
        self,
        xysize: Union[List[float], Tuple[float, float]],
        verbose: Union[None, bool, int, str] = None,
    ) -> None:
        """
        Parameters
        ----------
        xysize : list or tuple of float
            XY dimensions of the simulation box in Å, e.g. ``[15.0, 15.0]``.
        verbose : None, bool, int, or str, optional
            Controls package-wide log verbosity via :func:`set_verbosity`.
            ``None`` (default) leaves the current level unchanged.

            Integer scale:

            - ``0`` / ``False`` — WARNING (quiet)
            - ``1`` / ``True``  — INFO  (normal)
            - ``2``             — DEBUG (detailed)
            - ``3``, ``4``, …  — DEBUG (same as 2 for values < 10)

            Integers ≥ 10 are treated as raw ``logging`` level constants.
            Strings are resolved by name (``"DEBUG"``, ``"INFO"``, …).
        """
        if verbose is not None:
            set_verbosity(verbose)

        xsize, ysize = self._validate_xysize(xysize)
        self._xsize = xsize
        self._ysize = ysize
        self._layers: List[Compartment] = []
        self._all_species: List[Any] = []
        self._universe: Optional[mda.Universe] = None

        log_banner(logger, "mdinterface :: SimCell", f"version {_get_mdi_version()}")
        logger.info("  >> xysize: %.2f x %.2f Å", xsize, ysize)

    # ------------------------------------------------------------------
    # Layer-adding methods (chainable)
    # ------------------------------------------------------------------

    def add_slab(self, species: Any, nlayers: int = 1) -> None:
        """
        Add a solid-interface layer (slab).

        Parameters
        ----------
        species : Specie
            The unit-cell species to tile into a slab.
        nlayers : int
            Number of unit-cell layers to stack along Z.
        """
        slab_sp = species.copy()
        slab_idx = sum(1 for lay in self._layers if isinstance(lay, SlabCompartment))
        suffix = f"_s{slab_idx}"
        for atom in slab_sp._stype:
            atom.set_label(atom.label + suffix)
        self._register(slab_sp)
        self._layers.append(SlabCompartment(species=slab_sp, nlayers=nlayers, label="slab"))
        logger.info("  + slab     %s,  %d layer(s)",
                    getattr(slab_sp, "resname", "?"), nlayers)

    def add_prebuilt(self, species: Any, nlayers: int = 1) -> None:
        """
        Add a pre-built layer whose atomic positions come from a prior MD run.

        Semantically equivalent to :meth:`add_slab` but intended for species
        whose positions have already been set (e.g. via
        :meth:`~mdinterface.core.specie.Specie.update_positions`).  No XY
        tiling is performed; the species cell is used as-is.

        Any pre-processing of the positions (centering, trimming, …) should
        be done on the species before passing it here.

        Parameters
        ----------
        species : Specie or Polymer
            Species with positions already set from a prior trajectory.
        nlayers : int
            Number of repeat units to stack along Z (usually 1).
        """
        slab_sp = species.copy()
        slab_idx = sum(1 for lay in self._layers if isinstance(lay, SlabCompartment))
        suffix = f"_s{slab_idx}"
        for atom in slab_sp._stype:
            atom.set_label(atom.label + suffix)
        self._register(slab_sp)
        self._layers.append(SlabCompartment(species=slab_sp, nlayers=nlayers, label="prebuilt"))
        logger.info("  + prebuilt %s,  %d layer(s)",
                    getattr(slab_sp, "resname", "?"), nlayers)

    def add_solvent(
        self,
        solvent: Optional[Any] = None,
        solute: Optional[List[Any]] = None,
        nsolute: Optional[Union[int, List[int]]] = None,
        zdim: Optional[float] = None,
        density: Optional[float] = None,
        nsolvent: Optional[Union[int, List[int]]] = None,
        concentration: Optional[float] = None,
        conmodel: Optional[dict] = None,
        solute_pos: Optional[str] = None,
        dilate: float = 1.0,
        packmol_tolerance: float = 2.0,
        ratio: Optional[List[float]] = None,
        regions: Optional[List[FilledRegion]] = None,
        seed: Optional[int] = None,
    ) -> None:
        """
        Add a solvent (liquid) layer, optionally with dissolved species.

        Parameters
        ----------
        solvent : Specie or list of Specie or None
            Solvent molecule(s). Pass ``None`` for a solute-only region.
        solute : list of Specie, optional
            Species to dissolve (ions, neutral molecules, …).
        nsolute : int or list of int, optional
            Number of each solute species (alternative to *concentration*).
        zdim : float
            Thickness of this region in Angstroms.
        density : float, optional
            Solvent density in g/cm³ (alternative to *nsolvent*).
        nsolvent : int or list of int, optional
            Explicit number of solvent molecules (alternative to *density*).
        concentration : float, optional
            Solute concentration in Molar (alternative to *nsolute*).
        conmodel : dict, optional
            Spatially varying concentration model.
        solute_pos : str, optional
            Solute placement strategy:

            - ``None`` / ``"packmol"`` — PACKMOL random placement in the full
              box (default).
            - ``"left"``   — PACKMOL random placement in the left half
              (z ≤ zdim/2).
            - ``"right"``  — PACKMOL random placement in the right half
              (z ≥ zdim/2).
            - ``"center"`` — each molecule fixed at the box centre.
        dilate : float, optional
            Dilation factor > 1 for concentrated systems. PACKMOL will pack
            into a box ``dilate`` times taller at ``density / dilate``,
            keeping the number of solvent molecules identical. The larger box
            gives PACKMOL more breathing room; subsequent NpT equilibration
            compresses the system to the correct density. Default is ``1.0``
            (no dilation).
        packmol_tolerance : float, optional
            Minimum distance (Å) between atoms of different molecules during
            PACKMOL packing. Reduce from the default of ``2.0`` for very
            concentrated systems if PACKMOL cannot find a valid packing.
        regions : list of FilledRegion, optional
            Spatially-heterogeneous sub-regions within this layer (e.g. a gas
            pocket embedded in the bulk solvent). Build with
            ``SomeRegion(...).fill(solvent=..., density=...)`` — see
            :mod:`mdinterface.build.regions`. A region's ``center`` may be
            ``"random"`` to have a non-overlapping placement chosen
            automatically within its parent volume.
            Bulk solvent and solute exclude these regions; density and
            concentration use the remaining volume. Cannot be combined with
            ``conmodel`` or ``solute_pos="center"``.
        seed : int, optional
            Seed for the RNG used to resolve ``center="random"`` regions.
            Pass the same seed to reproduce an identical placement.
        """
        if zdim is None:
            raise ValueError("'zdim' is required for add_solvent()")
        if dilate <= 0:
            raise ValueError(f"'dilate' must be positive, got {dilate}")
        if packmol_tolerance <= 0:
            raise ValueError(f"'packmol_tolerance' must be positive, got {packmol_tolerance}")

        # Normalise solvent to a list; copy each species.
        if solvent is None:
            solv_copies = []
        elif isinstance(solvent, (list, tuple)):
            solv_copies = [s.copy() for s in solvent]
        else:
            solv_copies = [solvent.copy()]

        solute_copies = [sp.copy() for sp in (solute or [])]

        region_copies = [self._copy_filled_region(fr) for fr in (regions or [])]
        region_species: List[Any] = []
        for fr in region_copies:
            region_species.extend(_collect_region_species(fr))

        self._register(*solv_copies, *solute_copies, *region_species)

        self._layers.append(SolventCompartment(
            solvent=solv_copies,
            solute=solute_copies,
            nsolute=nsolute,
            zdim=zdim,
            density=density,
            nsolvent=nsolvent,
            concentration=concentration,
            conmodel=conmodel,
            solute_pos=solute_pos,
            dilate=dilate,
            packmol_tolerance=packmol_tolerance,
            ratio=ratio,
            regions=region_copies,
            seed=seed,
        ))
        solv_str = "+".join(getattr(s, "resname", "?") for s in solv_copies) or "ions"
        rho_str  = (f"ρ={density:.2f} g/cm³" if density is not None
                    else f"nsolvent={nsolvent}" if nsolvent is not None else "density=?")
        solute_info = ""
        if solute_copies:
            sol_names = "+".join(getattr(s, "resname", "?") for s in solute_copies)
            solute_info = f",  solute: {sol_names} (n={nsolute})"
        region_info = f",  {len(region_copies)} region(s)" if region_copies else ""
        logger.info("  + solvent  %s,  zdim=%.1f Å,  %s%s%s",
                    solv_str, zdim, rho_str, solute_info, region_info)

    def add_vacuum(self, zdim: float = 0.0) -> None:
        """
        Add an empty vacuum gap.

        Parameters
        ----------
        zdim : float
            Thickness of the vacuum gap in Angstroms.
        """
        self._layers.append(VacuumCompartment(zdim=zdim))
        logger.info("  + vacuum   zdim=%.1f Å", zdim)

    # ------------------------------------------------------------------
    # Build
    # ------------------------------------------------------------------

    def build(
        self,
        padding: float = 0.5,
        center: bool = False,
        layered: bool = False,
        match_cell: Union[bool, Any] = True,
        hijack: Optional[ase.Atoms] = None,
        stack_axis: str = "z",
    ) -> None:
        """
        Assemble all layers into a simulation box.

        Parameters
        ----------
        padding : float
            Spacing (Å) inserted between adjacent layers.
        center : bool
            If True, shift the system so the center of the first layer falls
            in the middle of the box (z=zdim/2), using its allocated thickness,
            then wrap individual atoms. The periodic seam falls in the
            opposite part of the system, which need not be vacuum.
            Before version 2.0.0, True centered the first layer on the
            periodic boundary instead. Default is False.
        layered : bool
            Assign distinct molecule indices to each slab layer for LAMMPS.
        match_cell : bool or Specie
            Controls XY cell matching:

            - ``True`` — scale all slabs to the largest fitted XY cell
              (default). Ensures the solid/liquid interface is well-defined
              with no XY mismatch between layers.
            - ``False`` — no scaling; each slab keeps its natural tiled XY.
            - *Specie* — lock XY to that species' cell and stretch all other
              slabs to match. The reference species is left unscaled. Useful
              when a polymer or pre-relaxed slab should define the cell and
              electrodes should conform to it.
        hijack : ase.Atoms, optional
            After assembly, override both positions and cell dimensions with
            this external ``ase.Atoms`` object. Useful when you have a
            pre-relaxed structure you want to map the topology onto.
        stack_axis : str, optional
            Axis along which layers are stacked in the output file.
            Assembly always happens along Z; a final coordinate permutation
            reorients the box. Accepted values: ``"x"``, ``"y"``, ``"z"``
            (default ``"z"``).

        """
        stack_axis = stack_axis.lower()
        if stack_axis not in ("x", "y", "z"):
            raise ValueError(
                f"stack_axis must be 'x', 'y', or 'z', got {stack_axis!r}"
            )

        log_header(logger, "Build")
        logger.info("  >> %d layer(s)  |  inter-layer padding: %.2f Å",
                    len(self._layers), padding)

        self._update_topology_indexes()

        xsize, ysize, cell_ref, do_match = self._resolve_xy(match_cell)
        xsize, ysize = self._fit_slabs(xsize, ysize, cell_ref)

        all_sp_univs = [sp.to_universe() for sp in self._all_species]
        system, zdim, first_layer_zdim = self._stack_layers(
            xsize, ysize, all_sp_univs, padding, layered, do_match
        )

        system.dimensions = [xsize, ysize, zdim] + [90, 90, 90]
        res_counts = Counter(res.resname for res in system.residues)

        log_header(logger, "Done")
        logger.info("  >> %d atoms  |  %.3f x %.3f x %.3f Å",
                    len(system.atoms), xsize, ysize, zdim)
        logger.info("  >> %d layers  |  %d residues",
                    len(self._layers), len(system.residues))
        parts = [f"{name} ({n} mol)" for name, n in res_counts.items()]
        for i in range(0, len(parts), 3):
            logger.info("  >> %s", ",  ".join(parts[i:i + 3]))

        if center:
            shift = (zdim - first_layer_zdim) / 2
            system.atoms.translate([0, 0, shift])
            _ = system.atoms.wrap()
            logger.info("  >> system centered on first layer (shift %.2f Å)", shift)

        if hijack is not None:
            system.dimensions = hijack.get_cell_lengths_and_angles()
            system.atoms.positions = hijack.get_positions()
            logger.info("  >> positions and cell overridden by hijack ase.Atoms")

        if stack_axis != "z":
            system, xsize, ysize = self._apply_stack_axis(system, xsize, ysize, zdim, stack_axis)

        self._universe = system
        self._xsize = xsize
        self._ysize = ysize

        return

    # Output
    def write_lammps(
        self,
        filename: str = "data.lammps",
        atom_style: str = "full",
        write_coeff: bool = True,
        metadata: Optional[str] = None,
        expected_charge: Optional[float] = None,
    ) -> None:
        """
        Write a LAMMPS data file (and optional force-field coefficients).

        .. note::
            This method is only needed for classical MD with LAMMPS.  If you
            are using the assembled structure for AIMD, ML-MD, or any other
            workflow, use :meth:`to_ase` or :attr:`universe` instead — no
            force-field parameters are required for those paths.

        Parameters
        ----------
        filename : str
            Output file path.
        atom_style : str
            LAMMPS atom style (``"full"`` or ``"atomic"``).
        write_coeff : bool
            Whether to write force-field coefficient blocks.
        metadata : str, optional
            Path for a schema-versioned JSON description of the final export.
            Includes exported IDs, species, topology, coefficients and checksum.
            Requires elements and, for full style, charges. Does not select MD settings.
        expected_charge : float, optional
            Required total charge in elementary-charge units, checked with an
            absolute tolerance of 1e-5. Omit for systems with no specified target.
            Nonfinite charges and incomplete coefficients are always rejected
            when writing force-field coefficients.

        """
        if self._universe is None:
            raise RuntimeError("Call build() before write_lammps().")
        if metadata is not None:
            if Path(metadata).resolve() == Path(filename).resolve():
                raise ValueError("Metadata and LAMMPS data paths must differ.")
            if Path(metadata).exists():
                raise FileExistsError(metadata)
            if atom_style not in {"full", "atomic"}:
                raise ValueError("Metadata supports full or atomic exports.")
            if not hasattr(self._universe.atoms, "elements"):
                raise ValueError("Metadata requires explicit atom elements.")
            if atom_style == "full" and not hasattr(self._universe.atoms, "charges"):
                raise ValueError("Full-style metadata requires atom charges.")

        if hasattr(self._universe.atoms, "charges"):
            charges = self._universe.atoms.charges
            if not np.isfinite(charges).all():
                raise ValueError("Cannot export nonfinite partial charges.")
            if expected_charge is not None:
                if not np.isfinite(expected_charge) or not np.isclose(charges.sum(), expected_charge, atol=1e-5, rtol=0):
                    raise ValueError(f"Total charge {charges.sum():.8f} does not match expected_charge={expected_charge}.")
        elif expected_charge is not None:
            raise ValueError("expected_charge requires atomic charge data.")
        if write_coeff:
            active_types = set(self._universe.atoms.types)
            for specie in self._all_species:
                if active_types.intersection(specie.get_atom_types()):
                    specie.validate_force_field()

        log_header(logger, "Output")
        logger.info("  >> LAMMPS data file: %s  (style=%s,  coeff=%s)", filename, atom_style, write_coeff)
        system = self._universe.copy()

        if not write_coeff:
            for attr in ["bonds", "angles", "dihedrals", "impropers"]:
                try:
                    system.del_TopologyAttr(attr)
                except Exception:
                    pass

        with DATAWriter(filename) as dt:
            dt.write(system.atoms, atom_style=atom_style)

        if write_coeff:
            sorted_attrs = {
                "atoms":     self.get_sorted_attribute("atoms"),
                "bonds":     self.get_sorted_attribute("bonds"),
                "angles":    self.get_sorted_attribute("angles"),
                "dihedrals": self.get_sorted_attribute("dihedrals"),
                "impropers": self.get_sorted_attribute("impropers"),
            }
            temp_file = None
            try:
                with tempfile.NamedTemporaryFile(mode="w", prefix=".mdinterface-", suffix=".lammps",
                                                 dir=Path(filename).resolve().parent, delete=False) as tfile:
                    temp_file = Path(tfile.name)
                    with open(filename, "r") as ffile:
                        for fl in ffile:
                            if fl.startswith("Atoms"):
                                write_lammps_coefficients(system, sorted_attrs, fout=tfile)
                            tfile.write(fl)
                temp_file.replace(filename)
            finally:
                if temp_file is not None:
                    temp_file.unlink(missing_ok=True)

        if metadata is not None:
            from mdinterface.io.structure_metadata import lammps_metadata

            description = lammps_metadata(filename, system, atom_style, _get_mdi_version())
            Path(metadata).write_text(json.dumps(description, indent=2, allow_nan=False) + "\n")

        try:
            nbonds = len(system.atoms.bonds)
        except Exception:
            nbonds = 0
        logger.info("  >> Written: %d atoms,  %d bonds", len(system.atoms), nbonds)
        return

    def write_gromacs(
        self,
        prefix: str = "system",
        outdir: str = ".",
        system_name: str = "MD System",
    ) -> None:
        """
        Write GROMACS input files for the assembled system.

        Produces (all written to *outdir*):

        - ``{prefix}.gro`` - atomic coordinates and box vectors.
        - ``{resname}.itp`` - per-species include topology (one per unique
          species), with ``[ moleculetype ]``, ``[ atoms ]``,
          ``[ bonds ]``, ``[ pairs ]``, ``[ angles ]``, and ``[ dihedrals ]`` sections.
        - ``{prefix}.top`` - system topology that defines shared atom types, includes the ITP files, and
          lists molecule counts matching the GRO file.

        Parameters
        ----------
        prefix : str
            Base name for the ``.gro`` and ``.top`` files. Default ``"system"``.
        outdir : str
            Directory in which all files are written. Created if it does not
            exist. Default ``"."`` (current working directory).
        system_name : str
            Title written in the ``[ system ]`` section of the ``.top`` file.

        Raises
        ------
        RuntimeError
            If called before :meth:`build`.

        ValueError
            If parameters are incomplete or unsupported, species names collide
            with different definitions, or assembled residues do not match.

        .. warning::
            Experimental export for harmonic bonds/angles, four-term OPLS
            torsions, CVFF impropers, geometric LJ mixing and 0.5 LJ/Coulomb
            1-4 scaling. Constraints and simulation settings are not generated.
        """
        if self._universe is None:
            raise RuntimeError("Call build() before write_gromacs().")

        species = _prepare_species(self._all_species, self._universe)

        import os
        os.makedirs(outdir, exist_ok=True)

        log_header(logger, "Output")
        logger.warning(
            "write_gromacs is experimental -- verify output before production use."
        )

        itp_basenames = []
        for sp in species:
            itp_name = f"{sp.resname}.itp"
            itp_path = os.path.join(outdir, itp_name)
            write_gromacs_itp(sp, filename=itp_path, include_atomtypes=False)
            itp_basenames.append(itp_name)
            logger.info("  >> Species topology: %s", itp_path)

        # -- coordinates (.gro) --------------------------------------------
        gro_file = os.path.join(outdir, f"{prefix}.gro")
        self._universe.atoms.write(gro_file)
        logger.info("  >> Coordinates:      %s  (%d atoms)", gro_file, len(self._universe.atoms))

        # -- system topology (.top) ----------------------------------------
        top_file = os.path.join(outdir, f"{prefix}.top")
        write_gromacs_top(self._universe, itp_basenames,
                          filename=top_file, system_name=system_name, species=species)
        logger.info("  >> System topology:  %s", top_file)

    def to_ase(self) -> ase.Atoms:
        """
        Convert the assembled system to an ``ase.Atoms`` object.

        Returns
        -------
        ase.Atoms
            The simulation box as an ASE Atoms object with cell and PBC set.

        Raises
        ------
        RuntimeError
            If called before :meth:`build`.
        """
        if self._universe is None:
            raise RuntimeError("Call build() before to_ase().")

        system = self._universe
        positions = system.atoms.positions
        masses    = system.atoms.masses
        labels    = system.atoms.types

        symbols = [label_to_element(lab, mas) for lab, mas in zip(labels, masses)]
        atoms = ase.Atoms(symbols=symbols, positions=positions)

        if system.dimensions is not None:
            atoms.set_cell(system.dimensions)
            atoms.set_pbc(True)

        try:
            if system.atoms.charges is not None:
                atoms.set_initial_charges(system.atoms.charges)
        except Exception:
            pass

        return atoms

    def to_universe(self):
        return self.universe

    # ------------------------------------------------------------------
    # Properties
    # ------------------------------------------------------------------

    @property
    def universe(self) -> Optional[mda.Universe]:
        """The assembled MDAnalysis Universe (available after build())."""
        return self._universe

    # ------------------------------------------------------------------
    # Internal helpers
    # ------------------------------------------------------------------

    def _resolve_xy(self, match_cell):
        """
        Determine the starting XY cell dimensions and match-cell settings.

        Returns
        -------
        xsize, ysize : float
            Starting XY dimensions in Å.
        cell_ref : Specie or None
            Reference species that locks the XY cell, or ``None``.
        do_match : bool
            Whether slab conversion should apply XY scaling.
        """
        cell_ref, do_match = self._resolve_match_cell(match_cell)

        if cell_ref is not None:
            # Pre-fit through make_interface_slab to get *tiled* dimensions.
            # For a polymer (xrep=yrep=1) this is a no-op; for a crystalline
            # slab it tiles to xysize and returns the correct fitted cell.
            ref_fitted = make_interface_slab(cell_ref, self._xsize, self._ysize)
            if ref_fitted is not None:
                xsize = np.dot([1, 0, 0], ref_fitted.atoms.cell @ [1, 0, 0])
                ysize = np.dot([0, 1, 0], ref_fitted.atoms.cell @ [0, 1, 0])
            else:
                xsize = cell_ref.atoms.get_cell()[0][0]
                ysize = cell_ref.atoms.get_cell()[1][1]
            logger.info("  >> Reference cell: %s:  %.3f x %.3f Å",
                        getattr(cell_ref, "resname", "?"), xsize, ysize)
            logger.info("  >> All slabs will be stretched to match it")
        else:
            xsize, ysize = self._xsize, self._ysize
            if do_match:
                logger.info("  >> No reference cell: using largest XY slab")

        return xsize, ysize, cell_ref, do_match

    def _fit_slabs(self, xsize: float, ysize: float, cell_ref) -> Tuple[float, float]:
        """
        First pass: tile every slab layer and finalise the XY cell size.

        Tiled slabs are cached in ``layer["_slab"]``.  When no *cell_ref* is
        set, *xsize*/*ysize* grow to accommodate the largest slab footprint.

        Returns
        -------
        xsize, ysize : float
            Final XY dimensions in Å.
        """
        slab_layers = [l for l in self._layers if isinstance(l, SlabCompartment)]
        if slab_layers:
            log_subheader(logger, "Slab tiling")

        slab_dims = []
        for layer in self._layers:
            if not isinstance(layer, SlabCompartment):
                continue
            tslab = make_interface_slab(
                layer.species, xsize, ysize, layers=layer.nlayers
            )
            layer._slab = tslab
            if tslab is None:
                continue
            xi = np.dot([1, 0, 0], tslab.atoms.cell @ [1, 0, 0])
            yi = np.dot([0, 1, 0], tslab.atoms.cell @ [0, 1, 0])
            zi = np.dot([0, 0, 1], tslab.atoms.cell @ [0, 0, 1])
            slab_dims.append((xi, yi))
            layer._native_xy = (xi, yi)
            logger.info("  >> %s:  %.3f x %.3f x %.3f Å,  %d atoms",
                        getattr(layer.species, "resname", "?"),
                        xi, yi, zi, len(tslab.atoms))

        if slab_dims:
            if cell_ref is None:
                # Cell XY is defined by the largest fitted slab, not the
                # requested xysize (which was only used as a tiling target).
                xsize = max(d[0] for d in slab_dims)
                ysize = max(d[1] for d in slab_dims)
            logger.info("  >> Final XY size: %.3f x %.3f Å", xsize, ysize)
            self._check_xy_mismatch(slab_dims, xsize, ysize)

        return xsize, ysize

    def _stack_layers(
        self,
        xsize: float,
        ysize: float,
        all_sp_univs: list,
        padding: float,
        layered: bool,
        do_match: bool,
    ) -> Tuple[mda.Universe, float, float]:
        """
        Second pass: assemble all layers into a single Universe.

        Returns
        -------
        system : mda.Universe
        zdim : float
            Total Z height in Å (excluding final cell padding).
        first_layer_zdim : float
            Z height after the first layer only, used for centering.
        """
        system = None
        zdim = 0.0
        first_layer_zdim = 0.0
        n_layers = len(self._layers)

        for ii, layer in enumerate(self._layers):
            tag = f"[{ii + 1}/{n_layers}]"

            if isinstance(layer, SlabCompartment):
                name  = getattr(layer.species, "resname", "?")
                log_subheader(logger, f"Layer {tag}")
                logger.info("  >> %s: %s,  %d layer(s)", layer.label, name, layer.nlayers)
                slab_u = layer.build(xsize, ysize, all_sp_univs,
                                     layered=layered, do_match=do_match)
                zdim_before = zdim
                system, zdim = add_component(system, slab_u, zdim, padding=padding)
                sx, sy, sz = slab_u.dimensions[:3]
                native_xi, native_yi = layer._native_xy or (sx, sy)
                layer_zdim = zdim - zdim_before
                stretched = (do_match and not
                             (np.isclose(native_xi, sx, rtol=1e-2) and
                              np.isclose(native_yi, sy, rtol=1e-2)))
                extras = ([f"native XY: {native_xi:.3f} x {native_yi:.3f} Å"]
                          if stretched else [])
                _log_layer_result(len(slab_u.atoms), (sx, sy, sz),
                                  zdim, layer_zdim, extra_lines=extras)

            elif isinstance(layer, SolventCompartment):
                solv_names = " + ".join(
                    getattr(s, "resname", "?") for s in layer.solvent
                ) or "ions"
                log_subheader(logger, f"Layer {tag}")
                logger.info("  >> solvent: %s", solv_names)
                if layer.density is not None:
                    logger.info("  >> density: %.2f g/cm³", layer.density)
                solv_box = layer.build(xsize, ysize, all_sp_univs)
                if solv_box is not None:
                    zdim_before = zdim
                    system, zdim = add_component(system, solv_box, zdim, padding=padding)
                    layer_zdim = zdim - zdim_before
                    svx, svy, svz = solv_box.dimensions[:3]
                    mol_counts = Counter(res.resname for res in solv_box.residues)
                    # one line listing every species (solvent + solute) with counts
                    all_sp = list(layer.solvent) + list(layer.solute)
                    mol_line = ",  ".join(
                        f"{getattr(s, 'resname', '?')} ({mol_counts.get(getattr(s, 'resname', '?'), 0)} mol)"
                        for s in all_sp
                    )
                    if mol_line:
                        logger.info("  >> %s", mol_line)
                    _log_layer_result(len(solv_box.atoms), (svx, svy, svz),
                                      zdim, layer_zdim)
                else:
                    logger.warning("  >> empty; check packmol.log")
                    system, zdim = add_component(system, solv_box, zdim, padding=padding)

            elif isinstance(layer, VacuumCompartment):
                layer_zdim = layer.zdim
                zdim += layer_zdim
                log_subheader(logger, f"Layer {tag}")
                logger.info("  >> vacuum: %.1f Å", layer_zdim)
                _log_layer_result(None, None, zdim, layer_zdim)

            if ii == 0:
                first_layer_zdim = zdim

        return system, zdim, first_layer_zdim

    @staticmethod
    def _resolve_match_cell(match_cell):
        """
        Parse the *match_cell* parameter used in :meth:`build`.

        Returns
        -------
        cell_ref : Specie or None
            The reference species whose XY cell locks *xsize/ysize*, or
            ``None`` when a plain bool was supplied.
        do_match : bool
            Whether ``to_universe`` should apply XY scaling.
        """
        if match_cell is False or match_cell is True:
            return None, bool(match_cell)
        return match_cell, True

    @staticmethod
    def _validate_xysize(
        xysize: Union[List[float], Tuple[float, float]]
    ) -> Tuple[float, float]:
        if not isinstance(xysize, (list, tuple, np.ndarray)):
            raise TypeError(
                f"xysize must be a list, tuple, or numpy array, "
                f"got {type(xysize).__name__}."
            )
        if len(xysize) != 2:
            raise ValueError(
                f"xysize must have exactly 2 elements [x, y], "
                f"got {len(xysize)}: {xysize}."
            )
        try:
            xsize, ysize = float(xysize[0]), float(xysize[1])
        except (ValueError, TypeError) as e:
            raise ValueError(f"xysize elements must be numeric: {e}") from e
        if xsize <= 0 or ysize <= 0:
            raise ValueError(
                f"xysize values must be positive, got [{xsize}, {ysize}]."
            )
        return xsize, ysize

    def _check_xy_mismatch(
        self, slab_dims: list, xsize: float, ysize: float
    ) -> None:
        """Warn if slabs have inconsistent XY sizes or differ much from requested."""
        threshold = 0.1   # 10 % relative difference triggers a warning

        # Check consistency across slabs.
        if len(slab_dims) > 1:
            xs = [d[0] for d in slab_dims]
            ys = [d[1] for d in slab_dims]
            x_spread = (max(xs) - min(xs)) / xsize
            y_spread = (max(ys) - min(ys)) / ysize
            if x_spread > threshold or y_spread > threshold:
                logger.warning(
                    "Slab XY dimensions are inconsistent: "
                    "x range [%.2f – %.2f Å], y range [%.2f – %.2f Å]. "
                    "Consider using match_cell=True in build().",
                    min(xs), max(xs), min(ys), max(ys),
                )

        # Check against original requested xysize.
        req_x, req_y = self._xsize, self._ysize
        dx = abs(xsize - req_x) / req_x
        dy = abs(ysize - req_y) / req_y
        if dx > threshold or dy > threshold:
            logger.warning(
                "Actual cell size (%.2f × %.2f Å) differs from requested "
                "(%.2f × %.2f Å) by %.1f%% / %.1f%%. "
                "The slab periodicity determines the final XY size.",
                xsize, ysize, req_x, req_y, dx * 100, dy * 100,
            )

    @staticmethod
    def _apply_stack_axis(
        system: mda.Universe,
        xsize: float,
        ysize: float,
        zdim: float,
        stack_axis: str,
    ):
        """
        Permute atomic coordinates so that the stacking direction is *stack_axis*.

        Assembly always builds along Z.  This method applies a lossless
        coordinate permutation and updates the cell dimensions accordingly:

        - ``"x"`` : (x, y, z) → (z, y, x)  |  cell [zdim, ysize, xsize]
        - ``"y"`` : (x, y, z) → (x, z, y)  |  cell [xsize, zdim, ysize]

        Returns the updated (system, new_xsize, new_ysize).
        """
        pos = system.atoms.positions
        if stack_axis == "x":
            system.atoms.positions = pos[:, [2, 1, 0]]
            system.dimensions = [zdim, ysize, xsize, 90, 90, 90]
            new_xsize, new_ysize = zdim, ysize
        else:  # "y"
            system.atoms.positions = pos[:, [0, 2, 1]]
            system.dimensions = [xsize, zdim, ysize, 90, 90, 90]
            new_xsize, new_ysize = xsize, zdim
        logger.info("  >> axis permuted Z -> %s", stack_axis.upper())
        return system, new_xsize, new_ysize

    def _register(self, *species: Any) -> None:
        """Add species to the global registry if not already present."""
        for sp in species:
            if sp is not None and not any(s is sp for s in self._all_species):
                self._all_species.append(sp)

    @staticmethod
    def _copy_filled_region(fr: FilledRegion) -> FilledRegion:
        """Deep-copy a FilledRegion's species and nested regions (region geometry is immutable, shared as-is)."""
        if fr.solvent is None:
            solvent_copy = None
        elif isinstance(fr.solvent, (list, tuple)):
            solvent_copy = [s.copy() for s in fr.solvent]
        else:
            solvent_copy = fr.solvent.copy()
        solute_copy = [sp.copy() for sp in fr.solute] if fr.solute else None
        return FilledRegion(
            region=fr.region,
            solvent=solvent_copy,
            solute=solute_copy,
            nsolute=fr.nsolute,
            density=fr.density,
            nsolvent=fr.nsolvent,
            concentration=fr.concentration,
            conmodel=fr.conmodel,
            ratio=fr.ratio,
            regions=[SimCell._copy_filled_region(child) for child in fr.regions],
        )

    def _update_topology_indexes(self) -> None:
        n_species = len(self._all_species)
        nitems: dict = {
            "_btype": [],
            "_atype": [],
            "_dtype": [],
            "_itype": [],
        }
        for attribute in nitems:
            for specie in self._all_species:
                for attr in specie.__getattribute__(attribute):
                    if attr.id not in nitems[attribute]:
                        nitems[attribute].append(attr.id)
                    else:
                        idx = find_smallest_missing(nitems[attribute], start=1)
                        attr.set_id(idx)
                        nitems[attribute].append(attr.id)

        # Sort atom types alphabetically by extended label.
        atom_types = []
        for specie in self._all_species:
            atom_types.extend([stype.extended_label for stype in specie._stype])
        atom_types.sort()

        for specie in self._all_species:
            for stype in specie._stype:
                idx = np.argwhere(stype.extended_label == np.array(atom_types))[0][0]
                stype.set_id(idx + 1)

        logger.debug(
            "  >> %d species,  %d atom types,  %d bond,  %d angle,  %d dihedral,  %d improper",
            n_species, len(atom_types), len(nitems["_btype"]),
            len(nitems["_atype"]), len(nitems["_dtype"]), len(nitems["_itype"]),
        )

    def get_sorted_attribute(self, attribute: str) -> list:
        attr_map = {
            "bonds":     "_btype",
            "angles":    "_atype",
            "dihedrals": "_dtype",
            "impropers": "_itype",
            "atoms":     "_stype",
        }
        key = attr_map.get(attribute.lower(), attribute)

        indexes = []
        attributes = []
        for specie in self._all_species:
            for attr in specie.__getattribute__(key):
                indexes.append(attr.id)
                attributes.append(attr)

        return [attributes[ii] for ii in np.argsort(indexes)]

universe property

The assembled MDAnalysis Universe (available after build()).

__init__(xysize, verbose=None)

Parameters:

Name	Type	Description	Default
xysize	list or tuple of float	XY dimensions of the simulation box in Å, e.g. [15.0, 15.0].	required
verbose	None, bool, int, or str	Controls package-wide log verbosity via :func:set_verbosity. None (default) leaves the current level unchanged. Integer scale: 0 / False — WARNING (quiet) 1 / True — INFO (normal) 2 — DEBUG (detailed) 3, 4, … — DEBUG (same as 2 for values < 10) Integers ≥ 10 are treated as raw logging level constants. Strings are resolved by name ("DEBUG", "INFO", …).	None

Source code in mdinterface/build/builder.py

def __init__(
    self,
    xysize: Union[List[float], Tuple[float, float]],
    verbose: Union[None, bool, int, str] = None,
) -> None:
    """
    Parameters
    ----------
    xysize : list or tuple of float
        XY dimensions of the simulation box in Å, e.g. ``[15.0, 15.0]``.
    verbose : None, bool, int, or str, optional
        Controls package-wide log verbosity via :func:`set_verbosity`.
        ``None`` (default) leaves the current level unchanged.

        Integer scale:

        - ``0`` / ``False`` — WARNING (quiet)
        - ``1`` / ``True``  — INFO  (normal)
        - ``2``             — DEBUG (detailed)
        - ``3``, ``4``, …  — DEBUG (same as 2 for values < 10)

        Integers ≥ 10 are treated as raw ``logging`` level constants.
        Strings are resolved by name (``"DEBUG"``, ``"INFO"``, …).
    """
    if verbose is not None:
        set_verbosity(verbose)

    xsize, ysize = self._validate_xysize(xysize)
    self._xsize = xsize
    self._ysize = ysize
    self._layers: List[Compartment] = []
    self._all_species: List[Any] = []
    self._universe: Optional[mda.Universe] = None

    log_banner(logger, "mdinterface :: SimCell", f"version {_get_mdi_version()}")
    logger.info("  >> xysize: %.2f x %.2f Å", xsize, ysize)

add_slab(species, nlayers=1)

Add a solid-interface layer (slab).

Parameters:

Name	Type	Description	Default
species	Specie	The unit-cell species to tile into a slab.	required
nlayers	int	Number of unit-cell layers to stack along Z.	1

Source code in mdinterface/build/builder.py

def add_slab(self, species: Any, nlayers: int = 1) -> None:
    """
    Add a solid-interface layer (slab).

    Parameters
    ----------
    species : Specie
        The unit-cell species to tile into a slab.
    nlayers : int
        Number of unit-cell layers to stack along Z.
    """
    slab_sp = species.copy()
    slab_idx = sum(1 for lay in self._layers if isinstance(lay, SlabCompartment))
    suffix = f"_s{slab_idx}"
    for atom in slab_sp._stype:
        atom.set_label(atom.label + suffix)
    self._register(slab_sp)
    self._layers.append(SlabCompartment(species=slab_sp, nlayers=nlayers, label="slab"))
    logger.info("  + slab     %s,  %d layer(s)",
                getattr(slab_sp, "resname", "?"), nlayers)

add_prebuilt(species, nlayers=1)

Add a pre-built layer whose atomic positions come from a prior MD run.

Semantically equivalent to :meth:add_slab but intended for species whose positions have already been set (e.g. via :meth:~mdinterface.core.specie.Specie.update_positions). No XY tiling is performed; the species cell is used as-is.

Any pre-processing of the positions (centering, trimming, …) should be done on the species before passing it here.

Parameters:

Name	Type	Description	Default
species	Specie or Polymer	Species with positions already set from a prior trajectory.	required
nlayers	int	Number of repeat units to stack along Z (usually 1).	1

Source code in mdinterface/build/builder.py

def add_prebuilt(self, species: Any, nlayers: int = 1) -> None:
    """
    Add a pre-built layer whose atomic positions come from a prior MD run.

    Semantically equivalent to :meth:`add_slab` but intended for species
    whose positions have already been set (e.g. via
    :meth:`~mdinterface.core.specie.Specie.update_positions`).  No XY
    tiling is performed; the species cell is used as-is.

    Any pre-processing of the positions (centering, trimming, …) should
    be done on the species before passing it here.

    Parameters
    ----------
    species : Specie or Polymer
        Species with positions already set from a prior trajectory.
    nlayers : int
        Number of repeat units to stack along Z (usually 1).
    """
    slab_sp = species.copy()
    slab_idx = sum(1 for lay in self._layers if isinstance(lay, SlabCompartment))
    suffix = f"_s{slab_idx}"
    for atom in slab_sp._stype:
        atom.set_label(atom.label + suffix)
    self._register(slab_sp)
    self._layers.append(SlabCompartment(species=slab_sp, nlayers=nlayers, label="prebuilt"))
    logger.info("  + prebuilt %s,  %d layer(s)",
                getattr(slab_sp, "resname", "?"), nlayers)

add_solvent(solvent=None, solute=None, nsolute=None, zdim=None, density=None, nsolvent=None, concentration=None, conmodel=None, solute_pos=None, dilate=1.0, packmol_tolerance=2.0, ratio=None, regions=None, seed=None)

Add a solvent (liquid) layer, optionally with dissolved species.

Parameters:

Name	Type	Description	Default
solvent	Specie or list of Specie or None	Solvent molecule(s). Pass None for a solute-only region.	None
solute	list of Specie	Species to dissolve (ions, neutral molecules, …).	None
nsolute	int or list of int	Number of each solute species (alternative to concentration).	None
zdim	float	Thickness of this region in Angstroms.	None
density	float	Solvent density in g/cm³ (alternative to nsolvent).	None
nsolvent	int or list of int	Explicit number of solvent molecules (alternative to density).	None
concentration	float	Solute concentration in Molar (alternative to nsolute).	None
conmodel	dict	Spatially varying concentration model.	None
solute_pos	str	Solute placement strategy: None / "packmol" — PACKMOL random placement in the full box (default). "left" — PACKMOL random placement in the left half (z ≤ zdim/2). "right" — PACKMOL random placement in the right half (z ≥ zdim/2). "center" — each molecule fixed at the box centre.	None
dilate	float	Dilation factor > 1 for concentrated systems. PACKMOL will pack into a box dilate times taller at density / dilate, keeping the number of solvent molecules identical. The larger box gives PACKMOL more breathing room; subsequent NpT equilibration compresses the system to the correct density. Default is 1.0 (no dilation).	1.0
packmol_tolerance	float	Minimum distance (Å) between atoms of different molecules during PACKMOL packing. Reduce from the default of 2.0 for very concentrated systems if PACKMOL cannot find a valid packing.	2.0
regions	list of FilledRegion	Spatially-heterogeneous sub-regions within this layer (e.g. a gas pocket embedded in the bulk solvent). Build with SomeRegion(...).fill(solvent=..., density=...) — see :mod:mdinterface.build.regions. A region's center may be "random" to have a non-overlapping placement chosen automatically within its parent volume. Bulk solvent and solute exclude these regions; density and concentration use the remaining volume. Cannot be combined with conmodel or solute_pos="center".	None
seed	int	Seed for the RNG used to resolve center="random" regions. Pass the same seed to reproduce an identical placement.	None

Source code in mdinterface/build/builder.py

def add_solvent(
    self,
    solvent: Optional[Any] = None,
    solute: Optional[List[Any]] = None,
    nsolute: Optional[Union[int, List[int]]] = None,
    zdim: Optional[float] = None,
    density: Optional[float] = None,
    nsolvent: Optional[Union[int, List[int]]] = None,
    concentration: Optional[float] = None,
    conmodel: Optional[dict] = None,
    solute_pos: Optional[str] = None,
    dilate: float = 1.0,
    packmol_tolerance: float = 2.0,
    ratio: Optional[List[float]] = None,
    regions: Optional[List[FilledRegion]] = None,
    seed: Optional[int] = None,
) -> None:
    """
    Add a solvent (liquid) layer, optionally with dissolved species.

    Parameters
    ----------
    solvent : Specie or list of Specie or None
        Solvent molecule(s). Pass ``None`` for a solute-only region.
    solute : list of Specie, optional
        Species to dissolve (ions, neutral molecules, …).
    nsolute : int or list of int, optional
        Number of each solute species (alternative to *concentration*).
    zdim : float
        Thickness of this region in Angstroms.
    density : float, optional
        Solvent density in g/cm³ (alternative to *nsolvent*).
    nsolvent : int or list of int, optional
        Explicit number of solvent molecules (alternative to *density*).
    concentration : float, optional
        Solute concentration in Molar (alternative to *nsolute*).
    conmodel : dict, optional
        Spatially varying concentration model.
    solute_pos : str, optional
        Solute placement strategy:

        - ``None`` / ``"packmol"`` — PACKMOL random placement in the full
          box (default).
        - ``"left"``   — PACKMOL random placement in the left half
          (z ≤ zdim/2).
        - ``"right"``  — PACKMOL random placement in the right half
          (z ≥ zdim/2).
        - ``"center"`` — each molecule fixed at the box centre.
    dilate : float, optional
        Dilation factor > 1 for concentrated systems. PACKMOL will pack
        into a box ``dilate`` times taller at ``density / dilate``,
        keeping the number of solvent molecules identical. The larger box
        gives PACKMOL more breathing room; subsequent NpT equilibration
        compresses the system to the correct density. Default is ``1.0``
        (no dilation).
    packmol_tolerance : float, optional
        Minimum distance (Å) between atoms of different molecules during
        PACKMOL packing. Reduce from the default of ``2.0`` for very
        concentrated systems if PACKMOL cannot find a valid packing.
    regions : list of FilledRegion, optional
        Spatially-heterogeneous sub-regions within this layer (e.g. a gas
        pocket embedded in the bulk solvent). Build with
        ``SomeRegion(...).fill(solvent=..., density=...)`` — see
        :mod:`mdinterface.build.regions`. A region's ``center`` may be
        ``"random"`` to have a non-overlapping placement chosen
        automatically within its parent volume.
        Bulk solvent and solute exclude these regions; density and
        concentration use the remaining volume. Cannot be combined with
        ``conmodel`` or ``solute_pos="center"``.
    seed : int, optional
        Seed for the RNG used to resolve ``center="random"`` regions.
        Pass the same seed to reproduce an identical placement.
    """
    if zdim is None:
        raise ValueError("'zdim' is required for add_solvent()")
    if dilate <= 0:
        raise ValueError(f"'dilate' must be positive, got {dilate}")
    if packmol_tolerance <= 0:
        raise ValueError(f"'packmol_tolerance' must be positive, got {packmol_tolerance}")

    # Normalise solvent to a list; copy each species.
    if solvent is None:
        solv_copies = []
    elif isinstance(solvent, (list, tuple)):
        solv_copies = [s.copy() for s in solvent]
    else:
        solv_copies = [solvent.copy()]

    solute_copies = [sp.copy() for sp in (solute or [])]

    region_copies = [self._copy_filled_region(fr) for fr in (regions or [])]
    region_species: List[Any] = []
    for fr in region_copies:
        region_species.extend(_collect_region_species(fr))

    self._register(*solv_copies, *solute_copies, *region_species)

    self._layers.append(SolventCompartment(
        solvent=solv_copies,
        solute=solute_copies,
        nsolute=nsolute,
        zdim=zdim,
        density=density,
        nsolvent=nsolvent,
        concentration=concentration,
        conmodel=conmodel,
        solute_pos=solute_pos,
        dilate=dilate,
        packmol_tolerance=packmol_tolerance,
        ratio=ratio,
        regions=region_copies,
        seed=seed,
    ))
    solv_str = "+".join(getattr(s, "resname", "?") for s in solv_copies) or "ions"
    rho_str  = (f"ρ={density:.2f} g/cm³" if density is not None
                else f"nsolvent={nsolvent}" if nsolvent is not None else "density=?")
    solute_info = ""
    if solute_copies:
        sol_names = "+".join(getattr(s, "resname", "?") for s in solute_copies)
        solute_info = f",  solute: {sol_names} (n={nsolute})"
    region_info = f",  {len(region_copies)} region(s)" if region_copies else ""
    logger.info("  + solvent  %s,  zdim=%.1f Å,  %s%s%s",
                solv_str, zdim, rho_str, solute_info, region_info)

add_vacuum(zdim=0.0)

Add an empty vacuum gap.

Parameters:

Name	Type	Description	Default
zdim	float	Thickness of the vacuum gap in Angstroms.	0.0

Source code in mdinterface/build/builder.py

def add_vacuum(self, zdim: float = 0.0) -> None:
    """
    Add an empty vacuum gap.

    Parameters
    ----------
    zdim : float
        Thickness of the vacuum gap in Angstroms.
    """
    self._layers.append(VacuumCompartment(zdim=zdim))
    logger.info("  + vacuum   zdim=%.1f Å", zdim)

build(padding=0.5, center=False, layered=False, match_cell=True, hijack=None, stack_axis='z')

Assemble all layers into a simulation box.

Parameters:

Name	Type	Description	Default
padding	float	Spacing (Å) inserted between adjacent layers.	0.5
center	bool	If True, shift the system so the center of the first layer falls in the middle of the box (z=zdim/2), using its allocated thickness, then wrap individual atoms. The periodic seam falls in the opposite part of the system, which need not be vacuum. Before version 2.0.0, True centered the first layer on the periodic boundary instead. Default is False.	False
layered	bool	Assign distinct molecule indices to each slab layer for LAMMPS.	False
match_cell	bool or Specie	Controls XY cell matching: True — scale all slabs to the largest fitted XY cell (default). Ensures the solid/liquid interface is well-defined with no XY mismatch between layers. False — no scaling; each slab keeps its natural tiled XY. Specie — lock XY to that species' cell and stretch all other slabs to match. The reference species is left unscaled. Useful when a polymer or pre-relaxed slab should define the cell and electrodes should conform to it.	True
hijack	Atoms	After assembly, override both positions and cell dimensions with this external ase.Atoms object. Useful when you have a pre-relaxed structure you want to map the topology onto.	None
stack_axis	str	Axis along which layers are stacked in the output file. Assembly always happens along Z; a final coordinate permutation reorients the box. Accepted values: "x", "y", "z" (default "z").	'z'

Source code in mdinterface/build/builder.py

def build(
    self,
    padding: float = 0.5,
    center: bool = False,
    layered: bool = False,
    match_cell: Union[bool, Any] = True,
    hijack: Optional[ase.Atoms] = None,
    stack_axis: str = "z",
) -> None:
    """
    Assemble all layers into a simulation box.

    Parameters
    ----------
    padding : float
        Spacing (Å) inserted between adjacent layers.
    center : bool
        If True, shift the system so the center of the first layer falls
        in the middle of the box (z=zdim/2), using its allocated thickness,
        then wrap individual atoms. The periodic seam falls in the
        opposite part of the system, which need not be vacuum.
        Before version 2.0.0, True centered the first layer on the
        periodic boundary instead. Default is False.
    layered : bool
        Assign distinct molecule indices to each slab layer for LAMMPS.
    match_cell : bool or Specie
        Controls XY cell matching:

        - ``True`` — scale all slabs to the largest fitted XY cell
          (default). Ensures the solid/liquid interface is well-defined
          with no XY mismatch between layers.
        - ``False`` — no scaling; each slab keeps its natural tiled XY.
        - *Specie* — lock XY to that species' cell and stretch all other
          slabs to match. The reference species is left unscaled. Useful
          when a polymer or pre-relaxed slab should define the cell and
          electrodes should conform to it.
    hijack : ase.Atoms, optional
        After assembly, override both positions and cell dimensions with
        this external ``ase.Atoms`` object. Useful when you have a
        pre-relaxed structure you want to map the topology onto.
    stack_axis : str, optional
        Axis along which layers are stacked in the output file.
        Assembly always happens along Z; a final coordinate permutation
        reorients the box. Accepted values: ``"x"``, ``"y"``, ``"z"``
        (default ``"z"``).

    """
    stack_axis = stack_axis.lower()
    if stack_axis not in ("x", "y", "z"):
        raise ValueError(
            f"stack_axis must be 'x', 'y', or 'z', got {stack_axis!r}"
        )

    log_header(logger, "Build")
    logger.info("  >> %d layer(s)  |  inter-layer padding: %.2f Å",
                len(self._layers), padding)

    self._update_topology_indexes()

    xsize, ysize, cell_ref, do_match = self._resolve_xy(match_cell)
    xsize, ysize = self._fit_slabs(xsize, ysize, cell_ref)

    all_sp_univs = [sp.to_universe() for sp in self._all_species]
    system, zdim, first_layer_zdim = self._stack_layers(
        xsize, ysize, all_sp_univs, padding, layered, do_match
    )

    system.dimensions = [xsize, ysize, zdim] + [90, 90, 90]
    res_counts = Counter(res.resname for res in system.residues)

    log_header(logger, "Done")
    logger.info("  >> %d atoms  |  %.3f x %.3f x %.3f Å",
                len(system.atoms), xsize, ysize, zdim)
    logger.info("  >> %d layers  |  %d residues",
                len(self._layers), len(system.residues))
    parts = [f"{name} ({n} mol)" for name, n in res_counts.items()]
    for i in range(0, len(parts), 3):
        logger.info("  >> %s", ",  ".join(parts[i:i + 3]))

    if center:
        shift = (zdim - first_layer_zdim) / 2
        system.atoms.translate([0, 0, shift])
        _ = system.atoms.wrap()
        logger.info("  >> system centered on first layer (shift %.2f Å)", shift)

    if hijack is not None:
        system.dimensions = hijack.get_cell_lengths_and_angles()
        system.atoms.positions = hijack.get_positions()
        logger.info("  >> positions and cell overridden by hijack ase.Atoms")

    if stack_axis != "z":
        system, xsize, ysize = self._apply_stack_axis(system, xsize, ysize, zdim, stack_axis)

    self._universe = system
    self._xsize = xsize
    self._ysize = ysize

    return

write_lammps(filename='data.lammps', atom_style='full', write_coeff=True, metadata=None, expected_charge=None)

Write a LAMMPS data file (and optional force-field coefficients).

.. note:: This method is only needed for classical MD with LAMMPS. If you are using the assembled structure for AIMD, ML-MD, or any other workflow, use :meth:to_ase or :attr:universe instead — no force-field parameters are required for those paths.

Parameters:

Name	Type	Description	Default
filename	str	Output file path.	'data.lammps'
atom_style	str	LAMMPS atom style ("full" or "atomic").	'full'
write_coeff	bool	Whether to write force-field coefficient blocks.	True
metadata	str	Path for a schema-versioned JSON description of the final export. Includes exported IDs, species, topology, coefficients and checksum. Requires elements and, for full style, charges. Does not select MD settings.	None
expected_charge	float	Required total charge in elementary-charge units, checked with an absolute tolerance of 1e-5. Omit for systems with no specified target. Nonfinite charges and incomplete coefficients are always rejected when writing force-field coefficients.	None

Source code in mdinterface/build/builder.py

def write_lammps(
    self,
    filename: str = "data.lammps",
    atom_style: str = "full",
    write_coeff: bool = True,
    metadata: Optional[str] = None,
    expected_charge: Optional[float] = None,
) -> None:
    """
    Write a LAMMPS data file (and optional force-field coefficients).

    .. note::
        This method is only needed for classical MD with LAMMPS.  If you
        are using the assembled structure for AIMD, ML-MD, or any other
        workflow, use :meth:`to_ase` or :attr:`universe` instead — no
        force-field parameters are required for those paths.

    Parameters
    ----------
    filename : str
        Output file path.
    atom_style : str
        LAMMPS atom style (``"full"`` or ``"atomic"``).
    write_coeff : bool
        Whether to write force-field coefficient blocks.
    metadata : str, optional
        Path for a schema-versioned JSON description of the final export.
        Includes exported IDs, species, topology, coefficients and checksum.
        Requires elements and, for full style, charges. Does not select MD settings.
    expected_charge : float, optional
        Required total charge in elementary-charge units, checked with an
        absolute tolerance of 1e-5. Omit for systems with no specified target.
        Nonfinite charges and incomplete coefficients are always rejected
        when writing force-field coefficients.

    """
    if self._universe is None:
        raise RuntimeError("Call build() before write_lammps().")
    if metadata is not None:
        if Path(metadata).resolve() == Path(filename).resolve():
            raise ValueError("Metadata and LAMMPS data paths must differ.")
        if Path(metadata).exists():
            raise FileExistsError(metadata)
        if atom_style not in {"full", "atomic"}:
            raise ValueError("Metadata supports full or atomic exports.")
        if not hasattr(self._universe.atoms, "elements"):
            raise ValueError("Metadata requires explicit atom elements.")
        if atom_style == "full" and not hasattr(self._universe.atoms, "charges"):
            raise ValueError("Full-style metadata requires atom charges.")

    if hasattr(self._universe.atoms, "charges"):
        charges = self._universe.atoms.charges
        if not np.isfinite(charges).all():
            raise ValueError("Cannot export nonfinite partial charges.")
        if expected_charge is not None:
            if not np.isfinite(expected_charge) or not np.isclose(charges.sum(), expected_charge, atol=1e-5, rtol=0):
                raise ValueError(f"Total charge {charges.sum():.8f} does not match expected_charge={expected_charge}.")
    elif expected_charge is not None:
        raise ValueError("expected_charge requires atomic charge data.")
    if write_coeff:
        active_types = set(self._universe.atoms.types)
        for specie in self._all_species:
            if active_types.intersection(specie.get_atom_types()):
                specie.validate_force_field()

    log_header(logger, "Output")
    logger.info("  >> LAMMPS data file: %s  (style=%s,  coeff=%s)", filename, atom_style, write_coeff)
    system = self._universe.copy()

    if not write_coeff:
        for attr in ["bonds", "angles", "dihedrals", "impropers"]:
            try:
                system.del_TopologyAttr(attr)
            except Exception:
                pass

    with DATAWriter(filename) as dt:
        dt.write(system.atoms, atom_style=atom_style)

    if write_coeff:
        sorted_attrs = {
            "atoms":     self.get_sorted_attribute("atoms"),
            "bonds":     self.get_sorted_attribute("bonds"),
            "angles":    self.get_sorted_attribute("angles"),
            "dihedrals": self.get_sorted_attribute("dihedrals"),
            "impropers": self.get_sorted_attribute("impropers"),
        }
        temp_file = None
        try:
            with tempfile.NamedTemporaryFile(mode="w", prefix=".mdinterface-", suffix=".lammps",
                                             dir=Path(filename).resolve().parent, delete=False) as tfile:
                temp_file = Path(tfile.name)
                with open(filename, "r") as ffile:
                    for fl in ffile:
                        if fl.startswith("Atoms"):
                            write_lammps_coefficients(system, sorted_attrs, fout=tfile)
                        tfile.write(fl)
            temp_file.replace(filename)
        finally:
            if temp_file is not None:
                temp_file.unlink(missing_ok=True)

    if metadata is not None:
        from mdinterface.io.structure_metadata import lammps_metadata

        description = lammps_metadata(filename, system, atom_style, _get_mdi_version())
        Path(metadata).write_text(json.dumps(description, indent=2, allow_nan=False) + "\n")

    try:
        nbonds = len(system.atoms.bonds)
    except Exception:
        nbonds = 0
    logger.info("  >> Written: %d atoms,  %d bonds", len(system.atoms), nbonds)
    return

write_gromacs(prefix='system', outdir='.', system_name='MD System')

Write GROMACS input files for the assembled system.

Produces (all written to outdir):

• {prefix}.gro - atomic coordinates and box vectors.
• {resname}.itp - per-species include topology (one per unique species), with [ moleculetype ], [ atoms ], [ bonds ], [ pairs ], [ angles ], and [ dihedrals ] sections.
• {prefix}.top - system topology that defines shared atom types, includes the ITP files, and lists molecule counts matching the GRO file.

Parameters:

Name	Type	Description	Default
prefix	str	Base name for the .gro and .top files. Default "system".	'system'
outdir	str	Directory in which all files are written. Created if it does not exist. Default "." (current working directory).	'.'
system_name	str	Title written in the [ system ] section of the .top file.	'MD System'

Raises:

Type	Description
RuntimeError	If called before :meth:build.
ValueError	If parameters are incomplete or unsupported, species names collide with different definitions, or assembled residues do not match.
.. warning::	Experimental export for harmonic bonds/angles, four-term OPLS torsions, CVFF impropers, geometric LJ mixing and 0.5 LJ/Coulomb 1-4 scaling. Constraints and simulation settings are not generated.

Source code in mdinterface/build/builder.py

def write_gromacs(
    self,
    prefix: str = "system",
    outdir: str = ".",
    system_name: str = "MD System",
) -> None:
    """
    Write GROMACS input files for the assembled system.

    Produces (all written to *outdir*):

    - ``{prefix}.gro`` - atomic coordinates and box vectors.
    - ``{resname}.itp`` - per-species include topology (one per unique
      species), with ``[ moleculetype ]``, ``[ atoms ]``,
      ``[ bonds ]``, ``[ pairs ]``, ``[ angles ]``, and ``[ dihedrals ]`` sections.
    - ``{prefix}.top`` - system topology that defines shared atom types, includes the ITP files, and
      lists molecule counts matching the GRO file.

    Parameters
    ----------
    prefix : str
        Base name for the ``.gro`` and ``.top`` files. Default ``"system"``.
    outdir : str
        Directory in which all files are written. Created if it does not
        exist. Default ``"."`` (current working directory).
    system_name : str
        Title written in the ``[ system ]`` section of the ``.top`` file.

    Raises
    ------
    RuntimeError
        If called before :meth:`build`.

    ValueError
        If parameters are incomplete or unsupported, species names collide
        with different definitions, or assembled residues do not match.

    .. warning::
        Experimental export for harmonic bonds/angles, four-term OPLS
        torsions, CVFF impropers, geometric LJ mixing and 0.5 LJ/Coulomb
        1-4 scaling. Constraints and simulation settings are not generated.
    """
    if self._universe is None:
        raise RuntimeError("Call build() before write_gromacs().")

    species = _prepare_species(self._all_species, self._universe)

    import os
    os.makedirs(outdir, exist_ok=True)

    log_header(logger, "Output")
    logger.warning(
        "write_gromacs is experimental -- verify output before production use."
    )

    itp_basenames = []
    for sp in species:
        itp_name = f"{sp.resname}.itp"
        itp_path = os.path.join(outdir, itp_name)
        write_gromacs_itp(sp, filename=itp_path, include_atomtypes=False)
        itp_basenames.append(itp_name)
        logger.info("  >> Species topology: %s", itp_path)

    # -- coordinates (.gro) --------------------------------------------
    gro_file = os.path.join(outdir, f"{prefix}.gro")
    self._universe.atoms.write(gro_file)
    logger.info("  >> Coordinates:      %s  (%d atoms)", gro_file, len(self._universe.atoms))

    # -- system topology (.top) ----------------------------------------
    top_file = os.path.join(outdir, f"{prefix}.top")
    write_gromacs_top(self._universe, itp_basenames,
                      filename=top_file, system_name=system_name, species=species)
    logger.info("  >> System topology:  %s", top_file)

to_ase()

Convert the assembled system to an ase.Atoms object.

Returns:

Type	Description
Atoms	The simulation box as an ASE Atoms object with cell and PBC set.

Raises:

Type	Description
RuntimeError	If called before :meth:build.

Source code in mdinterface/build/builder.py

def to_ase(self) -> ase.Atoms:
    """
    Convert the assembled system to an ``ase.Atoms`` object.

    Returns
    -------
    ase.Atoms
        The simulation box as an ASE Atoms object with cell and PBC set.

    Raises
    ------
    RuntimeError
        If called before :meth:`build`.
    """
    if self._universe is None:
        raise RuntimeError("Call build() before to_ase().")

    system = self._universe
    positions = system.atoms.positions
    masses    = system.atoms.masses
    labels    = system.atoms.types

    symbols = [label_to_element(lab, mas) for lab, mas in zip(labels, masses)]
    atoms = ase.Atoms(symbols=symbols, positions=positions)

    if system.dimensions is not None:
        atoms.set_cell(system.dimensions)
        atoms.set_pbc(True)

    try:
        if system.atoms.charges is not None:
            atoms.set_initial_charges(system.atoms.charges)
    except Exception:
        pass

    return atoms