Improvements to take_elementary_step.py

I've improved take_elementary_step.py that I wrote about previously.  

First of all I have fixed a bug in the bond order energy function so the program now suggests significantly fewer structures. Second I have written some code that creates input files for a QM program (GFN-xTB) and analyses the output. I choose GFN-xTB because it's a super fast semiempirical QM method that is also reasonably accurate. I have tried to keep things fairly modular so that it's easy to interface it with other QM programs (more details in the code comments). 

Using the Diels-Alder reaction as an example, staring with ethene and butadiene as the reactant the program now suggests 28 other structures with energies no more than 100 kcal/mol higher than the reactants. I choose 100 kcal/mol because the energy function is very crude.

All structures are then optimized with GFN-xTB and I now select structures with energies no more than 20 kcal/mol higher than the reactants. I choose a lower value here because I trust the energy function a bit more. I use the electronic energy plus an estimate of the translational entropy, but the code can easily be modified to use free energies. This narrows it down to 21 structures (ranked by energy below) in addition to the reactants, so the bond order energy function is actually doing quite well.

You'll notice that some of the structures appears to have two H atoms, but this is in fact an H₂ molecule. The reason for this is that the code that translates the Cartesian coordinates back to SMILES uses a distance criterion to assign bonds and that apparently needs some tweaking for H₂.  It's important to recompute the SMILES because the bonds can rearrange during the geometry optimization.

I also need to fix the code for calculations on atoms where GFN-xTB doesn't print out coordinates.

The next step is to find TSs connecting all these structures to the reactants to find the lowest barrier.

This work is licensed under a Creative Commons Attribution 4.0

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

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

xyz2mol: converting an xyz file to an RDKit mol object

I have written a program called xyzmol that converts an xyz file to an RDKit mol object, based on this paper.  For ions the molecular charge has to be specified in the xyz file.

In principle this could be accomplished by using OpenBabel to convert an xyz file to and sdf file and reading the sdf file with RDKit. However, in my experience OpenBabel occasionally mis-assigns bond orders and I believe this code does this less often.

The code works by constructing an atom connectivity (AC) matrix from the xyz file, converts the AC matrix to a bond order (BO) matrix, computes the formal atomic charges, and then constructs a molecule object using all this information, including the 3D coordinates.

The code that constructs the AC matrix from the xyz file (xyz2AC) uses scaled covalent radii to find bonds and the current scale factor is 1.25, but this may need adjusting.

The code currently works only for organic molecules, i.e. molecules containing the following elements: H, B, C-F, Si-Cl, Br, and I.

The current version assumes a close shell molecule, i.e. it does not work with radicals yet.

I hope you find it useful.

This work is licensed under a Creative Commons Attribution 4.0

TS_conf_search: Conformer search for transition states

Here's som prototype code for performing a conformer search for a transition state.  The idea is that you have found a TS for small model of the system and now you want to attach floppy ligands, generate different conformers for these ligands, and MMFF minimise them while keeping the core (model) frozen.  Each conformer is then used as an initial guess for a TS search, perhaps preceded by a constrained minimisation using the QM method.

Usage
To use the program type "python TS_conf_search ts.sdf small_model.sdf  5".

ts.sdf is an sdf file that contains a guess structure of the TS you want to do a conformational search for.  The coordinates of the small model-part of the structures does not have to match that of the small model and it doesn't even have to be very TS-like but it has reflect the bonding in the small model. Alternatively you can provide a SMILES string of the molecule.

small_model.sdf is an sdf file with with the optimized coordinates for the small model TS.  You can use any name for the file as long as the extension is sdf. TS_conf will only use the coordinates of the non-hydrogen atoms, unless the word "keepHs" is part of the file name.  If you use keepHs then the hydrogen atoms have to have corresponding  hydrogen atoms in the big model (see examples below)

5 is the number conformers you want to generate and you can use any integer you want here.  The default is 20.

What TS_conf_search does
TS_conf_search uses the ConstrainedEmbed function implemented in RDKit, which generates a conformer of a molecule, finds the atoms in the molecule corresponding to the small model, and then performs a FF minimization in which the Cartesian coordinates of these atoms are constrained to match the Cartesian coordinates of the small model.  If the small model contains two separate molecules then intramolecular energy terms are ignored.

sdf files
sdf files contain the Cartesian coordinates of the atoms but also explicitly defines the bond orders (single, double, etc) and any formal charges.  It is important the the bond orders in the ts.sdf and small_model.sdf files match. sdf files are basically the same as mol files, but if your program generates mol files be sure to change the file extension to sdf.  You can use OpenBabel to generate sdf files from nearly any other format (however, see next section).

Examples (and when to use keepHs)
The GitHub repository contains two Diels-Alder examples, where the small models are the PM3 optimized transition states for R = H, and the big model is R = CH2CH2CH2CH3.

I created the small model sdf files by converting the GAMESS output files to sdf using Open Babel.  However, the bond order didn't correspond to that shown in the picture above, so I manually edited the bond orders in the sdf file and remved the "RAD line" that OpenBabel created.

For Diels_Alder1 we type 'python TS_conf_search.py "Diels_Alder1;ClC=C.C=C(CCCC)C=C" Diels_Alder1_small.sdf 5' i.e. we build the TS guess from the SMILES string of the molecule.

For Diels_Alder2 we type "python TS_conf_search.py Diels_Alder2.sdf Diels_Alder2_small_keepHs.sdf 5". Here I got Diels_Alder2.sdf by loading the small model in to Avogadro, adding the CH2CH2CH2CH3 chain in some arbitrary conformation, and minimising the  structure.

In the second example, I use keepHs because the dienophile (ethene) only has two heavy atoms, which is not enough to define the orientation of the plane of the molecule relative to the diene. Because I use keepHs I have removed the hydrogen corresponding to the R group, so that all the hydrogens in the small model can be matched to hydrogens in the large model.

While Diels-Alder2 provides one example where keepHs is useful.  Another is where hydrogens atoms are an integral part of the mechanism, such as in hydrogen or proton transfer TSs.

Results

I used each of the geometries as a starting guess for a TS search using GAMESS where I calculated the Hessian to start with (hess=calc) and after every fifth step (ihrep=5).

In the first example none of the five runs found the correct TS, while for the latter example all 5 worked. When I re-ran the first example with python TS_conf_search.py "Diels_Alder1;ClC=C.C=C(CCCC)C=C" Diels_Alder1_small_keepHs.sdf 5 it worked fine.

Not using keepHs has the advantage that one doesn't have to make different small model sdf files when different hydrogens are replaced by substituents. But in this case an extra constrained minimisation step appears to be needed to get the hydrogen atoms placed correctly.



This work is licensed under a Creative Commons Attribution 3.0 Unported License.

Atom_mapper: determining atom correspondence between reactants and products

Many TS search algorithms interpolate between reactants and products and require that the atom order is identical. This is time consuming to enforce by hand and hard to automate, but this paper by the folks at Schrödinger describe a solution that seems to work pretty well. The general idea is also fairly easy to implement if you are familiar with RDKit or similar cheminformatics tools and here I describe my implementation so far.

The basic idea
If two molecules are identical but have different atom numbering then RDKit can map the atoms using "GetSubstructMatch". So (pseudocode) "OCC.GetSubstructMatch(CCO)" gives "(2,1,0)" meaning that atom 2 in CCO corresponds to atom 0 in OCC and so forth.  And "RenumberAtoms(CCO, (2,1,0))" gives OCC.

Let's say our reaction is CH3-CH=O -> CH2=CH-OH.  "product.GetSubstructMatch(reactant)" will not give a match because the molecules are different. To focus on connectivity we remove bond orders: CH3-CH-O -> CH2-CH-OH, but the molecules are still different. The next step is to break each bond in turn and compare. For the reactant you get [H + CH2-CH-O, CH3 + CH-OH, H + CH3-C-O, O + CH3-CH2] and the first entry will match the last entry in the corresponding list for the product [H + CH-CH-OH, CH2 + CH-OH, H + CH2-C-OH, H + CH2-CH-O].

It is very easy to check to for equivalence by converting the list of fragmented molecules to a set of corresponding SMILES strings and looking for identical elements in the two sets: "matches = list(set(react_frag_smiles) & set(prod_frag_smiles))". If there is a match then the atoms can be matched as before because they are identical (note that the bond breaking does not alter the atom numbering of the molecule). If no match is found then the procedure is repeated for the fragments.

One problem is that equivalent atoms are not matched in 3D. For example, if you label (i.e. number) the equivalent H atoms in a methyl group then the methyl group is chiral and RDKit will create a random chirality.  So, the Cartesian coordinated of methyl hydrogens in the CCO example above are not necessarily superimposable, which creates problems in TS search algorithms that use interpolation.

Following the ideas presented in the Schrödinger paper, I try to solve this issue by switching pairs of equivalent atoms and comparing the RMSDs between reactant and product before and after.  However, this will only work if the reactants and products have the same conformation, so I use a modified version of ConstrainedEmbed (ConstrainedEmbedRP) which uses either Cartesian or distance constraints to match two structures using a UFF force field minimization of one of the structures.

Using ConstrainedEmbed has the added advantage of producing reactant and product structures that are in the same orientation and, therefore, often good starting points for interpolation.

Atom_mapper
* I have implemented these ideas in Atom_mapper and made it available on GitHub under the MIT license.
Usage  
* The input is SMILES strings of reactant and products but the code can be adapted to work with user-supplied coordinates via sdf files.
* The output is two sdf files labelled "template" and "mol" with (hopefully) identical atom ordering. Template is the most rigid of the two and the coordinates come from an unconstrained UFF minimization. Mol is superimposed on template using ConstrainedEmbedRP.
* Both geometries should probably be minimized using the energy function you plan to use for the TS search and, depending on the optimization algorithm, the structures should probably be rigidly superimposed before interpolation
Current limitations
* I haven't tested the code extensively. If you run into cases where it doesn't work, file a bug report on GitHub
* No conformation search is done so the conformations chosen will be random. Solution: start with coordinates rather than SMILES.
* If reactants or products consists of two molecules their relative orientation resulting from the UFF minimization will be somewhat arbitrary. This is not a problem if only the reactants or products consists of two molecules, since the coordinates will be superimposed on a more rigid molecule. However, it is a problem if the reactants and products both consist of two molecules. Solution: start with coordinates rather than SMILES, but not sure how to determine best intermolecular orientation automatically.
* The code only permutes pairs of equivalent atoms so it won't work for equivalent groups, e.g. methyl groups in CC(C)(C)CO because we need to switch several pairs (corresponding to C and 3 H atoms) simultaneously.  Not sure how to fix that easily. Suggestions welcome.

There are also some general limitations to this "active bonds" approach as described in the article.


This work is licensed under a Creative Commons Attribution 4.0

Automatic generation of a set of molecules

Many quantum chemistry projects have reached a point where setup and analysis consumes more human time than CPU time, e.g. it takes all day to set-up enough input files to keep the computer busy overnight. Many people use scripts to automatically extract and analyse the data, but few use scripts to generate the coordinates.

Here I show how this can easily be done using the RDKit toolkit.  The following python script adds F and OH substituents to benzene making all possible 91 combinations of mono-, di-, hexa-substituted molecules.

However, the script can easily be changed to do something else by changing parent_smiles and rxn_smarts_list in line 4 and 5.  If you are not familiar with SMILES start here and there are plenty of GUIs, such as this, that generate SMILES.

To use the Reaction SMARTS you have to learn SMARTS, which can be a bit tricky, but it is a very powerful tool. For example, if you change [cX3;H1:1]>>[*:1]F to [*;H1:1]>>[*:1]F then the program will add H to any atom with one H, i.e. also the OH group to create the OF substituent.  So it you set substitutions = 2, you'll get mono-substituted Ph-OF in addition to mono- and di-substituted Ph-F and Ph-OH.

Similarly, ff you add [cX3;H1:1][cX3;H1:2]>>[c:1]1[c:2]cccc1 to the list (and use substitutions = 2) you'll get un- and mono-substituted napthalene as well as un-substituted anthracene and phenanthrene.

In my experience, the only thing that limits what I can build with this approach is my understanding of SMARTS.  Hope this is of some use to you.



This work is licensed under a Creative Commons Attribution 4.0

Finding the reference molecule for a pKa calculation using RDKit

This is prototype code related to this project.  I use Histamine as an example which has two ionizable sites: a primary amine and an imidazole ring.

The code figures out that imidazole is the best reference for the imidazole ring, while ethylamine is the best reference for the primary amine.  The code does this by figuring out which atom is being deprotonated, computes the Morgan fingerprint around this atom, and compares it to the Morgan fingerprints of imizadole and ethylamine.

Conformer search with RDKit

I'm teaching myself how to use RDKit.  Here is code for conformer search using RDKit that also computes the energy of each conformer using the MMFF94 force field.

Comments welcome


This work is licensed under a Creative Commons Attribution 3.0 Unported License.  

A cheminformatics problem: protonate with SMILES and InChI

I am teaching a python programming course and this is one of the projects I want to try this year.

The over all goal of the project is to write one or more programs that generates protonation states for a list of nitrogen containing molecules specified by name. The project uses SMILES, and maybe InChI, which you can read more about here.

Getting started
1. Write code that from this list molecules = ["CCN", "CNC", "CN(C)C"] generates this output

CCN
CC[NH3+]

C[N@H]C
C[N@@H]C
C[NH2+]C

C[N@](C)C
C[N@@](C)C
C[N@H+](C)C
C[N@@H+](C)C

2. Write code that from this list  molecules = ["C(C(=O)O)N"] generates this output (the order is not important)

C(C(=O)O)N
C(C(=O)O)[NH3+]
C(C(=O)[O-])N
C(C(=O)[O-])[NH3+]

3. C(C(=O)O)N is the amino acid glycine.  Extend this program to work for alanine, asparagine, aspartate, and lysine.  Use this site to get SMILES strings for these amino acids. Find a picture of asparagine and make sure you're treating the side-chain correctly.

The project
4. (optional) Figure out how to generate a file containing SMILES strings from a file containing names. The best way is probably bash.  Get inspiration here, here, and here.

5. Generate all possible protonation state SMILES for the molecules in Table 2 in this paper. If you completed step 4 you can use tools like https://pdftables.com/ to generate a file with the names.

6. Repeat for Table 1 and 3

7. (optional) The neutral form of the amino acids histidine and arginine side chain groups have tautomers. Generate SMILES for all tautomers (InChI might help you identify tautomers).

8. (optional) Do any of the molecules in step 5 and 6 have tautomers? If so generate SMILES for all tautomers.

Some code snippets to get you started



This work is licensed under a Creative Commons Attribution 4.0