from simtk.openmm.app import PDBFile

from openff.toolkit.utils import get_data_file_path
from openff.toolkit.topology import Molecule, Topology
from openff.toolkit.typing.engines.smirnoff import ForceField

# Create an OpenFF Topology of toluene from a pdb file.
toluene_pdb_file_path = get_data_file_path('molecules/toluene.pdb')
toluene_pdbfile = PDBFile(toluene_pdb_file_path)
toluene = Molecule.from_smiles('Cc1ccccc1')
off_topology = Topology.from_openmm(openmm_topology=toluene_pdbfile.topology,
                                    unique_molecules=[toluene])

# Load the Parsley force field from disk.
force_field = ForceField('openff_unconstrained-1.0.0.offxml')

# Parametrize the toluene molecule.
toluene_system = force_field.create_openmm_system(off_topology)

"""
Explanation: Create a system mixing SMIRNOFF and non-SMIRNOFF-formatted force fields
This example shows how to create a receptor-ligand OpenMM System where the ligand (toluene) is parametrized with a SMIRNOFF force field and the protein (T4 Lysozyme) solvated in water is assigned AMBER and TIP3P-FB parameters through the ParmEd library.
We'll need two PDB files. One for the ligand in vacuum, and one for the solvated protein without ligand. The coordinates of the protein-ligand complex will be determined by these PDB files, so their positions needs to be consistent with the ligand being positioned in the binding pocket if this is the desired initial configuration of your simulation.
Parametrize a molecule with the Parsley force field
First, we parametrize the ligand (toluene) with the SMIRNOFF-format Parsley force field through the usual route to create an OpenMM System.
End of explanation
"""


import parmed

# Convert OpenMM System into a ParmEd Structure.
toluene_structure = parmed.openmm.load_topology(toluene_pdbfile.topology,
                                                toluene_system,
                                                xyz=toluene_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
"""


# Load the AMBER protein force field parameters through OpenMM.
from simtk.openmm import app
omm_forcefield = app.ForceField('amber99sbildn.xml', 'tip3pfb.xml')

# Load the solvated receptor from a PDB file.
t4_pdb_file_path = get_data_file_path('systems/test_systems/T4_lysozyme_water_ions.pdb')
t4_pdb_file = PDBFile(t4_pdb_file_path)

# Obtain the updated OpenMM Topology and positions.
omm_topology = t4_pdb_file.getTopology()
positions = t4_pdb_file.getPositions()

"""
Explanation: Create a ParmEd Structure of an AMBER-parametrized receptor in TIP3P-FB water
We have to create a ParmEd Structure of the receptor (T4 Lysozyme) to combine to the toluene Structure. Here we assign amber99sbildn to a T4 Lysozyme receptor solvated in TIP3P-FB water using OpenMM.
First, we load the parameters and the PDB file including positions for the protein, water, and ion atoms.
<div class="alert alert-block alert-info">
    <b>Note:</b> If you already have AMBER (prmtop/inpcrd), GROMACS (top/gro), or any other file supported by ParmEd specifying the parameters for the solvated protein, 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
"""


# Parameterize the protein.
t4_system = omm_forcefield.createSystem(omm_topology, rigidWater=False)

# Convert the protein System into a ParmEd Structure.
t4_structure = parmed.openmm.load_topology(omm_topology,
                                           t4_system,
                                           xyz=positions)

"""
Explanation: We then create a parametrized OpenMM System and convert it to a Structure.
Note the rigidWater=False argument in ForceField.createSystem(). This is necessary to work around a problem araising with ParmEd in reading the parameters of constrained bonds from an OpenMM System (see https://github.com/openforcefield/openff-toolkit/issues/259 for more details). We'll re-add the hydrogen bond constraints when we'll create the System for the complex.
<div class="alert alert-block alert-info">
    <b>Note:</b> If you don't solvate the system or if you load it directly from AMBER, GROMACS, or other files directly to ParmEd, you won't need extra precautions.
</div>
End of explanation
"""


complex_structure = toluene_structure + t4_structure

"""
Explanation: Combine receptor and ligand structures
We can then merge the receptor and ligand Structure objects to form the complex. Note that the coordinates of protein and ligand in this example are determined by the PDB file, and they are already consistent with the ligand being positioned in the binding pocket.
<div class="alert alert-block alert-info">
    <b>Note:</b> If you want to include water molecules in the binding site, you will have to be careful to place them so that they won't create steric clashes once the ligand is inserted.
</div>
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)

"""
Explanation: Convert back the structure into an OpenMM System
Once we have the Structure of the complex, we can chose to create a System object that we can simulate with OpenMM.
End of explanation
"""


