Molecular and chemical dynamics describes the motions of the atoms within the molecule and the surrounding solvent The collisions among molecules and resulting energy transfers among translational, vibrational, rotational, and electronic modes, as well as chemical reactions that occur intramolecularly or in bi-molecular encounters lie within the realm of molecular and chemical dynamics theory. The vibrational and rotational motions that a molecule's nuclei undergo on any one of the potential energy surfaces EkQ) is also a subject for molecular dynamics and provides a logical bridge to the subject of molecular vibration-rotation spectroscopy. 1. Classical Newtonian Dynamics Can Often Be Used For any particular molecule with its electrons occupying one of its particular electronic states, the atomic centers (i.e., the N nuclei) undergo translational, rotational, and vibrational movements. The translational and rotational motions do not experience any forces and are thus "free motions" unless (1) surrounding solvent or lattice species are present or (2) an external electric or magnetic field is applied. Either of the latter influences will cause the translations and rotations to experience potential energies that depend on the location of the molecule's center of mass and the molecule's orientation in space, respectively. In contrast, the vibrational coordinates of a molecule experience forces that result from the dependence of the electronic state energy Ek on the internal coordinates Fi = - dEk({Q})/dQi. By expressing the kinetic energy T for internal vibrational motions and the corresponding potential energy V= Ek({Q}) in terms of 3N-6 internal coordinates {Qk}, classical equations of motion based on the so-called Lagrangian L = T - V: d/dt d[T-V]/di = d[T-V]/dQi can be developed. If surrounding solvent species or external fields are present, equations of motion can also be developed for the three center of mass coordinates R and the three orientational coordinates W . To do so one must express the molecule-solvent intermolecular potential energy in terms of R and W, and the translational and orientational kinetic energies must also be written in terms of the time rates of change of R and W. The equations of motion discussed above are classical. Most molecular dynamics theory and simulations are performed in this manner. When light atoms such as H, D, or He appear, it is often essential to treat the equations of motion that describe their motions (vibrations as well as translations and rotations in the presence of solvent or lattice surroundings) using the Schrödinger equation instead. This is a much more difficult task. a. A collision between an atom and a diatomic molecule Let us consider an example to illustrate the classical and quantum treatments of dynamics. In particular, consider the collision of an atom A with a diatomic molecule BC with all three atoms constrained to lie in a plane. I. The coordinates in which the kinetic energy has no cross terms Here the three atoms have a total of 2N = 6 coordinates because they are constrained to lie in the X,Y plane. The center of mass of the three atoms x = (mA xA + mB xB + mC xC )/M y = (mA yA + mB yB + mC yC )/M requires two coordinates to specify, which leaves 2N-2 = 4 coordinates to describe the internal and overall orientational arrangement of the three atoms. It is common (the reason will be explained below) in such triatomic systems to choose the following specific set of internal coordinates: 1. r, the distance between two of the atoms (usually the two that are bound in the A + BC collision); 2. R, the distance of the third atom (A) from the center of mass of the first two atoms (BC); 3. a and b , two angles between the r and R vectors and the laboratory-fixed X axis, respectively. II. The Hamiltonian or total energy In terms of these coordinates, the total energy (i.e., the Hamiltonian) can be expressed as follows: H = 1/2 m' 2 + 1/2 m 2 + 1/2 m'R22 + 1/.2 m r22 + 1/2 M (2 + 2 ) + V. Here, m is the reduced mass of the BC molecule (m = mBmC/(mB +mc )) M is the mass of the entire molecule M = mA + mB + mC , and m' is the reduced mass of A relative to BC (m' = mAmBC/(mA +mBC )). The potential energy V is a function of R, r, and q = a + b the angle between the r and R vectors. III. The conjugate momenta To eventually make a connection to the quantum mechanical Hamiltonian, it is necessary to rewrite the above H in terms of the coordinates and their so-called conjugate momenta. For any coordinate Q, the conjugate momentum PQ is defined by PQ = L/ . For example, PR = m', Pr = m , Pa = m' R2, Pb = m r2 allow H to be rewritten as H = PR2/2m' + Pr2 /2m + Pa2 /2m'R2 + Pb2 /2mr2 + (PX2 + PY2 )/2M +V. Because the potential V contains no dependence on X or Y, the center of mass motion (i.e., the kinetic energy (PX 2 + PY 2 )/2M) will time evolve in a straight-line trajectory with no change in its energy. As such, it can be removed from further consideration. The total angular momentum of the three atoms about the Z- axis can be expressed in terms of the four remaining coordinates as follows: LZ = m' R2 - m r2 = Pa - Pb . Since this total angular momentum is a conserved quantity (i.e., for each trajectory, the value of LZ remains unchanged throughout the trajectory), one can substitute this expression for Pa and rewrite H in terms of only three coordinates and their momenta: H = PR2 /2m' + pr2 /2m + (LZ -Pb )2 /2m'R2 + Pb2 /2mr2 +V. Notice that Pb is the angular momentum of the BC molecule (Pb = m r2 db/dt). IV. The equations of motion From this Hamiltonian (or the Lagrangian L = T - V), one can obtain, using d/dt d[T-V]/di = d[T-V]/dQi , the equations of motion for the three coordinates and three momenta: dPR /dt = L/R = -V/R - (LZ -Pb )2 /m'R3 , dPr /dt = -V/r - Pb2 /mr3 , dPb /dt = -V/b = -V/q , dR/dt = PR /m', dr/dt = Pr /m , db/dt = Pb /mr2 . V. The initial conditions These six classical equations of motion that describe the time evolution of r, R, b, PR , pr , and Pb can then be solved numerically as discussed earlier. The choice of initial values of the coordinates and momenta will depend on the conditions in the A + BC collision that one wishes to simulate. For example, R will be chosen as very large, and PR will be negative (reflecting inward rather than outward relative motion) and with a magnitude determined by the energy of the collision PR2 /2m' = Ecoll.. If the diatomic BC is initially in a particular vibrational state, say v, the coordinate r will be chosen from a distribution of values given as the square of the v-state's vibrational wavefunction |yv (r)|2 . The value of pr can then be determined within a sign from the energy ev of the v state: pr2 /2m + VBC = ev , where VBC is the BC molecule's potential energy as a function of r. Finally, the angle b can be selected from a random distribution within 02 ]-1/2 exp (-(R-R0 )2 /(42 )). Professor Rick Heller Here, the parameter 2 gives the uncertainty or "spread" along the R- degree of freedom for this wavefunction, defined as the mean squared displacement away from the average coordinate R0 , since it can be shown for the above function that Ú (R-R0)2 yR* (R) yR (R) R2 dR = 2 and that R0 = Ú R yR *(R) yR (R) R2 dR. It can also be shown that the parameter PR0 is equal to the average value of the momentum along the R coordinate: Ú yR (R)*{- i h yR/R } R2 dR = PR0, and that the uncertainty in the momentum along the R coordinate: = Ú yR (R)*{- i h yRR + PR0 yR (R) }2 R2 dR is related to the uncertainty in the R-coordinate by = h2/4. Note that for these coherent state wavefunctions, the uncertainties fulfill the general Heisenberg uncertainty criterion for any coordinate Q and its conjugate momentum P: > h2/4 but have the smallest possible product of uncertainties. In this sense, the coherent state function is the closest possible to a classical description (where uncertainties are absent). II. Time propagation of the wavefunction must be effected- the basis expansion approach The difficult part of the propagation process involves applying the time evolution operator exp (- iHt/h) to this initial wavefunction. The most powerful tools for applying this operator fall under the class of methods termed Feynman path-integral methods (see the links to Professors Greg Voth, Nancy Makri, and Jim Doll). However, let us first examine an approach that has been used for a longer time, for example in the hands of Professors George Schatz and John Light. a. Expanding the wavefunction in a basis If, alternatively, one prefers to attempt to solve the time-independent Schrödinger equation Hyk = Ek yk the most widely used approaches involve expanding the unknown y(R,r,a,b) in a basis consisting of products of functions of R, functions of r, and functions of a and of b: yk = Sn,v,L,M Ck (n,v,L,M) Fn (R) yv(r) yL (b) yM (a). This expansion is what one would use in the case of a non-reactive collision. If a chemical reaction (e.g., A + BC Æ AB + C) is to be examined, then one must add to the sum of terms given above another sum which contains products of functions of the so-called exit-channel coordinates (r, t, g, and f). In this case, one uses yk = Sn,v,L,M Ck (n,v,L,M) Fn (R) yv(r) yL (b) yM (a) + Sn,v,L,M Ck (n,v,L,M) Fn (r) yv(t) yL (g) yM (f). Of course, the yv , yL, Fn, and yM functions relating to the exit-channel are different than those for the entrance channel (e.g., yv (r) would be a AB-molecule vibrational wavefunction, but yv (r) would be a BC-molecule vibrational function). As will be seen from the discussion to follow shortly, expanding y in this manner and substituting into the Schrödinger equation produces a matrix eigenvalue equation whose eigenvalues are the Ek and the eigenvectors are the expansion coefficients Ck (n,v,L,M). Given any wavefunction y(R,r,a,b;t=0) at t=0 (e.g., the one shown above may be appropriate for a collision of A with a BC molecule in a specified vibrational/rotational state), one can then express this wavefunction at later times as follows: y(R,r,a,b;t) = Sk yk exp(-i Ekt/h) yk * y(R,r,a,b;t=0) dr . This expression arises when one uses the facts that: 1. the {yk } form a complete set of functions so one can expand y(R,r,a,b;t=0) in this set, and 2. the time evolution operator exp (-i Ht/h) can be applied to any eigenfunction yk and produces exp (-i Ht/h) yk = exp (-i Ekt/h) yk . So, it is possible to follow the time development of an initial quantum wavefunction by first solving the time-independent Schrödinger equation for all of the {yk } and then expressing the initial wavefunction in terms of the eigenfunctions (i.e, the integral Ú yk * y(R,r,a,b;t=0) dr is nothing but the expansion coefficient of y(R,r,a,b;t=0) in terms of the yk ) after which the wavefunction at a later time can be evaluated by the expression y(R,r,a,b;t) = Sk yk exp(-i Ekt/h) Ú yk * y(R,r,a,b;t=0) dr. b. The resulting matrix eigenvalue problem Let us now return to see how the expansion of y in the manner yk = Sn,v,L,M Ck (n,v,L,M) Fn (R) yv(r) yL (b) yM (a) (or the generalization in which both entrance- and exit-channel product functions are used) produces a matrix eigenvalue problem. Substituting this expansion into Hy =Ey and subsequently multiplying through by the complex conjugate of one particular Fn' (R) yv' (r) yL' (b) yM' (a)and integrating over R, r, b, and a one obtains a matrix eigenvalue equation: H Ck = Ek Ck in which the eigenvector is the vector of expansion coefficients Ck (n.v,L,M) and the H matrix has elements H(n,v,L,M|n',v',L',M') = < Fn' (R) yv' (r) yL' (b) yM' (a)| H | Fn (R) yv(r) yL (b) yM (a)> given in terms of integrals over the basis functions with the above Hamiltonian operator in the middle. Notice that in writing the integral on the right-hand side of the above equation, the so-called Dirac notation has been used. In this notation, any integral with a wavefunction on the right, a complex conjugated wavefunction on the left, and an operator in the middle is written as follows: ÚY*Opfdq = < Y|Op|f> c. The dimension of the matrix scales as the number of basis functions to number of coordinates power What makes the solution of the four-variable Schrödinger equation much more difficult than the propagation of the four-coordinate (plus four-momenta) Newtonian equations? The Hamiltonian differential operator is non-separable because V depends on R, r, and q (which is related to b and a by q = a + b ) in a manner that can not be broken apart into an R-dependent piece plus an r-dependent piece plus a q-dependent piece. As a result, the four dimensional second-order partial differential equation can not be separated into four second-order ordinary differential equations. This is what makes its solution especially difficult. To make progress solving the four-dimensional Schrödinger equation on a computer, one is thus forced either to 1. expand y (R,r,a,b) as S C(n,v,L,M) Fn (R) yv(r) yL (b) yM (a), and then solve the resultant matrix eigenvalue problem, or to 2. represent y (R,r,a,b) on a grid of points in R, in r, and in a and b space and use finite-difference expressions (such as [F(R+dR) + F(R-dR) -2F(R)]/d2 to represent d2F(R)/dR2 ) to also express the derivatives on this same grid to ultimately produce a sparse matrix whose eigenvalues must be found. In either case, one is left with the problem of finding eigenvalues of a large matrix. In, the former case, the dimension of the H matrix is equal to the product of the number of Fn (R) functions used in the expansion times the number of yv(r) functions used times the number of yL (b) used and times the number of yM (a) functions. In the latter case, the dimension of H is equal to the number of grid points along the R coordinate times the number of grid points along r times the number along b and times the number along a. It is the fact that the dimension of H grows quartically with the size of the basis or grid set used (because there are four coordinates in this problem; for a molecule with N atoms, the dimension of H would grow as the number of grid points or basis functions to the 3N-6 power!) and that the computer time needed to find all of the eigenvalues of a matrix grows as the cube of the dimension of the matrix (to find one eigenvalue requires time proportional to the square of the dimension) that makes the solution of this quantum problem very difficult. Let us consider an example. Suppose that a grid of 100 points along the R-coordinate and 100 points along the r-coordinate were used along with say only 10 points along the a and b-coordinates. This would produce a (albeit very sparse) H matrix of dimension 106 . To find just one eigenvalue of such a large matrix on a 100 Mflop desktop workstation would require of the order of several times (106)2 /(100x106 ) sec, or at least 100 minutes. To find all 106 eigenvalues would require 106 times as long. III. Time propagation of the wavefunction can be effected using Feynman path integrals instead Given a wavefunction at t=0, y(R,r,a,b;t=0), it is possible to propagate it forward in time using so-called path integral techniques that Richard Feynman pioneered and which have become more popular and commonly used in recent years (see links to Professors Greg Voth, Jim Doll, and Nancy Makri). In these approaches, one divides the time interval between t=0 and t=t into N small increments of length t = t/N, and then expresses the time evolution operator as a product of N short-time evolution operators: exp[-i Ht/h ] = {exp[-i Ht/h ]}N . For each of the short time steps t , one then approximates the propagator as exp[-i Ht/h] = exp[-i Vt/2h] exp[-i Tt/h] exp[-i Vt/2h], where V and T are the potential and kinetic energy operators that appear in H = T + V. It is possible to show that the above approximation is valid up to terms of order t4 . So, for short times (i.e., small t ), these so-called symmetric split operator approximations to the propagator should be accurate. The time evolved wavefunction y (t) can then be expressed as y (t) = { exp[-i Vt/2h] exp[-i Tt/h] exp[-i Vt/2h]}N y (t=0). The potential V is (except when external magnetic fields are present) a function only of the coordinates {qj } of the system, while the kinetic term T is a function of the momenta {pj } (assuming cartesian coordinates are used). By making use of the completeness relations for eigenstates of the coordinate operator 1 = dqj| qj> < qj| and inserting this identity N times (once between each power of exp[-i Vt/2h] exp[-i Tt/h] exp[-i Vt/2h]), the expression given above for y (t) can be rewritten as follows: y (qN ,t)= dqN-1 dqN-2 . . . dq1 dq0 exp{(-it/2h)[V(qj) + V(qj-1)]} < qj| exp(-itT / h ) |qj-1>y (q0 ,0). Then, by using the analogous completeness identity for the momentum operator 1 = 1 / h dpj| pj>< pj | one can write < qj| exp(-itT / h ) |qj-1> = 1 / hdp < qj|p > exp(-ip2t / 2mh ) < p|qj-1 >. Finally, by using the fact that the momentum eigenfunctions |p>, when expressed as functions of coordinates q are given by < qj|p > = (1/2p)1/2 exp(ipq/h), the above integral becomes < qj | exp(-itT / h) |qj-1> = 1/2 ph dp exp(-ip2 t / 2mh) exp[ip(qj - qj - 1)]. This integral over p can be carried out analytically to give < qj | exp(-itT / h) |qj-1> =exp[im(qj - qj - 1)2 / 2ht]. When substituted back into the multidimensional integral for y (qN ,t), we obtain y (qN ,t)= dqN-1 dqN-2 . . . dq1 dq0 exp{(-it/2h)[V(qj) + V(qj-1)]} exp[im(qj - qj-1)2 /2ht] y (qo,0) or y (qN ,t) = dqN-1 dqN-2 . . . dq1 dq0 exp{i/h [m(qj - qj-1)2 / 2t - t(V(qj) + V(qj-1))/2]} y (q0 , t=0). Why are such quantities called path integrals? The sequence of positions q0 , q2 , ... qN describes a "path" connecting q0 to qN . By integrating over all of the intermediate positions q1 , q2 ,... qN-1 one is integrating over all such paths. Further insight into the meaning of the above is gained by first realizing that (qj - qj-1)2 = (qj - qj-1)2 t = Tdt is the representation, within the N discrete time steps of length t, of the integral of Tdt over the jth time step, and that [V(qj) + V(qj-1)] = V(q)dt is the representation of the integral of Vdt over the jth time step. So, for any particular path (i.e., any specific set of q0 , q1, , ... qN-1 , qN values), the sum over all N such terms [m(qj - qj-1)2 / 2t - t(V(qj) + V(qj-1))/2] is the integral over all time from t=0 until t=t of the so-called Lagrangian L = T - V: [m(qj - qj-1)2 / 2t - t(V(qj) + V(qj-1))/2] = Ldt. This time integral of the Lagrangian is called the "action" S in classical mechanics. Hence, the N-dimensional integral in terms of which y (qN ,t) is expressed can be written as y (qN ,t) = exp{i / h dt L } y (q0 ,t=0). Here, the notation "all paths" is realized in the earlier version of this equation by dividing the time axis from t=0 to t=t into N equal divisions, and denoting the coordinates of the system at the jth time step by qj . By then allowing each qj to assume all possible values (i.e., integrating over all possible values of qj ), one visits all possible paths that begin at q0 at t=0 and end at qN at t=t. By forming the classical action S S = dtL for each path and then summing exp(iS/h) y ( q0 , t=0) over all paths and multiplying by , one is able to form y (qN ,t). dqN-1 dqN-2 . . . dq1 dq0 exp{i/h [m(qj - qj-1)2 / 2t - t(V(qj) + V(qj-1))/2]} The difficult step in implementing this Feynman path integral method in practice involves how one identifies all paths connecting q0 , t=0 to qN , t. Each path contributes an additive term involving the complex exponential of the quantity [m(qj - qj-1)2 / 2t - t(V(qj) + V(qj-1))/2] This sum of many, many (actually, an infinite number) oscillatory exp(iS/h) = cos (S/h) + i sin(S/h) terms is extremely difficult to evaluate because of the tendency of contributions from one path to cancel those of another path. The evaluation of such sums remains a very active research subject. The most commonly employed approximation to this sum involves finding the path(s) for which the action S= [m(qj - qj-1)2 / 2t - t(V(qj) + V(qj-1))/2] is smallest because such paths produce the lowest frequency oscillations in exp(iS/h), and thus are not subject to cancellation by contributions from other paths. The path(s) that minimize the action S are, in fact the classical paths. That is, they are the paths that the system whose quantum wavefunction is being propagated in time would follow if the system were undergoing classical Newtonian mechanics subject to the conditions that the system be at q0 at t=0 and at qN at t=t. In this so-called semi-classical approximation to the propagation of the initial wavefunction using Feynman path integrals, one finds all classical paths that connect q0 at t=0 and at qN at t=t, and one evaluates the action S for each such path. One then applies the formula y (qN ,t) = exp{i / hdt L } y (q0 ,t=0) but includes in the sum only the contribution from the classical path(s). In this way, one obtains an approximate quantum propagated wavefunction via a procedure that requires knowledge of only classical propagation paths. b. How is information extracted from a quantum dynamics simulation? Once an initial quantum wavefunction has been propagated for a time long enough for the event of interest to have occurred (e.g., long enough for an A + BC Æ AB + C reactive collision to take place in the above example problem), one needs to interpret the wavefunction in terms of functions that are applicable to the final-state. In the example considered here, this means interpreting y(R, r, a, b; t) in terms of exit-channel basis states that describe an intact AB molecule with a C atom moving away from it (i.e., the terms Sn,v,L,M Ck (n,v,L,M) Fn (r) yv(t) yL (g) yM (f)). As explained earlier when the classical trajectory approach was treated, to describe this final-state configuration of the three atoms, one uses different coordinates than R,r,a,b that were used for the initial state. The appropriate coordinates are shown below in the figure. That is, one projects the long-time wavefunction y(R,r,a,b;t) = Sk yk exp(-i Ekt/h) Ú yk * y(R,r,a,b;t=0) dr onto a particular exit-channel product function Fn (r) yv(t) yL (g) yM (f) to obtain the amplitude A of that exit channel in the final wavefunction: A = < Fn (r) yv(t) yL (g) yM (f)| y(R,r,a,b;t)> = Sk < Fn (r) yv(t) yL (g) yM (f)| yk> exp(-iEkt/h) Ú yk * y(R,r,a,b;t=0) dr. Using the earlier expansion for yk , one sees that < Fn (r) yv(t) yL (g) yM (f)| yk> = Ck (n,v,L,M) for the exit channel, so the modulus squared of A, which gives the probability of finding the AB + C system in the particular final state Fn (r) yv(t) yL (g) yM (f), is given by: |A|2 = | Sk Ck (n,v,L,M) exp(-i Ekt/h) Ú yk * y(R,r,a,b;t=0) dr |2 . This result can be interpreted as saying that the probability of finding the state Fn (r) yv(t) yL (g) yM (f) is computed by (1) first, evaluating the projection of the initial wavefunction along each eigenstate (these components are the Ú yk * y(R,r,a,b;t=0) dr), (2) multiplying by the projection of the specified final wavefunction along each eigenstate (the < Fn (r) yv(t) yL (g) yM (f)| yk> = Ck (n,v,L,M)), (3) multiplying by a complex phase factor exp(-i h Ek t) that details how each eigenstate evolves in time, and (4) summing all such products after which (5) taking the modulus squared of the entire sum. 3. Present Day Challenges in Chemical Dynamics In addition to the people specifically mentioned in this text for whom I have already provided web links, I encourage you to look at the following web pages for further information: a. Large Biomolecules and Polymers Following the motions of large molecules even using classical Newton equations involves special challenges. First, there is the fact that large molecules contain more atoms, so the Newtonian equations that must be solved involve many coordinates (N) and momenta. Because the potential energy functions,and thus the Newtonian forces, depend on the relative positions of the atoms, of which there are N(N-1)/2 , the computer time needed to carry out a classical trajectory varies as (at least) the square of the number of atoms in the molecule. Moreover, to adequately sample an ensemble of initial coordinates and momenta appropriate to a large molecule, one needs to run many trajectories. That is, one needs to choose a range initial values for each of the molecule's internal coordinates and momenta; the number of such choices is proportional to N. The net result is that computer time proportional to N3 (or worse) is needed,and this becomes a major challenge for large biomolecules or polymers. One way this problem is being addressed is to use parallel computers which allow one to run trajectories with different initial conditions on different processors. Another problem that is special to large molecules involves the force field itself (see, for example, the web pages of Professors Andy McCammon, Charles Brooks and the late Peter Kollman, as well as of Biosym and CHARMm). Usually, the potential energy is expressed as a sum of terms involving: 1. bond stretching potentials of either harmonic V(r) = 1/2 k (r-re)2 or Morse form V(r) = De (1-exp(-b(r-re)))2 2. bond angle bending potentials of harmonic form depending on the angles between any three atoms A, B, C that are bonded to one another (e.g., the H-C-H angles in -CH2 - groups) 3. dihedral angles (e.g., the H-C-C-H angle in -H2C-CH2- groups 4. Van der Waals interactions V(r) = Ar-12 - Br-6 among all pairs of atoms (these potentials allow for repulsions among, for example, the H atoms of a methyl group in H3C-CH2-CH2-CH2-CH2-C H3 to repel the other H atoms as this "floppy" molecule moves into geometries that bring such groups into strong steric interaction) 5. electrostatic attraction (for opposite charges) and repulsions (for like charges) among charged or highly polar groups (e.g., between pairs of phosphate groups in a biopolymer, between a Na+ ion and a phosphate group, or between a charged group and a neighboring H2O molecule) 6. polarization of the electronic clouds of the large solute molecule or of surrounding solvent molecules induced by ionic or highly polar groups approaching highly polarizable parts of the molecule. One of the earliest pioneers in applying fundamental statistical mechanics, molecular mechanics, and solvation modelling theories to bio-molecules is Professor Arieh Warshel (below). His groups more recent interests include the dynamics of photobiological processes, enzyme catalysis and protein action, as well as the study of basic chemical reactions in solutions. The primary difficulties in using such a multi-parameterized force field are 1. that the values of the parameters characterizing each kind of potential have to be obtained either by requiring the results of a simulation to "fit" some experimental data or by extracting these parameters from ab initio quantum chemical calculations of the interaction potentials on model systems; 2. these developments and calibrations of the parameters in such force fields are still undergoing modifications and improvements, so they are not yet well established for a wide range of functional groups, solvents, ionic strengths, etc. One more difficulty that troubles most large molecule simulations is the wide range of time scales that must be treated when solving the Newton equations. For example, in attempting to monitor the uncoiling of a large biomolecule, which may occur over a time scale of ca. 10-5 to 1 sec, one is usually forced to "freeze" the much faster motions that occur (i.e., to treat the C-H, O-H, and N-H bonds, which oscillate over ca.10-14 sec, as rigid). If one were to attempt to follow the faster motions, one would have to take time steps in solving the Newton equations of ca. 10-14 sec, as a result of which the uncoiling event would require 109 to 1014 time steps. Even with a computer that could perform 109 floating point operations per second (i.e., a one gigaflop computer), and assuming that a single time step would use at least 1000 floating point operations (which is very optimistic), a single such trajectory would use 103 to 108 seconds of computer time. This is simply unrealistic. As a result, one is forced to freeze (i.e., by constraining to fixed geometries) the faster motions so that the Newton equations can be solved only for the slower motions. It is an area of current research development to invent new methods that allow such multiple time scale issues to be handled in a more efficient and accurate manner. b. Mixed Classical and Quantum Dynamics One of the most pressing issues in chemical reaction dynamics involves how to treat a very large molecular system that is undergoing a chemical reaction or a photochemical event (i.e., a change that affects its electronic structure). Because the system's bonding and other attributes of its electronic wavefunction are changing during such a process, one is required to use quantum mechanical methods. However, such methods are simply impractical (because their computer time, memory, and disk space needs scale as the number of electrons in the system to the fourth (or higher) power) to use on very large molecules. If the chemical reaction and/or photon excitation is localized within a small portion of the molecular system (that may involve solvent molecules too), there are so-called quantum mechanics/molecular mechanics (QMMM) methods currently under much development that can be employed. In these methods, one uses an explicit quantum chemical (ab initio or semi-empirical) method to handle that part of the system that is involved in the electronic process. For example, in treating the tautomerization of formamide H2 N-CHO to formamidic acid HN=CHOH in aqueous solution, one must treat the six atoms in the formamide explicitly using quantum chemistry tools. The surrounding water molecules (H2O)n could then be handled using a classical force field (i.e., the formamide-water and water-water interaction potentials as well as the internal potential energy function of each H2O molecule could be expressed as a classical potential like those described earlier in this section). If, on the other hand, one wished to treat one or more of the H2O solvent molecules as actively involved in promoting the tautomerization (as shown in the figure below) one would use quantum chemistry tools to handle the formamide plus the actively involved water molecule(s), and use classical potentials to describe the interactions of the remaining solvent molecules with these active species and among one another. At present, much effort is being devoted to developing efficient and accurate ways to combine quantum treatment of part of the system with classical treatment of the remainder of the system. Professor Chris Cramer has been very active developing new ways to model the influence of surrounding solvent molecules on chemical species whose spectroscopy or reactions one wishes to study. c. The Car-Parrinello Method As discussed earlier, it is often important to treat the electronic degrees of freedom, at least for the electrons involved in a chemical reaction or in a photon absorption or emission event, at a quantum level while still, at least for practical reasons, treating the time evolution of the nuclei classically. A novel approach developed by Car and Parrinello, and reviewed very nicely by Remler and Madden, incorporates the task of optimizing the electronic wavefunction into the Newtonian dynamics equations used to follow the nuclear movements. The LCAO-MO coefficients Ci,k relating molecular orbital fi to atomic orbital ck are assumed to minimize an energy function E[{Ci,k}] which can be a Hartree-Fock or DFT-type function, for example. At each time step during which the Newtonian equations ma d2Xa /dt2 = - E/Xa are propagated, a corresponding propagation equation is used to follow the time evolution of the {Ci,k } coefficients. The latter equation was put forth by Car and Parrinello by defining a (fictitious) kinetic energy T = Si, k 1/2 m |dCi,k/dt|2 in which the time derivatives of the Ci,k coefficients are viewed as analogs to the time rate of changes of the molecule's Cartesian coordinates {Xa}, and m is a (fictitious mass introduced to give this expression the units of energy. Then, using the dependence of the electronic energy E on the {Ci,k } to define the potential energy V = E[{Ci,k}], and then introducing a (fictitious) Lagrangian L = T-V whose classical action S = Ú L dt is made stationary with respect to variations in the {Ci,k } coefficients but subject to the constraint that the orbitals {fi} remain orthonormal: di,j = S k,l C*i,k Sk,l Cj,l the following dynamical equations are obtained for the Ci,k coefficients: m d2Ci,k /dt2 = -E/Ci,k - Sj,l li,j Sk,. Cj,l . The left hand side of this equation gives the mass m multiplied by the "acceleration" of the Ci,k variable. The right side gives the "force" -E/Ci,k acting on this variable, and - Sj,l li,j Sk,. Cj,l is the coupling force that connects Ci.k to the other Cj,l variables (this term arises, with its Lagrange multipliers li,j , from the constraints di,j = S k,l C*i,k Sk,l Cj,l relating to the orthonormality of the molecular orbitals. The result is that the time evolution of the Ci,k coefficients can be followed by solving an equation that is exactly the same in form as the Newtonian equations ma d2Xa /dt2 = - E/Xa used to follow the nuclear positions. Of course, one wonders whether the idea used by Car and Parrinello to obtain the classical-like equations for the time development of the Ci,k coefficients is valid. These equations were obtained by making stationary, subject to orthonormality constraints, an action. However, if one were to set the mass m equal to zero, the Lagrangian L would reduce to L = T-V = -V, since the (fictitious) kinetic energy would vanish. Since the Ci,k coefficients, in a conventional quantum chemical study, would be chose to make E[{Ci.k }] stationary, it can be shown that making the action stationary when T vanishes is equivalent to making V stationary, which, in turn, makes E[{Ci,k}] stationary. So, the Car-Parrinello approach is valid, but one has to be aware of the restriction to the ground-state (since E[{Ci,k }] is made stationary with no constraint of orthogonality to lower-energy states) and one must realize that finite-time-step propagations with non-zero m render the {Ci,k} amplitudes obtained via classical propagation not equivalent to those obtained by minimizing E at each geometry. Before closing this section on dynamics, I wish to bring to the readers' attention several other leading researchers in this exciting area of theoretical chemistry by showing photos of them below.