## Thursday, March 31, 2016

### Reviews for Towards a barrier height benchmark set for biologically relevant systems

2016.04.06: our rebuttal can be found here

The reviews our latest PeerJ submission (editor assigned February 25) came back last evening with the alway welcome verdict "minor revision".  The preprint has already received over 400 views and 85 downloads.

## Sunday, March 27, 2016

### Generating protonation states and conformations

We're working on pKa prediction using semiempirical methods and need to generate coordinates for all possible protonation states of a large number of amines and perform conformer search.  Here are my current thoughts so far.

We have the amine names so we can use cactus to generate SMILES strings and then Babel to generate coordinates and do the conformer search.  However, Babel only (1) generates neutral protonation states, (2) does not consider C-NH2 bonds in primary amines rotateable, (3) does not consider conformers relates to inversion of N.  Babel also has some additional problems that I'll describe at the end of the post.

Some easy fixes
Most of the problems can be addressed fairly easily. It was easy to write a program that reads a mol2 file (where atom typing is mostly done) and protonates all appropriate N atoms.

Similarly, it also seems easy to write a program that generates all possible protonation states by removing select protons from the fully protonated coordinate file.  Furthermore, by removing all protons on primary and secondary amines one at a time one will help solve the problems with C-NH2 rotation in primary amines and the inversion of secondary amines.

The main problem is the inversion of tertiary amines
Actually, it is not heard to generate a crappy guess at the inverted structure simply by placing the H atom on the opposite side of the tertiary N.  When do this manually in Avogadro by moving the H atom and optimising the structure it works fine.  But generating such coordinates in an xyz file and then passing them to Babel doesn't work because the atom typing is screwed up.  I tried generating mol2 files instead but for some reason Babel perform a conformer search starting from a mol2 file.

Another option could be to generate a z-matrix (gzmat) and then write code that does the inversion by changing the sign of select dihedral angles.

Yet another option could be to do it directly from the SMILES string, e.g. C[N@H+](F)Cl vs C[N@@H+](F)Cl.  The challenge here for me would be to only do this for select N atoms. In general writing a SMILES parser seems like a big job.

I am open to suggestions here

Babel does not generate the correct stereoisomer for molecules with chiral centers involving several rings.  It does print out an error message, but it remains to be seen how many of the amines we are interesting in are affected by this.

Babel also doesn't consider C-OH bonds rotateable so we have to somehow generate the start conformers separately.  This also applies to COOH groups.

Babel generates -COOH groups, so we also have to generate the COO- form.

Other programs
Perhaps we should look at CDK, RDkit, or PubChemPy instead.

## Thursday, March 17, 2016

Here are some very useful spreadsheet commands that I have been using a lot lately.  The syntax is for Google Sheets.  To use in Excel replace the ","s by ":"s

SUMIFS Returns a sum based on one or more criteria (if you only have one criteria you can use SUMIF).  Say you have

A C 2
B C 1
A D 4

You can use SUMIFS to add the "A-numbers", i.e. 6.  If you want the sum for  "A C" then use SUMIFS.

One of my side jobs is to keep track of teaching hours reported using a Google Form.  I use SUMIF to add all the reported hours together for a given person.

VLOOKUP Searches down the first column of a range for a key and returns the value of a specified cell in the row found.

You can use VLOOKUP to find the value corresponding to "B", i.e. 1.  If you use the same command to find the A value it will return the first one it found, i.e. 2.  VLOOKUP does not support multiple criteria, but if you want to look up the value for "A and D" you can use CONCATENATE to combine the two cells and lookup the value for "AD".

A more concrete example: let's say you have a long list of conformer names and energies.  You can use MIN to find the lowest energy and combine it with VLOOKUP to give you the conformer name.

CONTATENATE I've already mentioned this command, but here's another example.  In some of the papers I make tables with RMSD and $r$ values in the same column, e.g.: 0.2 (0.92).  This is easy to make from separate tables with RMSD and $r$ values with  CONCATENATE(TEXT(I4,"0.0")," (",W4,")")

SPLIT is the opposite of CONTATENATE.  SPLIT can be combined with INDEX to do some pretty nifty things.  For example INDEX(SPLIT(INDEX(SPLIT(B5,"("),0,2),"%"),0,1) will extract the number 92 from the text "B8FX10 Buried chi1: 23 correct out of 25 (92%)"

Finally the IF command is incredibly useful.  IF is actually an "if-then-else" statement.  For example IF(B15>3,"the number is bigger than 3","the number is smaller than 3") returns "the number is smaller than 3" if B15 contains the number 2.

Notice that IF commands can be combined: IF(B15>3,"the number is bigger than 3",if(B15>1,"the number is smaller than 3 but bigger than 1","the number is smaller than 1"))returns "the number is smaller than 3 but bigger than 1" B15 contains the number 2.