Externals

Optional integrations with third-party tools.

LigParGen (OPLS-AA parameters)

mdinterface.externals.ligpargen.LigParGenError

Bases: RuntimeError

LigParGen configuration, execution, or output-processing failure.

Source code in mdinterface/externals/ligpargen.py

class LigParGenError(RuntimeError):
    """LigParGen configuration, execution, or output-processing failure."""

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

mdinterface.externals.ligpargen.run_ligpargen(system, charge=None, is_snippet=False)

Generate OPLS-AA parameters by running LigParGen.

Parameters:

Name	Type	Description	Default
system	Atoms	Atomic system to parameterize. Stored RDKit chemistry is transferred through MOL files; coordinate-only inputs use XYZ and Open Babel.	required
charge	int or None	Total molecular charge. LigParGen detects it when omitted.	None
is_snippet	bool	Whether the system is a capped molecular snippet.	False

Returns:

Type	Description
tuple	Parameterized system, atom types, bonds, angles, dihedrals, and impropers.

Raises:

Type	Description
LigParGenError	If LigParGen or BOSS is not configured, execution fails, or the output is missing or unreadable.

Source code in mdinterface/externals/ligpargen.py

def run_ligpargen(system, charge=None, is_snippet=False):
    """Generate OPLS-AA parameters by running LigParGen.

    Parameters
    ----------
    system : ase.Atoms
        Atomic system to parameterize. Stored RDKit chemistry is transferred
        through MOL files; coordinate-only inputs use XYZ and Open Babel.
    charge : int or None, default None
        Total molecular charge. LigParGen detects it when omitted.
    is_snippet : bool, default False
        Whether the system is a capped molecular snippet.

    Returns
    -------
    tuple
        Parameterized system, atom types, bonds, angles, dihedrals, and impropers.

    Raises
    ------
    LigParGenError
        If LigParGen or BOSS is not configured, execution fails, or the output
        is missing or unreadable.
    """

    if len(system) > 200:
        raise ValueError(f"LigParGen accepts at most 200 atoms including caps; received {len(system)}.")

    ligpargen_executable = shutil.which("ligpargen")
    if ligpargen_executable is None:
        raise LigParGenError(
            "LigParGen executable was not found on PATH. Install it in the active "
            "environment with `python -m pip install "
            "\"git+https://github.com/roncofaber/ligpargen.git@ad78036842318f166531be41cfcbc3563d7c5476\"` and verify "
            "the installation with `ligpargen -h`"
        )

    from rdkit import Chem
    from mdinterface.core.chemistry import stored_molecule

    mol = stored_molecule(system)
    if mol is not None:
        formal_charge = Chem.GetFormalCharge(mol)
        if charge is not None and charge != formal_charge:
            raise ValueError("LigParGen charge conflicts with the molecular formal charge.")
        charge = formal_charge
    if mol is None and shutil.which("obabel") is None:
        raise LigParGenError(
            "Open Babel executable `obabel` was not found on PATH. LigParGen "
            "requires it to read mdinterface's XYZ input. Install it with "
            "`conda install -c conda-forge openbabel` and verify the installation "
            "with `obabel -V`"
        )

    if not os.environ.get("BOSSdir"):
        from mdinterface.config import load_config

        load_config()

    if not os.environ.get("BOSSdir"):
        mdint = os.environ.get("MDINT_CONFIG_DIR", "~/.config/mdinterface")
        config_file = os.path.join(os.path.expanduser(mdint), "config.ini")
        raise LigParGenError(
            "BOSSdir is not configured. Export BOSSdir or set it under [settings] "
            f"in {config_file}"
        )

    # all ligpargen files go in a temp dir; kept on failure for inspection
    tmpdir   = tempfile.mkdtemp(prefix="ligpargen_")
    mol_name = os.path.basename(tmpdir)
    extension = "mol" if mol is not None else "xyz"
    input_file = os.path.join(tmpdir, f"{mol_name}.{extension}")
    log_file = os.path.join(tmpdir, "ligpargen.log")

    if mol is None:
        ase.io.write(input_file, system)
    else:
        Chem.MolToMolFile(mol, input_file)

    # use relative filenames and cwd=tmpdir -- ligpargen does not accept
    # absolute paths for -i
    ligpargen_command = [ligpargen_executable, "-i", f"{mol_name}.{extension}", "-p", tmpdir,
                         "-debug", "-o", "0", "-cgen", "CM1A"]
    if charge is not None:
        ligpargen_command.extend(["-c", str(charge)])

    try:
        result = subprocess.run(ligpargen_command, check=True,
                                stdout=subprocess.PIPE, stderr=subprocess.PIPE,
                                cwd=tmpdir, text=True, encoding="utf-8",
                                errors="replace")
        with open(log_file, "w") as fh:
            fh.write("STDOUT:\n" + result.stdout + "\n")
            fh.write("STDERR:\n" + result.stderr + "\n")
        logger.debug("ligpargen completed successfully")
        logger.debug("ligpargen stdout:\n%s", result.stdout)

    except subprocess.CalledProcessError as e:
        with open(log_file, "w") as fh:
            fh.write("STDOUT:\n" + (e.stdout or "") + "\n")
            fh.write("STDERR:\n" + (e.stderr or "") + "\n")
        logger.error("ligpargen failed; temp files kept at: %s", tmpdir)
        logger.debug("ligpargen stderr:\n%s", e.stderr)
        raise LigParGenError(
            "LigParGen exited unsuccessfully",
            tmpdir,
            log_file,
            returncode=e.returncode,
        ) from e
    except OSError as e:
        raise LigParGenError(
            f"LigParGen could not be started: {e}",
            tmpdir,
            log_file,
        ) from e

    output_file = os.path.join(tmpdir, f"{mol_name}.lammps.lmp")
    if not os.path.isfile(output_file):
        diagnostic = (result.stderr or result.stdout).strip()
        message = f"LigParGen did not create the expected output file {output_file}"
        if diagnostic:
            message += f". LigParGen output: {diagnostic[-500:]}"
        raise LigParGenError(
            message,
            tmpdir,
            log_file,
            returncode=result.returncode,
        )

    try:
        original_numbers = system.numbers.copy()
        system, atoms, bonds, angles, dihedrals, impropers = read_lammps_data_file(
            output_file, is_snippet=is_snippet)
        if not np.array_equal(original_numbers, system.numbers):
            raise ValueError("LigParGen output atom ordering differs from the input.")
    except Exception as e:
        raise LigParGenError(
            f"LigParGen output could not be read from {output_file}",
            tmpdir,
            log_file,
            returncode=result.returncode,
        ) from e

    # success -- clean up
    shutil.rmtree(tmpdir, ignore_errors=True)

    return system, atoms, bonds, angles, dihedrals, impropers

mdinterface.externals.ligpargen.refine_large_specie_topology(specie, snippet_radius=12, cap_element='H', charge_correction='none', segment_size=200)

Assign LigParGen parameters atomically; see Specie.parameterize.

Parameters:

Name	Type	Description	Default
specie	Specie	Species to parameterize without changing its coordinates.	required
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 molecular charge.	"none"
segment_size	int	Maximum segment size including caps, at most 200.	200

Returns:

Type	Description
dict	Charge audit. No changes are applied if any calculation fails.

Source code in mdinterface/externals/ligpargen.py

def refine_large_specie_topology(specie, snippet_radius=12, cap_element="H",
                                 charge_correction="none", segment_size=200):
    """Assign LigParGen parameters atomically; see ``Specie.parameterize``.

    Parameters
    ----------
    specie : Specie
        Species to parameterize without changing its coordinates.
    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 molecular charge.
    segment_size : int, default 200
        Maximum segment size including caps, at most 200.

    Returns
    -------
    dict
        Charge audit. No changes are applied if any calculation fails.
    """
    if not isinstance(segment_size, (int, np.integer)) or not 4 <= segment_size <= 200:
        raise ValueError("segment_size must be an integer between 4 and 200, including caps.")
    if not isinstance(snippet_radius, (int, np.integer)) or snippet_radius < 4:
        raise ValueError("snippet_radius must be an integer of at least 4.")
    if charge_correction not in {"none", "uniform"}:
        raise ValueError("charge_correction must be 'none' or 'uniform'.")
    if cap_element not in {"H", "F", "Cl", "Br", "I"}:
        raise ValueError("cap_element must be a neutral monovalent element.")
    staged, attributes = specie._parameterization_copy()
    initial = _clean_charge(staged.charges.sum())
    target = staged._resolve_charge(None)
    if len(staged.atoms) <= segment_size:
        result, atom_types, bonds, angles, dihedrals, impropers = run_ligpargen(staged.atoms, charge=target)
        if not np.array_equal(result.numbers, staged.atoms.numbers) or len(atom_types) != len(result):
            raise ValueError("Parameterization changed atom ordering or omitted atom types.")
        staged._setup_topology(atom_types, bonds, angles, dihedrals, impropers)
        staged.atoms.set_initial_charges(result.get_initial_charges())
        junctions = 0
    else:
        junctions = _refine_large_specie_topology(staged, snippet_radius, cap_element, segment_size)
    charges = staged.charges
    if not np.isfinite(charges).all():
        raise ValueError("Parameterization returned nonfinite partial charges.")
    refined = _clean_charge(charges.sum())
    residual = _clean_charge(refined - target)
    correction = -residual / len(charges) if charge_correction == "uniform" and residual != 0 else 0.0
    staged.atoms.set_initial_charges(charges + correction)
    report = dict(formal_charge=target, initial_charge=initial, refined_charge=refined,
                  residual=residual, correction_per_atom=correction,
                  final_charge=_clean_charge(staged.charges.sum()), junctions=junctions)
    staged.validate_force_field()
    specie._apply_parameterization(staged, attributes)
    logger.info("Parameterization charge audit: %s", report)
    return report

RESP charges (PySCF)

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

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, charge=None):

    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, charge=specie._resolve_charge(charge))

    # 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