prepare

Preparing systems ready for simulation.

Functions:

• load_ligand –

  Load a ligand from disk.

• load_ligands –

  Load the first, and optionally second, ligand from disk.

• load_receptor –

  Loads a receptor from disk.

• apply_hmr –

  Apply hydrogen mass repartitioning to a system.

• prepare_system –

  Solvates and parameterizes a system.

load_ligand

load_ligand(
    path: Path, residue_name: str = "LIG"
) -> Topology

Load a ligand from disk.

Parameters:

• path (Path) –

  The path to the ligand file (.sdf, .mol2)

• residue_name (str, default: 'LIG' ) –

  The residue name to assign to the ligand.

Returns:

• Topology –

  The loaded ligand

Source code in femto/md/prepare.py

def load_ligand(path: pathlib.Path, residue_name: str = "LIG") -> mdtop.Topology:
    """Load a ligand from disk.

    Args:
        path: The path to the ligand file (.sdf, .mol2)
        residue_name: The residue name to assign to the ligand.

    Returns:
        The loaded ligand
    """
    residue_name = "LIG" if residue_name is None else residue_name

    if path.suffix.lower() == ".sdf":
        mol = Chem.AddHs(Chem.MolFromMolFile(str(path), removeHs=False))
    elif path.suffix.lower() == ".mol2":
        mol = Chem.AddHs(Chem.MolFromMol2File(str(path), removeHs=False))
    else:
        raise NotImplementedError(f"unsupported file format: {path.suffix}")

    top = mdtop.Topology.from_rdkit(mol, residue_name)

    return top

load_ligands

load_ligands(
    ligand_1_path: Path, ligand_2_path: Path | None
) -> tuple[Topology, Topology | None]

Load the first, and optionally second, ligand from disk.

Parameters:

• ligand_1_path (Path) –

  The path to the first ligand.

• ligand_2_path (Path | None) –

  The (optional) path of the second ligand.

Returns:

• tuple[Topology, Topology | None] –

  The loaded ligands.

Source code in femto/md/prepare.py

def load_ligands(
    ligand_1_path: pathlib.Path,
    ligand_2_path: pathlib.Path | None,
) -> tuple[mdtop.Topology, mdtop.Topology | None]:
    """Load the first, and optionally second, ligand from disk.

    Args:
        ligand_1_path: The path to the first ligand.
        ligand_2_path: The (optional) path of the second ligand.

    Returns:
        The loaded ligands.
    """

    ligand_1 = load_ligand(ligand_1_path, femto.md.constants.LIGAND_1_RESIDUE_NAME)

    if ligand_2_path is None:
        return ligand_1, None

    ligand_2 = load_ligand(ligand_2_path, femto.md.constants.LIGAND_2_RESIDUE_NAME)

    return ligand_1, ligand_2

load_receptor

load_receptor(path: Path) -> Topology

Loads a receptor from disk.

Parameters:

• path (Path) –

  The path to the receptor (.pdb, .mol2, .sdf).

Returns:

• Topology –

  The loaded receptor.

Source code in femto/md/prepare.py

def load_receptor(path: pathlib.Path) -> mdtop.Topology:
    """Loads a receptor from disk.

    Args:
        path: The path to the receptor (.pdb, .mol2, .sdf).

    Returns:
        The loaded receptor.
    """
    if path.suffix.lower() == ".pdb":
        return mdtop.Topology.from_file(path)
    elif path.suffix.lower() in {".mol2", ".sdf"}:
        return femto.md.prepare.load_ligand(path, "REC")

    raise NotImplementedError(f"unsupported file format: {path.suffix}")

apply_hmr

apply_hmr(
    system: System,
    topology: Topology,
    hydrogen_mass: Quantity = 1.5 * amu,
)

Apply hydrogen mass repartitioning to a system.

Parameters:

• system (System) –

  The system to modify in-place.

• topology (Topology) –

  The topology of the system.

• hydrogen_mass (Quantity, default: 1.5 * amu ) –

  The mass to use ofr hydrogen atoms.

Source code in femto/md/prepare.py

def apply_hmr(
    system: openmm.System,
    topology: mdtop.Topology,
    hydrogen_mass: openmm.unit.Quantity = 1.5 * openmm.unit.amu,
):
    """Apply hydrogen mass repartitioning to a system.

    Args:
        system: The system to modify in-place.
        topology: The topology of the system.
        hydrogen_mass: The mass to use ofr hydrogen atoms.
    """

    atoms = topology.atoms

    for bond in topology.bonds:
        atom_1: mdtop.Atom = atoms[bond.idx_1]
        atom_2: mdtop.Atom = atoms[bond.idx_2]

        if atom_1.atomic_num == 1:
            (atom_1, atom_2) = (atom_2, atom_1)

        if atom_2.atomic_num != 1:
            continue
        if atom_1.atomic_num == 1:
            continue

        elements = sorted(a.atomic_num for a in atom_2.residue.atoms)

        if elements == [1, 1, 8]:
            continue

        mass_delta = hydrogen_mass - system.getParticleMass(atom_2.index)

        system.setParticleMass(atom_2.index, hydrogen_mass)
        system.setParticleMass(
            atom_1.index, system.getParticleMass(atom_1.index) - mass_delta
        )

prepare_system

prepare_system(
    receptor: Topology | None,
    ligand_1: Topology | None,
    ligand_2: Topology | None,
    cofactors: list[Topology] | None = None,
    config: Prepare | None = None,
    ligand_1_offset: Quantity | None = None,
    ligand_2_offset: Quantity | None = None,
    cavity_formers: list[Topology] | None = None,
    extra_params: list[Path] | None = None,
) -> tuple[Topology, System]

Solvates and parameterizes a system.

Parameters:

• receptor (Topology | None) –

  A receptor to include in the system.

• ligand_1 (Topology | None) –

  A primary ligand to include in the system.

• ligand_2 (Topology | None) –

  A secondary ligand to include in the system if doing RBFE.

• cofactors (list[Topology] | None, default: None ) –

  Any cofactors to include in the system.

• config (Prepare | None, default: None ) –

  Configuration options for the preparation.

• ligand_1_offset (Quantity | None, default: None ) –

  The amount to offset the first ligand by before computing the box size if using a padded box.

• ligand_2_offset (Quantity | None, default: None ) –

  The amount to offset the second ligand by before computing the box size if using a padded box.

• cavity_formers (list[Topology] | None, default: None ) –

  A (optional) list of topologies that should be considered 'present' when placing the solvent molecules such that they leave cavities, but are not added to the final topology themselves.

  Note that cavity formers will be considered when determining the box size.

• extra_params (list[Path] | None, default: None ) –

  The paths to any extra parameter files (.xml, .parm) to use when parameterizing the system.

Returns:

• tuple[Topology, System] –

  The solvated and parameterized topology and system, containing the ligands, the receptor, cofactors, ions and solvent.

Source code in femto/md/prepare.py

def prepare_system(
    receptor: mdtop.Topology | None,
    ligand_1: mdtop.Topology | None,
    ligand_2: mdtop.Topology | None,
    cofactors: list[mdtop.Topology] | None = None,
    config: femto.md.config.Prepare | None = None,
    ligand_1_offset: openmm.unit.Quantity | None = None,
    ligand_2_offset: openmm.unit.Quantity | None = None,
    cavity_formers: list[mdtop.Topology] | None = None,
    extra_params: list[pathlib.Path] | None = None,
) -> tuple[mdtop.Topology, openmm.System]:
    """Solvates and parameterizes a system.

    Args:
        receptor: A receptor to include in the system.
        ligand_1: A primary ligand to include in the system.
        ligand_2: A secondary ligand to include in the system if doing RBFE.
        cofactors: Any cofactors to include in the system.
        config: Configuration options for the preparation.
        ligand_1_offset: The amount to offset the first ligand by before computing the
            box size if using a padded box.
        ligand_2_offset: The amount to offset the second ligand by before computing the
            box size if using a padded box.
        cavity_formers: A (optional) list of topologies that should be considered
            'present' when placing the solvent molecules such that they leave cavities,
            but are not added to the final topology themselves.

            Note that cavity formers will be considered when determining the box size.
        extra_params: The paths to any extra parameter files (.xml, .parm) to use
            when parameterizing the system.

    Returns:
        The solvated and parameterized topology and system, containing the ligands, the
        receptor, cofactors, ions and solvent.
    """
    cofactors = [] if cofactors is None else cofactors
    cavity_formers = [] if cavity_formers is None else copy.deepcopy(cavity_formers)

    config = config if config is not None else femto.md.config.Prepare()

    force_field = _load_force_field(
        *config.default_protein_ff, *([] if extra_params is None else extra_params)
    )

    if config.default_ligand_ff is not None:
        _register_openff_generator(ligand_1, ligand_2, cofactors, force_field, config)

    topology = mdtop.Topology.merge(
        *([] if ligand_1 is None else [ligand_1]),
        *([] if ligand_2 is None else [ligand_2]),
        *([] if receptor is None else [receptor]),
        *cofactors,
    )

    box_size = None

    if config.box_padding is not None:
        box_size = _compute_box_size(
            receptor,
            ligand_1,
            ligand_2,
            cofactors,
            config.box_padding,
            ligand_1_offset,
            ligand_2_offset,
            cavity_formers,
        )

        if config.box_shape.lower() == "cube":
            box_size = (
                numpy.array([max(box_size.value_in_unit(openmm.unit.angstrom))] * 3)
                * openmm.unit.angstrom
            )

        _LOGGER.info(f"using a box size of {box_size}")

    for former in cavity_formers:
        for residue in former.residues:
            residue.name = "CAV"

    cavity = mdtop.Topology.merge(topology, *cavity_formers)

    _LOGGER.info("adding solvent and ions")
    modeller = openmm.app.Modeller(cavity.to_openmm(), cavity.xyz)
    modeller.addExtraParticles(force_field)
    modeller.addSolvent(
        force_field,
        model=config.water_model.lower(),
        boxSize=box_size,
        numAdded=None if box_size is not None else config.n_waters,
        boxShape="cube",
        positiveIon=config.cation,
        negativeIon=config.anion,
        neutralize=config.neutralize,
        ionicStrength=config.ionic_strength,
    )

    topology = mdtop.Topology.from_openmm(modeller.topology)
    topology.xyz = (
        numpy.array(modeller.positions.value_in_unit(openmm.unit.angstrom))
        * openmm.unit.angstrom
    )
    topology = topology["not r. CAV"]  # strip cavity formers

    _LOGGER.info("parameterizing the system")

    system = force_field.createSystem(
        topology.to_openmm(),
        nonbondedMethod=openmm.app.PME,
        nonbondedCutoff=0.9 * openmm.unit.nanometer,
        constraints=openmm.app.HBonds,
        rigidWater=True,
    )

    # TODO: is this still needed??
    bound = mdtop.Topology.merge(
        *([] if ligand_1 is None else [ligand_1]),
        *([] if ligand_2 is None else [ligand_2]),
        *([] if receptor is None else [receptor]),
        *cofactors,
    )

    center_offset = (
        bound.xyz.value_in_unit(openmm.unit.angstrom).mean(axis=0)
        * openmm.unit.angstrom
    )
    box_center = (
        numpy.diag(topology.box.value_in_unit(openmm.unit.angstrom)) * 0.5
    ) * openmm.unit.angstrom

    topology.xyz = topology.xyz - center_offset + box_center

    return topology, system