Specie

mdinterface.core.specie.Specie

Bases: object

A single molecular species with geometry and force-field parameters.

Combines an ASE :class:ase.Atoms object with the topology (atom types, bonds, angles, dihedrals, impropers) and Lennard-Jones parameters needed to write a LAMMPS data file. All database entries (Water, Ion, Metal111, …) and Polymer inherit from this class.

.. note:: All topology parameters (bonds, angles, dihedrals, impropers, lj, atom_types, charges) are optional. A bare Specie(atoms) containing only positions is perfectly valid and sufficient for geometry tasks, AIMD, or ML-MD workflows. Force-field parameters are only needed when writing LAMMPS data files.

Parameters:

Name	Type	Description	Default
atoms	(Atoms, Mol or str)	Molecular geometry. A chemical formula string (e.g. "H2O") is accepted and converted to an ase.Atoms object.	None
charges	float, list, or array-like	Partial charges in elementary charge units. A scalar broadcasts to all atoms.	None
atom_types	list of Atom	Force-field atom-type objects. If None and lj is provided, types are inferred from element symbols.	None
bonds	Bond or list of Bond	Bond interaction definitions.	None
angles	Angle or list of Angle	Angle interaction definitions.	None
dihedrals	Dihedral or list of Dihedral	Dihedral interaction definitions.	None
impropers	Improper or list of Improper	Improper dihedral definitions.	None
lj	dict	Lennard-Jones parameters keyed by element symbol: {element: [epsilon (kcal/mol), sigma (Å)]}.	None
cutoff	float	Scale factor applied to covalent radii when building the bond graph. Increase slightly for flexible molecules with long bonds.	1.0
name	str	Residue name (truncated to 4 characters for PACKMOL compatibility). Defaults to the chemical formula.	None
lammps_data	str	Path to a LAMMPS data file. When provided, geometry and topology are read from this file and all other structural arguments are ignored.	None
fix_missing	bool	If True, attempt to auto-generate missing bond/angle/dihedral entries based on the bond graph.	False
chg_scaling	float	Multiplicative scale factor applied to charges.	1.0
pbc	bool	Whether to treat the atomic cell as periodic.	False
ligpargen	bool	If True, call LigParGen to generate OPLS-AA force-field parameters. Requires a working LigParGen installation and BOSSdir in config.	False
tot_charge	int	Total molecular charge (integer). Inferred from nominal_charge arrays if not supplied.	None
prune_z	bool	Remove the Z-coordinate component from atom positions.	False
calc	Calculator	ASE calculator to attach to the ase.Atoms object.	None
keep_ids	bool	Preserve atom IDs as read from the input rather than re-indexing.	False
smiles	str	Explicit SMILES input, including formal charges and stereochemistry. Cannot be combined with atoms or lammps_data. Hydrogens and a three-dimensional conformer are generated by RDKit. An RDKit molecule may also be passed as atoms.	None
seed	int	Random seed for RDKit conformer generation when coordinates are absent.	0

Examples:

Load from the built-in database::

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

Define a custom molecule::

from ase import Atoms
from mdinterface import Specie
mol = Atoms("CO2", positions=[[0,0,0],[1.16,0,0],[-1.16,0,0]])
co2 = Specie(mol, charges=[-0.3298, 0.6596, -0.3298])

Generate OPLS-AA parameters with LigParGen::

methanol = Specie("CH3OH", ligpargen=True)

Source code in mdinterface/core/specie.py

class Specie(object):
    """
    A single molecular species with geometry and force-field parameters.

    Combines an ASE :class:`ase.Atoms` object with the topology (atom types,
    bonds, angles, dihedrals, impropers) and Lennard-Jones parameters needed
    to write a LAMMPS data file.  All database entries (Water, Ion, Metal111,
    …) and Polymer inherit from this class.

    .. note::
        All topology parameters (``bonds``, ``angles``, ``dihedrals``,
        ``impropers``, ``lj``, ``atom_types``, ``charges``) are **optional**.
        A bare ``Specie(atoms)`` containing only positions is perfectly valid
        and sufficient for geometry tasks, AIMD, or ML-MD workflows.
        Force-field parameters are only needed when writing LAMMPS data files.

    Parameters
    ----------
    atoms : ase.Atoms, rdkit.Chem.Mol or str, optional
        Molecular geometry.  A chemical formula string (e.g. ``"H2O"``) is
        accepted and converted to an ``ase.Atoms`` object.
    charges : float, list, or array-like, optional
        Partial charges in elementary charge units.  A scalar broadcasts to
        all atoms.
    atom_types : list of Atom, optional
        Force-field atom-type objects.  If *None* and *lj* is provided, types
        are inferred from element symbols.
    bonds : Bond or list of Bond, optional
        Bond interaction definitions.
    angles : Angle or list of Angle, optional
        Angle interaction definitions.
    dihedrals : Dihedral or list of Dihedral, optional
        Dihedral interaction definitions.
    impropers : Improper or list of Improper, optional
        Improper dihedral definitions.
    lj : dict, optional
        Lennard-Jones parameters keyed by element symbol:
        ``{element: [epsilon (kcal/mol), sigma (Å)]}``.
    cutoff : float, default 1.0
        Scale factor applied to covalent radii when building the bond graph.
        Increase slightly for flexible molecules with long bonds.
    name : str, optional
        Residue name (truncated to 4 characters for PACKMOL compatibility).
        Defaults to the chemical formula.
    lammps_data : str, optional
        Path to a LAMMPS data file.  When provided, geometry and topology are
        read from this file and all other structural arguments are ignored.
    fix_missing : bool, default False
        If True, attempt to auto-generate missing bond/angle/dihedral entries
        based on the bond graph.
    chg_scaling : float, default 1.0
        Multiplicative scale factor applied to *charges*.
    pbc : bool, default False
        Whether to treat the atomic cell as periodic.
    ligpargen : bool, default False
        If True, call LigParGen to generate OPLS-AA force-field parameters.
        Requires a working LigParGen installation and ``BOSSdir`` in config.
    tot_charge : int, optional
        Total molecular charge (integer).  Inferred from ``nominal_charge``
        arrays if not supplied.
    prune_z : bool, default False
        Remove the Z-coordinate component from atom positions.
    calc : ase.Calculator, optional
        ASE calculator to attach to the ``ase.Atoms`` object.
    keep_ids : bool, default False
        Preserve atom IDs as read from the input rather than re-indexing.
    smiles : str, optional
        Explicit SMILES input, including formal charges and stereochemistry.
        Cannot be combined with ``atoms`` or ``lammps_data``. Hydrogens and
        a three-dimensional conformer are generated by RDKit. An RDKit
        molecule may also be passed as ``atoms``.
    seed : int, default 0
        Random seed for RDKit conformer generation when coordinates are absent.

    Examples
    --------
    Load from the built-in database::

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

    Define a custom molecule::

        from ase import Atoms
        from mdinterface import Specie
        mol = Atoms("CO2", positions=[[0,0,0],[1.16,0,0],[-1.16,0,0]])
        co2 = Specie(mol, charges=[-0.3298, 0.6596, -0.3298])

    Generate OPLS-AA parameters with LigParGen::

        methanol = Specie("CH3OH", ligpargen=True)
    """

    def __init__(self, atoms=None, charges=None, atom_types=None, bonds=None,
                 angles=None, dihedrals=None, impropers=None, lj=None, cutoff=1.0,
                 name=None, lammps_data=None, fix_missing=False, chg_scaling=1.0,
                 pbc=False, ligpargen=False, tot_charge=None, prune_z=False,
                 calc=None, keep_ids=False, smiles=None, seed=0):

        # store int. variables
        self.cutoff = cutoff
        if smiles is not None:
            if atoms is not None or lammps_data is not None:
                raise ValueError("smiles cannot be combined with atoms or lammps_data.")
            parser = Chem.SmilesParserParams()
            parser.removeHs = False
            mol = Chem.MolFromSmiles(smiles, parser)
            if mol is None:
                raise ValueError(f"Invalid SMILES: {smiles!r}")
            atoms = atoms_from_molecule(mol, seed=seed)
        elif isinstance(atoms, Chem.Mol):
            atoms = atoms_from_molecule(atoms, seed=seed)

        # if file provided, read it
        if lammps_data is not None:
            atoms, atom_types, bonds, angles, dihedrals, impropers = read_lammps_data_file(lammps_data)
        else:
            # read/update atoms atoms
            atoms, atom_types_tmp = self._read_atoms(atoms, charges,
                                                 chg_scaling=chg_scaling, pbc=pbc)
            if atom_types is None:
                atom_types = atom_types_tmp

        chemical_mol = stored_molecule(atoms)
        if chemical_mol is not None and tot_charge is not None:
            if Chem.GetFormalCharge(chemical_mol) != tot_charge:
                raise ValueError("tot_charge conflicts with the SMILES/RDKit formal charge.")

        # setup nominal charge
        if tot_charge is None:
            if not "nominal_charge" in atoms.arrays:
                atoms.set_array("nominal_charge", np.array(len(atoms)*[0]))
            tot_charge = int(np.sum(atoms.arrays["nominal_charge"]))

        # assign name
        if name is None:
            name = atoms.get_chemical_formula()
        if len(name) > 4:
            logger.warning(
                "Resname '%s' is longer than 4 characters and will be "
                "truncated to '%s'. PACKMOL requires resnames <= 4 chars.",
                name, name[:4],
            )
        self.resname = name[:4]

        # set up atoms and generate graph of specie
        self.set_atoms(atoms, cutoff, prune_z=prune_z)

        # read atom_types from LJ
        if atom_types is None:
            atom_types = self._atom_types_from_lj(lj if lj is not None else {})

        # set up internal topology attributes
        self._setup_topology(atom_types, bonds, angles, dihedrals, impropers,
                             keep_ids=keep_ids)

        if fix_missing:
            self._fix_missing_interactions()

        # initialize topology indexing
        self._update_topology_indexing()

        # run ligpargen to calculate parameters; for large molecules use the
        # segment-and-junction strategy automatically
        self._tot_charge = tot_charge
        if ligpargen:
            self.parameterize()
            self.atoms.set_initial_charges(self.charges * chg_scaling)

        # add calculator if provided
        self.assign_calculator(calc)

        return

    @property
    def atoms(self):
        return self._atoms

    @property
    def graph(self):
        return self._graph

    def to_rdkit(self):
        """Return an independent RDKit molecule with current coordinates.

        Returns
        -------
        rdkit.Chem.Mol
            Explicit-hydrogen molecular graph in ASE atom order. For inputs
            without stored chemistry, bonds are perceived from coordinates
            using the specified total charge.

        Raises
        ------
        ValueError
            If coordinate-based perception is requested for a periodic
            structure, or chemistry conflicts with formal-charge annotations.
        """
        return perceive_molecule(self.atoms, charge=self._tot_charge)

    def generate_conformer(self, seed=0, minimize=False, max_iterations=1000):
        """Generate molecular coordinates with RDKit ETKDGv3.

        Parameters
        ----------
        seed : int, default 0
            Nonnegative 32-bit embedding seed.
        minimize : bool, default False
            Also minimize the new coordinates with MMFF94.
        max_iterations : int, default 1000
            Maximum MMFF94 iterations when ``minimize=True``.

        Returns
        -------
        float or None
            MMFF94 energy in kcal/mol when minimized, otherwise None.

        Raises
        ------
        ValueError
            If options or chemistry are invalid, MMFF parameters are missing,
            or the structure is periodic, disconnected, or constrained.
        RuntimeError
            If embedding or minimization fails. Coordinates are unchanged.

        Notes
        -----
        Preserves atom ordering, topology and charges. A conformer is an
        isolated-molecule starting geometry, not an equilibrated ensemble.
        """
        from rdkit.Chem import rdDistGeom
        from mdinterface.core.chemistry import minimize_molecule

        if not isinstance(seed, (int, np.integer)) or not 0 <= seed < 2**31:
            raise ValueError("seed must be a nonnegative 32-bit integer.")
        if not isinstance(max_iterations, (int, np.integer)) or max_iterations < 1:
            raise ValueError("max_iterations must be a positive integer.")
        mol = self._geometry_molecule()
        params = rdDistGeom.ETKDGv3()
        params.randomSeed = int(seed)
        if rdDistGeom.EmbedMolecule(mol, params) < 0:
            raise RuntimeError("RDKit could not embed the molecule; try another seed.")
        energy = minimize_molecule(mol, max_iterations) if minimize else None
        positions = mol.GetConformer().GetPositions()
        if not np.isfinite(positions).all():
            raise RuntimeError("RDKit produced nonfinite coordinates.")
        self.atoms.set_positions(positions)
        return energy

    def minimize_geometry(self, max_iterations=1000):
        """Minimize existing molecular coordinates with RDKit MMFF94.

        Parameters
        ----------
        max_iterations : int, default 1000
            Maximum minimization iterations.

        Returns
        -------
        float
            MMFF94 energy in kcal/mol, not an OPLS-AA energy.

        Raises
        ------
        ValueError
            If arguments or chemistry are invalid, parameters are unavailable,
            or the structure is periodic, disconnected, or constrained.
        RuntimeError
            If minimization fails. Original coordinates remain unchanged.

        Notes
        -----
        Does not embed a new conformer or change simulation parameters.
        """
        from mdinterface.core.chemistry import minimize_molecule

        mol = self._geometry_molecule()
        energy = minimize_molecule(mol, max_iterations)
        self.atoms.set_positions(mol.GetConformer().GetPositions())
        return energy

    def _geometry_molecule(self):
        if any(self.atoms.pbc) or self.atoms.constraints:
            raise ValueError("Geometry preparation requires an unconstrained, nonperiodic molecule.")
        mol = self.to_rdkit()
        if len(Chem.GetMolFrags(mol)) != 1:
            raise ValueError("Geometry preparation requires one connected molecule.")
        return mol

    def to_smiles(self, mapped=False):
        """Return canonical isomeric SMILES for the molecular structure.

        Parameters
        ----------
        mapped : bool, default False
            Include persistent atom-map numbers.

        Returns
        -------
        str
            SMILES with formal charges and stereochemistry.
        """
        mol = self.to_rdkit()
        if not mapped:
            for atom in mol.GetAtoms():
                atom.SetAtomMapNum(0)
        return Chem.MolToSmiles(Chem.RemoveHs(mol), isomericSmiles=True)

    @property
    def bonds(self):
        if not self._bmap: return [[], []]
        return self._find_interactions(1, self._bmap)

    @property
    def angles(self):
        if not self._amap: return [[], []]
        return self._find_interactions(2, self._amap)

    @property
    def dihedrals(self):
        if not self._dmap: return [[], []]
        return self._find_interactions(3, self._dmap)

    @property
    def impropers(self):
        if not self._imap: return [[], []]
        return self._find_interactions(3, self._imap, impropers=True)

    @property
    def _sids(self):
        return self.atoms.arrays["sids"]

    @_sids.setter
    def _sids(self, atom_ids):
        self.atoms.arrays["sids"] = atom_ids

    def copy(self):
        return copy.deepcopy(self)

    @property
    def charges(self):
        return self.atoms.get_initial_charges()

    @property
    def calc(self):
        return self.atoms.calc

    # read atoms to return ase.Atoms
    @staticmethod
    def _read_atoms(atoms, charges, chg_scaling=1.0, pbc=False):

        # initialize atoms obj
        if isinstance(atoms, str):
            try:
                atoms = ase.io.read(atoms)
            except Exception:
                try:
                    atoms = ase.build.molecule(atoms)
                except Exception:
                    atoms = ase.Atoms(atoms)
        elif isinstance(atoms, ase.Atoms):
            atoms = atoms.copy()

        # assign charges
        if charges is not None:
            charges = as_list(charges)
            if len(charges) == 1:
                charges = len(atoms)*charges
        else:
            charges = atoms.get_initial_charges()

        # rescale charges if needed
        charges = chg_scaling*np.asarray(charges)
        atoms.set_initial_charges(charges)

        # see if stype is already present and contains Atom objects (not raw
        # strings from a serialized xyz -- those can't be used without LJ data)
        if "stype" in atoms.arrays:
            stype = atoms.arrays["stype"]
            if len(stype) == 0 or not hasattr(stype[0], "label"):
                stype = None
        else:
            stype = None

        return atoms, stype

    def _setup_topology(self, atoms, bonds, angles, dihedrals, impropers,
                        keep_ids=False):
        """
        Setup topology from scratch - this is run in the beginning
        """

        # define atoms
        atoms_list, atom_map, atom_ids = pmap.map_atoms(as_list(atoms), keep_ids=keep_ids)

        self._stype = atoms_list
        self._smap = atom_map
        self._sids = atom_ids

        # initialize topology values
        self._old_bonds = copy.deepcopy(as_list(bonds))
        self._old_angles = copy.deepcopy(as_list(angles))
        self._old_dihedrals = copy.deepcopy(as_list(dihedrals))
        self._old_impropers = copy.deepcopy(as_list(impropers))

        # map list of inputs
        self._update_topology_mappings()

        return

    def _update_topology_mappings(self):
        """
        Update topology mappings.
        """
        # Remap topology lists
        bonds_list, bond_map, bond_ids = pmap.map_bonds(self._old_bonds)
        angles_list, angle_map, angle_ids = pmap.map_angles(self._old_angles)
        dihedrals_list, dihedral_map, dihedral_ids = pmap.map_dihedrals(self._old_dihedrals)
        impropers_list, improper_map, improper_ids = pmap.map_impropers(self._old_impropers)

        self._btype = bonds_list
        self._atype = angles_list
        self._dtype = dihedrals_list
        self._itype = impropers_list

        self._bmap = bond_map
        self._amap = angle_map
        self._dmap = dihedral_map
        self._imap = improper_map

        self._bids = bond_ids
        self._aids = angle_ids
        self._dids = dihedral_ids
        self._iids = improper_ids

        # Update topology to ensure consistency
        self._update_topology_indexing()

        return

    def _update_topology_indexing(self):
        formula = self.atoms.get_chemical_formula()

        for attributes in ["_btype", "_atype", "_dtype", "_itype", "_stype"]:
            attr_type = []
            for attr in self.__getattribute__(attributes):
                attr.set_formula(formula)
                attr.set_resname(self.resname)
                if attr.id is None or attr.id in attr_type:
                    idx = find_smallest_missing(attr_type, start=1)
                    attr.set_id(idx)
                else:
                    idx = attr.id
                attr_type.append(idx)

        return

    def _update_junction_lj_types(self, cut_idxs, snippet_idxs, sn_atypes):
        """Update LJ atom types for junction atoms using a snippet LigParGen result.

        Called after a junction snippet run to correct the OPLS types of atoms
        that had H caps substituting their real bonded partner during the
        segment/monomer LigParGen run.  Only the atoms in *cut_idxs* are
        updated; atoms further from the cut already have correct context.

        Parameters
        ----------
        cut_idxs : iterable of int
            Original atom indices whose types need correcting (typically the
            two atoms on either side of the cut bond).
        snippet_idxs : array-like of int
            Mapping from snippet positions to original atom indices.
        sn_atypes : list of Atom
            Atom-type objects returned by ``run_ligpargen`` for the snippet.
        """
        for cut_idx in cut_idxs:
            sn_pos     = int(np.argwhere(snippet_idxs == cut_idx)[0][0])
            new_atom   = sn_atypes[sn_pos]
            atom_label = str(self._sids[cut_idx])
            existing_idx = next(
                (idx for idx, a in enumerate(self._stype) if a == new_atom),
                None)
            if existing_idx is None:
                added = copy.deepcopy(new_atom)
                added.set_label(atom_label)
                self._stype.append(added)
                existing_idx = len(self._stype) - 1
            self._smap[atom_label] = existing_idx

    def _add_to_topology(self, bonds=None, angles=None, dihedrals=None, impropers=None):
        if bonds is None:
            bonds = []
        if angles is None:
            angles = []
        if dihedrals is None:
            dihedrals = []
        if impropers is None:
            impropers = []

        # check for uniqueness
        new_bonds = []
        new_angles = []
        new_dihedrals = []
        new_impropers = []
        for bond in copy.deepcopy(bonds):
            if not any(bond.__eq_strict__(bb) for bb in self._old_bonds):
                new_bonds.append(bond)

        for angle in copy.deepcopy(angles):
            if not any(angle.__eq_strict__(aa) for aa in self._old_angles):
                new_angles.append(angle)

        for dihedral in copy.deepcopy(dihedrals):
            if not any(dihedral.__eq_strict__(dd) for dd in self._old_dihedrals):
                new_dihedrals.append(dihedral)

        for improper in copy.deepcopy(impropers):
            if not any(improper.__eq_strict__(ii) for ii in self._old_impropers):
                new_impropers.append(improper)

        # add to topology
        self._old_bonds     += new_bonds
        self._old_angles    += new_angles
        self._old_dihedrals += new_dihedrals
        self._old_impropers += new_impropers

        # update topology mapping
        self._update_topology_mappings()

        return

    def _cleanup_topology(self):
        """
        Remove topology interactions that reference deleted atoms.

        This method should be called after atoms are deleted from the system
        to ensure topology consistency.

        Parameters:
        deleted_atom_ids (list): List of atom IDs (from self._sids) that were deleted
        """

        atom_ids = self._sids

        # Clean bonds
        valid_bonds = []
        for bond in self._old_bonds:
            if all(sim in atom_ids for sim in bond.symbols):
                valid_bonds.append(bond)
        self._old_bonds = valid_bonds

        # Clean angles
        valid_angles = []
        for angle in self._old_angles:
            if all(sim in atom_ids for sim in angle.symbols):
                valid_angles.append(angle)
        self._old_angles = valid_angles

        # Clean dihedrals
        valid_dihedrals = []
        for dihedral in self._old_dihedrals:
            if all(sim in atom_ids for sim in dihedral.symbols):
                valid_dihedrals.append(dihedral)
        self._old_dihedrals = valid_dihedrals

        # Clean impropers
        valid_impropers = []
        for improper in self._old_impropers:
            if all(sim in atom_ids for sim in improper.symbols):
                valid_impropers.append(improper)
        self._old_impropers = valid_impropers

        # Update topology mappings
        self._update_topology_mappings()

        return

    def _fix_missing_interactions(self):

        mss_bnd = pmap.generate_missing_interactions(self, "bonds")
        mss_ang = pmap.generate_missing_interactions(self, "angles")
        mss_dih = pmap.generate_missing_interactions(self, "dihedrals")
        mss_imp = pmap.generate_missing_interactions(self, "impropers")

        self._add_to_topology(mss_bnd, mss_ang, mss_dih, mss_imp)

        return

    # function to setup atom types
    def _atom_types_from_lj(self, lj):

        # use function to retrieve IDs
        atom_type_ids, types_map = find_atom_types(self.graph, max_depth=1)

        atom_types = []
        for atom_id in atom_type_ids:

            atom_symbol = types_map[atom_id][0]
            atom_neighs = "".join(types_map[atom_id][1])

            label = "{}_{}".format(atom_symbol, atom_neighs)
            # types_map[atom_type] = label

            if label in lj:
                eps, sig = lj[label]
            elif atom_symbol in lj:
                eps, sig = lj[atom_symbol]
            else:
                eps, sig = None,None

            atom = Atom(atom_symbol, label=label, eps=eps, sig=sig)
            atom_types.append(atom)

        return atom_types

    def round_charges(self, nround=7):

        rcharges = round_list_to_sum(self.charges, round(sum(self.charges), nround), nround)
        self.atoms.set_initial_charges(rcharges)

        return

    def set_atoms(self, atoms, cutoff=1.0, prune_z=False):

        # make sure no mess
        atoms = atoms.copy()

        if prune_z:
            zmin = atoms.get_positions(wrap=True)[:,2].min()
            zmax = atoms.get_positions(wrap=True)[:,2].max()
            atoms.translate([0,0,-zmin])
            atoms.cell[2][2] = zmax-zmin
            atoms.pbc[2] = False

        if not any(atoms.pbc):
            atoms.set_center_of_mass([0,0,0])

        self._atoms = atoms
        mol = stored_molecule(atoms)
        if mol is not None:
            self._graph = graph_from_molecule(mol)
        elif "mdinterface_bonds" in atoms.info:
            ids = atoms.arrays["topology_atom_id"]
            if len(set(ids)) != len(ids):
                raise ValueError("Imported topology atom identifiers must be unique.")
            index = {int(number): position for position, number in enumerate(ids)}
            self._graph = nx.Graph()
            self._graph.add_nodes_from((i, {"element": symbol}) for i, symbol in enumerate(atoms.get_chemical_symbols()))
            self._graph.add_edges_from((index[a], index[b]) for a, b in atoms.info["mdinterface_bonds"] if a in index and b in index)
        else:
            self._graph = molecule_to_graph(atoms, cutoff_scale=cutoff)

        return

    def assign_calculator(self, calc=None):
        self.atoms.calc = calc
        return

    # covnert to mda.Universe
    def to_universe(self, charges=True, layered=False, match_cell=False, xydim=None):

        # empty top object
        top = mda.core.topology.Topology(n_atoms=len(self.atoms))

        # empty universe
        uni = mda.Universe(top, self.atoms.get_positions(), to_guess=())

        # add some stuff
        uni.add_TopologyAttr("masses", self.atoms.get_masses())
        uni.add_TopologyAttr("resnames", [self.resname])
        # add atom symbols
        atom_symbols = self.atoms.get_chemical_symbols()
        uni.add_TopologyAttr("names", atom_symbols)
        uni.add_TopologyAttr("elements", atom_symbols)

        # generate type
        types, indexes = self.get_atom_types(return_index=True)
        uni.add_TopologyAttr("types", types)
        uni.add_TopologyAttr("type_index", indexes)

        # populate with bonds and angles
        for att in ["bonds", "angles", "dihedrals", "impropers"]:
            attribute, types = self.__getattribute__(att)
            att_types = self._type2id(att, types)
            uni._add_topology_objects(att, attribute, types=att_types)

        # add charges
        if charges:
            uni.add_TopologyAttr("charges", self.atoms.get_initial_charges())

        # layer it nicely
        if layered:
            layer_idxs, __ = ase.geometry.get_layers(self.atoms, [0,0,1],
                                                      tolerance=0.01)
            groups = []
            for idx in np.unique(layer_idxs):
                idxs = np.where(layer_idxs == idx)[0]

                groups.append(uni.atoms[idxs])

            uni = mda.Merge(*groups)

            # uni.residues[0].resids = layer_idxs

        # match the cell!
        if match_cell:
            assert xydim is not None
            zz = self.atoms.get_cell()[2][2]
            tatoms = self.atoms.copy()

            tatoms.set_cell([xydim[0], xydim[1], zz], scale_atoms=True)

            uni.dimensions = tatoms.cell.cellpar()
            uni.atoms.positions = tatoms.get_positions()

        # if has cell info, pass them along
        elif self.atoms.get_cell():
            uni.dimensions = self.atoms.cell.cellpar()

        return uni

    def repeat(self, rep, make_cubic=False):

        atoms = self.atoms.repeat(rep)
        multiplier = int(np.prod(rep))
        mol = stored_molecule(self.atoms)
        if mol is not None:
            combined = Chem.Mol(mol)
            map_stride = max(a.GetAtomMapNum() for a in mol.GetAtoms())
            for repetition in range(1, multiplier):
                extra = Chem.Mol(mol)
                for atom in extra.GetAtoms():
                    atom.SetAtomMapNum(atom.GetAtomMapNum() + repetition * map_stride)
                combined = Chem.CombineMols(combined, extra)
            store_molecule(atoms, combined)

        if "mdinterface_bonds" in atoms.info:
            count = len(self.atoms)
            edges = list(self.graph.edges)
            atoms.info["mdinterface_bonds"] = [(a + offset, b + offset) for offset in range(0, len(atoms), count) for a, b in edges]
            atoms.set_array("topology_atom_id", np.arange(len(atoms), dtype=int))

        if make_cubic:  # Note: assumes orthogonal cell; off-diagonal components are discarded

            xsize = [1,0,0]@atoms.cell@[1,0,0]
            ysize = [0,1,0]@atoms.cell@[0,1,0]
            zsize = [0,0,1]@atoms.cell@[0,0,1]

            atoms.set_cell([xsize, ysize, zsize, 90, 90, 90])
            atoms.wrap()

        self.set_atoms(atoms)
        self._tot_charge *= multiplier

        return

    def _find_interactions(self, path_length, tag_map, impropers=False):

        if not impropers:
            paths = find_unique_paths_of_length(self.graph, path_length)
        else:
            paths = find_improper_idxs(self.graph)

        interaction_list = []
        interaction_type = []

        for indices in paths:
            atoms = [self._sids[idx] for idx in indices]
            symbols = [atom.split("_")[0] for atom in atoms]

            if tuple(atoms) in tag_map:
                interaction_list.append(indices)
                interaction_type.append(tag_map[tuple(atoms)])
            elif tuple(atoms[::-1]) in tag_map:
                interaction_list.append(indices[::-1])
                interaction_type.append(tag_map[tuple(atoms[::-1])])
            elif tuple(symbols) in tag_map:
                interaction_list.append(indices)
                interaction_type.append(tag_map[tuple(symbols)])
            elif tuple(symbols[::-1]) in tag_map:
                interaction_list.append(indices[::-1])
                interaction_type.append(tag_map[tuple(symbols[::-1])])

        interaction_list = np.array(interaction_list, dtype=int)
        interaction_type = np.array(interaction_type, dtype=int)

        return [interaction_list.tolist(), interaction_type.tolist()]

    def get_atom_types(self, return_index=False):

        type_indexes = np.array([self._smap[ii] for ii in self._sids])
        atom_types   = np.array([self._stype[ii].extended_label for ii in type_indexes])

        if return_index:
            return atom_types, type_indexes

        return atom_types

    # method to estimate the volume of the specie
    def estimate_specie_volume(self, probe_radius=0):

        try:
            from libarvo import molecular_vs
        except ImportError:
            raise ImportError(
                "libarvo is required for volume estimation but is not installed. "
                "Install it with: pip install libarvo"
            )

        centers = self.atoms.get_positions()
        radii = [vdw_radii[ii] for ii in self.atoms.get_atomic_numbers()]

        volume, surface = molecular_vs(centers, radii, probe_radius)

        return volume

    # method to estimate the radius of the specie if it were a sphere
    def estimate_specie_radius(self, probe_radius=0):
        volume = self.estimate_specie_volume(probe_radius=probe_radius)
        return (3 * volume / (4 * np.pi)) ** (1/3)

    def view(self):
        from ase.visualize import view

        view(self.atoms)
        return

    def _type2id(self, attribute, types):

        # Get the corresponding attribute list
        if attribute == "bonds":
            attribute_list = self._btype
        elif attribute == "angles":
            attribute_list = self._atype
        elif attribute == "dihedrals":
            attribute_list = self._dtype
        elif attribute == "impropers":
            attribute_list = self._itype
        elif attribute == "atoms":
            attribute_list = self._stype
        else:
            raise ValueError("Invalid topology attribute type.")

        return [attribute_list[idx].id for idx in types]

    def suggest_missing_interactions(self, stype="all"):

        # Check for missing interactions
        missing_bonds = pmap.find_missing_bonds(self)
        missing_angles = pmap.find_missing_angles(self)
        missing_dihedrals = pmap.find_missing_dihedrals(self)
        missing_impropers = pmap.find_missing_impropers(self)

        suggestions = {
            "bonds": missing_bonds,
            "angles": missing_angles,
            "dihedrals": missing_dihedrals,
            "impropers": missing_impropers
        }

        if stype == "all":
            return suggestions
        return suggestions[stype]

    def __repr__(self):
        return f"{self.__class__.__name__}({self.resname})"

    def find_relevant_distances(self, Nmax, Nmin=0, centers=None, Ninv=0):

        unique_pairs_list = find_relevant_distances(self.graph, Nmax, Nmin=Nmin,
                                                    centers=centers, Ninv=Ninv)
        return unique_pairs_list

    def plot_graph(self, show_bonds=False):
        import matplotlib.pyplot as plt

        from mdinterface.utils.draw import draw_bond_markers


        colors = [jmol_colors[a.number] for a in self.atoms]

        # node_pos = nx.spring_layout(self.graph, k=1.5/np.sqrt(self.graph.order()))
        node_pos = nx.kamada_kawai_layout(self.graph)

        fig, ax = plt.subplots()
        nx.draw(self.graph, pos=node_pos, with_labels=True, node_color=colors,
                node_size=1000, edge_color='black', linewidths=2, font_size=15,
                edgecolors="black", ax=ax, width=2)

        if show_bonds:
            draw_bond_markers(ax, self, node_pos, jmol_colors)

        plt.show()

        return

    def update_positions(self, positions=None, cellpar=None, atoms=None, prune_z=False):

        if atoms is not None:
            positions = atoms.get_positions()
            cellpar   = atoms.get_cell()

        atoms = self._atoms.copy()

        if positions is not None:
            atoms.set_positions(positions)
        if cellpar is not None:
            atoms.set_cell(cellpar)

        self.set_atoms(atoms, cutoff=self.cutoff, prune_z=prune_z)

        return

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

        See :func:`~mdinterface.io.gromacswriter.write_gromacs_itp` for full
        documentation.

        Parameters
        ----------
        filename : str, optional
            Output filename. Defaults to ``{resname}.itp``.
        include_atomtypes : bool, default True
            Include atom types in this file. For multi-species exports, use
            False and supply ``species`` to ``write_gromacs_top`` to define
            all atom types before all molecule definitions.
        """
        logger.warning(
            "write_gromacs_itp is experimental -- verify output before production use."
        )
        from mdinterface.io.gromacswriter import write_gromacs_itp
        write_gromacs_itp(self, filename=filename, include_atomtypes=include_atomtypes)

    def validate_force_field(self):
        """Raise if molecular parameters required for classical export are missing.

        Returns
        -------
        None
            Successful validation returns None.

        Raises
        ------
        ValueError
            If charges or coefficients are nonfinite, pair parameters are
            missing, or molecular bonds, angles or proper torsions are unassigned.

        Notes
        -----
        Completeness uses explicit chemistry or an existing bonded topology.
        Improper coefficients are checked when assigned; missing impropers
        cannot be inferred universally from connectivity. Constraints and
        force-field compatibility remain the caller's responsibility.
        """
        errors = []
        if not np.isfinite(self.charges).all():
            errors.append("nonfinite partial charges")
        for index in set(self._smap[sid] for sid in self._sids):
            atom = self._stype[index]
            if any(value is None or not np.isfinite(value) or value < 0 for value in (atom.eps, atom.sig)):
                errors.append(f"missing or invalid pair parameters for {atom.label}")
        molecular = stored_molecule(self.atoms) is not None or any((self._old_bonds, self._old_angles, self._old_dihedrals))
        for kind, type_key, finder in (("bonds", "_btype", pmap.find_missing_bonds),
                                       ("angles", "_atype", pmap.find_missing_angles),
                                       ("dihedrals", "_dtype", pmap.find_missing_dihedrals),
                                       ("impropers", "_itype", None)):
            if molecular and finder is not None:
                missing = finder(self)
                if missing:
                    errors.append(f"{len(missing)} unassigned {kind} (first: {missing[0]})")
            _, indices = getattr(self, kind)
            for index in set(indices):
                parameter = getattr(self, type_key)[index]
                values = parameter.values[:4] if kind == "dihedrals" and parameter.values[-1] is None else parameter.values
                if any(value is None or not np.isfinite(value) for value in values):
                    errors.append(f"missing or nonfinite {kind} coefficients for {parameter.symbols}")
        if errors:
            raise ValueError(f"Incomplete force field for {self.resname}: " + "; ".join(errors))

    def _parameterization_copy(self):
        attributes = ["_stype", "_smap"]
        attributes += [f"_{kind}{suffix}" for kind in "badi" for suffix in ("type", "map", "ids")]
        attributes += [f"_old_{kind}" for kind in ("bonds", "angles", "dihedrals", "impropers")]
        staged = copy.copy(self)
        staged._atoms = self.atoms.copy()
        staged.__dict__.update(copy.deepcopy({key: getattr(self, key) for key in attributes}))
        return staged, attributes

    def _apply_parameterization(self, staged, attributes):
        self.atoms.set_initial_charges(staged.charges)
        self._sids = staged._sids.copy()
        self._graph = staged._graph
        self.atoms.info.update(staged.atoms.info)
        for name in ("nominal_charge", "atom_map"):
            if name in staged.atoms.arrays:
                self.atoms.set_array(name, staged.atoms.arrays[name].copy())
        self.__dict__.update({key: getattr(staged, key) for key in attributes})

    def parameterize(self, snippet_radius=12, charge_correction="none",
                     cap_element="H", segment_size=200):
        """Assign LigParGen OPLS-AA parameters atomically.

        Parameters
        ----------
        snippet_radius : int, default 12
            Graph radius of junction snippets when segmentation is needed.
        charge_correction : {"none", "uniform"}, default "none"
            Optional uniform correction to the stored molecular charge.
        cap_element : str, default "H"
            Neutral monovalent cap used at segment boundaries.
        segment_size : int, default 200
            Maximum atoms per segment, including caps, at most 200.

        Returns
        -------
        dict
            Charge audit with initial, refined and final charges, formal
            target, residual, correction per atom, and junction count.

        Raises
        ------
        ValueError
            If options, returned parameters, or capped fragment sizes are invalid.
        RuntimeError
            If fragmentation or parameterization fails. Original parameters
            and coordinates remain unchanged.

        Notes
        -----
        Small molecules are parameterized directly. Larger molecules use
        capped segments and junction refinement. Coordinates are preserved.
        """
        from mdinterface.externals.ligpargen import refine_large_specie_topology
        return refine_large_specie_topology(
            self, snippet_radius=snippet_radius, charge_correction=charge_correction,
            cap_element=cap_element, segment_size=segment_size)

    def refine_large_topology(self, snippet_radius=12, cap_element="H",
                              charge_correction="none", segment_size=200):
        """Parameterize using the same atomic workflow as :meth:`parameterize`.

        Parameters
        ----------
        snippet_radius : int, default 12
            Graph radius for junction snippets.
        cap_element : str, default "H"
            Neutral monovalent capping element.
        charge_correction : {"none", "uniform"}, default "none"
            Optional uniform correction to the stored total charge.
        segment_size : int, default 200
            Maximum segment size including caps.

        Returns
        -------
        dict
            Parameterization charge audit.
        """
        return self.parameterize(snippet_radius, charge_correction, cap_element, segment_size)

    def mark_attachment_sites(self, head_map, tail_map, leaving_element="H"):
        """Mark neutral terminal leaving atoms using RDKit atom-map identifiers.

        Parameters
        ----------
        head_map, tail_map : int
            Maps of the head and tail leaving atoms, or their anchor atoms.
            Anchors with multiple candidate leaving atoms are accepted only
            for equivalent hydrogens on a non-stereogenic center.
        leaving_element : str, default "H"
            Element of the neutral singly bonded leaving atoms.

        Raises
        ------
        ValueError
            If maps are absent, sites overlap, or a leaving atom is ambiguous,
            charged, or not terminal. Existing annotations remain unchanged.
        """
        mol = self.to_rdkit()
        by_map = {a.GetAtomMapNum(): a for a in mol.GetAtoms()}
        selected = []
        for number in (head_map, tail_map):
            if number not in by_map:
                raise ValueError(f"Atom map {number} is not present in the molecule.")
            atom = by_map[number]
            candidates = ([atom] if atom.GetSymbol() == leaving_element else
                          [a for a in atom.GetNeighbors() if a.GetSymbol() == leaving_element])
            candidates = [a for a in candidates if a.GetDegree() == 1 and a.GetFormalCharge() == 0
                          and mol.GetBondBetweenAtoms(a.GetIdx(), a.GetNeighbors()[0].GetIdx()).GetBondType() == Chem.BondType.SINGLE]
            if len(candidates) > 1:
                equivalent_h = leaving_element == "H" and len({a.GetIsotope() for a in candidates}) == 1
                if not equivalent_h or atom.GetChiralTag() != Chem.ChiralType.CHI_UNSPECIFIED:
                    raise ValueError(f"Ambiguous leaving atoms at map {number}; specify the leaving atom's map.")
            if not candidates:
                raise ValueError(f"Map {number} has no neutral terminal {leaving_element} leaving atom.")
            selected.append(candidates[0].GetIdx())
        if selected[0] == selected[1]:
            raise ValueError("Head and tail must select different leaving atoms.")
        marks = np.zeros(len(self.atoms), dtype=int)
        marks[selected] = [1, 2]
        self.atoms.set_array("polymerize", marks)

    def write_gro(self, filename=None):
        """
        Write a GROMACS structure (.gro) file for this species.

        Parameters
        ----------
        filename : str, optional
            Output filename. Defaults to ``{resname}.gro``.
        """
        if filename is None:
            filename = f"{self.resname}.gro"
        universe = self.to_universe()
        universe.atoms.write(filename)

    def _resolve_charge(self, charge):
        target = self._tot_charge
        mol = stored_molecule(self.atoms)
        if mol is not None and Chem.GetFormalCharge(mol) != target:
            raise ValueError("Stored charge conflicts with the chemical graph.")
        if charge is not None and charge != target:
            raise ValueError("charge conflicts with the stored molecular charge; set tot_charge when constructing the species.")
        return target

    def estimate_charges(self, method="obabel", charge=None, assign=False, **respargs):
        """Estimate partial charges using the stored molecular charge.

        Parameters
        ----------
        method : {"obabel", "ligpargen", "resp"}, default "obabel"
            Charge estimation backend.
        charge : int, optional
            Defaults to the stored total charge. Conflicting values raise.
        assign : bool, default False
            Assign the result, including optimized RESP coordinates if used.
        **respargs
            Additional RESP calculation options.

        Returns
        -------
        numpy.ndarray
            Estimated partial charges in elementary-charge units.
        """
        charge = self._resolve_charge(charge)
        if method == "obabel":
            charges = run_OBChargeModel(self.atoms, charge=charge)
        elif method == "ligpargen":
            system, _, _, _, _, _ = run_ligpargen(self.atoms, charge=charge)
            charges = system.get_initial_charges()
        elif method == "resp":
            charges, atoms = calculate_RESP_charges(self, charge=charge, **respargs)
        else:
            raise ValueError("method must be 'obabel', 'ligpargen', or 'resp'.")
        if assign:
            self.atoms.set_initial_charges(charges)
            if method == "resp":
                self.atoms.set_positions(atoms.get_positions())
        return charges

    def estimate_OPLSAA_parameters(self, charge=None):
        """Return LigParGen parameters using the stored molecular charge.

        Parameters
        ----------
        charge : int, optional
            Defaults to the stored total charge. Conflicting values raise.

        Returns
        -------
        tuple
            Parameterized atoms, atom types, bonds, angles, dihedrals and
            impropers. The species is not modified.
        """
        return run_ligpargen(self.atoms, charge=self._resolve_charge(charge))

    def relax_structure(self, optimizer='FIRE', fmax=0.05,
                       steps=200, update_positions=True, trajectory=None,
                       logfile=None, charge=None, spin=0, **kwargs):
        """
        Relax the Specie structure using ASE or UMA optimization.

        Parameters
        ----------
        optimizer : str, default 'FIRE'
            Optimizer for ASE method ('BFGS', 'LBFGS', 'FIRE')
        fmax : float, default 0.05
            Maximum force threshold for convergence (eV/Å)
        steps : int, default 200
            Maximum number of optimization steps
        update_positions : bool, default True
            Whether to update the Specie's atomic positions after relaxation
        trajectory : str, optional
            Path to save optimization trajectory
        logfile : str, optional
            Path to save optimization log
        charge : int or None, default None
            Total charge of the system. Defaults to the stored tot_charge.
        spin : int, default 0
            Total spin multiplicity (2S). Use >0 for open-shell systems.
        **kwargs
            Additional arguments passed to the relaxation function

        Returns
        -------
        relaxed_atoms : ase.Atoms
            The relaxed atomic structure
        converged : bool
            Whether the optimization converged
        """

        if charge is None:
            charge = self._tot_charge

        atoms_to_relax = self.atoms
        atoms_to_relax.info["charge"] = charge
        atoms_to_relax.info["spin"]   = spin

        relaxed_atoms = relax_structure(
            atoms_to_relax, optimizer=optimizer, fmax=fmax, steps=steps,
            trajectory=trajectory, logfile=logfile, **kwargs)

        # Update positions if requested
        if update_positions:
            self.update_positions(atoms=relaxed_atoms, prune_z=False)

        return

    def run_aimd(self, timestep=0.5, temperature_K=300, friction=0.001,
              steps=1000, update_positions=True, trajectory=None,
              logfile=None, **kwargs):
        """
        Run AIMD for the Specie structure using ASE and FAIRChem.

        Parameters
        ----------
        timestep : float, default 0.1
            Time step for the simulation in fs.
        temperature_K : float, default 300
            Target temperature for the Langevin dynamics in Kelvin.
        friction : float, default 0.001
            Frictional damping coefficient in 1/fs.
        steps : int, default 1000
            Number of time steps to run the AIMD.
        update_positions : bool, default True
            Whether to update the Specie's atomic positions after AIMD.
        trajectory : str, optional
            Path to save the MD trajectory.
        logfile : str, optional
            Path to save the MD log.
        **kwargs
            Additional arguments passed to the Langevin integrator.

        Returns
        -------
        None
        """

        atoms_to_simulate = self.atoms
        atoms_to_simulate.info["charge"] = self._tot_charge
        atoms_to_simulate.info["spin"] = 0  # Update this if needed

        # Call the previously defined function to perform AIMD
        run_aimd(atoms_to_simulate, timestep=timestep, temperature_K=temperature_K,
                 friction=friction, steps=steps, trajectory=trajectory,
                 logfile=logfile, **kwargs)

        # Update positions if requested
        if update_positions:
            self.update_positions(atoms=atoms_to_simulate, prune_z=False)

        return

    def _find_rings(self, max_ring_size=8):
        return find_rings(self.graph, max_ring_size=max_ring_size)

    def _get_rings_containing_atoms(self, atom_indices, max_ring_size=8):
        return get_rings_containing_atoms(self.graph, atom_indices, max_ring_size=max_ring_size)

to_rdkit()

Return an independent RDKit molecule with current coordinates.

Returns:

Type	Description
Mol	Explicit-hydrogen molecular graph in ASE atom order. For inputs without stored chemistry, bonds are perceived from coordinates using the specified total charge.

Raises:

Type	Description
ValueError	If coordinate-based perception is requested for a periodic structure, or chemistry conflicts with formal-charge annotations.

Source code in mdinterface/core/specie.py

def to_rdkit(self):
    """Return an independent RDKit molecule with current coordinates.

    Returns
    -------
    rdkit.Chem.Mol
        Explicit-hydrogen molecular graph in ASE atom order. For inputs
        without stored chemistry, bonds are perceived from coordinates
        using the specified total charge.

    Raises
    ------
    ValueError
        If coordinate-based perception is requested for a periodic
        structure, or chemistry conflicts with formal-charge annotations.
    """
    return perceive_molecule(self.atoms, charge=self._tot_charge)

generate_conformer(seed=0, minimize=False, max_iterations=1000)

Generate molecular coordinates with RDKit ETKDGv3.

Parameters:

Name	Type	Description	Default
seed	int	Nonnegative 32-bit embedding seed.	0
minimize	bool	Also minimize the new coordinates with MMFF94.	False
max_iterations	int	Maximum MMFF94 iterations when minimize=True.	1000

Returns:

Type	Description
float or None	MMFF94 energy in kcal/mol when minimized, otherwise None.

Raises:

Type	Description
ValueError	If options or chemistry are invalid, MMFF parameters are missing, or the structure is periodic, disconnected, or constrained.
RuntimeError	If embedding or minimization fails. Coordinates are unchanged.

Notes

Preserves atom ordering, topology and charges. A conformer is an isolated-molecule starting geometry, not an equilibrated ensemble.

Source code in mdinterface/core/specie.py

def generate_conformer(self, seed=0, minimize=False, max_iterations=1000):
    """Generate molecular coordinates with RDKit ETKDGv3.

    Parameters
    ----------
    seed : int, default 0
        Nonnegative 32-bit embedding seed.
    minimize : bool, default False
        Also minimize the new coordinates with MMFF94.
    max_iterations : int, default 1000
        Maximum MMFF94 iterations when ``minimize=True``.

    Returns
    -------
    float or None
        MMFF94 energy in kcal/mol when minimized, otherwise None.

    Raises
    ------
    ValueError
        If options or chemistry are invalid, MMFF parameters are missing,
        or the structure is periodic, disconnected, or constrained.
    RuntimeError
        If embedding or minimization fails. Coordinates are unchanged.

    Notes
    -----
    Preserves atom ordering, topology and charges. A conformer is an
    isolated-molecule starting geometry, not an equilibrated ensemble.
    """
    from rdkit.Chem import rdDistGeom
    from mdinterface.core.chemistry import minimize_molecule

    if not isinstance(seed, (int, np.integer)) or not 0 <= seed < 2**31:
        raise ValueError("seed must be a nonnegative 32-bit integer.")
    if not isinstance(max_iterations, (int, np.integer)) or max_iterations < 1:
        raise ValueError("max_iterations must be a positive integer.")
    mol = self._geometry_molecule()
    params = rdDistGeom.ETKDGv3()
    params.randomSeed = int(seed)
    if rdDistGeom.EmbedMolecule(mol, params) < 0:
        raise RuntimeError("RDKit could not embed the molecule; try another seed.")
    energy = minimize_molecule(mol, max_iterations) if minimize else None
    positions = mol.GetConformer().GetPositions()
    if not np.isfinite(positions).all():
        raise RuntimeError("RDKit produced nonfinite coordinates.")
    self.atoms.set_positions(positions)
    return energy

minimize_geometry(max_iterations=1000)

Minimize existing molecular coordinates with RDKit MMFF94.

Parameters:

Name	Type	Description	Default
max_iterations	int	Maximum minimization iterations.	1000

Returns:

Type	Description
float	MMFF94 energy in kcal/mol, not an OPLS-AA energy.

Raises:

Type	Description
ValueError	If arguments or chemistry are invalid, parameters are unavailable, or the structure is periodic, disconnected, or constrained.
RuntimeError	If minimization fails. Original coordinates remain unchanged.

Notes

Does not embed a new conformer or change simulation parameters.

Source code in mdinterface/core/specie.py

def minimize_geometry(self, max_iterations=1000):
    """Minimize existing molecular coordinates with RDKit MMFF94.

    Parameters
    ----------
    max_iterations : int, default 1000
        Maximum minimization iterations.

    Returns
    -------
    float
        MMFF94 energy in kcal/mol, not an OPLS-AA energy.

    Raises
    ------
    ValueError
        If arguments or chemistry are invalid, parameters are unavailable,
        or the structure is periodic, disconnected, or constrained.
    RuntimeError
        If minimization fails. Original coordinates remain unchanged.

    Notes
    -----
    Does not embed a new conformer or change simulation parameters.
    """
    from mdinterface.core.chemistry import minimize_molecule

    mol = self._geometry_molecule()
    energy = minimize_molecule(mol, max_iterations)
    self.atoms.set_positions(mol.GetConformer().GetPositions())
    return energy

to_smiles(mapped=False)

Return canonical isomeric SMILES for the molecular structure.

Parameters:

Name	Type	Description	Default
mapped	bool	Include persistent atom-map numbers.	False

Returns:

Type	Description
str	SMILES with formal charges and stereochemistry.

Source code in mdinterface/core/specie.py

def to_smiles(self, mapped=False):
    """Return canonical isomeric SMILES for the molecular structure.

    Parameters
    ----------
    mapped : bool, default False
        Include persistent atom-map numbers.

    Returns
    -------
    str
        SMILES with formal charges and stereochemistry.
    """
    mol = self.to_rdkit()
    if not mapped:
        for atom in mol.GetAtoms():
            atom.SetAtomMapNum(0)
    return Chem.MolToSmiles(Chem.RemoveHs(mol), isomericSmiles=True)

write_gromacs_itp(filename=None, *, include_atomtypes=True)

Write a GROMACS include topology (.itp) file for this species.

See :func:~mdinterface.io.gromacswriter.write_gromacs_itp for full documentation.

Parameters:

Name	Type	Description	Default
filename	str	Output filename. Defaults to {resname}.itp.	None
include_atomtypes	bool	Include atom types in this file. For multi-species exports, use False and supply species to write_gromacs_top to define all atom types before all molecule definitions.	True

Source code in mdinterface/core/specie.py

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

    See :func:`~mdinterface.io.gromacswriter.write_gromacs_itp` for full
    documentation.

    Parameters
    ----------
    filename : str, optional
        Output filename. Defaults to ``{resname}.itp``.
    include_atomtypes : bool, default True
        Include atom types in this file. For multi-species exports, use
        False and supply ``species`` to ``write_gromacs_top`` to define
        all atom types before all molecule definitions.
    """
    logger.warning(
        "write_gromacs_itp is experimental -- verify output before production use."
    )
    from mdinterface.io.gromacswriter import write_gromacs_itp
    write_gromacs_itp(self, filename=filename, include_atomtypes=include_atomtypes)

validate_force_field()

Raise if molecular parameters required for classical export are missing.

Returns:

Type	Description
None	Successful validation returns None.

Raises:

Type	Description
ValueError	If charges or coefficients are nonfinite, pair parameters are missing, or molecular bonds, angles or proper torsions are unassigned.

Notes

Completeness uses explicit chemistry or an existing bonded topology. Improper coefficients are checked when assigned; missing impropers cannot be inferred universally from connectivity. Constraints and force-field compatibility remain the caller's responsibility.

Source code in mdinterface/core/specie.py

def validate_force_field(self):
    """Raise if molecular parameters required for classical export are missing.

    Returns
    -------
    None
        Successful validation returns None.

    Raises
    ------
    ValueError
        If charges or coefficients are nonfinite, pair parameters are
        missing, or molecular bonds, angles or proper torsions are unassigned.

    Notes
    -----
    Completeness uses explicit chemistry or an existing bonded topology.
    Improper coefficients are checked when assigned; missing impropers
    cannot be inferred universally from connectivity. Constraints and
    force-field compatibility remain the caller's responsibility.
    """
    errors = []
    if not np.isfinite(self.charges).all():
        errors.append("nonfinite partial charges")
    for index in set(self._smap[sid] for sid in self._sids):
        atom = self._stype[index]
        if any(value is None or not np.isfinite(value) or value < 0 for value in (atom.eps, atom.sig)):
            errors.append(f"missing or invalid pair parameters for {atom.label}")
    molecular = stored_molecule(self.atoms) is not None or any((self._old_bonds, self._old_angles, self._old_dihedrals))
    for kind, type_key, finder in (("bonds", "_btype", pmap.find_missing_bonds),
                                   ("angles", "_atype", pmap.find_missing_angles),
                                   ("dihedrals", "_dtype", pmap.find_missing_dihedrals),
                                   ("impropers", "_itype", None)):
        if molecular and finder is not None:
            missing = finder(self)
            if missing:
                errors.append(f"{len(missing)} unassigned {kind} (first: {missing[0]})")
        _, indices = getattr(self, kind)
        for index in set(indices):
            parameter = getattr(self, type_key)[index]
            values = parameter.values[:4] if kind == "dihedrals" and parameter.values[-1] is None else parameter.values
            if any(value is None or not np.isfinite(value) for value in values):
                errors.append(f"missing or nonfinite {kind} coefficients for {parameter.symbols}")
    if errors:
        raise ValueError(f"Incomplete force field for {self.resname}: " + "; ".join(errors))

parameterize(snippet_radius=12, charge_correction='none', cap_element='H', segment_size=200)

Assign LigParGen OPLS-AA parameters atomically.

Parameters:

Name	Type	Description	Default
snippet_radius	int	Graph radius of junction snippets when segmentation is needed.	12
charge_correction	(none, uniform)	Optional uniform correction to the stored molecular charge.	"none"
cap_element	str	Neutral monovalent cap used at segment boundaries.	"H"
segment_size	int	Maximum atoms per segment, including caps, at most 200.	200

Returns:

Type	Description
dict	Charge audit with initial, refined and final charges, formal target, residual, correction per atom, and junction count.

Raises:

Type	Description
ValueError	If options, returned parameters, or capped fragment sizes are invalid.
RuntimeError	If fragmentation or parameterization fails. Original parameters and coordinates remain unchanged.

Notes

Small molecules are parameterized directly. Larger molecules use capped segments and junction refinement. Coordinates are preserved.

Source code in mdinterface/core/specie.py

def parameterize(self, snippet_radius=12, charge_correction="none",
                 cap_element="H", segment_size=200):
    """Assign LigParGen OPLS-AA parameters atomically.

    Parameters
    ----------
    snippet_radius : int, default 12
        Graph radius of junction snippets when segmentation is needed.
    charge_correction : {"none", "uniform"}, default "none"
        Optional uniform correction to the stored molecular charge.
    cap_element : str, default "H"
        Neutral monovalent cap used at segment boundaries.
    segment_size : int, default 200
        Maximum atoms per segment, including caps, at most 200.

    Returns
    -------
    dict
        Charge audit with initial, refined and final charges, formal
        target, residual, correction per atom, and junction count.

    Raises
    ------
    ValueError
        If options, returned parameters, or capped fragment sizes are invalid.
    RuntimeError
        If fragmentation or parameterization fails. Original parameters
        and coordinates remain unchanged.

    Notes
    -----
    Small molecules are parameterized directly. Larger molecules use
    capped segments and junction refinement. Coordinates are preserved.
    """
    from mdinterface.externals.ligpargen import refine_large_specie_topology
    return refine_large_specie_topology(
        self, snippet_radius=snippet_radius, charge_correction=charge_correction,
        cap_element=cap_element, segment_size=segment_size)

refine_large_topology(snippet_radius=12, cap_element='H', charge_correction='none', segment_size=200)

Parameterize using the same atomic workflow as :meth:parameterize.

Parameters:

Name	Type	Description	Default
snippet_radius	int	Graph radius for junction snippets.	12
cap_element	str	Neutral monovalent capping element.	"H"
charge_correction	(none, uniform)	Optional uniform correction to the stored total charge.	"none"
segment_size	int	Maximum segment size including caps.	200

Returns:

Type	Description
dict	Parameterization charge audit.

Source code in mdinterface/core/specie.py

def refine_large_topology(self, snippet_radius=12, cap_element="H",
                          charge_correction="none", segment_size=200):
    """Parameterize using the same atomic workflow as :meth:`parameterize`.

    Parameters
    ----------
    snippet_radius : int, default 12
        Graph radius for junction snippets.
    cap_element : str, default "H"
        Neutral monovalent capping element.
    charge_correction : {"none", "uniform"}, default "none"
        Optional uniform correction to the stored total charge.
    segment_size : int, default 200
        Maximum segment size including caps.

    Returns
    -------
    dict
        Parameterization charge audit.
    """
    return self.parameterize(snippet_radius, charge_correction, cap_element, segment_size)

mark_attachment_sites(head_map, tail_map, leaving_element='H')

Mark neutral terminal leaving atoms using RDKit atom-map identifiers.

Parameters:

Name	Type	Description	Default
head_map	int	Maps of the head and tail leaving atoms, or their anchor atoms. Anchors with multiple candidate leaving atoms are accepted only for equivalent hydrogens on a non-stereogenic center.	required
tail_map	int	Maps of the head and tail leaving atoms, or their anchor atoms. Anchors with multiple candidate leaving atoms are accepted only for equivalent hydrogens on a non-stereogenic center.	required
leaving_element	str	Element of the neutral singly bonded leaving atoms.	"H"

Raises:

Type	Description
ValueError	If maps are absent, sites overlap, or a leaving atom is ambiguous, charged, or not terminal. Existing annotations remain unchanged.

Source code in mdinterface/core/specie.py

def mark_attachment_sites(self, head_map, tail_map, leaving_element="H"):
    """Mark neutral terminal leaving atoms using RDKit atom-map identifiers.

    Parameters
    ----------
    head_map, tail_map : int
        Maps of the head and tail leaving atoms, or their anchor atoms.
        Anchors with multiple candidate leaving atoms are accepted only
        for equivalent hydrogens on a non-stereogenic center.
    leaving_element : str, default "H"
        Element of the neutral singly bonded leaving atoms.

    Raises
    ------
    ValueError
        If maps are absent, sites overlap, or a leaving atom is ambiguous,
        charged, or not terminal. Existing annotations remain unchanged.
    """
    mol = self.to_rdkit()
    by_map = {a.GetAtomMapNum(): a for a in mol.GetAtoms()}
    selected = []
    for number in (head_map, tail_map):
        if number not in by_map:
            raise ValueError(f"Atom map {number} is not present in the molecule.")
        atom = by_map[number]
        candidates = ([atom] if atom.GetSymbol() == leaving_element else
                      [a for a in atom.GetNeighbors() if a.GetSymbol() == leaving_element])
        candidates = [a for a in candidates if a.GetDegree() == 1 and a.GetFormalCharge() == 0
                      and mol.GetBondBetweenAtoms(a.GetIdx(), a.GetNeighbors()[0].GetIdx()).GetBondType() == Chem.BondType.SINGLE]
        if len(candidates) > 1:
            equivalent_h = leaving_element == "H" and len({a.GetIsotope() for a in candidates}) == 1
            if not equivalent_h or atom.GetChiralTag() != Chem.ChiralType.CHI_UNSPECIFIED:
                raise ValueError(f"Ambiguous leaving atoms at map {number}; specify the leaving atom's map.")
        if not candidates:
            raise ValueError(f"Map {number} has no neutral terminal {leaving_element} leaving atom.")
        selected.append(candidates[0].GetIdx())
    if selected[0] == selected[1]:
        raise ValueError("Head and tail must select different leaving atoms.")
    marks = np.zeros(len(self.atoms), dtype=int)
    marks[selected] = [1, 2]
    self.atoms.set_array("polymerize", marks)

write_gro(filename=None)

Write a GROMACS structure (.gro) file for this species.

Parameters:

Name	Type	Description	Default
filename	str	Output filename. Defaults to {resname}.gro.	None

Source code in mdinterface/core/specie.py

def write_gro(self, filename=None):
    """
    Write a GROMACS structure (.gro) file for this species.

    Parameters
    ----------
    filename : str, optional
        Output filename. Defaults to ``{resname}.gro``.
    """
    if filename is None:
        filename = f"{self.resname}.gro"
    universe = self.to_universe()
    universe.atoms.write(filename)

estimate_charges(method='obabel', charge=None, assign=False, **respargs)

Estimate partial charges using the stored molecular charge.

Parameters:

Name	Type	Description	Default
method	(obabel, ligpargen, resp)	Charge estimation backend.	"obabel"
charge	int	Defaults to the stored total charge. Conflicting values raise.	None
assign	bool	Assign the result, including optimized RESP coordinates if used.	False
**respargs		Additional RESP calculation options.	{}

Returns:

Type	Description
ndarray	Estimated partial charges in elementary-charge units.

Source code in mdinterface/core/specie.py

def estimate_charges(self, method="obabel", charge=None, assign=False, **respargs):
    """Estimate partial charges using the stored molecular charge.

    Parameters
    ----------
    method : {"obabel", "ligpargen", "resp"}, default "obabel"
        Charge estimation backend.
    charge : int, optional
        Defaults to the stored total charge. Conflicting values raise.
    assign : bool, default False
        Assign the result, including optimized RESP coordinates if used.
    **respargs
        Additional RESP calculation options.

    Returns
    -------
    numpy.ndarray
        Estimated partial charges in elementary-charge units.
    """
    charge = self._resolve_charge(charge)
    if method == "obabel":
        charges = run_OBChargeModel(self.atoms, charge=charge)
    elif method == "ligpargen":
        system, _, _, _, _, _ = run_ligpargen(self.atoms, charge=charge)
        charges = system.get_initial_charges()
    elif method == "resp":
        charges, atoms = calculate_RESP_charges(self, charge=charge, **respargs)
    else:
        raise ValueError("method must be 'obabel', 'ligpargen', or 'resp'.")
    if assign:
        self.atoms.set_initial_charges(charges)
        if method == "resp":
            self.atoms.set_positions(atoms.get_positions())
    return charges

estimate_OPLSAA_parameters(charge=None)

Return LigParGen parameters using the stored molecular charge.

Parameters:

Name	Type	Description	Default
charge	int	Defaults to the stored total charge. Conflicting values raise.	None

Returns:

Type	Description
tuple	Parameterized atoms, atom types, bonds, angles, dihedrals and impropers. The species is not modified.

Source code in mdinterface/core/specie.py

def estimate_OPLSAA_parameters(self, charge=None):
    """Return LigParGen parameters using the stored molecular charge.

    Parameters
    ----------
    charge : int, optional
        Defaults to the stored total charge. Conflicting values raise.

    Returns
    -------
    tuple
        Parameterized atoms, atom types, bonds, angles, dihedrals and
        impropers. The species is not modified.
    """
    return run_ligpargen(self.atoms, charge=self._resolve_charge(charge))

relax_structure(optimizer='FIRE', fmax=0.05, steps=200, update_positions=True, trajectory=None, logfile=None, charge=None, spin=0, **kwargs)

Relax the Specie structure using ASE or UMA optimization.

Parameters:

Name	Type	Description	Default
optimizer	str	Optimizer for ASE method ('BFGS', 'LBFGS', 'FIRE')	'FIRE'
fmax	float	Maximum force threshold for convergence (eV/Å)	0.05
steps	int	Maximum number of optimization steps	200
update_positions	bool	Whether to update the Specie's atomic positions after relaxation	True
trajectory	str	Path to save optimization trajectory	None
logfile	str	Path to save optimization log	None
charge	int or None	Total charge of the system. Defaults to the stored tot_charge.	None
spin	int	Total spin multiplicity (2S). Use >0 for open-shell systems.	0
**kwargs		Additional arguments passed to the relaxation function	{}

Returns:

Name	Type	Description
relaxed_atoms	Atoms	The relaxed atomic structure
converged	bool	Whether the optimization converged

Source code in mdinterface/core/specie.py

def relax_structure(self, optimizer='FIRE', fmax=0.05,
                   steps=200, update_positions=True, trajectory=None,
                   logfile=None, charge=None, spin=0, **kwargs):
    """
    Relax the Specie structure using ASE or UMA optimization.

    Parameters
    ----------
    optimizer : str, default 'FIRE'
        Optimizer for ASE method ('BFGS', 'LBFGS', 'FIRE')
    fmax : float, default 0.05
        Maximum force threshold for convergence (eV/Å)
    steps : int, default 200
        Maximum number of optimization steps
    update_positions : bool, default True
        Whether to update the Specie's atomic positions after relaxation
    trajectory : str, optional
        Path to save optimization trajectory
    logfile : str, optional
        Path to save optimization log
    charge : int or None, default None
        Total charge of the system. Defaults to the stored tot_charge.
    spin : int, default 0
        Total spin multiplicity (2S). Use >0 for open-shell systems.
    **kwargs
        Additional arguments passed to the relaxation function

    Returns
    -------
    relaxed_atoms : ase.Atoms
        The relaxed atomic structure
    converged : bool
        Whether the optimization converged
    """

    if charge is None:
        charge = self._tot_charge

    atoms_to_relax = self.atoms
    atoms_to_relax.info["charge"] = charge
    atoms_to_relax.info["spin"]   = spin

    relaxed_atoms = relax_structure(
        atoms_to_relax, optimizer=optimizer, fmax=fmax, steps=steps,
        trajectory=trajectory, logfile=logfile, **kwargs)

    # Update positions if requested
    if update_positions:
        self.update_positions(atoms=relaxed_atoms, prune_z=False)

    return

run_aimd(timestep=0.5, temperature_K=300, friction=0.001, steps=1000, update_positions=True, trajectory=None, logfile=None, **kwargs)

Run AIMD for the Specie structure using ASE and FAIRChem.

Parameters:

Name	Type	Description	Default
timestep	float	Time step for the simulation in fs.	0.1
temperature_K	float	Target temperature for the Langevin dynamics in Kelvin.	300
friction	float	Frictional damping coefficient in 1/fs.	0.001
steps	int	Number of time steps to run the AIMD.	1000
update_positions	bool	Whether to update the Specie's atomic positions after AIMD.	True
trajectory	str	Path to save the MD trajectory.	None
logfile	str	Path to save the MD log.	None
**kwargs		Additional arguments passed to the Langevin integrator.	{}

Returns:

Type	Description
None

Source code in mdinterface/core/specie.py

def run_aimd(self, timestep=0.5, temperature_K=300, friction=0.001,
          steps=1000, update_positions=True, trajectory=None,
          logfile=None, **kwargs):
    """
    Run AIMD for the Specie structure using ASE and FAIRChem.

    Parameters
    ----------
    timestep : float, default 0.1
        Time step for the simulation in fs.
    temperature_K : float, default 300
        Target temperature for the Langevin dynamics in Kelvin.
    friction : float, default 0.001
        Frictional damping coefficient in 1/fs.
    steps : int, default 1000
        Number of time steps to run the AIMD.
    update_positions : bool, default True
        Whether to update the Specie's atomic positions after AIMD.
    trajectory : str, optional
        Path to save the MD trajectory.
    logfile : str, optional
        Path to save the MD log.
    **kwargs
        Additional arguments passed to the Langevin integrator.

    Returns
    -------
    None
    """

    atoms_to_simulate = self.atoms
    atoms_to_simulate.info["charge"] = self._tot_charge
    atoms_to_simulate.info["spin"] = 0  # Update this if needed

    # Call the previously defined function to perform AIMD
    run_aimd(atoms_to_simulate, timestep=timestep, temperature_K=temperature_K,
             friction=friction, steps=steps, trajectory=trajectory,
             logfile=logfile, **kwargs)

    # Update positions if requested
    if update_positions:
        self.update_positions(atoms=atoms_to_simulate, prune_z=False)

    return