get_conformations.py: a simple script to generate rotamers

Here's a simple script for generating rotamers. It assumes the input geometry is reasonable and perturbs the dihedral angles of $n$ rotateable bonds by $\pm 120^\circ$. It makes all $3^n$ combinations of dihedral angles but if that is larger than a user-defined cutoff it will pick a random subset. Rotateable bonds include alcohol, phenols, and primary amines, but also amide bonds.  It doesn't know about equivalent atoms, so it will happily rotate a tert-butyl or trifluoromethyl group.

This work is licensed under a Creative Commons Attribution 4.0

Predicting pKa values using PM3 - conformer search part 3

This is a follow up to this post, but now we (that is Mads) use up to 1000 structures. To summarise: we use RDKit to generate $n$ starting conformations of neutral ("0") and protonated ("+") acebutolol and minimise them in aqueous solution using PM3/COSMO (MOPAC) and PM3/SMD /GAMESS) for a maximum of 200 steps.

The plot shows the lowest energy structure relative to the lowest energy found among all values of $n$. For example, the lowest PM3/COSMO value for the neutral form is found when generating 400 starting geometries and using 1000 starting structure, the lowest energy found is 1.1 kcal/mol higher.

I was a bit worried about the numerical stability of the COSMO and SMD methods, i.e. that the $n$ = 400 structure was generated in the $n$ = 1000 run, but had a higher energy for some reason. The figure above shows the lowest energy structure found ($n$ = 400) in purple and the structure most closely resembling it in the $n$ = 1000 search in green. So the $n$ = 400 structure was not generated in the $n$ = 1000 run.

Acebutolol has 12 rotateable bonds and roughly $3^{12}$ = 500,000 possible rotamers and RDKit samples these randomly. So, on hindsight, it is not surprising that $n$ = 400 and $n$ = 1000 randomly chosen samples do not both contain the lowest energy structure.

Incredibly (at least to me), the green structure is 1.7 kcal/mol higher in energy than the green. But, further optimisation of the green structure does not change energy, nor does changing the optimiser to EF instead of the default. I even rotated the molecule and restarted to the optimisation to test the robustness of the COSMO grid. But, no, this energy difference due to what I think is a minor strutural difference appears to be real. This means that the conformational PES is quite rugged, i.e. small structural changes can result in relatively big energy changes.

As it stands, pKa-predictions using PM3/COMSO or SMD (where I use $n$ > 100) will have a ca 1 pKa unit random error due to the conformational search and increasing $n$ is not really going to decrease that very much.

I'm not sure what to do about that.  We might try combining several low-$n$ conformational searches, with different random seeds. Alternatively, we may have to implement some sort of directed search, e.g. a genetic algorithm but that is going to be quite expensive for an SQM based method.

Alternatively, we have to pick the structures we use based on an MM screen where we can do a more thorough search but then we really need a MM/continuum solvation method interface.  If you have other ideas I'd be happy to hear about it.

Let me end with a challenge. Can you make a robust method that can find a PM3/COSMO energy equal to or lower than -183.01 kcal/mol for neutral acebutolol? And that is converged with respect to starting geometry and number of conformers tested? The MOPAC keywords are

pm3 eps=78.4 charge=0 cycles=200



This work is licensed under a Creative Commons Attribution 4.0

take_elementary_step.py: exploring chemical space one elementary step at a time

I have written a program called take_elementary_step.py. It is basically an implementation of the idea described in this paper by Zimmerman, but using the atom connectivity matrix approach described in this paper by Woo Youn Kim.

Zimmerman defines an elementary reaction step as a chemical rearrangement "with no more than two connections breaking and two connections forming simultaneously, while maintaining the upper and lower limits for coordination number of each atom." Two structures that are related in this way are very likely to be connected by a single transition state.

The code generates all such states for a given molecule and then selects the "best ones" based on a crude energy function based on bond dissociation energies. The idea is then to optimize the molecules using a QM method to refine the selection and then locate the TS connecting the selected structures to the input structures.  The atom ordering should be preserved (haven't checked yet) so they are suitable for interpolation-based methods such as NEB or GSM.

The default is to lists all structures with energies no higher than 200 kJ/mol compared to the input structure and to use homolytic bond cleavage, i.e. a C-H bonds is broken into a C and H radical, rather than, say, C^- and H⁺, but this can be controlled by the "heterolytic" keyword.

Using the heterolytic option will result in incorrectly charged fragments in some cases, but these should be weeded out by the energy criterion. For example, if you input acetylene the code will generate HC^3- + HC^3- instead of HC^3- + HC³⁺, but I didn't feel like writing complicated code to fix this since these possibilities will no be presented to the user because they are too high in energy.

The figure at the bottom shows the results using 1,3-butadiene + ethene as input. For some reason, RDKit has trouble displaying H₂ molecules so they have been removed and acetylene looks very weird (e.g. for comp64).  There are 88 structures in all, but that number could be reduced significantly be removing molecules with unsaturated atoms and/or three- and four-membered rings.

In his paper, Zimmerman shows 5 selected structures obtained using his implementation:

Taken from P.M. Zimmerman, J. Comp. Chem. 2013, 34, 1385. Copyright John Wiley & Sons. 

The first three also appear below as comp8, comp21, and comp54, while the last two are seemingly missing. However, they actually don't correspond to elementary steps since more than 2 bonds are broken or formed, but result from the subsequent QM optimization of comp7 and comp52 where additional bonds are formed between unsaturated C atoms.

This work is licensed under a Creative Commons Attribution 4.0

Simple transformation rules lead to complex and beautiful molecules

I'm messing around with some code for a new project. The code starts with the first molecule shown below (NCCN) and repeatedly applies two transformation: it inserts a C atom between 2 C atoms (CC>>CCC) and cyclizes 3 C atoms (CCC>>C1CC1). The image below shows the result for 7 cycles of this, i.e. cycle 1 creates NCCCN, cycle 2 creates  NC1CC1N and NCCCCN, etc.

After 7 cycles these two simples transformations lead to 121 molecules, some of which are quite beautiful in my opinion. I wonder how many of these have been made.

This work is licensed under a Creative Commons Attribution 4.0