⇐ Tutorial 4 Kimball Tutorial 5 Tutorial 6 ⇒

Content: How do we get structure in chemistry (first part)

ES 26 Sept 2017

This might be an eye opener for inquisitive minds who are not satisfied with "of course" or "stupid !" answers.

It takes some time, we hope it is well spent!

1) Suppose, you had a good Quantum Chemistry Program, advertised for "ab initio" results, e.g. GAMESS-US. You expect it to predict the stable atoms, molecules, and crystals of everyday matter by just feeding it with the atomic nuclei of interest in arbitrary positions of space for a start, and the necessary number of electrons. An easy case: Select e.g. one C(+6) and 4 H(+) nuclei, add 10 electrons for electrical neutrality, and trap these 15 particles in a cage of 10x10x10 bohr^3. Let the program start with random locations of the nuclei and trust it with automatically handling the electrons. Here is the input to the program:

$CONTRL SCFTYP=RHF RUNTYP=optimize COORD=unique

MULT=1 ICHARG=0 $END not a free radical but a neutral molecule expected

$BASIS GBASIS=PM3 $END this defines the electron handling scheme, a MOPAC hamiltonian PM3

$statpt nstep=400 $end you give the program 400 optimization steps

$DATA

CH4 verzerrt

C1

C 6.0 0.0 0.0 0.0 these are the nuclei, their charge and x,y,z coordinates in Å

H 1.0 2.1 0.0 0.0

H 1.0 -0.6 3.1 -1.7

H 1.0 1.0 2.0 3.0

H 1.0 -4.0 2.3 0.7

$END

This read, it deploys its arsenal of interactions, looking for the most stable arrangement of the given atoms. You might expect a stable, well constructed CH4 molecule being produced with the blackboard chemistry

C(+6) + 4 H(+) + 10 e- → [ C + 4H ] or other intermediates → CH4 + energy

and properties near to those known from spectroscopic observation. But, perhaps, two molecules are coming out, CH2 and H2, who knows? or CH(-) and H3(+), C and H2 and 2 H atoms? or many more possibilities, such as CH + H + H + H ... CH2 and CH are stable molecules of which billions of tons exist - in interstellar space - and H2 you know; CH(-) (discovered 1976) and H3(+) (well known to any mass spectroscopist). Since earth side chemistry is of more immediate interest, we hope for CH4, a stable molecule in our environment, whereas CH2, CH, CH(-), C, H3(+), H, and even H2, will be involved in several chemical reactions in our atmosphere and vanish, eventually, as a consequence of them, to be recreated e.g. as CO2, H2O.

2) But, nothing of this sort is happening:

There is a fundamental problem, about which the program does not tell you anything. Connected to the expected process of formation of CH4 the selected particles have to get rid of quite a large chunk of energy, which the program would be telling you if it had found a solution. How can the particles do this? We help by giving the cage, within which we have enclosed the particles, impenetrable but translucent walls. Thus, the process of atom and molecule formation may start by radiating light energy into the universe as a vast sink. There are some restrictions to this, and a conundrum with spin states of collision pairs as well, but let's not fret about those. Neither does the program which continues unimpressed (you had given the expected spin multiplicity MULT=1 and ICHARG=0 in the input, meaning the expected CH4 to be electrically neutral and having only up-down spin pairs):

Almost in no time, with a modern desk computer, you get the sad answer:

OPTIMIZATION ABORTED.

-- GRADIENT OUT OF RANGE

-- MAXIMUM ALLOWED FORCE (FMAXT) = 10.000

EXECUTION OF GAMESS TERMINATED -ABNORMALLY- AT Tue Feb 03 21:07:55 2015

580000 WORDS OF DYNAMIC MEMORY USED

STEP CPU TIME = 0.04 TOTAL CPU TIME = 1.6 ( 0.0 MIN)

TOTAL WALL CLOCK TIME= 1.6 SECONDS, CPU UTILIZATION IS 100.00%

A bit thunderstruck, eh? Repeat the computation, using a better electron handling Basis-Set, e.g. 6-31G, with the rest of the input the same:

$CONTRL SCFTYP=RHF RUNTYP=optimize COORD=unique

MULT=1 ICHARG=0 $END

$BASIS GBASIS=N31 NGAUSS=6 $END

$statpt nstep=400 $end

$DATA

CH4 verzerrt

C1

C 6.0 0.0 0.0 0.0

H 1.0 2.1 0.0 0.0

H 1.0 -0.6 3.1 -1.7

H 1.0 1.0 2.0 3.0

H 1.0 -4.0 2.3 0.7

$END

After 132 optimization steps in 1.6 sec the search is converged and Gamess produces stable CH2 and H2 molecules, well separated in space, not CH4! I told you so! (BtW: Gaussian03™ aborts ungraciously with the given coordinates and the same model chemistry without a result!).

You react again, and offer the program starting coordinates of the nuclei with successively better coordinates of a CH4 molecule, e.g.: in Å

B) C) D) Result coordinates of D) Td not assumed! C-H distance

C 6 0.0 0.0 0.0 C 6 0.0 0.0 0.0 C 6 0.0 0.0 0.0 C 6.0 0.020074 0.059996 -0.239989

H 1 2.1 0.0 0.0 H 1 2.1 0.0 0.0 H 1 1.1 1.0 1.0 H 1.0 0.594972 0.630022 0.478015 1.0821125

H 1 -0.6 1.5 -1.7 H 1 -0.6 1.5 -1.7 H 1 -1.0 1.5 -1.7 H 1.0 -0.535935 0.736553 -0.875715 1.0821382

H 1 1.0 2.0 1.2 H 1 1.0 -0.9 1.2 H 1 1.0 -0.9 -1.2 H 1.0 0.689186 -0.536438 -0.846170 1.0820820

H 1 -1.0 2.3 0.7 H 1 -1.0 1.3 0.7 H 1 -1.0 -1.3 0.7 H 1.0 -0.668297 -0.590132 0.283859 1.0820996

Put the protons nearer to the expected tetrahedral positions around a central C nucleus, C-H distances about 0.9 to 1.5 Å. The coordinates in B) and C) still lead to CH2 + H2. When you arrange the nuclei quite near to the observed structure in D), the computer runs gracefully to an end, finds CH4, and proudly announces: " ***** EQUILIBRIUM GEOMETRY LOCATED *****! This is an occasion to celebrate" (this last remark is uttered by the outstanding program "Turbomole" of R. Ahlrichs, Karlsruhe).

It means that the program has found a stable molecule within small error limits for total energy, bond lengths and symmetry of CH4. The energy of formation can be extracted from the quantities determined, if you give it data of the standard states of Carbon and Hydrogen elements, although the program does not know how it got rid of that. None of these QC programs is a simulation of the physical processes by which molecules are formed.

From here on, you may use the converged results for several levels of sophistication, taking more and more computer time with better and better approximations. Finally, a very good structure is definitely found, and the vibrational/rotational spectrum, observed in the Infrared and microwave region, and by Raman scattering, is quite accurately simulated. The enthalpy of formation is within less than 4 kJ/mol of the experiment (using the G3MP2 sequence of optimizations). The quantum chemical code may even produce useful UV-visible spectra of electronic excitations, and NMR chemical shifts and transition frequencies/intensities between spin states.

3) Is this prediction? Is this "ab initio"? Not by any margin, in regard to what we set out to do in 1)! Unless you know the answer pretty well, you don't get an answer! So, the program is much more "a posteriori", good at refining what you already know (and explaining this, if you are an expert at interpreting sheets of numbers), but mostly unsuitable for predicting what you don't know yet. Linus Pauling once said: "Usually, I am not very much interested in calculations because I know what the answers are going to be anyway!". He was an exceptionally gifted, intuitive chemist. Meanwhile, we expect much more guidance from theory than available at his time, around 1939 when he wrote his "Nature of the Chemical Bond".

You can easily reproduce the above runs with GAMESS-US, since this program is freely downloadable for non commercial applications and runs perfectly under Windows, using all the cores your CPU provides (mine has 4). Or, you can directly run your own (deformed) methane computation free on a WebMO demo server in USA (just DO it! you learn it by doing) without installing anything!

Let's try to find out what's amiss, or is the expectation in 1) naive?

Yes, in several ways: This problem would have been unsolvable already in “classical” physics. The equations of motion for three “bodies” with known forces on each other, e.g. gravitation, cannot be solved exactly. That was well known to the astronomers of the 18th century. Our approach with the 15 particles of building CH4 would need the solution of a system of 15*(15-1)/2 = 105 hamiltonian differential equations for some million of timesteps of femtoseconds, giving just the first billionth of a second of their trajectories. In quantum mechanics we would only move the nuclei and describe the electrons by Schrödingers equation, and the Pauli exclusion principle. The forces are the wellknown electrostatic interactions between charged particles, and new quantum mechanical forces connected to the change of the kinetic energy of the electrons. It is not easier! So, we actually leave the nuclei at rest in a plausible arrangement based on experience and try to find the distribution of electrons compatible with this, using a much simplified “Schrödinger” equation. Then we change the nuclear arrangement in small steps, get the electrons adjusted to them, and repeat this for many, perhaps a few 100 steps, to find a (local) minimum of energy. If we do this carefully, we might get millions of answers to our first question - meaning, there is a large number of local energy minima representing (slightly) different "chemistry". In fact, already one H atom has infinitly many of those. It is almost impossible for any program to find the one we hope for, if we do not reduce the search space for its discovery. It is not predictable where the global minimum energy system may be located. There is no direct way to find this. The program might get stuck in a local minimum, like CH2 + H2 or any other of the varieties presented in 1), and not find its way out to the lowest minimum energy CH4 molecule. That is the reason, quantum chemical program packages usually start from a reasonably good (experimental) structure to determine the energy and other properties of the target molecule.

We want to use the model of this site to see a bit clearer what is involved. We construct a very simple system:

We put a number of electrons into one sphere with radius R and 2 to 8 protons nearby (we enter nuclei of different charges in the next tutorial). Then we let the interactions begin and look for a stationary energy minimum, if it exists. What “shape” do the nuclei form? The following program can be used for this exercise. Only n, the number of protons, needs to be changed for the series to be computed, the statement in red, see below.

This is an occasion when you should use Mathematica for actually doing the manipulations described. Since this is a very expensive program, which might not be available at your desk, you can download a free trial version for 15 days. For a cheap permanent install buy a Raspberry Pi 2 computer for $ 35 and use Mathematica, Version 10, installed on the freely available Raspbian operating system. With the exception of very large 3D-graphics the “notebooks” on this site will run perfectly on the Raspberry Pi (1B to 3).

Input and Definitions

The large orange sphere is the equilibrium Kimball free cloud. The dark spheres are the final equilibrium positions of the protons, while the small red dots are their initial coordinates. The stereo image which you can see by squinting, shows the 8 protons at the vertices of a regular square antiprism, i.e.on a cubus whose bottom face has been turned 45° from the top.

Substitute n=4,5,6,7 to produce the following final results:

We summarize:

n=4 gives a regular tetrahedron

n=5 a trigonal bipyramid

n=6 a regular octahedron

n=7 a pentagonal bipyramid, and

n=8 a regular square antiprism, and more:

n=10 an axially doubly capped square antiprism

n=12 an icosahedron

n=13 an icosahedron with center, and still more, see below.

If you run the program for these different n (and you might include n=2,3, yielding a stick and an equilateral triangle), you can check on these structures by looking at the proton distances from the origin which is printed below the results of the FindMinimum function call.

This exercise shows, what kind of shapes are produced by the attraction of the protons by the electron sphere, and repulsion of each proton from all the others.

Some colleagues may remember that several of these structures are observed with ligand free Na-metal clusters Na_{n}. No wonder: n Na+ in a drop of n 3s electrons may have similar properties to those computed for n H+ in a drop of n electrons, above. However, for n < 6 sodium clusters prefer flat structures, not predicted by electrostatics alone.

Except for metal clusters, this part of the tutorial may chemically be somewhat a "dry run". We seldom have more than two nuclei in an electron sphere (e.g. H3(+)). As soon as we have to deal with nuclei of different charges the situation may change dramatically. Already one higher charged, different nucleus can destroy the symmetry of the proton locations. It tends to go nearer to the negative center and thus displaces the protons. If it sits in the center, however, the shape does not change much: It just reduces the negative charge of the cloud and we have the symmetry of the polyhedron with one nucleus less. However, such a structure is unstable for small n, unless we divide the electron density into 2-electron domains in order to obey the Pauli exclusion principle. See next tutorial.

The remarkable result here is, that all the structures produced above by electrostatic interactions only (in a charge sphere spanned by its zero point kinetic energy), are found and known in molecular and coordination chemistry, see the interesting series from H,H2,...,H24,H25.