That's Classic!

I've very recently had the joy of setting up classical simulations for use in my research developing the next generation of force-fields. Since I've dealt with quantum mechanics and classical mechanics derived from quantum mechanics, my need for classical simulations on proteins and ligands has been void. Until now.

I only know of two tools to actually calculate Molecular Dynamics simulations on entire proteins: The Tinker package and the GROMACS package. I've tinkered(!) around with the first package and kept a good distance to the latter, but times have changed since I needed to include a ligand. GROMACS (and Tinker for that matter) are tuned to either proteins (in solution) or small molecules (in solution)*. If you want another molecule to use for the MD, you must first obtain the force-field parameters somehow. Since the most of us have no time for that, or that the next generation force-field is on its way but not quite there yet, what do you do?

To start from nothing and finish two days later, I did what any (in)sane scientist would do if no-one in the department has apparent experience with a particular piece of software: I Googled and I browsed the GROMACS "guides". The first is good if you know what to look for, and the latter is good to get inspiration, but it does not quite get you there since information always seem to be outdated just a bit.

Here is how I did it:

1. Find the structure of interest to you on the Protein Data Bank. Build your ligand into that system using your favorite build tool. You should now have a PDB file with a protein and a ligand.
2. Use SwissParam (more info below) to make a topology file (actually it is a force-field) of your ligand.
3. Follow the SwissParam Gromacs guide until step 6 (you must include this step!)
4. Follow the Protein-Ligand Complex guide by Justin Lemkul from the solvation step (it is good to actually read the entire guide). This will get you to make a production run MD of 1 ns of your own system.
5. Analyze your MD! This is out of scope for this post, but man was I happy when the MD actually successfully finished.

Note on SwissParam
The SwissParam tool takes a small molecule of interest and generates what is know in GROMACS-speak as a topology file, but in regular office-chatter, this is known as a force-field. Here is an exert from the main page of SwissParam

The data are derived from the Merck Molecular ForceField (MMFF). Dihedral angle terms are taken as is, while only the harmonic part of the bond, angle and improper terms are retained. Charges are taken from MMFF. Van der Waals parameters are taken from the closest atom type in CHARMM22.

The paper of SwissParam(doi:10.1002/jcc.21816) has been published so give it a read if you're really feeling geeky.

* I should note that if one chooses the MMFF 94 force field in Tinker, one should potentially be able to use it for everything, but it is a bit messy to actually set up and run not to mention that I found absolutely no way to automate the atom type-setting.