Posts

Showing posts from 2015

Native contact (Q) biasing in GROMACS and plumed2

Fraction of native contacts is a pretty darn good coordinate for looking at protein folding, definitely when it comes to small proteins. Just have a look:

http://www.pnas.org/content/110/44/17874.abstract

Fraction of native contacts can be used to analyze long, unbiased MD runs – like those from DESRES – or in biasing simulations. But how does it fare when used in biasing? Well I don't know – but you can now find out by yourself! Either in umbrella sampling or metadynamics – or a method of your choice that doesn't use energy as your collective variable.

Here, I'm happy to report that I incorporated an implementation of native contact CV into plumed2 – an enhanced sampling package that can be used with a number of popular MD simulation engines. The pull request is now closed https://github.com/plumed/plumed2/pull/177 and it's likely that the next release of plumed2 will include this.

"Tutorial" is a regression test in plumed2 sources  regtest/basic/rt29/ which …

Scripting paramchem.org access

Below is a simple python script that allows one to programatically access paramchem.org, a small-molecule parametrization service. Autocorrect said "problematic" instead of "programatic", maybe I'm noting getting the hint...

The script has minimal dependencies and is very "light-weight", meaning no error handling or error messages. A full working example for a benzene (duh!) is attached.

Long time ago, I had to implement something similar to this script, so this was a fairly easy job. It relies on mechanize and beautiful soup to do it's thing.

python paramchem.py -u "YOUR USERNAME" -p "YOUR PASSWORD" -c benzene.pdb

# manually correct the "RESI" line in benzene.str to be 
# RESI BNZ 0.000 ! param...

python cgenff_charmm2gmx.py BNZ benzene.mol2 benzene.str charmm36-jun2015.ff/

Source code and examples
git clone https://github.com/jandom/paramchem
Example on google drive

Edited: note on errors, link to github repo

Feeling the burn: DPPC lipid CHARMM36 with gromacs

The "oh shit" moment came when I started running a short DPPC run on trusty gromacs 4.6 with CHARMM36. The bilayer switched from liquid crystalline to gel phase and I knew I was in trouble. Played with the system size, increased the temperature from 323 to 333K -- no luck. Turns out there is some pretty strong sensitivity of the lipid properties to electrostatics treatment in gromacs.

The key paper is probably Molecular Dynamics Simulations of Phosphatidylcholine Membranes: A Comparative Force Field Study (paywall) which demonstrates the sensitivity to electrostatic parameters and proposes some solutions.

The gromacs authors have suggested their own parameters for running CHARMM with gromacs

The topic has attracted a lot of attention on the gromacs mailing list, two useful threads are below:
https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-users/2013-August/083370.htmlhttps://www.mail-archive.com/gmx-users@gromacs.org/msg64050.html

Publication-quality matplotlib figures

Image
We're going to go from this...
... to this


Install prettyplotlib, a useful tool for graph styling


sudo pip install prettyplotlib 

Add a new method to style the axes to prettyplotlib/util.py 

def styling(ax, xpad=10.0, ypad=10.0, axiswidth=2.0, axistickwidth=2.0, axiscolor="#333333"):

    for axis in ['bottom','left']: ax.spines[axis].set_linewidth(axiswidth)

    for tick in ax.get_xaxis().get_major_ticks():
        tick.set_pad(xpad)
        tick.label1 = tick._get_text1()

    for tick in ax.get_yaxis().get_major_ticks():
        tick.set_pad(ypad)
        tick.label1 = tick._get_text1()

    ax.get_yaxis().set_tick_params(direction='out', width=axistickwidth)
    ax.get_xaxis().set_tick_params(direction='out', width=axistickwidth)

    ax.spines['bottom'].set_color(axiscolor)
    ax.spines['left'].set_color(axiscolor)

    ax.tick_params(axis='x', colors=axiscolor)
    ax.tick_params(axis='y', colors=axiscolor)

Set figure siz…

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. Atta…

Using gdb with gromacs

Download and compile gromacs 5.x

cd gromacs; mkdir build ; cd build 
cmake .. -DCMAKE_BUILD_TYPE=Debug

Then, go to the methanol example

cd share/gromacs/tutor/methanol/
grompp 

Then enter gdb console and issue some commands

$ gdb
# that's where the gmx command lives
exec ../../build/bin/gmx
# read in the symbol table
file ../../build/bin/gmx 
# break on errors
break _gmx_error
# break on do_force function
break do_force
run mdrun -debug 1 -v -s

This is not very practical - you can put the gdb commands into a file, gdb_cmds, and execute

$ gdb -x gdb_cmds

Note gdb has a ton of useful commands, google around

References:

- add pretty command to gdb

http://stackoverflow.com/questions/23578312/gdb-pretty-printing-with-python-a-recursive-structure

- gdb tutorial handout
http://www.cs.umd.edu/~srhuang/teaching/cmsc212/gdb-tutorial-handout.pdf

Gromacs games

Additional energy terms in custom Hamiltonians, exchanges The standard energy function in looks something like this:

E = E_bonds + E_angles + E_dihedrals + E_LJ + E_Culombic

This energy term, at some temperature T then becomes:

U = exp(-E/kb*T)

Let's say one has a n such Hamiltonians, U1...Un. These can then exchange, if the energies overlap. The details of the exchange scheme are skipped here, for gromacs implementation if the exchange see src/kernel/repl_ex.c

In temperature replica exchange one uses a ladder of T values, hoping that U1...Un overlap. In solute tampering replica exchange, one modifies the energy function E, keeping the T constant, hoping again to see exchange.

In deciding if a pair of replicas on a ladder will exchange, one compares their potential energy (Epot).

When including extra energy terms into the hamiltonian (custom biases, for example) one has to make sure they will be included in the Epot.

For instance, let's see if the gromacs pull code adds its e…

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 my gu…

complie plumed 2 with gromacs 4.6 on ubuntu 14.04

Should be simple right? Just you wait... Because I do it for the 5th time now, this will deserve a blog entry, if anything for my own sanity. 


sudo apt-get install libblas-dev liblapack-dev
tar xvf plumed-2.1.1.tgz
cd plumed-2.1.1
export PLUMED_PREFIX=~/Programs/plumed/2.1.1/
./configure.sh  << EOF
1
EOF
make
make install


[Edit]For OpenMPI plumed installation do the following

sudo apt-get install libcr-dev mpich2 mpich2-doc./configure.sh  << EOF
3
EOF
sed -i s#'-lblas'#'-lblas -L/usr/lib/openmpi/lib'#g Makefile.conf sed -i s#'PLUMED_INCLUDE)'#'PLUMED_INCLUDE) -I/usr/lib/openmpi/include'#g Makefile.conf
The default ./configure command doesn't do the job for me. So let's ignore that for now.

cmake .. -DCMAKE_INSTALL_PREFIX=/path/ -DGMX_MPI=OFF -DGMX_GPU=OFF
make 
make install 

And then if all goes well

source /path/
mdrun -h 2>&1 >/dev/null | grep plumed
-plumed  plumed.dat  Input, Opt.  Generic data file

Sadly, this is just a serial vers…

gromacs 5.0 - first look

Image
Oh, gromacs, old friend let's see what you have to offer in your new release


wget ftp://ftp.gromacs.org/pub/gromacs/gromacs-5.0.4.tar.gz
tar xvf gromacs-5.0.4.tar.gz
cd gromacs-5.0.4
mkdir build
cd build
cmake .. 

Problem #1

-- Could NOT find LibXml2 (missing:  LIBXML2_LIBRARIES LIBXML2_INCLUDE_DIR)
CMake Warning at CMakeLists.txt:601 (message):
  libxml2 not found.  Will build GROMACS without unit-tests.  This is not
  recommended, because the unit-tests help to verify that GROMACS functions
  correctly.  Most likely you are missing the libxml2-dev(el) package.  After
  you installed it, set GMX_BUILD_UNITTESTS=ON.


Problem #2

CMake Warning at cmake/gmxManageFFTLibraries.cmake:94 (message):
  The FFTW library was compiled with --enable-avx to enable AVX SIMD
  instructions.  That might sound like a good idea for your processor, but
  for FFTW versions up to 3.3.3, these are slower than the SSE/SSE2 SIMD
  instructions for the way GROMACS uses FFTs.  Limitations in the way FFTW
  allows GROMACS to…