solvate

Solvate protein ligand systems

Modules:

• amber_utils –

  Utilities for manipulating AMBER data and tools

• femto –

  A comprehensive toolkit for predicting free energies

Functions:

• solvate_system –

  Solvates a system.

solvate_system

solvate_system(
    receptor: Structure | None,
    ligand_1: Structure,
    ligand_2: Structure | None,
    solvent: Solvent,
    ligand_1_offset: Quantity | None = None,
    ligand_2_offset: Quantity | None = None,
    cavity_formers: list[Structure] | None = None,
) -> Structure

Solvates a system.

Parameters:

• receptor (Structure | None) –

  The (optional) receptor.

• ligand_1 (Structure) –

  The first ligand.

• ligand_2 (Structure | None) –

  The (optional) second ligand.

• solvent (Solvent) –

  The solvent configuration.

• 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[Structure] | None, default: None ) –

  The (optional) list of structures 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.

Returns:

• Structure –

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

Source code in femto/md/solvate.py

def solvate_system(
    receptor: parmed.Structure | None,
    ligand_1: parmed.Structure,
    ligand_2: parmed.Structure | None,
    solvent: femto.md.config.Solvent,
    ligand_1_offset: openmm.unit.Quantity | None = None,
    ligand_2_offset: openmm.unit.Quantity | None = None,
    cavity_formers: list[parmed.Structure] | None = None,
) -> parmed.Structure:
    """Solvates a system.

    Args:
        receptor: The (optional) receptor.
        ligand_1: The first ligand.
        ligand_2: The (optional) second ligand.
        solvent: The solvent configuration.
        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: The (optional) list of structures 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.

    Returns:
        The solvated system containing the ligands, the receptor, ions and the solvent.
    """

    bound = ligand_1

    if ligand_2 is not None:
        bound = bound + ligand_2
    if receptor is not None:
        bound = bound + receptor

    cavity = bound

    if cavity_formers is not None:
        for former in cavity_formers:
            cavity = cavity + former

    modeller = openmm.app.Modeller(cavity.topology, cavity.positions)

    box_size = None

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

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

    modeller.addSolvent(
        _MockForceField(cavity),
        model=solvent.water_model.lower(),
        boxSize=box_size,
        numAdded=None if box_size is not None else solvent.n_waters,
        boxShape="cube",
        positiveIon=solvent.cation,
        negativeIon=solvent.anion,
        neutralize=solvent.neutralize,
        ionicStrength=solvent.ionic_strength,
    )

    solution: parmed.Structure = parmed.openmm.load_topology(
        modeller.topology, None, modeller.positions
    )
    solution_box = solution.box

    _LOGGER.info("parameterizing solvent")

    solvent_mask = amber_utils.extract_water_and_ions_mask(bound)
    bound.strip(solvent_mask)

    complex_mask = [1 - i for i in amber_utils.extract_water_and_ions_mask(solution)]
    solution.strip(complex_mask)

    solution = amber_utils.parameterize_structure(solution, solvent.tleap_sources)

    system = bound + solution
    system.box = solution_box

    center_offset = numpy.array(bound.coordinates).mean(axis=0)

    coordinates = numpy.array(system.coordinates)
    coordinates -= center_offset
    coordinates += system.box[:3] * 0.5

    system.coordinates = coordinates.tolist()

    return system