Using gromacs with two force-fields are for ligand one for protein? Wanting to do that but not knowing how? Look no further.

The technical problem is as follows: let’s say you have your favourite custom force-field. For me it resides in charmm.ff and it’s my hybrid charmm36+charmm22* set of parameters.

To simulate a ligand inside of a system parametrized using the ff-above, I went to paramchem.org and obtained small molecule parameters for my ligand. These were in CHARMM format, so (courtesy of Alex Mackerell) I converted those to gromacs itp files. That left me with a dependency on charmm36-jun2015.ff.tgz

There is only a handful of lines I want from charmm36-jun2015.ff, the ligand is very simple. There is no easy way to get those out, or to merge my custom charmm with Alex’s.

Prerequisites

Install networkx, a python package for dealing with graphs:

pip install networkx

Let’s take some example ligand molecules, such as the TPP - the ligand for EmrE drug transporter. Attached here, after protonating it in maestro.

First iteration, the fast way

On paramchem.org, select the checkbox “Include parameters that are already in CGenFF” in paramchem to include copies of the parameters in the paramchem output. This removes the dependency on charmm36-jun2015 - everything needed is explicitly copied into the .str and then the .prm files. Then directly include the .prm into your force-filed and done.

Second iteration, the clever way

Instead of downloading an .str with all the parameters explicitly included, we’ll take it the default way. Then we’re going to use our own code to extract the parameters for the isolated component in the system. The end result is the same: a stand-alone itp with all the parameters for a given molecule of interest.

python cgenff_charmm2gmx.py Ligand Ligand.mol2 Ligand.str charmm36-jun2015.ff
grompp -pp -p ligand.top
# exctract only the parameters needed from processed.top
python forcefield.py processed.top processed.itp

The script forcefield.py can be downloaded from github.

processed.itp contains only the parameters needed to run this ligand molecule: append all the [ *types ] sections to your ffbonded.itp and ffnonbond.itp, then #include the reminder of processed.itp in your topol.top and fly off with your simulations. With any force-field.

Advantages of the slow way - separation of parameters

Using forcefield.py you could obtain molecule.itp file for your lipid, protein and small molecule parameters - potentially using 3 different force-fields. These can then be sanely combined together and used in simulation. How to combine will be the topic of the next entry.