Showing posts from August, 2014

ROCS-like shape overlap in rdkit

One year ago, I gave a brief talk at the RDKit user group meeting in Cambridge. As usual, I ranted about how it would be fantastic to have robust open-source computational chemistry tools. In particular, a commonly used tool from the 80s – gaussian-based shape overlap for pairs of molecules – was only covered by a single open-source package. So, together with help from some friends (big shout out to Greg Landrum, rdkit) I set out to change that. Took almost a year but we're getting there!

The only open-source package doing shape overlap was shape-it from silicos-it. We basically ripped the code out of it and incorporated it into rdkit. The aim was to provide python-scriptable shape-overlap tool. The silicos-it is run by Hans de Winter – big thanks to him – I'm not sure what's going on now but the silicos-it webpage down.

Greg did most (or certainly a _lot_) of work cleaning up the code. Now, we're finally moving towards testing against standard benchmarking sets. Her…

Conformer generation

Edited: update the figure after working around the RMSD alignment (rdkit tries all the atom order permutations to find the best alignment, rather than keep the atom order fixed, i think); removed the technical problems paragraph. New figure. 

"Conformer generation is an art". That statement from a more senior computational chemist initially made me chuckle – I mean, how hard could that be? If computational methods are claiming to bring a lot to the table, surely conformer generation is a problem solved long ago... Oh the naive young mind ;) 

Let's start with a simple question here is: if we start from the native (experimentally observed) conformation of a ligand and minimize it using two popular ligand force-fields (UFF and MMFF), we should stay close to the crystal structure. Is that indeed the case?

Figure 1
So that doesn't look so bad! Let's hydrogenate the structures – right now they're mixed bag.

Figure 2
Woops, adding hydrogens shouldn't incre…

Make a molecule dissapear in gromacs

This is a free energy perturbation basics – making molecules appear and disappear and using that to estimate free energies.

This code can also be used ror other, imaginative purposes. To make a molecule 'Ligand' disappear, use the following

integrator = ld
dt = 0.001
nsteps = 10000
free_energy          = yes
init_lambda          = 0.00
delta_lambda         = 1e-4
sc-alpha             = 1.5 #
sc-power             = 1.0
couple-moltype       = Pore
couple-lambda0       = vdw-q
couple-lambda1       = none

Now let's break it down:
free_energy = yes  is an on-switch for this feature of gromacs
init_lambda = 0.0 defines how invisible is the molecule at the start of the run, a value of 0 indicates that it's fully visible

delta_lambda = 1e-4 will define by how much the molecule is disappeared at every step of the simulation. Here this has to match the nsteps value - after 10000 steps the molecule will be fully gone.

couple-lambda0 = vdw-q shows that at the st…

High-frequency lunch booking

It's summer time and a 750-year-old tradition is that Merton college SCR (the senior common room) opens up its lunch sittings to the MCR members (middle common room). It's quite a treat – the quality of the meal greatly surpasses what MCR members get in other times of the year.

There is a catch obviously: the number of seats available to the MCR members is limited. The seats become available some time between Fri-Monday each week. Nobody has the time to sit there and keep refreshing the page so why not use a simple python web-crawler, executed periodically from crontab to maximize your odds of getting a spot?

Source code is attached! Enjoy ;)

EDIT: Due to popular demand from college folks, the tool has been converted to web app. It's called, HackAHall and it's available (invite only) here ;)



Acetylated lysine and negatively charged cysteine for CHARMM protein force-field

The CHARMM protein force-fields in gromacs provides topologies for many commonly found protein residues. However, when one wants to include a slightly more exotic residue - say a modified amino acid - the process can be a little tedious.

Below, two example topologies are shown: the acetylated lysine (ALY) and negatively charged Cys (CYN); treat those as an extension to CHARMM protein force-field. ALY is can be used to model modifications on histone proteins, while CYN is supposed to model Cys in zinc-finger motifs where it coordinates Zn2+ ions.

The two topologies below are simple hacks. At no point am I claiming that they are 'correct' in any way – but they provide a starting point for any improvements one wishes to make.

Edit: I recently bumped into a paper called Force field parameters for the simulation of modified histone tails where quantum calculations are done to obtain CHARMM parameters for the modified aminoacid residues. This could be a more reliable set of paramet…