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/rt

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 gi

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.html https://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'

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 tra

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 pul

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

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. 

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 t