<jupyter_start><jupyter_text>Replacing ligand parameters in an already-parametrized system using `openmm-forcefields`This example applies SMIRNOFF-format parameters to a BRD4 inhibitor from the [living review on binding free energy benchmark systems](https://www.annualreviews.org/doi/abs/10.1146/annurev-biophys-070816-033654) by <NAME> Gilson. The BRD4 system comes from the [accompanying GitHub repository](https://github.com/MobleyLab/benchmarksets/tree/master/input_files/BRD4).This example uses the [`openmm-forcefields` package](http://github.com/openmm/openmm-forcefields) to add a [residue template generator](http://docs.openmm.org/latest/userguide/application.htmladding-residue-template-generators) to the [`openmm.app.ForceField`](http://docs.openmm.org/latest/api-python/generated/openmm.app.forcefield.ForceField.htmlopenmm.app.forcefield.ForceField) to allow Open Force Field small molecule parameters to be generated on the fly when parameterizing a system containing protein, small molecules, ions, and water. This example is meant to illustrate how to apply parameters to a single ligand, but it's also easy to process many ligands. Loading the already-parametrized system<jupyter_code># Retrieve protein and ligand files for BRD4 and a docked inhibitor from the benchmark systems GitHub repository
# https://github.com/MobleyLab/benchmarksets
import requests

repo_url = (
    "https://raw.githubusercontent.com/MobleyLab/benchmarksets/master/input_files/"
)
sources = {
    "system.pdb": repo_url
    + "BRD4/prmtop-coords/BRD4-1.pdb",  # complete system (protein+ligand+solvent+ions)
    "ligand.sdf": repo_url + "BRD4/sdf/ligand-1.sdf",  # ligand molecular identity
}
for (filename, url) in sources.items():
    r = requests.get(url)
    open(filename, "w").write(r.text)
try:
    from openmm.app import PDBFile
except ImportError:
    from simtk.openmm.app import PDBFile

# Read complete system in OpenMM PDBFile
system_pdb = "system.pdb"
pdbfile = PDBFile(system_pdb)

# We have to remove H1-H2 bonds in waters if they are present
# AMBER's 'ambpdb -conect' adds these H1-H2 bonds, so we must remove them
def fix_water_bonds(topology):
    # TODO: We should create a simpler way to do this within OpenMM's Topology object
    n_bonds_before = sum(1 for bond in topology.bonds())
    try:
        from openmm.app.element import hydrogen
    except ImportError:
        from simtk.openmm.app.element import hydrogen

    bonds_to_delete = [
        index
        for (index, bond) in enumerate(topology.bonds())
        if ((bond.atom1.element == hydrogen) and (bond.atom2.element == hydrogen))
    ]
    bonds_to_delete.reverse()
    for index in bonds_to_delete:
        topology._bonds.pop(index)
    n_bonds_after = sum(1 for bond in topology.bonds())
    print(f"{n_bonds_before - n_bonds_after} H-H bonds removed")


fix_water_bonds(pdbfile.topology)
import numpy as np

try:
    from openmm import unit
except ImportError:
    from simtk import unit

# Load the definition of the small molecule in the system from an SDF file
from openff.toolkit.topology import Molecule

ligand = Molecule.from_file("ligand.sdf")
# Patch until openmmforcefield 0.7.3 is released
if ligand.partial_charges is None:
    ligand.partial_charges = np.zeros([ligand.n_particles]) * unit.elementary_charge
# Create an OpenMM ForceField object with AMBER ff14SB and TIP3P with compatible ions
try:
    from openmm import app
except ImportError:
    from simtk.openmm import app

forcefield = app.ForceField(
    "amber/protein.ff14SB.xml",
    "amber/tip3p_standard.xml",
    "amber/tip3p_HFE_multivalent.xml",
)
# Use the SMIRNOFF residue template generator to load the openff-2.0.0 ("Sage") that knows about the ligand
from openmmforcefields.generators import SMIRNOFFTemplateGenerator

smirnoff = SMIRNOFFTemplateGenerator(forcefield="openff-2.0.0", molecules=ligand)
# Register the SMIRNOFF template generator
forcefield.registerTemplateGenerator(smirnoff.generator)
# Create a parameterized OpenMM System from the PDB topology without bond constraints so we can convert to other packages
system = forcefield.createSystem(
    pdbfile.topology,
    nonbondedMethod=app.PME,
    constraints=None,
    rigidWater=False,
    removeCMMotion=False,
)<jupyter_output><empty_output><jupyter_text>Create a ParmEd Structure object to export to other formats<jupyter_code># Create the complex Structure
import parmed

complex_structure = parmed.openmm.load_topology(pdbfile.topology, system=system)

# Copy over the original coordinates and box vectors
complex_structure.coordinates = pdbfile.positions
complex_structure.box_vectors = pdbfile.topology.getPeriodicBoxVectors()
# Save the final PDB file to make sure the conversion worked
complex_structure.save("new-system.pdb", overwrite=True)<jupyter_output><empty_output><jupyter_text>Export to AMBER and GROMACS formatsWe started off in AMBER format, and presumably may want to continue in that format -- so let's write out to AMBER and GROMACS format:<jupyter_code># Export the Structure to AMBER files
complex_structure.save("complex-openff.prmtop", overwrite=True)
complex_structure.save("complex-openff.inpcrd", overwrite=True)

# Export the Structure to Gromacs files
complex_structure.save("complex-openff.gro", overwrite=True)
complex_structure.save("complex-openff.top", overwrite=True)<jupyter_output><empty_output><jupyter_text>That should conclude our work in this example. However, perhaps we should just doublecheck by ensuring we can actually run some dynamics on the combined system without any trouble. As a test, run some dynamics on the combined systemFirst, we create an OpenMM system, as we've done in other examples here. We can do this, in this case, using ParmEd's built-in `createSystem` functionality already attached to the combined `Structure`. We ask for a reasonable cutoff, constrained hydrogen bonds (note that **this keyword argument overrides the fact that we use the `unconstrained` force field above**; the ligand (and all other molecules in the system) **will** have covalent bonds to hydrogen constrainted), PME, and rigid water:<jupyter_code>from parmed.openmm import NetCDFReporter

try:
    from openmm import LangevinIntegrator, app
except ImportError:
    from simtk.openmm import LangevinIntegrator, app

system = complex_structure.createSystem(
    nonbondedMethod=app.PME,
    nonbondedCutoff=9 * unit.angstrom,
    constraints=app.HBonds,
    rigidWater=True,
)<jupyter_output><empty_output><jupyter_text>Next we'll set up the integrator, a reporter to write the trajectory, pick the timestep, and then go on to minimize the energy and run a very short amount of dynamics after setting the temperature to 300K:<jupyter_code>integrator = LangevinIntegrator(
    300 * unit.kelvin, 1 / unit.picosecond, 0.001 * unit.picoseconds
)
simulation = app.Simulation(complex_structure.topology, system, integrator)

# Depending on where your system came from, you may want to
# add something like (30, 30, 30)*Angstrom to center the protein
# (no functional effect, just visualizes better)
# simulation.context.setPositions(complex_structure.positions + np.array([30, 30, 30])*unit.angstrom)

simulation.context.setPositions(complex_structure.positions)

nc_reporter = NetCDFReporter("trajectory.nc", 10)
simulation.reporters.append(nc_reporter)
# Show the initial potential energy
potential_energy = simulation.context.getState(getEnergy=True).getPotentialEnergy()
print(potential_energy)
# Minimize the energy
simulation.minimizeEnergy()
minimized_coords = simulation.context.getState(getPositions=True).getPositions()
# Run some dynamics
simulation.context.setVelocitiesToTemperature(300 * unit.kelvin)
simulation.step(1000)
# Show the final potential energy
potential_energy = simulation.context.getState(getEnergy=True).getPotentialEnergy()
print(potential_energy)<jupyter_output>-500762.9999424317 kJ/mol
{'fraction_non_alphanumeric': 0.08254209393594125, 'fraction_numerical': 0.0126598303582732, 'mean_word_length': 5.36583400483481, 'pattern_counts': {'":': 2, '<': 21, '<?xml version=': 0, '>': 21, 'https://': 4, 'lorem ipsum': 0, 'www.': 1, 'xml': 3}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '508899', 'license': 'MIT', 'n_tokens_mistral': 2505, 'n_tokens_neox': 2340, 'n_words': 730, 'path': 'examples/swap_amber_parameters/swap_existing_ligand_parameters_with_openmmforcefields.ipynb', 'repo_name': 'andrew-abimansour/openff-toolkit', 'size': '7903'}