Merge force-fields in gromacs

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 interaction, 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 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.




Comments

Post a Comment

Popular posts from this blog

React.js – edit and delete comments in CommentBox

It appears customary for high-intelligence, high-skills groups of engineers to "roll their own" whenever they have a chance to. Counter-intuitively to most managerial/MBA types out there, re-inventing the wheel can give you... a better wheel. Success is not guaranteed though – but these people may leave sooner rather than later if you don't let them... Having looked at number of JS frameworks (ember, backbone, angular) there was no clear winner. They were all a mess, at least to me. So I approached React.js with no hopes – how much better could it possible be? It turns out, a lot! This post will be an extension of the official React.js tutorial . The tutorial stops at the stage where you can add comments to a dynamic lists, that will grow as you add comments. This falls short of the functionality displayed on facebook: once you add a comment you can edit or delete it, so that's what're going to do. Full disclosure, I only started using React this weekendo so m...
Read more

Example slurm cluster on your laptop (multiple VMs via vagrant)

To anybody who's considering a switch to modern queuing systems, like slurm, this will be an useful guide. Rather than doing a roll out on production nodes, we'll try things out in vagrant. This text assumes you have vagrant installed on your local machine (it's very easy). The ultimate aim will be to setup 2 servers (worker nodes) and one master/controller, and run a simple job on them. There are multiple great guides out there: most of them are out of date, and won't work without some tweaking on Ubuntu 16.04 Xenial. The entirety of the code is hosted at https://github.com/jandom/gromacs-slurm-openmpi-vagrant The README.md file will show you how to set it up, with checks on every step. References https://mussolblog.wordpress.com/2013/07/17/setting-up-a-testing-slurm-cluster/ https://github.com/gabrieleiannetti/slurm_cluster_wiki/wiki/Installing-a-Slurm-Cluster http://philipwfowler.me/2016/04/14/how-to-setup-a-gramble/ https://github.com/dakl/slurm-clus...
Read more

Wrapping openbabel in python - using cython

The previous post showed a contrived example of how one could access the openbabel functionality in python using the boost libraries. There is alternative to boost to wrap around c/c++ code, it's called cython. Here is an example openbabel wrapper, exposing some very simple functionality. Just before diving in, you should know that there already exists a wrapper around openbabel in python (it is called pybel). What I'm showing here doesn't even come close to pybel in terms of usefulness, the idea here is to show a demo of how easy wrapping things in python is. #! /usr/bin/python from ez_setup import use_setuptools use_setuptools() from setuptools import setup, Extension, find_packages from Cython.Distutils import build_ext import sys, os import glob import subprocess as sub def get_include(name="openbabel"): p = sub.Popen('locate %s' % name ,stdout=sub.PIPE,stderr=sub.PIPE, shell=True) output, errors = p.communicate() ...
Read more