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