## Sunday, November 12, 2017

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

## Sunday, September 3, 2017

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

## Thursday, August 10, 2017

### Predicting most labile CH bond using semiempirical methods

Yesterday, I mentioned on Twitter that I had written some prototype code to predict labile proton using PM3/COSMO.  A few people expressed interest so I put the code on GitHub.

The code is based on some previous work and removes all CH protons/hydrides/H atoms and finds the position with the lowest free energy.  I've just finished the code so I have no idea how well it works or if PM3/COSMO is the best choice.  Also, there are very few comments and some aspects of the setup (paths, etc) are particular to my machine and needs to be changed.

In its current form the method doesn't consider different reference molecules, so the prediction is "pure PM3/COSMO".  The basic code to do it is there but I haven't extended it yet. Also, the code only considers CH groups, so if have other groups (OH, NH, etc) with lower pKa they won't be considered.  This is very much a work in progress but if you have suggestions, test cases, ideas, or any other kind of feedback please let me know.  I hope you find it useful.

## Saturday, July 15, 2017

### Planned papers for 2017 - six months in

In January I wrote about the papers I plan to publish and made this list:

Accepted
1. Protein structure refinement using a quantum mechanics-based chemical shielding predictor
2. Prediction of pKa values for drug-like molecules using semiempirical quantum chemical methods

Probable
3. Intermolecular Interactions in the Condensed Phase: Evaluation of Semi-empirical Quantum Mechanical Methods
4. Fast Prediction of the Regioselectivity of Electrophilic Aromatic Substitution Reactions of Heteroaromatic Systems Using Semi-Empirical Quantum Chemical Methods
5. Benchmarking cost vs. accuracy for computation of NMR shielding constants by quantum mechanical methods
6. Improved prediction of chemical shifts using machine learning
7. PM6 for all elements in GAMESS, including PCM interface

Probably not in 2017
8. Protonator: an open source program for the rapid prediction of the dominant protonation states of organic molecules in aqueous solution
9. pKa prediction using semi-empirical methods: difficult cases
10. Prediction of C-H pKa values and homolytic bond strengths using semi-empirical methods
11. High throughput transition state determination using semi-empirical methods

The status is

Published
1. Protein structure refinement using a quantum mechanics-based chemical shielding predictor
2. Prediction of pKa values for drug-like molecules using semiempirical quantum chemical methods

Probable (at least submission)
4. Fast Prediction of the Regioselectivity of Electrophilic Aromatic Substitution Reactions of Heteroaromatic Systems Using Semi-Empirical Quantum Chemical Methods

We had a draft ready to be sent in in mid-April, but decided to include a "few" more types of heteroaromatics and the study has now ballooned to nearly 600 compounds (up from about 150 in the original draft)! The calculations and most of the analysis is done, but the paper basically has to be rewritten, and this really has to be done by my synthetic chemistry co-author, since the study is now much more about the chemistries of these heteroaromatic and much less about the method.  Even though its out of my hands I am still hopeful that we can submit it this year.

9. pKa prediction using semi-empirical methods: difficult cases automation

This paper is 2/3 written and presents a completely automated PM3-based pKa prediction protocol. The method works quite well, but most outliers turn out to be due to high energy conformations. The main remaining issue is to find a conformer-search protocol that consistently gives low-energy conformations. Depending on how much time I have to devote to paper 4 and the proposal mentioned below, I am still hopefull I can get this published this year.

7. PM6 for all elements in GAMESS, including PCM interface

This will basically be a paper on the accuracy of the SMD method using semiempirical methods.  I'm hopefull we will submit this year, but there is still a lot to do.  Among other things, it explains why PM3/COSMO works best for pKa predictions: the solvation energy errors happen to be smallest for this method.

10. Prediction of C-H pKa values and homolytic bond strengths using semi-empirical methods

I am planning to submit a proposal on prediction of pKa,  homolytic bond strengths, etc as predictors of CH activation-sites. The proposals is due September 25th, so much of my "spare" time until then will be spent on getting preliminary results and writing.

## Sunday, May 14, 2017

### Predicting pKa values using PM3 - conformer search using implicit solvation

Disclaimer: These are preliminary results and may contain errors

In a previous post I showed that minimising RDKit-generated conformers generated with MMFF led to slightly worse PM3/COSMO pKa predictions. One reason might be that the MMFF minimisation is done in the gas phase.  Anders Christensen reminded me that TINKER can do MMFF/GBSA minimisations so I used TINKER so minimise 20 and 50 RDKit-generated conformers in used these as initial guesses for the PM3/COSMO energy minimisations (20mmtk and 50mmtk

As you can see there is relatively little difference in the overall performance of Xmm and Xmmtk. The error distribution is a bit tighter for the tk variants and 50mmtk does better at finding structures closer to the "global" minimum, but 50nomm is still the best choice for avoiding >1 pKa errors.

So, either GBSA is not a good substitute for COSMO or MMFF is not a good substitute for PM3.  In light of this recent paper my money is on the latter.  So I'll go with 50nomm for now.

## Saturday, May 13, 2017

### Computing apparent pKa values using QM energies of isodesmic reactions

A few days ago I presented some preliminary pKa results.  Almost all the molecules in the set have more than one ionisable group. Here's how I compute the apparent pKa values

The microscopic pKa values are computed by

\label{eqn:pka}
\mathrm{pK_a}=\mathrm{pK_a^{ref}} + \frac{\Delta G^\circ}{RT\ln (10)}

where $\Delta G^\circ$ denotes the change in standard free energy for the isodesmic reaction

\mathrm{ BH^+ + B_{ref} \rightleftharpoons B + B_{ref}H^+ }

If there are more than one titrateable sites then several protonation states can contribute to the apparent pKa.  Here I follow the approach described by Bochevarov et al. If $\mathrm{D}$ and $\mathrm{P}$ differ by one proton and there are $M$ deprotonated sites and $N$ protonated sites then

\mathrm{K_{app}}  = \frac{([\mathrm{D}_1] + [\mathrm{D}_2] + ... [\mathrm{D}_M])[\mathrm{H}^+]}{[\mathrm{P}_1] + [\mathrm{P}_2] + ... [\mathrm{P}_N]}

which can be rewritten in terms of microscopic pKa values

\mathrm{K_{app}}  = \sum^M_j \frac{1}{\sum^N_i 10^{pK_{ij}} }

The sum contains contains microscopic pKa values for which more than one protonation state is changed. For example, molecule with three ionizable groups $(\mathrm{B_1H^+B_2H^+B_3H^+})$ will have the following microscopic pKa value

\label{eqn:pkapp}
K_{ij} = \frac{ [\mathrm{B_1B_2H^+B_3][H^+]}}{[\mathrm{B_1H^+B_2B_3H^+}]}

in the expression for the apparent pKa value for going from the doubly to the singly protonated state. However, the error in such a pKa value is considerably higher due to less error cancellation and such pKa values are therefore neglected in the sum.

In the Bochevarov et al. implementation the user defines the titrateable group for which the apparent pKa is computed.  However, in my approach all possible protonation states so the assignment of the apparent pKa value to a particular ionizable group is not immediately obvious.  Inspection of Eq $\ref{eqn:pkapp}$ shows that the largest microscopic pKa values will dominate the sum over $N$, while the smallest of these maximum pKa values will dominate the sum over $M$. Thus, the apparent pKa is assigned to the functional group corresponding to the microscopic pKa

\label{eqn:max}
pK^\prime_{ij} = \min_{j}(\{ \max_{i}(\{ pK_{ij} \} ) \} )

In some cases there are several microscopic pKa that contribute almost equally to the respective sums and for these cases the apparent pKa value cannot meaningfully be assigned to a single ionizable sites.  In my approach include all microscopic pKa values that are within 0.6 pKa units of the respective maximum and minimum values defined in Eq $\ref{eqn:max}$.   I choose 0.6 because that is the site-site interaction energy below which two groups titrate simultaneously.

You can find the code I wrote to do this here. It is not cleaned up yet.

## Thursday, May 11, 2017

### Small molecule energy minimisation and conformational search using Tinker

In this post I outlined the (possible) need for conformational analysis using a continuum model. Anders Christensen reminded me that Tinker can do this so I looked in to it. This blogpost and this GitHub repo were particularly helpful.  Here's what I found.

Generating atomic charges
Tinker has MMFF force field parameters except the atomic charges.  These can be generated from an sdf file using sdf2tinkerxyz (which relies on OpenBabel). When you download from SourceForge you get a bin file, which I didn't know what to do with so I downloaded the Linux executable from here instead.

./sdf2tinkerxyz -k default.key < xxx.sdf

creates yyy.xyz (coordinates) and yyy.key (charges), where yyy is the title in the sdf file (which I'll assume = xxx in the following). "default.key" is a file with Tinker keywords that is added to xxx.key. Here's mine

# Solvation Model
SOLVATE GBSA

#Force Field Parameters
PARAMETERS /opt/tinker/tinker/params/mmff.prm

Energy minimisation
To isolate the effect of continuum solvation on the pKa predictions, I want to generate conformers with RDKit and minimise them with MMFF/GBSA using Tinker.  This is done by

/opt/tinker/tinker/bin/minimize xxx.xyz -k xxx.key 0.01 > xxx.tout

"0.01" is the convergence criterium (in kcal/molÅ??). xxx.tout contains the energy information and a xxx.xyz_2 file is generated with the optimized coordinates.  This is in a format that OpenBabel can't handle, so I need to write a converter.  The main challenge is to translate the atom types into names. See the list called lookup in this code.

Conformational Search using SCAN
Tinker also has its own conformational search method. It looks like it's based on eigenvector following but I haven't looked closely at it yet.  This is done by

/opt/tinker/tinker/bin/scan xxx.xyz -k xxx.key 0 10 20 0.00001 > xxx.tout

Here I use the settings used in the DP4 GitHub repo.  "0" is automatic selection of torsional angles, "10" is the number of search directions (modes?), "20" is the energy threshold for local minima,and "0.00001" is the optimisation convergence criterium.  You get a bunch of coordinate files and xxx.tout holds the energy information.

I haven't really played with this option yet. If you have any tips, please leave comment.