Rather than make a post that is six miles long with updates, new posts with new content seem like a good idea. For most people this would be an easy idea!

I’ve continued to learn with GAMESS. The functionality that interests me the most is the capacity to use the background of automated quantum mechanical calculations to simulate the occurrence of chemical reactions. A constellation of eigenstates containing electrons as shaped by the positioning of atomic nuclei can be sculpted by the constraints of energy minimization to reveal which configurations are more likely to occur than others. You can sort of think about nuclei as tiny pearls trapped in a continuum of cotton candy: if you grab a pearl and try to move it, the candy clings to it like springy rubber to try to dampen the movement. The only trick is that the candy is not distributed in a uniform fashion and is actually redistributed by the movement of these pearls; it tends to spring and flex in preferred directions. As with most physics, this becomes essentially a system of masses on springs where you can calculate how the deflections occur in a potential energy “surface” that is a function of the positions of the nuclei. I put “surface” in quotes because it’s a hyper-surface with many more than simply two dimensions: it has three dimensions with every nucleus.

Figuring out a chemical reaction is very much about intuiting the shape of this hyper-surface from clues in the 3D structure of the atoms. You use the potential energy surface to generate a mathematical entity called a “hessian,” which is a 3Nx3N matrix formed from a function of the potential energy surface with second-derivatives taken with respect to pairs of the 3N coordinates. Derivatives of energies are forces and the hessian gives what are called force constants (as I understand it) of the potential energy surface and enable the identification of minima and maxima by use of inflection (second derivative!). Eigenvalues of the hessian matrix reveal vibration modes within the molecule and, if you find imaginary vibrations, you have stumbled over transits in the free energy surface between one configuration and another… literally non-vibratory motions inside the model corresponding to the occurrence of chemical reactions.

To find a chemical reaction with GAMESS, the first task is to find features in the potential energy surface that are called “saddle points.” These saddle points are locations in the space of the molecular coordinates where the vibratory modes of the hessian are all minimized while the one imaginary mode is maximized (again, as I understand it… these constructs are all such huge numbers of variables that you have the computer sitting between you and the calculation and I have not done the math myself directly). Saddle points can seem very arcane described in that way, but they mark transitory states in the geometric structure: literally the intermediate structural state of the chemicals positioned between the chemical products and reactants called a transition state. If you can guess accurately at the form of a transition state, the computer can use energy minimization to massage that configuration into the stationary configuration -the flat points along the potential energy surface when the derivative goes to zero- where the imaginary vibratory mode is maximized and the other modes are minimized in potential energy.

If you can find a saddle point, which turns out to be very tricky, the computer is then able to follow the curves of quickest descent in the potential energy surface forward to the reaction products and backward to the reactants. You can then stitch these half-paths together to produce the whole chemical reaction.

Here is a chemical reaction…

This is the entire reaction pathway between a water molecule and a molecule of singly protonated phosphate. The ball-and-stick model above is the usual modeling produced by the GAMESS satellite program wxMacMolPlt: oxygens are red balls, phosphorous is black and hydrogen is white. The wireframe surface is an electron equipotential surface and its coloration is for regions of negative electrostatic potential (blue) and positive potential (red). Bonds are drawn based on electron occupancy of orbitals in that region, with dashed lines depicting essentially hydrogen bonds, single sticks depicting single covalent bonds (~two electrons) and double covalent bonds (~four electrons). Orbitals are closed shells using restricted Hartree-Fock and the basis set was relatively small, only 3-21G. There is a polarized continuum model in the simulation in order to pretend to some extent that this takes place in a medium like continuous water.

In this reaction, phosphate with a charge of -2 sucks a proton (charge +1) off of water, creating a transition state of hydroxide (charge -1) which then sucks a proton back off of the phosphate to recover water.

Since the water is turned away from the camera, here’s another set of images of the same thing:

In performing this set of calculations, which took a ridiculous amount time, my initial objective was to try to find a proton exchange. Happily, I finally found that, though not without 5 iterations of the work. Even in this version, I was backward in my guess of the transition state: I was guessing the transition would be hydronium, with both protons present on the water… it ended up with both on the phosphate since the phosphorous apparently doesn’t withdraw electrons strongly enough from the oxygens to make unbonded water better at holding protons.

I found this in multiple steps. First, I modeled water with triply deprotonated phosphate and looked for an energy minimum. This ended up with water doubly hydrogen bonded to the phosphate (looked like a guy in stirrups). I then plunked a third proton onto the water, creating hydronium (+1 charge). Since this ended up with the third proton added in a funny direction, I used part of the saddle point search as if I was doing a regular optimization to get to the proper shape of hydronium. Finally, as I was looking for the actual saddle point, phosphate dragged the hydronium in by a pincer move and ripped off both of his arms.

This sort of surprised me since I wasn’t expecting such violence, but the transition state ended up with both exchanged protons on the phosphate and only one remaining in the water as hydroxide.

Following the intrinsic reaction coordinate using the imaginary hessian mode allowed me to slide along the potential energy surface to give the reaction depicted above. You can see the bonds shrink and swell as the electrons are redistributed in the molecule to accept the exchange with the water, electron density pulling out away from the phosphorous as the proton is accepted and filling back in on the other arm as the proton is reciprocally donated.

I have a strong feeling that this is generally similar to the state of acids and bases in water, where protons are passed around like bad kleenex along lines of hydrogen bonding. Likely, acidic protons are never free and hydronium and hydroxide constitute transition states. Since water has an overall concentration of about 55 molar, and aqueous strong acids and bases have concentrations of only about 12 molar, plenty of intact water exists even in strong acid or strong base where protons are constantly moving around. In base, the proton is scarce, so the transition state is probably hydroxide, while acid has the proton as common, giving hydronium as the transition state instead. It’s an interesting speculation.