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 benc

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 s

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 grompp.mdp ... integrator = ld dt = 0.001 nsteps = 10000 ... ; FREE ENERGY 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-lam

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 ;) Thanks, Jan

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 pa