Externals

Optional integrations with third-party tools.

LigParGen (OPLS-AA parameters)

mdinterface.externals.ligpargen.refine_large_specie_topology(specie, Nmax=12, ending='H', offset=True, segment_size=200)

Assign OPLS-AA force-field parameters to a Specie via LigParGen using a segment-and-junction strategy.

The molecule is split into segments of at most segment_size atoms along clean backbone bonds (avoiding rings, heteroatoms, and their immediate neighbours). Each segment is capped with ending atoms and passed to LigParGen independently. A local snippet centred on each cut bond is then refined to correct parameters for atoms adjacent to a cap.

Parameters:

Name	Type	Description	Default
specie	Specie	The molecule to parametrise. Modified in place.	required
Nmax	int	Neighbourhood radius (in bonds) used for junction snippet creation. Default 12.	12
ending	str	Element used to cap dangling bonds at segment boundaries. Default "H".	'H'
offset	bool	If True, redistribute any charge rounding error uniformly so the total partial charge matches nominal_charge.sum().	True
segment_size	int	Maximum number of atoms per segment. Default 200. Lower values force more splits and are useful for testing.	200

Raises:

Type	Description
ValueError	If nominal_charge is not set on specie.atoms.
RuntimeError	If not enough valid cut bonds can be found.

Source code in mdinterface/externals/ligpargen.py

def refine_large_specie_topology(specie, Nmax=12, ending="H", offset=True,
                                 segment_size=200):
    """
    Assign OPLS-AA force-field parameters to a ``Specie`` via LigParGen using
    a segment-and-junction strategy.

    The molecule is split into segments of at most *segment_size* atoms along
    clean backbone bonds (avoiding rings, heteroatoms, and their immediate
    neighbours).  Each segment is capped with *ending* atoms and passed to
    LigParGen independently.  A local snippet centred on each cut bond is then
    refined to correct parameters for atoms adjacent to a cap.

    Parameters
    ----------
    specie : Specie
        The molecule to parametrise.  Modified **in place**.
    Nmax : int
        Neighbourhood radius (in bonds) used for junction snippet creation.
        Default 12.
    ending : str
        Element used to cap dangling bonds at segment boundaries.  Default
        ``"H"``.
    offset : bool
        If ``True``, redistribute any charge rounding error uniformly so the
        total partial charge matches ``nominal_charge.sum()``.
    segment_size : int
        Maximum number of atoms per segment.  Default 200.  Lower values
        force more splits and are useful for testing.

    Raises
    ------
    ValueError
        If ``nominal_charge`` is not set on ``specie.atoms``.
    RuntimeError
        If not enough valid cut bonds can be found.
    """
    from mdinterface.build.snippets import make_snippet, remap_snippet_topology

    natoms = len(specie.atoms)

    if "nominal_charge" not in specie.atoms.arrays:
        raise ValueError(
            "nominal_charge array not found on specie.atoms.  "
            "Set it before calling refine_large_specie_topology() -- "
            "use 0 for uncharged atoms and the formal integer charge for "
            "charged ones (same convention as the Polymer class)."
        )

    n_segments = math.ceil(natoms / segment_size)
    n_cuts     = n_segments - 1
    logger.info("Refining topology for %d-atom molecule -- %d segment(s) via LigParGen",
                natoms, n_segments)

    # ------------------------------------------------------------------
    # 1. Find cut bonds
    # ------------------------------------------------------------------
    candidates = _candidate_cut_bonds(specie, n_needed=n_cuts)
    if len(candidates) < n_cuts:
        raise RuntimeError(
            f"Only {len(candidates)} valid cut bond(s) found but "
            f"{n_cuts} cut(s) are needed.  Consider relaxing the molecule "
            f"structure or cutting manually."
        )

    cut_edges, segments = _balanced_cuts(specie, n_cuts, candidates)
    for ii, seg in enumerate(segments):
        logger.info("  >> segment %d: %d atoms", ii, len(seg))
    logger.info("  >> cut bond(s): %s",
                ", ".join(f"{int(a)} -- {int(b)}" for a, b in cut_edges) or "none")

    # ------------------------------------------------------------------
    # 2. Run LigParGen on each segment, accumulate topology
    # ------------------------------------------------------------------
    charges            = specie.charges.copy()
    atom_types_ordered = [None] * natoms  # Atom-type objects in original order
    all_bonds, all_angles, all_dihs, all_imps = [], [], [], []

    for seg_idx, seg_set in enumerate(segments):
        seg_indices = sorted(seg_set)
        logger.info("  >> ligpargen on segment %d (%d atoms)...",
                    seg_idx, len(seg_indices))

        capped, n_real = _make_capped_segment(
            specie, seg_indices, cut_edges, ending)
        sn_charge = int(capped.arrays["nominal_charge"].sum())

        sn_sys, sn_atypes, sn_bonds, sn_angles, sn_dihs, sn_imps = \
            run_ligpargen(capped, charge=sn_charge, is_snippet=True)

        # Remap: real atoms -> original sids; caps -> placeholders
        real_sids     = list(specie._sids[seg_indices])
        cap_sids      = [f"__cap_{seg_idx}_{c}__"
                         for c in range(len(sn_atypes) - n_real)]
        original_idxs = np.array(real_sids + cap_sids)
        local_idxs    = np.array(real_sids)

        b, a, d, i = remap_snippet_topology(
            original_idxs, sn_sys, sn_atypes,
            sn_bonds, sn_angles, sn_dihs, sn_imps,
            local_idxs)
        all_bonds  += b
        all_angles += a
        all_dihs   += d
        all_imps   += i

        seg_charges = sn_sys.get_initial_charges()
        for pos, orig_idx in enumerate(seg_indices):
            charges[orig_idx]            = seg_charges[pos]
            atom_types_ordered[orig_idx] = sn_atypes[pos]

    # ------------------------------------------------------------------
    # 3. Rebuild topology from assembled segment results
    # ------------------------------------------------------------------
    specie._setup_topology(atom_types_ordered,
                           all_bonds, all_angles, all_dihs, all_imps)
    specie.atoms.set_initial_charges(charges)

    # ------------------------------------------------------------------
    # 4. Refine each junction with a local snippet run
    # ------------------------------------------------------------------
    for (ci, cj) in cut_edges:
        # One snippet per cut, centred on ci (it is bonded to cj so the
        # snippet naturally spans both sides of the cut)
        center = ci
        logger.info("  >> refining junction at atoms %d -- %d...", ci, cj)

        snippet, snippet_idxs = make_snippet(specie, center, Nmax,
                                             ending=ending)

        ldxs    = list(set(np.concatenate(
            specie.find_relevant_distances(4, centers=center))))
        mapping = [int(np.argwhere(snippet_idxs == ll)[0][0]) for ll in ldxs]

        if "nominal_charge" not in snippet.arrays:
            raise ValueError(
                f"nominal_charge missing in junction snippet at atom {center}.")
        sn_charge = int(snippet.arrays["nominal_charge"].sum())

        sn_sys, sn_atypes, sn_bonds, sn_angles, sn_dihs, sn_imps = \
            run_ligpargen(snippet, charge=sn_charge, is_snippet=True)

        original_idxs = specie._sids[snippet_idxs]
        local_idxs    = specie._sids[ldxs]

        new_bonds, new_angles, new_dihs, new_imps = remap_snippet_topology(
            original_idxs, sn_sys, sn_atypes,
            sn_bonds, sn_angles, sn_dihs, sn_imps,
            local_idxs)

        new_charges    = sn_sys.get_initial_charges()
        charges[ldxs]  = new_charges[mapping]

        # Correct OPLS types for the two cut-bond atoms (they had H caps in
        # their respective segments so their types may be wrong).
        specie._update_junction_lj_types([ci, cj], snippet_idxs, sn_atypes)

        specie._add_to_topology(bonds=new_bonds, angles=new_angles,
                                dihedrals=new_dihs, impropers=new_imps)

    # ------------------------------------------------------------------
    # 5. Final charge assignment (+ optional offset correction)
    # ------------------------------------------------------------------
    if offset:
        target   = int(specie.atoms.arrays["nominal_charge"].sum())
        charges -= (charges.sum() - target) / natoms

    specie.atoms.set_initial_charges(charges)
    logger.info("  >> done.  Total charge: %.4f  (target %d)",
                charges.sum(),
                int(specie.atoms.arrays["nominal_charge"].sum()))

RESP charges (PySCF)

mdinterface.externals.pyscf.calculate_RESP_charges(specie, basis='def2-svpd', xc='b3lyp', calc_type='RKS', gpu=True, optimize=False, maxit=250)

Source code in mdinterface/externals/pyscf.py

def calculate_RESP_charges(specie, basis='def2-svpd', xc="b3lyp", calc_type="RKS",
                           gpu=True, optimize=False, maxit=250):

    try:
        from gpu4pyscf.pop import esp
        from pymbxas.build.structure import ase_to_mole, mole_to_ase
        from pymbxas.build.input_pyscf import make_pyscf_calculator
        from pymbxas.md.solvers import Geometry_optimizer
    except ImportError:
        raise ImportError("You need to install gpu4pyscf AND pymbxas to run this.")

    logger.info("RESP charges: %d atoms,  basis=%s,  xc=%s,  optimize=%s",
                len(specie.atoms), basis, xc, optimize)

    # make pyscf mol
    mol = ase_to_mole(specie.atoms, basis=basis)

    # generate calculator
    mf = make_pyscf_calculator(mol, xc=xc, calc_type=calc_type, gpu=gpu)

    atoms = specie._atoms
    if optimize:
        logger.info("  >> geometry optimization (max 100 steps)")
        gopt = Geometry_optimizer(mf)
        gopt.optimize(100)
        mol = gopt.mol_eq
        atoms = mole_to_ase(mol)
        logger.info("  >> geometry optimized")

    mol.cart = True    # PySCF uses spherical basis by default

    # run calculation and calculate density matrix
    mf = make_pyscf_calculator(mol, xc=xc, calc_type=calc_type, gpu=gpu)
    logger.debug("  >> SCF kernel")
    mf.kernel()
    dm = mf.make_rdm1()

    # RESP charge // first stage fitting
    logger.debug("  >> RESP stage 1")
    q1 = esp.resp_solve(mol, dm, maxit=maxit)

    # Add constraint: fix those charges in the second stage
    sum_constraints   = []
    equal_constraints = []

    _, idxs = find_equivalent_atoms(specie.graph)
    for idx in set(idxs):
        tmp_idx = np.concatenate(np.argwhere(idxs == idx))

        if len(tmp_idx) == 1:
            sum_constraints.append([q1[tmp_idx[0]], tmp_idx[0]])
        else:
            equal_constraints.append(tmp_idx.tolist())

    # RESP charge // second stage fitting
    logger.debug("  >> RESP stage 2 (%d equal constraints)", len(equal_constraints))
    q2 = esp.resp_solve(mol, dm, resp_a=5e-4, resp_b=0.1, tol=1e-7,
                        sum_constraints=sum_constraints,
                        equal_constraints=equal_constraints, maxit=maxit)

    logger.info("  >> done: sum=%.4f,  min=%.4f,  max=%.4f", q2.sum(), q2.min(), q2.max())
    return q2, atoms

Structure relaxation (ASE)

mdinterface.externals.optimization.relax_structure(atoms, optimizer='FIRE', fmax=0.05, steps=200, trajectory=None, logfile=None, **kwargs)

Perform structure relaxation using ASE optimizers.

Parameters:

Name	Type	Description	Default
atoms	Atoms	The atomic structure to relax	required
optimizer	str	Optimizer to use ('BFGS', 'LBFGS', 'FIRE')	'BFGS'
fmax	float	Maximum force threshold for convergence (eV/Å)	0.05
steps	int	Maximum number of optimization steps	200
trajectory	str	Path to save optimization trajectory	None
logfile	str	Path to save optimization log	None
**kwargs		Additional arguments passed to the optimizer	{}

Returns:

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

Source code in mdinterface/externals/optimization.py

def relax_structure(atoms, optimizer='FIRE', fmax=0.05, steps=200,
                       trajectory=None, logfile=None, **kwargs):
    """
    Perform structure relaxation using ASE optimizers.

    Parameters
    ----------
    atoms : ase.Atoms
        The atomic structure to relax
    optimizer : str, default 'BFGS'
        Optimizer to use ('BFGS', 'LBFGS', 'FIRE')
    fmax : float, default 0.05
        Maximum force threshold for convergence (eV/Å)
    steps : int, default 200
        Maximum number of optimization steps
    trajectory : str, optional
        Path to save optimization trajectory
    logfile : str, optional
        Path to save optimization log
    **kwargs
        Additional arguments passed to the optimizer

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

    # Select optimizer
    optimizer_dict = {
        'BFGS': BFGS,
        'LBFGS': LBFGS,
        'FIRE': FIRE
    }

    if optimizer not in optimizer_dict:
        raise ValueError(f"Unknown optimizer: {optimizer}. Available: {list(optimizer_dict.keys())}")

    # Initialize optimizer
    opt_class = optimizer_dict[optimizer]

    logger.info("Relaxing: %d atoms,  optimizer=%s,  fmax=%.3f eV/Å,  max steps=%d",
                len(atoms), optimizer, fmax, steps)

    dyn = opt_class(atoms, trajectory=trajectory, logfile=logfile, **kwargs)
    converged = dyn.run(fmax, steps)

    if converged:
        logger.info("  >> converged in %d steps", dyn.get_number_of_steps())
    else:
        logger.warning("  >> NOT converged after %d steps", steps)

    return atoms

AIMD (FAIRChem)

mdinterface.externals.aimd.run_aimd(atoms, timestep=0.5, temperature_K=300, friction=0.1, steps=1000, trajectory=None, logfile=None, **kwargs)

Perform Ab Initio Molecular Dynamics (AIMD) using ASE and FAIRChem.

Parameters:

Name	Type	Description	Default
atoms	Atoms	The atomic structure to run the AIMD on.	required
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
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
Atoms	The atomic structure after AIMD simulation.

Raises:

Type	Description
ImportError	If fairchem is not installed.
ValueError	If there's an error loading FAIRChem models.

Source code in mdinterface/externals/aimd.py

def run_aimd(atoms, timestep=0.5, temperature_K=300, friction=0.1, steps=1000,
             trajectory=None, logfile=None, **kwargs):
    """
    Perform Ab Initio Molecular Dynamics (AIMD) using ASE and FAIRChem.

    Parameters
    ----------
    atoms : ase.Atoms
        The atomic structure to run the AIMD on.
    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.
    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
    -------
    ase.Atoms
        The atomic structure after AIMD simulation.

    Raises
    ------
    ImportError
        If fairchem is not installed.
    ValueError
        If there's an error loading FAIRChem models.
    """

    if not FAIRCHEM_AVAILABLE:
        raise ImportError(
            "FAIRChem is required for AIMD simulations but is not installed. "
            "Install it with: pip install fairchem-core"
        )

    logger.info("AIMD: %d atoms,  T=%d K,  dt=%.2f fs,  %d steps  (%.3f ps)",
                len(atoms), temperature_K, timestep, steps, timestep * steps / 1000)

    # Make a copy to avoid modifying the original
    aimd_atoms = atoms.copy()

    # Ensure the computation setup with a fairchem calculator
    try:
        predictor = pretrained_mlip.get_predict_unit("uma-s-1p1", device="cuda")
        calc = FAIRChemCalculator(predictor, task_name="omol")
        logger.debug("  >> UMA calculator loaded")
    except Exception as e:
        raise ValueError(f"Error loading UMA or OMol models: {e}")

    aimd_atoms.calc = calc

    # start velocities
    MaxwellBoltzmannDistribution(aimd_atoms, temperature_K=temperature_K)
    logger.debug("  >> Maxwell-Boltzmann velocities initialized")

    # Initialize Langevin dynamics
    dyn = Langevin(
            aimd_atoms,
            timestep=timestep * units.fs,
            temperature_K=temperature_K,
            friction=friction / units.fs,
            **kwargs
            )

    # define traj and logger
    if trajectory:
        npt_traj = Trajectory(trajectory, mode="w", atoms=aimd_atoms)
        dyn.attach(npt_traj.write, interval=1)
        logger.debug("  >> trajectory -> %s", trajectory)
    if logfile:
        npt_log  = MDLogger(dyn, aimd_atoms, logfile, stress=False, peratom=False, mode="w")
        dyn.attach(npt_log, interval=1)
        logger.debug("  >> log        -> %s", logfile)

    # Run the dynamics
    dyn.run(steps)
    logger.info("  >> done: %.3f ps simulated", timestep * steps / 1000)

    return aimd_atoms