# 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 = {
    'receptor.pdb' : repo_url + 'BRD4/pdb/BRD4.pdb',
    'ligand.pdb'   : repo_url + 'BRD4/pdb/ligand-1.pdb',
    'ligand.sdf'   : repo_url + 'BRD4/sdf/ligand-1.sdf',
}
for (filename, url) in sources.items():
    r = requests.get(url)
    open(filename, 'w').write(r.text)

"""
Explanation: Using SMIRNOFF with Amber on BRD4:inhibitor complexes: Exporting parameterized complexes to Amber, Gromacs, and CHARMM
This example applies SMIRNOFF-format parameters to BRD4 inhibitors from the living review on binding free energy benchmark systems by Mobley and Gilson. The BRD4 system comes from the accompanying GitHub repository.
This example uses ParmEd to combine a protein parameterized with Amber with a ligand parameterized with SMIRNOFF. This example is meant to illustrate how to apply parameters to a single ligand, but it's also easy to process many ligands.
End of explanation
"""


# Create an OpenFF Molecule object from the ligand SDF fiel
from openff.toolkit.topology import Molecule
ligand_off_molecule = Molecule('ligand.sdf')

# Load the SMIRNOFF-format Parsley force field
from openff.toolkit.typing.engines.smirnoff import ForceField
force_field = ForceField('openff_unconstrained-1.0.0.offxml')

# Parametrize the ligand molecule by creating a Topology object from it
ligand_system = force_field.create_openmm_system(ligand_off_molecule.to_topology())

"""
Explanation: Parametrize a molecule with the SMIRNOFF-format "Parsley" force field
First, we parametrize the ligand with the Parsley force field.
We use both a PDB file and an SDF file for the ligand.
End of explanation
"""


# Read in the coordinates of the ligand from the PDB file
from simtk.openmm.app import PDBFile
ligand_pdbfile = PDBFile('ligand.pdb')

# Convert OpenMM System object containing ligand parameters into a ParmEd Structure.
import parmed
ligand_structure = parmed.openmm.load_topology(ligand_pdbfile.topology,
                                                ligand_system,
                                                xyz=ligand_pdbfile.positions)

"""
Explanation: and we convert the OpenMM System to a ParmEd Structure that we'll be able to mix with the protein
<div class="alert alert-block alert-warning">
  <b>Warning:</b> ParmEd's Structure model is inspired by AMBER. Some information in an OpenMM System are not directly translatable into a Structure. In particular, long-range interaction treatment method (e.g., PME, CutoffPeriodic) and parameters (e.g., cutoff and cutoff switching distance, PME error tolerance) are known to be lost during the conversion.
</div>
End of explanation
"""


receptor_pdbfile = PDBFile('receptor.pdb')

# Load the AMBER protein force field through OpenMM.
from simtk.openmm import app
omm_forcefield = app.ForceField('amber14-all.xml')

# Parameterize the protein.
receptor_system = omm_forcefield.createSystem(receptor_pdbfile.topology)

# Convert the protein System into a ParmEd Structure.
receptor_structure = parmed.openmm.load_topology(receptor_pdbfile.topology,
                                                 receptor_system,
                                                 xyz=receptor_pdbfile.positions)

"""
Explanation: Create a ParmEd Structure of an AMBER-parametrized receptor
We have to create a ParmEd Structure containing positions and parameters of the receptor (BRD4) that we will then combine with the positions and parameters in the ligand Structure we created above. 
First, we use OpenMM to assign Amber14 parameters using OpenMM, but you can also create a Structure from Amber prmtop and inpcrd files.
<div class="alert alert-block alert-info">
    <b>Note:</b> If you already have AMBER (prmtop/inpcrd), GROMACS (top/gro), or any other file specifying the protein parameters supported by ParmEd, you can simply load the files directly into a Structure using ParmEd's functionalities. See https://parmed.github.io/ParmEd/html/readwrite.html .
</div>
End of explanation
"""


complex_structure = receptor_structure + ligand_structure

"""
Explanation: Combine receptor and ligand structures
We can then merge the receptor and ligand Structure objects into a single Structure containing both positions and parameters.
End of explanation
"""


from simtk.openmm.app import NoCutoff, HBonds
from simtk import unit

# Convert the Structure to an OpenMM System in vacuum.
complex_system = complex_structure.createSystem(nonbondedMethod=NoCutoff,
                                                nonbondedCutoff=9.0*unit.angstrom,
                                                constraints=HBonds,
                                                removeCMMotion=False)

# Export the System to an OpenMM System XML and PDB file
complex_structure.save('complex.pdb', overwrite=True)
from simtk.openmm import XmlSerializer
with open('complex.xml', 'w') as f:
    f.write(XmlSerializer.serialize(complex_system))

"""
Explanation: Export to OpenMM
Once we have the Structure for the complex, containing both positions and parameters, we can chose to create an OpenMM System object that we can simulate using OpenMM:
End of explanation
"""


# Export the Structure to AMBER files
complex_structure.save('complex.prmtop', overwrite=True)
complex_structure.save('complex.inpcrd', overwrite=True)

"""
Explanation: Export to Amber
We can also export the Structure to Amber prmtop and inpcrd files using ParmEd:
End of explanation
"""


# Export the Structure to Gromacs files
complex_structure.save('complex.gro', overwrite=True)
complex_structure.save('complex.top', overwrite=True)

"""
Explanation: Export to gromacs
We can export the Structure to gromacs gro and top files using ParmEd:
End of explanation
"""


