Classical Motion of the Nuclei in a Molecule: A Concept Without Alternatives

It is still true that chemical reactions can be explained by a quantum treatment of the electrons and a classical motion of the nuclei. The paper tries to make this finding more plausible by examining elementary ground and excited state reactions using Car-Parrinello and Born-Oppenheimer molecular dynamics. Several different reaction pathways leading to different products were observed for the chlor methane photoreaction, illustrating the classical chaos connected with molecular dynamics on-the-fly simulations. Kasha’s rule saying that all photochemistry occurs in the first excited state is preserved one more time.


Introduction
How to describe a many-particles molecular system in general? Conceptually, everything is still relatively easy as long as we have one atom only. The nucleus is fixed in space as if it were infinitely heavy, and the Schrödinger equation is solved for the surrounding electronic cloud using the electronic Hamilton operatorĤ el : [1] H el ¼T elec þV elecÀ nuc þV elecÀ eleĉ H el Y ¼ E el Ŷ H el is the kinetic energy of the electrons,V elecÀ nuc the electron-nucleus interaction,V elecÀ elec the electron-electron interaction. The electronic wave function is approximated by an anti-symmetrized product of the single orbitals, that is, it is represented by one or several Slater determinants in order to fulfill the Pauli principle and spin symmetry. The Schrödinger equation was introduced ad hoc in order to describe the hydrogen spectrum correctly, the Pauli principle was introduced ad hoc in order to describe multi-electron systems using the Schrödinger equation. While a more basic explanation is missing, this concept works excellently. It is clear that using the electronic Hamilton operator is an approximation (non-relativistic, time-independent, electrostatic interaction only). Nevertheless this approach describes most chemical systems with high accuracy and quantum chemists try to solve it as exactly as possible which led to the development of various quantum chemical ab-initio methods like Hartree-Fock, density functional theory, and many others. This conceptually relatively simple picture is destroyed if several nuclei shall be treated, that is, if we go from an atom to a molecule. What is normally applied here is the Born-Oppenheimer approximation. Theoretically, this approximation means separating electronic and nuclear motion after having made a Produktansatz. In practice, making the Born-Oppenheimer approximation normally means solving the electronic Schrödinger equation self-consistently and doing something about the nuclei later. This is also called the adiabatic approximation. Thanks to this approximation we have a well-defined connection of the electronic wave function with the molecular geometry in a multi-atom system: For a given geometry the wave function is fully optimized while the nuclei are fixed in space. As soon as we have several atoms, the nuclear-nuclear interaction,V nucÀ nuc starts to matter while before it was simply zero. It turns out thatV nucÀ nuc must be added in order to obtain what is called Born-Oppenheimer surfaces or potential energy surfaces by optimizing the wave function for given nuclear geometries. [2] H QC ¼T elec þV nucÀ nuc þV elecÀ nuc þV elecÀ elec WithoutV nucÀ nuc a purely attractive curve is obtained instead of a Morse potential which is repulsive at short distances. It goes without saying that the nuclear-nuclear interaction is easily computed under the assumption that the nuclei are classical point charges and this is done in every quantum chemical calculation for anything that is larger than a single atom. Note that while adding the nuclear-nuclear interaction term is numerically essential, it is not so clear if adding the vibrational zero-point energy to potential energy surfaces is an advancement. Today, most potential energy surfaces are computed without vibrational zero-point energy.
What aboutT nuc ? It is often argued that this term is small as the nuclei are heavy and hence slow. This argument is not helpful, because the kinetic energy raises linearly with mass. Indeed, if one computes the kinetic energy classically, it turns out that this term cannot be neglected at room temperature, i. e. ∼ 300 K. In this sense, traditional quantum chemical calculations can be viewed as being performed at zero Kelvin. For ab-initio molecular dynamics simulations (AIMD) using Newton dynamics for the nuclei,T nuc must be added in order to obtain conservation of energy during an MD run.
H AIMD ¼T nuc þT elec þV nucÀ nuc þV elecÀ nuc þV elecÀ elec Like this we arrive at an expression which contains all kinetic energy terms and all potential energy terms. No term is left for solving the nuclear Schrödinger equation. This happened because we made again and again the approximation that nuclear motion is classical. By making the singleconfiguration approximation and forming functional derivatives, Hartree-Fock or Kohn-Sham equations are obtained depending on the approach used for the electron-electron interaction. In this framework it is easy to compute the forces on the nuclei (see for example Ref. [3]): The alternative to treating the nuclei classically is obviously path integral molecular dynamics [4][5][6][7][8][9][10][11][12] , for a review focussing on the application to proton transfer and proton wires see Ref. [13]. We will compare both approaches in this study. In accordance with the literature, we will find that there can be significant differences between standard molecular dynamics and path integral molecular dynamics, for example when investigating rare events such as chemical reactions. However, we will interpret this difference as classical chaos and we will show that a deterministic interpretation gives a more convincing picture particularly of dissociation processes.

Results and Discussion
Let us have a look at a single hydrogen atom ( Figure 1). What do we expect? Well, the single atom should be stable and should not change with time as long as it is alone in empty space. Performing an ab-initio molecular dynamics simulation yields exactly this result. Does the single hydrogen atom have a zero-point energy? The zero-point energy is the difference between the bottom of the potential and the ground state. This difference is infinity for the bound electron in its 1/r potential. For the moment we have to live with this. A yet unknown, more complete theory must yield a finite size for the nucleus which prevents an infinite zero-point energy.
What happens if two hydrogen atoms are considered? They move linearly till hitting the other hydrogen atom? Not necessarily, if the atoms are slow, as Figure 1 shows: What is observed is a strongly inelastic interaction leading to the formation of a hydrogen molecule. The hydrogen atoms move towards each other like boats in a giant whirl.
As we use a program code which treats the nuclear motion classically, we can be certain that this observation is free of nuclear quantum effects. A quantum mechanical description of the electronic cloud is sufficient.
What changes if we use path integral molecular dynamics (PIMD)? The results are summarized in Figs. 2 and 3. Figure 2 b) illustrates what happens if the total initial rotational energy and the center of mass energy are put to zero. Now the hydrogen atoms approach linearly. In path integral calculations these quantities must be put to zero anyhow. Nevertheless we get a spiral motion (Figure 2 c)). This picture is more strongly smeared out if the number of Trotter dimensions is increased (Figure 2 d)). If the kinetic energy is increased (right column), dissociation can be observed. After such a dissociation event, the atoms may encounter each other in the next unit cell and form a molecule again. To conclude, the pattern observed in normal CPMD simulations and PIMD simulations is similar. However, the PIMD calculations have to be interpreted differently: We observe a motion in imaginary time. The stroboscopic picture is a static pattern and shows the structure of the nuclei. Particularly in the dissociation reactions it is obvious that the result obtained is wrong. No molecule looks like that.
What changes if the system is deuterated? The result is illustrated in Figure 3. The two deuterium atoms react slowlier than the hydrogen atoms (see also the Supplementary Material). In one case even a different reaction pattern is observed. A vibrational zero-point energy is not needed to get a different behaviour for hydrogen and deuterium.
Let us consider more complex problems. For example one might ask how things look if all this would occur in a solvent like water. Actually that is the situation encountered for the Figure 1. a) As trivial as essential: the lonely hydrogen atom surrounded by its electron density is well approximated by AIMD. Since we describe nuclear motion classically, there is no problem with the nucleus being at rest and perfectly localized. In order to make it possible to see the nucleus, it is shown several orders of magnitude bigger than in reality. On the Angstrom scale it is simply a point particle. b) The same approach for two atoms: Also here the nuclei are perfectly localized at any point of time. The electronic cloud causes bond formation. c) Like b), now monitoring the spiral motion of the nuclei by plotting the positions of the nuclei during the whole molecular dynamics run. A movie of the simulation, showing both nuclear and electronic motion, is deposited in the Supplementary Material. electrolysis of water. At the cathode protons are neutralized leading to neutral hydrogen atoms. These can easily dimerize even at a distance by a rearrangement of hydrogen bonds. [14] The mechanism is similar to the Grotthus mechanism (or proton wire mechanism), except that we do not deal with protons but with neutral hydrogen atoms. Also this switching of hydrogen bonds is perfectly obtained with a quantum mechanical description of the electronic cloud and a classical description of the nuclear motion. It is characteristic for the mechanism that the nuclei are not moving much, while the electronic cloud which determines the bond pattern, is  rearranging (see Figs. 7 and 8 in Ref. [14]). Note that such simulations work only if the solvent is described quantum mechanically just like the solute, atom by atom and molecule by molecule, without the application of an empirical solvent model. Density functional theory as a single-configuration method is doing a great job here: Since it is size consistent, the isolated molecules and the molecules in solution can be directly compared. For an application to proton wires see for example Ref. [15]. It is these huge simulations which lead one to the assumption that a classical picture for the motion of the nuclei is fully sufficient. Let us have a look at the opposite reaction direction, namely bond cleavage. Absorption of a photon often leads to dissociation as the orbitals which are occupied upon electronic excitation may have antibonding character.We are using restricted open-shell Kohn-Sham theory (ROKS) [16][17][18][19] together with Born-Oppenheimer molecular dynamics (BOMD) to describe electronically excited systems by forcing two electrons to be unpaired. The ROKS method is implemented in the CPMD code together with appropriate SCF convergence procedures to facilitate the simulation of low-spin excited states. (The highspin excited states are trivial.) The main drawback in the actual implementation is that GGA functionals only can be used, which leads to a red shift compared to experimental values. As a consequence, the systems have too low kinetic energies when leaving the Franck-Condon region. B3LYP should be able to heal this as it is a mixture of BLYP and Hartree-Fock with the latter method yielding too high values for the excitation energy. Due to the high cost of hybrid functionals when combined with a plane wave basis set, BLYP will stay for some more time the state-of-the-art approach for performing excited state on-the-fly simulations.
Let us start with rotationally excited H 2 and HCl. The stroboscopic picture (Figure 4) shows for H 2 initially the rotation of the two hydrogen atoms. In the very moment the system is put to the excited state, they dissociate and move linearly till, due to the periodic boundary conditions, they encounter again and collide essentially elastically. The velocity which the system gains while leaving the Franck-Cordon region, corresponds to a temperature above 6000 K which prevents immediate recombination. Figure 4 b) shows the same type of reaction for HCl. The rotation of the H atom around the heavier Cl atom is observed. Upon excitation the direction changes immediately and the hydrogen atom is flying away while the remaining chlorine atom hardly moves.
This is all not unexpected if one has accepted the fact that nuclear motion is classical.
A more interesting case is chlormethane, CH 3 Cl, as several different reaction pathways are observed leading to different products, depending on the excitation wave lengths. Figure 5 shows the potential energy curves for this molecule.
With our normal ROKS procedure (equilibration in the RKS ground state, vertical excitation to the excited singlet state, propagation of the excited state till it crosses or touches the ground state) we cannot describe the dissociation into ions unless the ionic state is stabilized by a solvent and is hence the lowest lying state. [20] However, we are able to describe the experimentally relevant cleavage of a CÀ H bond alternatively to breaking of the CÀ Cl bond in CH 3 Cl. [21][22][23] The latter predominates only at lower kinetic energies. At higher kinetic energies the CÀ H cleavage becomes relevant. In addition, at intermediate kinetic energies we observe HCl formation, see Figure 6 and the Supplementary Material.
All reactions are ultrafast. The dissociations take less than 0.02 ps, that is, they occur immediately from the Franck-Condon region. Obviously, it is possible to save Kasha's rule [24] which states that all photochemistry occurs in the first excited state as all other states are too short-lived for influencing ionic motion. Getting different excited state reaction channels by raising the temperature was observed before. [25] Using very high temperatures is allowed in photoreactions as the high absorption energy of roughly 1 to 10 eV has to be converted into kinetic energy unless a photon is ejected. As 1 eV corresponds to 11600 K, small molecules with  few degrees of freedom can get quite hot. In a liquid or solid environment the kinetic energy is quickly redistributed.
One remark may be made at this point, namely that the observation of different reaction channels is not necessarily always decided by the motion through conical intersections but rather by the motion in the Franck-Condon region. When the system reaches the conical intersection for the hydrogen abstraction or that of the chlorine abstraction, its fate is already decided: either hydrogen or chlorine abstraction. Nevertheless, conical intersections are important points on the potential energy surfaces and there may be several conical intersections close to the Franck-Condon region responsible for the efficient deactivation to the first excited state.

Conclusions
To conclude, in several reactive molecular dynamics simulations it was shown that the classical picture for the motion of the nuclei is sufficient to explain the observed phenomena. In principle this is in line with the existing literature on pathintegral MD applications. [26][27][28] There is no need to apply the Schrödinger equation to the nuclear motion. This does not attack the idea of an universal quantum theory. For bound electrons the Schrödinger equation is practically perfect since this is what it was developed for. For a more general treatment, one might want to start a derivation from something more similar to the wave equation or the Klein-Gordon equation. [29] Chemical reactions might not be the right field of science to investigate the inner structure of the nuclei which matters at completely different length scales (femtometer instead of Angstrom). However, the theoretical investigation of chemical reactions allows one to exclude erroneous views of quantum mechanics. The electron cloud moves continuously from educt to product without any collapse. Both the electronic wave function and the nuclei move according to their differential equations, there is no differential equation known which would describe a measurement or a collapse. Any differential equation which leads to diffusion of the complete system can be excluded as a potential candidate for describing the electrons, the nuclei, or both. For the motion of the nuclei, only the classical Newton dynamics gives a reasonable picture. The connection to experimental observation is made by a Maxwell-Boltzmann velocity distribution in combination with classical chaos as observed in molecular dynamics simulations.
Finally, the stroboscopic pictures lead to an idea how one could generate a molecular motor, namely by covering one side of each rotor blade in a propeller with a compound that upon irradiation undergoes dissociation. Like this it would be possible to convert light energy into directed motion in a very straight-forward way. One might go down to nano propellors with just one or a few dissociation sites.

Supporting Information Summary Paragraph
The Supplementary Material contains movies of the reactions, a comparison to path integral molecular dynamics, a comparison to the CD 3 Cl reaction and a more detailed description of the chlormethane photoreaction.
In short, the simulations were performed using the Car-Parrinello molecular dynamics (CPMD) code [3,30,31] with the BLYP functional. [32,33] Troullier-Martins pseudopotentials were employed for describing the core electrons. [34,35] Path integral molecular dynamics simulations have been carried out using the CPMD program package. [4,7] For the excited state ROKS simulations [16][17][18][19] Born-Oppenheimer molecular dynamics was used in connection with the modified Goedecker algorithm. [36] Figure 6. Complex motion, easy to understand: the three main reaction channels for thechlormethane photodissociation: a) Dissociation of Cl, both moieties are moving away from each other. b) Dissociation of HCl, the resulting HCl molecule is spiraling along the dissociation coordinate. c) Dissociation of H, the remaining CH 2 Cl radical stays essentially at rest. (Color code: white: hydrogen, black: carbon, red: chlorine). Orbital movies of the HCl formation are deposited in the Supplementary Material.