Build a parametreised protein ligand system using openmmforcefields ready for simulation.

protein_ligand_setup.py

#combine a protein and ligand to build an OpenMM system
# the molecules have been extracted from the PDB file and saved to sdf to work with openff, note they are still in the bound corrdinates.

from simtk import unit
from simtk.openmm import app
from openff.toolkit.topology import Molecule
import parmed
from openmmforcefields.generators import SystemGenerator
from simtk.openmm.app import modeller

pdbfilename = "bace_protein.pdb"
sdf_file_name = "bace_ligands_shifted.sdf"

# protein only pdb file
pdb_file = app.PDBFile(pdbfilename)
# an sdf of all ligands, returns a molecule or list of molecules depending on the number of molecules in the sdf
molecules = Molecule.from_file(sdf_file_name)

# create the paremd structure to combine the system
protein_structure = parmed.openmm.load_topology(pdb_file.topology, xyz=pdb_file.positions)

molecules_structure = parmed.load_file(sdf_file_name)
# just build the first system
complex_structure = protein_structure + molecules_structure[0]


# write out a pdb file of the whole system
with open("complex.pdb", "w") as outfile:
    app.PDBFile.writeFile(complex_structure.topology, complex_structure.positions, outfile)

# build a system generator using amber and openff
sysytem_gen = SystemGenerator(forcefields=['amber/protein.ff14SB.xml', 'amber/tip3p_standard.xml', 'amber/tip3p_HFE_multivalent.xml'], small_molecule_forcefield="openff-1.1.0", molecules=molecules)

modeller = app.Modeller(complex_structure.topology, complex_structure.positions)
modeller.addSolvent(sysytem_gen.forcefield, padding=0*unit.angstrom, ionicStrength=300*unit.millimolar)

# now make the system
system = sysytem_gen.create_system(modeller.topology)

# write out the full system
with open("solvated_system.pdb", "w") as outfile:
    app.PDBFile.writeFile(modeller.topology, modeller.positions, outfile)