Posts

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)  ...

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