Trajectory initial conditions

Procedures for selecting initial values of coordinates and momenta for an ensemble of trajectories has been described in detail in recent chapters entitled Monte Carlo Sampling for Classical Trajectory Simulations and Classical Trajectory Simulations Initial Conditions. In this section a brief review is given of methods for selecting initial conditions for trajectory simulations of unimolecular and bimolecular reactions and gas-surface collisions. 

---

Figure 3. The ensemble average of Che three normal mode energies vs. time. Total ensemble size Is 50 and the planar trajectory Initial conditions had most of the energy In the CO normal mode. Total energy of each trajectory Is 20000 cm and total angular momentum Is zero.

Figure 3. Definition of the coordinates used for sampling trajectory initial conditions for a gas-surface collision. (From Song et al. [69].)...

An ensemble of trajectories is calculated in a trajectory simulation, with each trajectory in the ensemble specified by the system s initial set of momenta p and coordinates q. The initial ensemble of p and q is chosen to represent the experiment under investigation or chosen so that a particular dynamical attribute of the system may be studied. Distribution functions are usually sampled randomly in choosing the ensemble of initial conditions and the methodology of sampling is often called Monte Carlo sampling. Procedures for choosing trajectory initial conditions to represent unimolecular and bimolecular reactions, and gas-surface collisions are described in the section on trajectory initial conditions. 

With the current development of quantum chemistry, it is routine to evaluate Eq. [17] for the QM + MM model and the application of QM + MM direct dynamics is described in the section on trajectory initial conditions. However, in many situations the QM/MM boundary must cut through a chemical bond in a molecule. In such a case, the total electronic Hamiltonian cannot be divided as for the QM + MM model. Different approaches have been developed to treat QM/MM interactions when the boundary cuts through a chemical bond. Gao et al. identified a criterion for treating a covalent bond at the QM/MM boundary. In general, a reasonable boundary method should be able to mimic the real physical properties of the model system as closely as possible. The obtained properties such as vibrational frequencies, energies, and electronegativities, etc. should be comparable to experiment or accurate ab initio calculations. 

Implementation of the FBSD schemes requires knowledge of the initial density matrix in the coherent state representation. Usually, the initial density corresponds either to the ground vibrational state of a polyatomic molecule or a Boltzmann distribution. Below we describe ways of obtaining the coherent state matrix element through closed form expressions or in terms of an imaginary time path integral evaluated along the same Monte Carlo random walk which samples the trajectory initial conditions. 

Even expression ( B3.4.31), altiiough numerically preferable, is not the end of the story as it does not fiilly account for the fact diat nearby classical trajectories (those with similar initial conditions) should be averaged over. One simple methodology for that averaging has been tln-ough the division of phase space into parts, each of which is covered by a set of Gaussians [159, 160]. This is done by recasting the initial wavefunction as... 

It is convenient to analyse tliese rate equations from a dynamical systems point of view similar to tliat used in classical mechanics where one follows tire trajectories of particles in phase space. For tire chemical rate law (C3.6.2) tire phase space , conventionally denoted by F, is -dimensional and tire chemical concentrations, CpC2,- are taken as ortliogonal coordinates of F, ratlier tlian tire particle positions and velocities used as tire coordinates in mechanics. In analogy to classical mechanical systems, as tire concentrations evolve in time tliey will trace out a trajectory in F. Since tire velocity functions in tire system of ODEs (C3.6.2) do not depend explicitly on time, a given initial condition in F will always produce tire same trajectory. The vector R of velocity functions in (C3.6.2) defines a phase-space (or trajectory) flow and in it is often convenient to tliink of tliese ODEs as describing tire motion of a fluid in F with velocity field/ (c p). 

As already mentioned, the motion of a chaotic flow is sensitive to initial conditions [H] points which initially he close together on the attractor follow paths that separate exponentially fast. This behaviour is shown in figure C3.6.3 for the WR chaotic attractor at /c 2=0.072. The instantaneous rate of separation depends on the position on the attractor. However, a chaotic orbit visits any region of the attractor in a recurrent way so that an infinite time average of this exponential separation taken along any trajectory in the attractor is an invariant quantity that characterizes the attractor. If y(t) is a trajectory for the rate law fc3.6.2] then we can linearize the motion in the neighbourhood of y to get... 

In applying minimal END to processes such as these, one finds that different initial conditions lead to different product channels. In Figure 1, we show a somewhat truncated time lapse picture of a typical trajectory that leads to abstraction. In this rendering, one of the hydrogens of NHaD" " is hidden. As an example of properties whose evolution can be depicted we display interatomic distances and atomic electronic charges. Obviously, one can similarly study the time dependence of various other properties during the reactive encounter. 

These charge-transfer structures have been studied [4] in terms a very limited number of END trajectories to model vibrational induced electron tiansfer. An electronic 3-21G-1- basis for Li [53] and 3-21G for FI [54] was used. The equilibrium structure has the geometry with a long Li(2)—FI bond (3.45561 a.u.) and a short Li(l)—H bond (3.09017 a.u.). It was first established that only the Li—H bond stietching modes will promote election transfer, and then initial conditions were chosen such that the long bond was stretched and the short bond compressed by the same (%) amount. The small ensemble of six trajectories with 5.6, 10, 13, 15, 18, and 20% initial change in equilibrium bond lengths are sufficient to illustrate the approach. 

In what is called BO MD, the nuclear wavepacket is simulated by a swarm of trajectories. We emphasize here that this does not necessarily mean that the nuclei are being treated classically. The difference is in the chosen initial conditions. A fully classical treatment takes the initial positions and momenta from a classical ensemble. The use of quantum mechanical distributions instead leads to a seraiclassical simulation. The important topic of choosing initial conditions is the subject of Section II.C. 

One way to overcome this problem is to start by setting up the ensemble of trajectories (or wavepacket) at the transition state. If these bajectories are then run back in time into the reactants region, they can be used to set up the distribution of initial conditions that reach the barrier. These can then be run forward to completion, that is, into the products, and by using transition state theory a reaction rate obtained [145]. These ideas have also been recently extended to non-adiabatic systems [146]. 

The chaotic nature of individual MD trajectories has been well appreciated. A small change in initial conditions (e.g., a fraction of an Angstrom difference in Cartesian coordinates) can lead to exponentially-diverging trajectories in a relatively short time. The larger the initial difference and/or the timestep, the more rapid this Lyapunov instability. Fig. 1 reports observed behavior for the dynamics of a butane molecule. The governing Newtonian model is the following set of two first-order differential equations ... 

It is appropriate to consider first the question of what kind of accuracy is expected from a simulation. In molecular dynamics (MD) very small perturbations to initial conditions grow exponentially in time until they completely overwhelm the trajectory itself. Hence, it is inappropriate to expect that accurate trajectories be computed for more than a short time interval. Rather it is expected only that the trajectories have the correct statistical properties, which is sensible if, for example, the initial velocities are randomly generated from a Maxwell distribution. 

Coordin ates of atom s can he set by n orm al translation orrotation of HyperCh cm molecules, fo set initial velocities, however, it is necessary to edit th e H l. file explicitly. The tin it o f velocity in the HIN file is. An gstrom s/picosecon d.. Areact.hin file and a script react.scr are in eluded with HyperChem to illustrate one simple reacting trajectory. In order to have these initial velocities used in a trajectory the Restart check box of the Molecular Dynamics Option s dialog box must he checked. If it is n ot, the in itial velocities in the HIN file will be ignored and a re-equilibration to the tern peratiire f of th e Molecular Dyn am ics Option s dialog box will occur. This destroys any imposed initial conditions on the molecular dynamics trajectory. 

If Lh c con Stan t Ictn pcratti re a Igori Lli rn is used in a trajectory analysis, then the initial conditions arc constantly being modified according to the sirn ulation of th c con stan t tern perattirc bath an d th e relaxation of th e m olecu lar system to that bath temperature, fhe effect of such a bath on a trajectory analysis is less studied than for th c sirn 11 lation of cqu i libriii m behavior. 

A classical molecular dynamics trajectory is simply a set of atoms with initial conditions consisting of the 3N Cartesian coordinates of N atoms A(X, Y, Z ) and the 3N Cartesian velocities (v a VyA v a) evolving according to Newton s equation of motion ... 

---