I/O

PACKMOL

mdinterface.build.box.PackmolError

Bases: RuntimeError

PACKMOL execution or output-processing failure.

Source code in mdinterface/build/box.py

class PackmolError(RuntimeError):
    """PACKMOL execution or output-processing failure."""

    def __init__(self, message, tempdir, log_path, returncode=None):
        self.tempdir = tempdir
        self.log_path = log_path
        self.returncode = returncode
        details = [message]
        if returncode is not None:
            details.append(f"return code: {returncode}")
        details.extend([f"temporary files: {tempdir}", f"log: {log_path}"])
        super().__init__("; ".join(details))

Reading structures

mdinterface.io.read.read_lammps_data_file(filename, pbc=False, ato_start_idx=0, is_snippet=False)

Source code in mdinterface/io/read.py

def read_lammps_data_file(filename, pbc=False, ato_start_idx=0, is_snippet=False):
    logger.debug("Reading LAMMPS data file: %s", filename)
    system = ase.io.lammpsdata.read_lammps_data(filename)

    if not pbc:
        system.set_pbc(False)
        system.set_cell(None)

    symbols = system.get_chemical_symbols()
    unique_symbols = set(symbols)

    pair_coeff     = []
    bond_coeff     = []
    angle_coeff    = []
    dihedral_coeff = []
    improper_coeff = []

    atoms     = []
    atosym    = []
    bonds     = []
    angles    = []
    dihedrals = []
    impropers = []

    section_map = {
        "masses"     :     "masses",
        'pair coeffs':     'pair_coeff',
        'bond coeffs':     'bond_coeff',
        'angle coeffs':    'angle_coeff',
        'dihedral coeffs': 'dihedral_coeff',
        'improper coeffs': 'improper_coeff',
        'atoms':            'atoms',
        'bonds':            'bonds',
        'angles':           'angles',
        'dihedrals':        'dihedrals',
        'impropers':        'impropers'
    }

    with open(filename, "r") as file:
        to_read = None

        for line in file:
            line = line.strip()

            if not line:
                continue

            section_key = line.lower()

            if section_key in section_map:
                to_read = section_map[section_key]
                continue

            line = line.split('#')[0].split()

            if to_read == "masses":
                mass = float(line[1])
                atosym.append(mass2symbol(mass, unique_symbols))

            elif to_read == "pair_coeff":
                eps = float(line[1])
                sig = float(line[2])
                pair_coeff.append([eps, sig])

            elif to_read == "bond_coeff":
                kr = float(line[1])
                r0 = float(line[2])
                bond_coeff.append([kr, r0])

            elif to_read == "angle_coeff":
                kr = float(line[1])
                theta0 = float(line[2])
                angle_coeff.append([kr, theta0])

            elif to_read == "dihedral_coeff":
                values = [float(ii) for ii in line[1:]]
                dihedral_coeff.append(values)

            elif to_read == "improper_coeff":
                values = [float(line[1]), int(line[2]), int(line[3])]
                improper_coeff.append(values)

            elif to_read == "atoms":
                atoidx = ato_start_idx + (int(line[0]) - 1)
                atotyp = int(line[2]) - 1
                eps, sig = pair_coeff[atotyp]
                symbol   = atosym[atotyp]

                if not is_snippet:
                    atom_name = f"{symbol}_{str(atoidx).zfill(3)}"
                else:
                    atom_name = f"{symbol}_{str(atoidx).zfill(3)}sn"

                atom = Atom(symbol, atom_name, eps, sig, atoid=atotyp)
                atoms.append(atom)

            elif to_read == "bonds":
                idx1 = int(line[2]) - 1
                idx2 = int(line[3]) - 1
                kr, r0 = bond_coeff[int(line[1]) - 1]
                bond = Bond(atoms[idx1].label, atoms[idx2].label, kr, r0)
                bonds.append(bond)

            elif to_read == "angles":
                idx1 = int(line[2]) - 1
                idx2 = int(line[3]) - 1
                idx3 = int(line[4]) - 1
                kr, theta0 = angle_coeff[int(line[1]) - 1]
                angle = Angle(atoms[idx1].label, atoms[idx2].label, atoms[idx3].label, kr, theta0)
                angles.append(angle)

            elif to_read == "dihedrals":
                idxs = [int(ii) - 1 for ii in line[2:]]
                values = dihedral_coeff[int(line[1]) - 1]
                dihedral = Dihedral(atoms[idxs[0]].label, atoms[idxs[1]].label,
                                    atoms[idxs[2]].label, atoms[idxs[3]].label, *values)
                dihedrals.append(dihedral)

            elif to_read == "impropers":
                idxs = [int(ii) - 1 for ii in line[2:]]
                values = improper_coeff[int(line[1]) - 1]
                improper = Improper(atoms[idxs[0]].label, atoms[idxs[1]].label,
                                    atoms[idxs[2]].label, atoms[idxs[3]].label,
                                    K=values[0], d=values[1], n=values[2])
                impropers.append(improper)

    system.new_array("stype", np.array(atoms))
    label_index = {atom.label: index for index, atom in enumerate(atoms)}
    system.set_array("topology_atom_id", np.arange(len(atoms), dtype=int))
    system.info["mdinterface_bonds"] = [(label_index[bond.symbols[0]], label_index[bond.symbols[1]]) for bond in bonds]
    logger.debug("  >> %d atoms, %d bonds, %d angles, %d dihedrals, %d impropers",
                 len(atoms), len(bonds), len(angles), len(dihedrals), len(impropers))
    return system, atoms, bonds, angles, dihedrals, impropers

mdinterface.io.read.read_lammps_nth_frame(filename, frame=-1)

Read a single frame from a LAMMPS dump file without loading the full trajectory.

Parameters:

Name	Type	Description	Default
filename	str	Path to the LAMMPS dump file.	required
frame	int	Frame index. 0 = first frame, -1 = last frame (default), -2 = second-to-last, etc. Reading the last frame (frame=-1) seeks from the end of the file and is O(frame size), making it fast even for multi-GB trajectories.	-1

Returns:

Type	Description
Atoms	The requested frame as an ASE Atoms object with cell and PBC set.

Raises:

Type	Description
IndexError	If the requested frame index is out of range.

Source code in mdinterface/io/read.py

def read_lammps_nth_frame(filename, frame=-1):
    """
    Read a single frame from a LAMMPS dump file without loading the full trajectory.

    Parameters
    ----------
    filename : str
        Path to the LAMMPS dump file.
    frame : int
        Frame index. ``0`` = first frame, ``-1`` = last frame (default),
        ``-2`` = second-to-last, etc.
        Reading the last frame (``frame=-1``) seeks from the end of the file
        and is O(frame size), making it fast even for multi-GB trajectories.

    Returns
    -------
    ase.Atoms
        The requested frame as an ASE Atoms object with cell and PBC set.

    Raises
    ------
    IndexError
        If the requested frame index is out of range.
    """
    logger.debug("Reading frame %d from dump: %s", frame, filename)
    if frame == -1:
        # Fast path: seek from EOF -- O(frame size) regardless of file size.
        return ase.io.lammpsrun.read_lammps_dump_text(
            StringIO(_read_last_frame_text(filename))
        )
    elif frame >= 0:
        for ii, chunk in enumerate(_iter_lammps_dump_frames(filename)):
            if ii == frame:
                return ase.io.lammpsrun.read_lammps_dump_text(StringIO(chunk))
        raise IndexError(
            f"Dump file '{filename}' has fewer than {frame + 1} frame(s)."
        )
    else:
        # Negative index: keep a rolling buffer of the last abs(frame) chunks.
        buf = collections.deque(maxlen=-frame)
        for chunk in _iter_lammps_dump_frames(filename):
            buf.append(chunk)
        if len(buf) < -frame:
            raise IndexError(
                f"Dump file '{filename}' has fewer than {-frame} frame(s)."
            )
        return ase.io.lammpsrun.read_lammps_dump_text(StringIO(buf[0]))

GROMACS writer

Warning

GROMACS output is experimental. Verify results against a reference before production use.

mdinterface.io.gromacswriter.write_gromacs_itp(specie, filename=None, *, include_atomtypes=True)

Write a GROMACS include topology (.itp) file for a Specie.

.. warning:: GROMACS output is experimental and assumes harmonic bonds and angles, four-term OPLS proper torsions, CVFF impropers, geometric LJ mixing, and 0.5 LJ/Coulomb 1-4 scaling. Constraints are not generated.

Parameters:

Name	Type	Description	Default
specie	Specie	The molecular species to write.	required
filename	str	Output filename. Defaults to {resname}.itp.	None
include_atomtypes	bool	Include atom-type definitions. For multi-species systems, set False and pass the species to :func:write_gromacs_top so all atom types precede all molecule definitions.	True

Raises:

Type	Description
ValueError	If parameters are incomplete, conflicting, or unsupported.

Source code in mdinterface/io/gromacswriter.py

def write_gromacs_itp(specie, filename=None, *, include_atomtypes=True):
    """
    Write a GROMACS include topology (.itp) file for a Specie.

    .. warning::
        GROMACS output is experimental and assumes harmonic bonds and angles,
        four-term OPLS proper torsions, CVFF impropers, geometric LJ mixing,
        and 0.5 LJ/Coulomb 1-4 scaling. Constraints are not generated.

    Parameters
    ----------
    specie : Specie
        The molecular species to write.
    filename : str, optional
        Output filename. Defaults to ``{resname}.itp``.
    include_atomtypes : bool, default True
        Include atom-type definitions. For multi-species systems, set False
        and pass the species to :func:`write_gromacs_top` so all atom types
        precede all molecule definitions.

    Raises
    ------
    ValueError
        If parameters are incomplete, conflicting, or unsupported.
    """
    _validate_specie(specie)
    records = _atomtype_records([specie])
    if filename is None:
        filename = f"{specie.resname}.itp"

    masses   = specie.atoms.get_masses()
    charges  = specie.atoms.get_initial_charges()
    atom_types = specie.get_atom_types()

    bond_idxs,  bond_tids  = specie.bonds
    angle_idxs, angle_tids = specie.angles
    dih_idxs,   dih_tids   = specie.dihedrals
    imp_idxs,   imp_tids   = specie.impropers

    # column width for name/type: fit the longest label, at least "name" (4)
    nw = max(max(len(t) for t in atom_types), 4)

    ATOMS_HDR = _atoms_hdr(nw)
    ATOMS_FMT = _atoms_fmt(nw)

    with open(filename, "w") as f:
        f.write("; GROMACS ITP file generated by mdinterface\n")
        f.write(f"; Molecule: {specie.resname}\n\n")

        if include_atomtypes:
            _write_atomtypes(f, records)

        # ------------------------------------------------------------------
        # [ moleculetype ]
        # ------------------------------------------------------------------
        f.write("[ moleculetype ]\n")
        f.write("; name      nrexcl\n")
        f.write(f"  {specie.resname:<10s}  3\n\n")

        # ------------------------------------------------------------------
        # [ atoms ]
        # ------------------------------------------------------------------
        f.write("[ atoms ]\n")
        f.write(ATOMS_HDR.format(
            "nr", "type", "resnr", "residu", "atom", "cgnr", "charge", "mass"))
        for ii, (atype, charge, mass, sid) in enumerate(
                zip(atom_types, charges, masses, specie._sids)):
            atom_name = sid.split("_")[0]
            f.write(ATOMS_FMT.format(
                ii+1, atype, 1, specie.resname, atom_name, ii+1, charge, mass))
        f.write("\n")

        # ------------------------------------------------------------------
        # [ bonds ]
        # ------------------------------------------------------------------
        if bond_idxs:
            f.write("[ bonds ]\n")
            f.write(_BONDS_HDR.format("ai", "aj", "funct", "b0(nm)", "kb(kJ/mol/nm2)"))
            for (ai, aj), tid in zip(bond_idxs, bond_tids):
                bond = specie._btype[tid]
                b0 = bond.r0 * _ANG_TO_NM
                # LAMMPS: E = K*(r-r0)^2  |  GROMACS: V = (kb/2)*(r-r0)^2  => kb = 2*K
                kb = 2 * bond.kr * _KCAL_TO_KJ / _ANG_TO_NM**2
                f.write(_BONDS_FMT.format(ai+1, aj+1, 1, b0, kb))
            f.write("\n")

        graph = nx.Graph()
        graph.add_edges_from(bond_idxs)
        pairs = sorted((a, b) for a, distances in nx.all_pairs_shortest_path_length(graph, cutoff=3)
                       for b, distance in distances.items() if a < b and distance == 3)
        if pairs:
            f.write("[ pairs ]\n; ai  aj  funct\n")
            for ai, aj in pairs:
                f.write(f"  {ai + 1}  {aj + 1}  1\n")
            f.write("\n")

        # ------------------------------------------------------------------
        # [ angles ]
        # ------------------------------------------------------------------
        if angle_idxs:
            f.write("[ angles ]\n")
            f.write(_ANGS_HDR.format("ai", "aj", "ak", "funct", "th0(deg)", "cth(kJ/mol/rad2)"))
            for (ai, aj, ak), tid in zip(angle_idxs, angle_tids):
                angle = specie._atype[tid]
                th0 = angle.theta0
                # LAMMPS: E = K*(θ-θ0)^2  |  GROMACS: V = (kt/2)*(θ-θ0)^2  => kt = 2*K
                cth = 2 * angle.kr * _KCAL_TO_KJ
                f.write(_ANGS_FMT.format(ai+1, aj+1, ak+1, 1, th0, cth))
            f.write("\n")

        # ------------------------------------------------------------------
        # [ dihedrals ] -- proper (OPLS Fourier -> funct=9)
        #
        # LAMMPS: E = (K1/2)(1+cos φ) + (K2/2)(1-cos 2φ)
        #           + (K3/2)(1+cos 3φ) + (K4/2)(1-cos 4φ)
        # GROMACS funct=9: V = kphi * (1 + cos(n*phi - phi0))
        #   K1 -> n=1, phi0=0,   kphi = K1/2
        #   K2 -> n=2, phi0=180, kphi = K2/2
        #   K3 -> n=3, phi0=0,   kphi = K3/2
        #   K4 -> n=4, phi0=180, kphi = K4/2
        # ------------------------------------------------------------------
        if dih_idxs:
            f.write("[ dihedrals ]\n")
            f.write("; proper dihedrals -- OPLS Fourier series, funct=9\n")
            f.write(_DIHS_HDR.format("ai", "aj", "ak", "al", "funct", "phi0(deg)", "kphi(kJ/mol)", "n"))
            opls_map = [(0, 1, 0.0), (1, 2, 180.0),
                        (2, 3, 0.0), (3, 4, 180.0)]
            for (ai, aj, ak, al), tid in zip(dih_idxs, dih_tids):
                dih  = specie._dtype[tid]
                vals = dih.values  # [K1, K2, K3, K4, K5?]
                for vi, n, phi0 in opls_map:
                    K = vals[vi]
                    if K == 0.0:
                        continue
                    kphi = K / 2.0 * _KCAL_TO_KJ
                    f.write(_DIHS_FMT.format(ai+1, aj+1, ak+1, al+1, 9, phi0, kphi, int(n)))
            f.write("\n")

        # ------------------------------------------------------------------
        # [ dihedrals ] -- improper (LAMMPS cvff -> funct=4)
        #
        # LAMMPS: E = K(1 + d*cos(n*phi))
        # GROMACS funct=4: V = kphi(1 + cos(n*psi - psi0))
        #   d= 1 -> psi0=0
        #   d=-1 -> psi0=180
        # ------------------------------------------------------------------
        if imp_idxs:
            f.write("[ dihedrals ]\n")
            f.write("; improper dihedrals, funct=4\n")
            f.write(_DIHS_HDR.format("ai", "aj", "ak", "al", "funct", "phi0(deg)", "kphi(kJ/mol)", "n"))
            for (ai, aj, ak, al), tid in zip(imp_idxs, imp_tids):
                imp      = specie._itype[tid]
                K, d, n  = imp.values
                phi0 = 0.0 if d == 1 else 180.0
                kphi = K * _KCAL_TO_KJ
                f.write(_DIHS_FMT.format(ai+1, aj+1, ak+1, al+1, 4, phi0, kphi, int(n)))
            f.write("\n")

mdinterface.io.gromacswriter.write_gromacs_top(universe, itp_files, filename='system.top', system_name='MD System', *, species=None)

Write a GROMACS system topology (.top) file.

Parameters:

Name	Type	Description	Default
universe	Universe	The assembled system (after SimCell.build()). Used to determine molecule counts from universe.residues.	required
itp_files	list of str	Basenames of the per-species ITP files to #include.	required
filename	str	Output filename. Default "system.top".	'system.top'
system_name	str	Title written in the [ system ] section.	'MD System'
species	sequence of Specie	Species whose atom types are written before the molecule includes. These ITP files must be written with include_atomtypes=False. Also validates residue identity against the assembled universe.	None

Source code in mdinterface/io/gromacswriter.py

def write_gromacs_top(universe, itp_files, filename="system.top",
                      system_name="MD System", *, species=None):
    """
    Write a GROMACS system topology (.top) file.

    Parameters
    ----------
    universe : mda.Universe
        The assembled system (after SimCell.build()).  Used to determine
        molecule counts from ``universe.residues``.
    itp_files : list of str
        Basenames of the per-species ITP files to ``#include``.
    filename : str, optional
        Output filename. Default ``"system.top"``.
    system_name : str, optional
        Title written in the ``[ system ]`` section.
    species : sequence of Specie, optional
        Species whose atom types are written before the molecule includes.
        These ITP files must be written with ``include_atomtypes=False``.
        Also validates residue identity against the assembled universe.
    """
    records = None
    if species is not None:
        records = _atomtype_records(_prepare_species(species, universe))
    runs = [(name, sum(1 for _ in group))
            for name, group in groupby(universe.residues.resnames)]

    with open(filename, "w") as f:
        f.write("; GROMACS topology file generated by mdinterface\n\n")

        # OPLS-AA defaults: LJ (nbfunc=1), geometric combining (comb-rule=3)
        f.write("[ defaults ]\n")
        f.write("; nbfunc  comb-rule  gen-pairs  fudgeLJ  fudgeQQ\n")
        f.write("  1       3          yes        0.5      0.5\n\n")

        if records:
            _write_atomtypes(f, records)

        for itp in itp_files:
            f.write(f'#include "{itp}"\n')
        f.write("\n")

        f.write("[ system ]\n")
        f.write(f"{system_name}\n\n")

        f.write("[ molecules ]\n")
        f.write("; {:<14}  {:>6}\n".format("molecule", "nmols"))
        for resname, count in runs:
            f.write(f"  {resname:<14s}  {count:>6d}\n")

Logging utilities

mdinterface.utils.logger.set_verbosity(level)

Set the log level for the entire mdinterface package.

A single StreamHandler is attached to the mdinterface root logger (at most once). All child loggers (mdinterface.build.builder, mdinterface.io.read, etc.) propagate to it automatically.

Parameters:

Name	Type	Description	Default
level	bool, int, or str	Small integers map to a simple verbosity scale: 0 / False -- WARNING (quiet) 1 / True -- INFO (normal) 2 -- DEBUG (detailed) Integers >= 10 are passed directly as Python logging level constants (e.g. logging.DEBUG, logging.WARNING). Strings are resolved by name ("DEBUG", "INFO", "WARNING" ...).	required

Source code in mdinterface/utils/logger.py

def set_verbosity(level) -> None:
    """
    Set the log level for the entire mdinterface package.

    A single StreamHandler is attached to the ``mdinterface`` root logger
    (at most once).  All child loggers (``mdinterface.build.builder``,
    ``mdinterface.io.read``, etc.) propagate to it automatically.

    Parameters
    ----------
    level : bool, int, or str
        Small integers map to a simple verbosity scale:

        - ``0`` / ``False`` -- WARNING (quiet)
        - ``1`` / ``True``  -- INFO (normal)
        - ``2``             -- DEBUG (detailed)

        Integers >= 10 are passed directly as Python logging level constants
        (e.g. ``logging.DEBUG``, ``logging.WARNING``).
        Strings are resolved by name (``"DEBUG"``, ``"INFO"``, ``"WARNING"`` ...).
    """
    if isinstance(level, bool):
        level = logging.INFO if level else logging.WARNING
    elif isinstance(level, int) and level < 10:
        level = _VERBOSITY_SCALE.get(level, logging.DEBUG)
    elif isinstance(level, str):
        level = getattr(logging, level.upper(), logging.INFO)

    parent = logging.getLogger(_PACKAGE)

    # Add a StreamHandler exactly once (NullHandler is not a StreamHandler).
    if not any(isinstance(h, logging.StreamHandler) for h in parent.handlers):
        handler = logging.StreamHandler()
        handler.setFormatter(_CompactFormatter(_FMT))
        parent.addHandler(handler)

    parent.setLevel(level)
    parent.propagate = False  # don't bubble up to the root logger

    # Make sure all existing child loggers propagate to the parent
    # and do not filter records themselves.
    for name, obj in logging.Logger.manager.loggerDict.items():
        if name.startswith(_PACKAGE + ".") and isinstance(obj, logging.Logger):
            obj.handlers.clear()
            obj.setLevel(logging.NOTSET)
            obj.propagate = True

mdinterface.utils.logger.log_header(log, title, level=logging.INFO)

Emit a fixed-width section header via log, preceded by a blank line.

The title is left-anchored after === and = fills the rest to :data:_HEADER_WIDTH characters, so headers align regardless of title length.

Example

::

from mdinterface.utils.logger import log_header
log_header(logger, "Build")
# [mdi] INFO |
# [mdi] INFO | ===  Build  ================================

Source code in mdinterface/utils/logger.py

def log_header(log: logging.Logger, title: str,
               level: int = logging.INFO) -> None:
    """
    Emit a fixed-width section header via *log*, preceded by a blank line.

    The title is left-anchored after ``===`` and ``=`` fills the rest to
    :data:`_HEADER_WIDTH` characters, so headers align regardless of title
    length.

    Example
    -------
    ::

        from mdinterface.utils.logger import log_header
        log_header(logger, "Build")
        # [mdi] INFO |
        # [mdi] INFO | ===  Build  ================================
    """
    log.log(level, "")
    left = f"=== {title} "
    fill = "=" * (_HEADER_WIDTH - len(left))
    log.log(level, "%s%s", left, fill)

mdinterface.utils.logger.log_subheader(log, title, level=logging.INFO)

Emit a fixed-width sub-section header via log.

Same width and left-anchor style as :func:log_header but uses - as the fill character to indicate a lower level of hierarchy.

Example

::

log_subheader(logger, "Layer [1/4]")
# [mdi] INFO | --  Layer [1/4]  ----------------------------

Source code in mdinterface/utils/logger.py

def log_subheader(log: logging.Logger, title: str,
                  level: int = logging.INFO) -> None:
    """
    Emit a fixed-width sub-section header via *log*.

    Same width and left-anchor style as :func:`log_header` but uses ``-``
    as the fill character to indicate a lower level of hierarchy.

    Example
    -------
    ::

        log_subheader(logger, "Layer [1/4]")
        # [mdi] INFO | --  Layer [1/4]  ----------------------------
    """
    log.log(level, "")
    left = f"-- {title} "
    fill = "-" * max(2, _HEADER_WIDTH - len(left))
    log.log(level, "%s%s", left, fill)

mdinterface.utils.logger.log_banner(log, *lines, level=logging.INFO)

Emit a prominent multi-line banner with = borders.

Each line is centred within :data:_HEADER_WIDTH characters. Intended for the top-level startup message of a major component (e.g. SimCell).

Example

::

log_banner(logger, "mdinterface :: SimCell", "version 1.5.0")
# [mdi] INFO | ============================================
# [mdi] INFO |        mdinterface :: SimCell
# [mdi] INFO |              version 1.5.0
# [mdi] INFO | ============================================

Source code in mdinterface/utils/logger.py

def log_banner(log: logging.Logger, *lines: str,
               level: int = logging.INFO) -> None:
    """
    Emit a prominent multi-line banner with ``=`` borders.

    Each line is centred within :data:`_HEADER_WIDTH` characters.  Intended
    for the top-level startup message of a major component (e.g. SimCell).

    Example
    -------
    ::

        log_banner(logger, "mdinterface :: SimCell", "version 1.5.0")
        # [mdi] INFO | ============================================
        # [mdi] INFO |        mdinterface :: SimCell
        # [mdi] INFO |              version 1.5.0
        # [mdi] INFO | ============================================
    """
    bar = "=" * _HEADER_WIDTH
    log.log(level, bar)
    for line in lines:
        log.log(level, line.center(_HEADER_WIDTH))
    log.log(level, bar)