Molecular dynamics

Molecular dynamic information is obtained from a study of the variation of the muon spin relaxation rate with temperature. Reorientation depolarizes the muons by causing anisotropic or dipolar terms in the electron-muon hyperhne interaction to fluctuate. Peaks in the relaxation rate (analogous to 2i minima in NMR) occur when the reorientation rate matches the frequency of the dominant transition between the coupled muon-electron spin states. The correlation time, T, at each temperature can be obtained using the derivations of Cox and Sivia [14]. The measured relaxation rate X is given by the following expression  

If the process under scrutiny shows Arrhenius behavior, a plot of ln(l/T) versus l/T should reveal the activation parameters. As the activation energy is determined from the gradient of such a plot, its value is independent of the chosen transition frequency, but the attempt frequency will vary. 

An approximate value of the hyperfine coupling constant and hence the type of radical formed can be deduced using repolarization curves. These are plots of the initial amplitude of the muon relaxation signal, which increases as the hyperfine interaction is decoupled by an applied magnetic field. However, these estimates can be distorted by anisotropic terms and motional effects and more accurate values of the hyperfine coupling constants would be given by TF- xSR (transverse field) and ALC- xSR (avoided level crossing) measurements. 

Molecular dynamics involves the calculation of the time-dependent movement of each atom in a molecule11121. This is achieved by solving Newton s equations of motion. For this process the energy surface and the derivative of the energy in terms of the nuclear coordinates are required (Eqs. 4.1 and 4.2 mass m, acceleration a, potential energy E, coordinates r, time t). 

The maximum time step depends on the highest frequencies (usually R-H bonds). The main problem in dynamics calculations is that the movements leading to appreciable conformational changes usually have lower frequencies (milliseconds e.g., torsions). The computation of such long time intervals is usually prohibitive due to the CPU time involved. 

Molecular dynamics uses classical mechanics to study the evolution of a system in time. At each point in time the classical equations of motion are solved for a system of particles (atoms), interacting via a set of predefined potential functions (force field), after which the solution obtained is applied to predict positions and velocities of the particles for a (short) step in time. This step-by-step process moves the system along a trajectory in phase space. Assuming that the trajectory has sampled a sufficiently large part of phase space and the ergodicity principle is obeyed, all properties of interest can then be computed by averaging along the trajectory. In contrast to the Monte Carlo method (see below), the MD method allows one to calculate both the structural and time-dependent characteristics of the system. An interested reader can find a comprehensive description of the MD method in the books by Allen and Tildesley or Frenkel and Smit.  

Application of the MD method to zeolite systems was pioneered by De-montis and coworkers in 1986. Since then, MD has been used in an increasing number of zeolite modeling studies. Those papers appearing before 1992 have been reviewed in Ref. 9, and a review of later studies can be found in Ref. 3 some recent examples are discussed in subsequent parts of this chapter. 

A typical MD simulation protocol consists of the following steps (1) assign the initial values to the atomic coordinates and velocities (2) integrate the equations of motion (this is the heart of the MD method) and (3) analyze the results of the simulation. These steps are discussed below in relation to zeolite modeling. 

An initial set of coordinates and velocities of all particles in the system must be known to start a simulation. For the atoms of a zeolite framework, those initial coordinates are usually taken from experimental data, for example, from X-ray or neutron diffraction measurements collected in the Atlas of Zeolite Structure Types A 

Having defined the atomic coordinates, initial velocities must next be assigned. The atomic velocity components may be chosen randomly from either a Gaussian distribution at the desired temperature or from a uniform distribution in the interval fmax) where can be chosen to be equal to the 

Molecular dynamics (MD) is a long standing molecular modelling technique [23], In molecular dynamics we solve by numerical integration the equations of motion of a system of interacting particles with or without periodic boundary conditions. 

Molecular dynamics requires the description of the interactions between particles, the force field, the validity and quality of the results depend critically on the accuracy of this parameterization. The force field approach fixes molecular connectivity and it is not possible to create (or break) chemical bonds or study chemical reactions. Recently an approach that involves the combination of classical mechanics with electronic structure caleulations allows the intemuclear forces to be calculated on the fly from electronic structure calculation as the molecular dynamics develops [19], 

We consider a system of iVatom interacting particles moving under the influence of the forces between them. Where the position of the atom, r t), and its velocity, v t), evolve over time, t. Since, from Newton s equations of motion, where mi is the atomic mass  

Molecular dynamics [MD] is a well-established technique for simulating the structure and properties of materials in the solid, liquid, and gas phases and will only be briefly overviewed here [1]. In these simulations, N atoms are placed in a simulation cell with an initial set of positions and interact via an interatomic potential 7[r]. The force on each particle is determined by its interaction with all other atoms to within an interaction cutoff specified by 7[r]. For a given set of initial particle positions, velocities, and a specification of the position- or time-dependent forces acting on the particles, MD simulations solve the classical Newton s equations of motions numerically via finite-difference methods to calculate the time evolution of the particle trajectories. 

The simulation of material in the bulk phase is achieved via the use of periodic boundary conditions [1], which eliminate any nonporous surfaces in a finite-volume simulation cell. Conceptually, for a system with rectangular periodic boundaries, a particle exiting through a boundary to the right would reemerge through the left boundary via a simple subtraction of the cell length from the particle position. 

When MD techniques are applied to simulate the structure of amorphous materials, a common starting point is to quench from the liquid state. However, due to the fact that an MD time step is in the order of femtoseconds and limitations on the total possible MD simulation time from computational resources, the quench rates used in MD are several orders of magnitude faster than those found experimentally, and this can lead to the generation of structures that are not found experimentally. In these cases Monte Carlo techniques can be beneficial in developing an initial structure of the amorphous material, which can be further refined using MD or geometry optimization methods. 

Molecular dynamics simulation is perhaps the most powerful computational technique available for obtaining information on time dependent properties of molecular or atomic motion in zeolite crystals. It is used to obtain thermodynamic quantities and detailed dynamical information on sorption and diffusion processes in zeolite systems. For instance, the extent to which intramolecular vibration and framework motion assist sorption and diffusion of molecules can be simulated. The major limitation is its inability to model diffusion of larger sorbed molecules and electronic polarisability due to the huge amount of computer time and memory requirements. However, with the improvement in supercomputers and improved computing facilities, the full application of M.D. simulation to zeolite studies is becoming feasible. 

A molecular dynamics simulation samples the phase space of a molecule (defined by the position of the atoms and their velocities) by integrating Newton s equations of motion. Because MD accounts for thermal motion, the molecules simulated may possess enough thermal energy to overcome potential barriers, which makes the technique suitable in principle for conformational analysis of especially large molecules. In the case of small molecules, other techniques such as systematic, random. Genetic Algorithm-based, or Monte Carlo searches may be better suited for effectively sampling conformational space. 

In molecular dynamics simulations of polymers, Newton s equation of motion for the nth monomer with mass (n = 1,2. AT) and forces F acting on this monomer, 

From a thermodynamic point of view, molecular dynamics samples the microcanon-ical ensemble. In order to incorporate canonical thermal fluctuations by coupling the polymer system to a heat bath of canonical temperature T by means of a thermostat, Newton s equations of motion must be modified. Not the system energy E has to be constant, but the canonical expectation value ( ) at the given temperature. Thus, the task of the thermostat is twofold to keep the temperature constant and to sample the thermal 

The exclusive-or (XOR or ) operation, applied to two bits, yields 1, if both bits are unlike otherwise the 

The Andersen thermostat is very simple. After each time step Si, each monomer experiences a random collision with a fictitious heat-bath particle with a collision probability / coll = vSt, where v is the collision frequency. If the collisions are assumed to be uncorrelated events, the collision probability at any time t is Poissonian,pcoll(v, f) = v exp(—vi). In the event of a collision, each component of the velocity of the hit particle is changed according to the Maxwell-Boltzmann distribution p(v,)= exp(—wv /2k T)/ /Inmk T (i = 1,2,3). The width of this Gaussian distribution is determined by the canonical temperature. Each monomer behaves like a Brownian particle under the influence of the forces exerted on it by other particles and external fields. In the limit i — oo, the phase-space trajectory will have covered the complete accessible phase-space, which is sampled in accordance with Boltzmann statistics. Andersen dynamics resembles Markovian dynamics described in the context of Monte Carlo methods and, in fact, from a statistical mechanics point of view, it reminds us of the Metropolis Monte Carlo method. 

Another popular thermostat used in molecular dynamics simulations is the Langevin thermostat. It covers the heat-bath coupling part of the Langevin equation by friction and Gaussian random forces f. The Langevin equation basically describes the dynamics of a Brownian particle in solvent under the influence of external forces F. Its simplest form therefore reads  

TaUe 4. Dynamic parameters of cholestane spin probes in liquid crystal side chain polymers and low molecular weight analogues Rotational activation energies and anisotropy ratios 

Computer simukitions of the experimental tectra, applying the tteoretk approach discussed above, provided detailed information about the in probe motions in the various sy ems [35]. In Fig. 15 the correlation times for long axis rotation 

Inspection of the logoithmic plots reveak several discontinuities, whidi occur at the phase trandtimis. Within a pmticular phase the plots are linear. From the dopes of the straight lii the rotational activation energies have been determined. 

Th are listed in TaUe 4. As one can see, decreases with increadng order of the 

In this coimection, the effect of the molecular wdght M on the dde drain motion has been studied [35]. Variation of M betwren 4500 M 21000 product onty minor changes of the correlation times, in contrast to detectable changes of the thermodynamic parameters (see Table 1) [103-10. Apraiently, the plateau valiKs for the sidechain motions are already reachedatalowde eeofpolymeri2ation(Mn 4500). 

Homonuclear (e.g. H H) Overhauser effect in laboratory (NOE) and rotating (ROE) frame as a function of molecular correlation time If the product 1.12 the 

NOE vanishes and dipolar interactions cannot be detected in this way. 

Tj 1/APy (see Fig. 9.5). Chemical shifts are typically averaged by processes in vhich the lifetimes are 10 s or shorter [14]. Since the major part of conformational motions is much faster, usually conformationally averaged NMR spectra are observed. It means that most often NMR derived structural information is conformationally averaged. It should be pointed out that other spectroscopic techniques, like IR or UV ones, are characterized by different time scales. 

In two-dimensional (2D) spectroscopy, after the preparation period the perturbed spin system is left to evolve at its characteristic frequencies (evolution). Owing to the mixing, these frequencies can modulate the oscillations detected during acqui- 

As previously mentioned, the highest theoretically possible sensitivity could be achieved when the spin states of large y nuclei are initially perturbed and their magnetization is detected. Therefore, among many possible heteronuclear correlation schemes, those exploiting excitation and detection of sensitive nuclei (e.g. H) are of great importance. Such methods are characterized by two-fold coherence 

Molecules are dynamic, undergoing vibrations and rotations continually. Therefore the static picture of molecular structure provided by MM is not realistic. Flexibility and motion are clearly important to the biological functioning of biomacromolecules. These molecules are not static structures, but exhibit a variety of complex motions both in solution and in the crystalline state. Energy minimization concerns only the potential energy term of the total energy and so it treats the biomacromolecule as a static entity. The dynamic properties of the atoms in a macromolecule or the momentum of the atoms in space requires the description of the kinetic term. The momentum (p) is related to the force exerted on the atom (Ft) and the potential energy (V) by 

Thus the kinetic energy is related to the potential energy functions of the FF. The kinetic energy of an atom is related to the temperature T of the system by a simplified form  

The change in velocity v, is equal to the integral of acceleration over time. One numerically and iteratively integrates the classical equations of motion for every explicit atom N in the system by marching forward in time via tiny time increments, At in MolD. A number of algorithms exist for this purpose (Brooks et al, 1988 McCammon and Harvey, 1987) and the simplest formulation is shown below  

The total energy of the system, called the Hamiltonian (H), is the sum of the kinetic (K) and potential (V) energies  

This process of raising the temperature of the system will cover a time interval of 10-50 ps. The period of heating to the temperature of interest is followed by a period of equilibration with no temperature changes. The stabilization period will cover another time interval of 10-50ps. The mean kinetic energy of the system is monitored, and when it remains constant, the system is ready for study. The structure is in an equilibrium state at the desired temperature. 

Frenkel, D. and Smit, B. Understanding Molecular Simulation Molecular dynamics Academic Press 2001 

and Tildesley, DJ. Computer Simulation of Liquids Clarendon Press, Oxford 1989 

and Toboshnik, J. Introduction to Computer Simulation Methods. Applications to Physical Systems, Part 1. Addison-Wesley 1988 

Molecular Dynamics Simulation. Elementary Methods Monte Carlo Wiley Sons 1992 

The structure of a liquid-liquid interface is difficult to define because, by definition, we deal with a dynamic molecular interface with thermal fluctuations. Our knowledge to date stems mainly from molecular dynamic calculations, from capacitance and surface tension measurements, and from some experimental spectroscopic investigations. 

Alternatively, the phonon frequencies can also be obtained from molecular dynamics mns [3]. The relative displacements ujy in Eq. 5.8 can be Fourier-transformed, leading to 

The spectral analysis of ek(co), i.e. finding sharp peaks in the power spectrum 

The advantage of this method is that it can be used for more complicated systems, where explicit calculation of the full dynamical matrix would be extremely expensive. Furthermore, we can calculate the temperature dependence of the phonon spectrum by simply performing molecular dynamics simulations at different temperatures. The temperature dependence of the phonon spectrum is due to anharmonic effects, i.e., at larger displacements when terms higher than second order contribute to the potential energy in Eq. 5.4. 

There are two main approaches used to simulate polymer materials molecular dynamics and Monte Carlo methods. The molecular dynamics approach is based on numerical integration of Newton s equations of motion for a system of particles (or monomers). Particles follow dctcr-ministic trajectories in space for a well-defined set of interaction potentials between them. In a qualitatively different simulation technique, called Monte Carlo, phase space is sampled randomly. Molecular dynamics and Monte Carlo simulation approaches are analogous to time and ensemble methods of averaging in statistical mechanics. Some modern computer simulation methods use a combination of the two approaches. 

Consider a molecular dynamics simulation of a system consisting of K particles in a cubic box of volume V=L. Periodic boundary conditions are typically used to minimize surface effects. Periodic boundary condi-tions correspond to densely filling space with identical copies of the simulation box (see Fig. 9.28 for a two-dimensional sketch of periodic boundary conditions). As particles leave the simulation box from one side, they automatically reenter it from the opposite side. 

In one of the simplest models, all particles are identical with mass m and size a interacting with each other via a Lennard-Jones potential [Eq. (3.96)1. A cutoff in the potential is introduced at in order to speed- 

A typical cutoff used for attractive interactions is rcui = 2.5o- and for purely repulsive interactions the cutoff is at the minimum of the potential 

A simulation starts with an initial set of positions r and velocities diy/dt of all particles (with i=, 2. n for a system of n particles). 

Consequently, it may be necessary to use a more sophisticated computer procedure, molecular dynamics (section 5.3), to obtain the lowest minimum energy value and as a result the best model for the molecule. This final structure may be moved around the screen and expanded or reduced in size. It can also be rotated about the x or y axis to view different elevations of the molecule. 

Molecular mechanics calculations are made at 0 K, that is on structures that are frozen in time and so do not show the natural motion of the atoms in those 

With the only exception of liquid helium, quantum effects are de facto negligible in liquids and gases because the thermal wavelength X = h/-s/2 7T mkB T (with the Planck constant h, the Boltzmann constant k, the molecular mass m, and the temperature T) is smaller than the average distance of the molecules. Therefore, one can describe the MD of the molecules in a classical fluid by Newton s equations of motion 

Here r (t) denotes the position of molecule n with mass m at time t, and F is the force exerted on this molecule by all other molecules in the fluid as well as by the molecules forming potential confinements, such as channel walls or substrates. These equations of motion are coupled ordinary differential equations, which in MD simulations are solved numerically.  

The main challenge for MD simulations is the large number of degrees of freedom. For each molecule and in each numerical time step, the interactions with all other molecules have to be calculated, which requires in the order of IV (IV —1)/2 force calculations in total. Even on modern supercomputers this would limit the simulations to a rather small number of molecules. Therefore, in such simulations, often the range of intermolecular interactions is truncated such that the molecules interact only with the small number of neighboring molecules. This implies that the number offeree calculations in each numerical time step is proportional to N rather than N. However, 

To date MD is the most versatile and best-developed simulation method for nonequilibrium systems in which the discrete size of the molecules is important As a bonus the smaller the size of the system under consideration, the less problematic are most of the challenging issues discussed earlier. 

Reven and coworkers showed that it is possible to tune the chain mobility and related properties of a PMAA-based complex with great precision by selecting the correct polymer partners. These may range from hydrophilic polymers (e.g., PEO or PAAM) with fast dynamics where a rapid isotropic motion is dominant, to partner polymers (e.g., PVPon or PVCL) with significant hydrophobic character and high Tg and slow dynamics 

This point being made, we have not yet provided a description of how to follow a phase-space trajectory. This is the subject of molecular dynamics, upon which we now focus. 

There is an extensive body of literature covering the nature of physical interactions between molecules. The field of study is known as molecular dynamics, and it relies heavily on theories of statistical mechanics. 

The purpose of this section is to present only those aspects of molecular dynamics which are useful in providing a theoretical rationalization for the semiempirical Redlich-Kwong equation of state and for the algebraic form of the combining rules for mixtures. No attempt is made to derive the relations presented in this section. For a rigorous treatment the reader should consult the published works on the siibject, such as Prausnitz (1969) or Hirschfelder et al. (1964), which have supplied most of the information summarized in this section. 

The forces between an isolated pair of molecules can be separated into the repulsive force which prevents the molecules from merging together, and the attractive force, which among other things, is responsible for holding molecules together in liquids. The total potential energy (F) in a system of two molecules can be written as  

This was first proposed by G. Mie in 1903 and is referred to as the Mie equation. A schematic representation of the potential function is shown in Fig. 2. The repulsive potential is not well understood but it is clear that it increases very rapidly as the intermolecular separation becomes small. The repulsive potential is commonly reresented as 

In the above equation, the repulsive potential is represented by the twelfth power term and the attractive potential by the sixth power term, consequently the equation is known as a 6-12 Lennard- 

Up to this point, all electronic-structure methods have been concerned with the movement of the electrons only. This approach was justified by the idea of light electrons moving in the field of the fixed nuclei with large masses, and this simplification is a very reasonable one for most solids, especially if the properties to be predicted do not depend on the atomic movement. For example, to theoretically characterize insulators, semiconductors, and metals, atomic dynamics is imimportant likewise, magnetic phenomena (linked to the electronic spin) are largely independent of atomic movement. 

Whenever dynamical properties are of interest, the motions of the nuclei must be taken care of from the very beginning, but the sheer computational complexity of describing both electrons and nuclei renders a full quantum-mechanical approach extremely difficult, even today. Historically, two differ- 

The other approach, called molecular dynamics [249], is essentially based on Newton s mechanics utilized to describe the nuclear motions. Thus, quantum mechanics is discarded (at least for the moment), and classical mechanics is bravely applied to a system of atoms, ions, or molecules pretending these are not quantum objects. For molecules, such classical parametrizations go by the name force field methods, and they make up an important class of computational chemistry [58]. Given the knowledge of an atom s mass m and the force F(f) acting on it at some time t, the atom s acceleration a(t) is calculated according to F(f) = ma i) the force itself is the derivative of the potential energy V(r(f)) such that we have 

An ingenious method of extending the molecular dynamics approach when incorporating quantum mechanics, is given by the scheme of Car and Parri-nello [255]. Their method combines classical MD with parameter-free quantum mechanics to result in an ab initio molecular dynamics method, sometimes also called dynamical simulated annealing. In this approach, the nuclei are treated as classical objects (the Bom-Oppenheimer approximation is still valid) but the electrons are understood from density-functional theory. As we have seen already, the interactions between the electrons and the nuclei can be described satisfactorily by pseudopotentials (see Section 2.15.2), together with plane-wave basis sets, supercells (see Section 2.18), and periodic boundary conditions. 

The essential trick of the Car-Parrinello method is not to parametrize the potential-energy surface, but to calculate it on the fly from first principles when generating the nuclear trajectories by Newton s mechanics. The electronic groimd state is calculated at the same time. To do so, one starts with an extended Lagrangian which corresponds to a fictitious dynamical system 

Once the mathematical models (i.e., force fields) for the internal structure of each molecule (i.e., the intramolecular potential) and the interaction between molecules (i.e., the intermolecular potential) are known through classical statistical mechanics, one can predict the properties of a macroscopic sample of such molecules based on statistical averaging over the possible microscopic states of the system as it evolves under the rules of classical mechanics (Chandler 1987 McQuarrie 1976 Wilde and Singh 1998). Two of the most common techniques used are molecular dynamics (MD) (Alder and Wainwright 1957, 1959 McCammon et al. 1977 Rahman 1964 Stillinger and Rahman 1974) and Monte Carlo (MC) (Milik and Skolnick 1993 Ojeda et al. 2009) simulations. 

In its simplest form, MD solves Newton s familiar equation of motion  

The classical equations of motion are deterministic therefore once the initial coordinates and velocities are known, the coordinates and velocities can be determined later. The coordinates and velocities for a complete dynamics run are called the trajectory. Thus solving the classical equations of motion as a function of time (typically over a period limited to tens of nanoseconds) generates the microscopic states of the system. The systan may thus relax to equilibrium (provided the time for the relaxation falls within the time accessible to MD simulation), leading to the extraction of transport properties, which at the macroscopic scale describe the relaxation of the system in response to inhomogeneities. 

MC simulations generate equilibrium configurations stochastically according to the probabilities rigorously known from statistical mechanics. Because it generates equilibrium states directly, one can use it to study the equilibrium configurations of systems that may be expensive or impossible to access via MD. The drawback of MC is that it cannot yield the kind of dynamic response information that leads directly to transport properties. 

In Chapter 2, a brief discussion of statistical mechanics was presented. Statistical mechanics provides, in theory, a means for determining physical properties that are associated with not one molecule at one geometry, but rather, a macroscopic sample of the bulk liquid, solid, and so on. This is the net result of the properties of many molecules in many conformations, energy states, and the like. In practice, the difficult part of this process is not the statistical mechanics, but obtaining all the information about possible energy levels, conformations, and so on. Molecular dynamics (MD) and Monte Carlo (MC) simulations are two methods for obtaining this information 

Molecular dynamics is a simulation of the time-dependent behavior of a molecular system, such as vibrational motion or Brownian motion. It requires a way to compute the energy of the system, most often using a molecular mechanics calculation. This energy expression is used to compute the forces on the atoms for any given geometry. The steps in a molecular dynamics simulation of an equilibrium system are as follows  

Choose initial positions for the atoms. For a molecule, this is whatever geometry is available, not necessarily an optimized geometry. For liquid simulations, the molecules are often started out on a lattice. For solvent-solute systems, the solute is often placed in the center of a collection of solvent molecules, with positions obtained from a simulation of the neat solvent. 

Choose an initial set of atom velocities. These are usually chosen to obey a Boltzmann distribution for some temperature, then normalized so that the net momentum for the entire system is zero (it is not a flowing system). 

Compute the momentum of each atom from its velocity and mass. 

In Chapter 1,1 discussed the concept of mutual potential energy and demonstrated its relationship to that of force. So, for example, the mutual potential energy of the diatomic molecule discussed in Section 1.1.2 is 

This is a general rule we differentiate U by the coordinates of the particle in question in order to recover the force on that particle, from the expression for the mutual potential energy. Knowing the force F, we can use Newton s second law 

If the position (r), velocity (v), acceleration (a) and time derivative of the acceleration (b) are known at time t, then these quantities can be obtained dX. t + 8t by a Taylor expansion  

One way to do this is afforded by the predictor-corrector method. We ignore terms higher than those shown explicitly, and calculate the predicted terms starting with bP(t). However, this procedure will not give the correct trajectory because we have not included the force law. This is done at the corrector step. We calculate from the new position rP the force at time t + St and hence the correct acceleration a (t -f 5t). This can be compared with the predicted acceleration aP(f -I- St) to estimate the size of the error in the prediction step 

This error, and the results from the predictor step, are fed into the corrector step to give 

Both MD and MC methods employ a temperature as a guiding parameter for generating new geometries. At sufficiently high temperature and long run time, all the 

MD simulation, we need a reliable numerical algorithm for integrating Newton s equations. 

Because time is explicitly present in the formulations of MD, this technique is the most straightforward way of computer simulating the motion of penetrant molecules in amorphous polymer matrices (97-99). The MD method allows one to look at a truly atomistic level within the system as it evolves in time. Recently, excellent reviews on the use of MD for simulating penetrant diffusion in polymers have been published (96-99). A summary of the basic concepts and some relevant results obtained so far with MD will be presented bellow. 

As announced above these findings are in astonishing agreement with the heuristic pictures of the diffusion mechanism discussed in the framework of some microscopic diffusion models. But, besides being free of the conceptual drawbacks (the ad hoc assumptions) of the classical diffusion models, the MD method of computer simulation of diffusion in polymers makes it possible to get an even closer look at the diffusion mechanism and explain from a true atomistic level well known experimental findings. For example the results reported in (119,120) on the hopping mechanism reveal the following additional features. 

In a rubbery polymer with flexible macromolecular chains (PDMS for example) the cavities forming the free-volume are clearly separated from each other. The detailed evaluation of the movement of a penetrant particle from cavity (1) to the neighboring (2), did not show any immediate back jumps (2) — (1). This is mainly do to the fact that the channel between (1) and (2) closes quiet quickly. In a polymer with stiff chains (glassy polyimide (PI) for example) the individual cavities are closer to each other and a rather large number of immediate back jumps ocurred during the time interval simulated (120). This indicates that once a channel between two adjacent cavities in a stiff chain polymer is formed it will stay open for some 100 ps. This makes the back jump (2) - (1) of the penetrant more probable than a jump to any other adjacent hole (3). This process seems to be one cause for the general tendency that the diffusion coefficient of small penetrants in stiff chain glassy polymers is smaller than in flexible chain rubbery polymers. 

The results of MD simulations will be useful if they are able to reproduce with sufficient accuracy diffusion coefficients measured experimentally. Given the scatter between the results of different experiments reported in the literature, a computational method can be considered accurate enough if, for absolute diffusion coefficients, it reproduces the experimental values within one order of magnitude. Such results are presented in Table 5-1. 

The results given in Table 5-1 show that the agreement between the diffusion coefficients predicted from MD simulations and experimental ones ranges from reasonable to excellent. At temperatures around 300 K this is found both for polymers which are above their glass transition temperature, Tg, (PDMS, PIB, PE and aPP) and for polymers which are below Tg (PET, PS, PTMSP, PI and PAI). As a trend one can notice, and this not only from Tab. 5-1 but also from other works published in the last six or seven years, that the agreement between MD simulations of diffusion and solvation of small penetrants in polymers and experiment steadily improved. These are encouraging developments, showing that modern softwares (some of them available for 

Taking the limit of the vibrational partition function as the frequency approaches zero gives us 

The product of the frequency of crossing the energy barrier, v, and the cotTiceix-tration of the activated complex isjust the rate of reaction [(mol/dm ) X (1/s)] 

The term kTth is the order of magnitude typically found for the frequency factor A. At 300 K, kTIh = 6.25 X 10 s . The task now is to evaluate the partition functions, and techniques for doing this evaluation can be found in Laidler and are beyond the scope of this discussion. 

The theory above does not fully take into consideration the energy states of the reacting molecules and the offset of the colhsion as measured by the impact parameter. To accormt for these parameters, Karplus carried out a number of dynamic simulations to calculate collision trajectories. 

Figine G-5 shows the distances from separation / ab ac and Rbc as a function of tipie as A and BC approach each other and then separate (e g. / ab is the distance of separation between species A and species B). One also notes in this figure the vibration of the B— molecule. 

Such studies are even more demanding of computing skill and time than Monte Carlo calculations. 

Molecular dynamics is open to all the possibilities listed for the Monte Carlo calculations. One very important feature of these techniques is that the limitations involved are not limitations forced onto the model or onto the theory they are not physical but are computer limitations. 

There has been a large amount of work published in this field, and there is a vast potential and excitement in computer simulation techniques, and they can equally well be applied to other areas of chemistry, physics and biology. 

However, because the mathematical aspects of these developments are so complex, only a qualitative description of each will be given. 

The fundamental requirement for any computer simulation study is a knowledge of the interactions which are being included in the physical model of the electrolyte solution. The expression for the potential energy for each interaction and its dependence on distance, or distances, between the particles chosen is set up. The simplest calculations involve interactions where the potential energy is calculated pairwise, i.e. only interactions between two 

The most simple description for the electronically diabadc processes was formulated independently by Landau and by Zener. Following the model of these authors, the probability of passing from one surface to the other, considered as unidimensional, in the avoided crossing region, is given by 

The first studies on molecular dynamics were made in the early 1930s by Michael Polanyi and co-workers. In these smdies, very dilute flames and diffiisional flames were used to smdy the reactions of alkali metals with halogens and halides in the vapour phase, hi many cases, chemiluminescence was observed. For example, the reaction 

It was only in the 1960s that the pioneering work of M.G. Evans, M. Polanyi and also of H. Eyring, was extended with the arrival of new techniques, which allowed experimental determination of energies of products, and the introduction of computers capable of effectively solving the equations of motion. Some experimental observations have been rationalised by J.C. Polanyi in terms of localisation of the barrier in 

Relative distributions of vibrational, rotational and translational energy of products in gas-phase atom abstraction reactions [9] 

Although these ideas of J. C. Polanyi were well received at the time, and he was awarded the Nobel Prize in chemistry in 1986, it does not appear possible to make generalisations with them for all triatomic systems. Eor example, the reaction H-1-F2 appears to have an anomalous behaviour being highly exothermic, but with particularly small translational energy (Table 5.1). 

---

Two simulation methods—Monte Carlo and molecular dynamics—allow calculation of the density profile and pressure difference of Eq. III-44 across the vapor-liquid interface [64, 65]. In the former method, the initial system consists of N molecules in assumed positions. An intermolecule potential function is chosen, such as the Lennard-Jones potential, and the positions are randomly varied until the energy of the system is at a minimum. The resulting configuration is taken to be the equilibrium one. In the molecular dynamics approach, the N molecules are given initial positions and velocities and the equations of motion are solved to follow the ensuing collisions until the set shows constant time-average thermodynamic properties. Both methods are computer intensive yet widely used. 

In Fig. III-7 we show a molecular dynamics computation for the density profile and pressure difference P - p across the interface of an argonlike system [66] (see also Refs. 67, 68 and citations therein). Similar calculations have been made of 5 in Eq. III-20 [69, 70]. Monte Carlo calculations of the density profile of the vapor-liquid interface of magnesium how stratification penetrating about three atomic diameters into the liquid [71]. Experimental measurement of the transverse structure of the vapor-liquid interface of mercury and gallium showed structures that were indistinguishable from that of the bulk fluids [72, 73]. 

There is, of course, a mass of rather direct evidence on orientation at the liquid-vapor interface, much of which is at least implicit in this chapter and in Chapter IV. The methods of statistical mechanics are applicable to the calculation of surface orientation of assymmetric molecules, usually by introducing an angular dependence to the inter-molecular potential function (see Refs. 67, 68, 77 as examples). Widom has applied a mean-held approximation to a lattice model to predict the tendency of AB molecules to adsorb and orient perpendicular to the interface between phases of AA and BB [78]. In the case of water, a molecular dynamics calculation concluded that the surface dipole density corresponded to a tendency for surface-OH groups to point toward the vapor phase [79]. 

The theoretical treatments of Section III-2B have been used to calculate interfacial tensions of solutions using suitable interaction potential functions. Thus Gubbins and co-workers [88] report a molecular dynamics calculation of the surface tension of a solution of A and B molecules obeying Eq. III-46 with o,bb/ o,aa = 0.4 and... 

Both the Monte Carlo and the molecular dynamics methods (see Section III-2B) have been used to obtain theoretical density-versus-depth profiles for a hypothetical liquid-vapor interface. Rice and co-workers (see Refs. 72 and 121) have found that density along the normal to the surface tends to be a... 

It was noted in connection with Eq. III-56 that molecular dynamics calculations can be made for a liquid mixture of rare gas-like atoms to obtain surface tension versus composition. The same calculation also gives the variation of density for each species across the interface [88], as illustrated in Fig. Ill-13b. The density profiles allow a calculation, of course, of the surface excess quantities. 

Since the development of grazing incidence x-ray diffraction, much of the convincing evidence for long-range positional order in layers has come from this technique. Structural relaxations from distorted hexagonal structure toward a relaxed array have been seen in heneicosanol [215]. Rice and co-workers combine grazing incidence x-ray diffraction with molecular dynamics simulations to understand several ordering transitions [178,215-219]. 

Molecular dynamics calculations have been made on atomic crystals using a Lennard-Jones potential. These have to be done near the melting point in order for the iterations not to be too lengthy and have yielded density functioi). as one passes through the solid-vapor interface (see Ref. 45). The calculations showed considerable mobility in the surface region, amounting to the presence of a... 

Molecular dynamics and density functional theory studies (see Section IX-2) of the Lennard-Jones 6-12 system determine the interfacial tension for the solid-liquid and solid-vapor interfaces [47-49]. The dimensionless interfacial tension ya /kT, where a is the Lennard-Jones molecular size, increases from about 0.83 for the solid-liquid interface to 2.38 for the solid-vapor at the triple point [49], reflecting the large energy associated with a solid-vapor interface. 

Small metal clusters are also of interest because of their importance in catalysis. Despite the fact that small clusters should consist of mostly surface atoms, measurement of the photon ionization threshold for Hg clusters suggest that a transition from van der Waals to metallic properties occurs in the range of 20-70 atoms per cluster [88] and near-bulk magnetic properties are expected for Ni, Pd, and Pt clusters of only 13 atoms [89] Theoretical calculations on Sin and other semiconductors predict that the stmcture reflects the bulk lattice for 1000 atoms but the bulk electronic wave functions are not obtained [90]. Bartell and co-workers [91] study beams of molecular clusters with electron dirfraction and molecular dynamics simulations and find new phases not observed in the bulk. Bulk models appear to be valid for their clusters of several thousand atoms (see Section IX-3). 

Bartell and co-workers have made significant progress by combining electron diffraction studies from beams of molecular clusters with molecular dynamics simulations [14, 51, 52]. Due to their small volumes, deep supercoolings can be attained in cluster beams however, the temperature is not easily controlled. The rapid nucleation that ensues can produce new phases not observed in the bulk [14]. Despite the concern about the appropriateness of the classic model for small clusters, its application appears to be valid in several cases [51]. 

The Langmuir-Hinshelwood picture is essentially that of Fig. XVIII-14. If the process is unimolecular, the species meanders around on the surface until it receives the activation energy to go over to product(s), which then desorb. If the process is bimolecular, two species diffuse around until a reactive encounter occurs. The reaction will be diffusion controlled if it occurs on every encounter (see Ref. 211) the theory of surface diffusional encounters has been treated (see Ref. 212) the subject may also be approached by means of Monte Carlo/molecular dynamics techniques [213]. In the case of activated bimolecular reactions, however, there will in general be many encounters before the reactive one, and the rate law for the surface reaction is generally written by analogy to the mass action law for solutions. That is, for a bimolecular process, the rate is taken to be proportional to the product of the two surface concentrations. It is interesting, however, that essentially the same rate law is obtained if the adsorption is strictly localized and species react only if they happen to adsorb on adjacent sites (note Ref. 214). (The apparent rate law, that is, the rate law in terms of gas pressures, depends on the form of the adsorption isotherm, as discussed in the next section.)... 

Because of limitations of space, this section concentrates very little on rotational motion and its interaction with the vibrations of a molecule. However, this is an extremely important aspect of molecular dynamics of long-standing interest, and with development of new methods it is the focus of mtense investigation [18, 19, 20. 21. 22 and 23]. One very interesting aspect of rotation-vibration dynamics involving geometric phases is addressed in section A1.2.20. 

We have alluded to the comrection between the molecular PES and the spectroscopic Hamiltonian. These are two very different representations of the molecular Hamiltonian, yet both are supposed to describe the same molecular dynamics. Furthemrore, the PES often is obtained via ab initio quairtum mechanical calculations while the spectroscopic Hamiltonian is most often obtained by an empirical fit to an experimental spectrum. Is there a direct link between these two seemingly very different ways of apprehending the molecular Hamiltonian and dynamics And if so, how consistent are these two distinct ways of viewing the molecule ... 

There are also approaches [, and M] to control that have had marked success and which do not rely on quantum mechanical coherence. These approaches typically rely explicitly on a knowledge of the internal molecular dynamics, both in the design of the experiment and in the achievement of control. So far, these approaches have exploited only implicitly the very simplest types of bifiircation phenomena, such as the transition from local to nonnal stretch modes. If fiittlier success is achieved along these lines m larger molecules, it seems likely that deliberate knowledge and exploitation of more complicated bifiircation phenomena will be a matter of necessity. 

Davis M J 1995 Trees from spectra generation, analysis, and energy transfer information Molecular Dynamics and Spectroscopy by Stimulated Emission Pumping ed H-L Dai and R W Field (Singapore World Scientific)... 

It is possible to use the quantum states to predict the electronic properties of the melt. A typical procedure is to implement molecular dynamics simulations for the liquid, which pemiit the wavefiinctions to be detemiined at each time step of the simulation. As an example, one can use the eigenpairs for a given atomic configuration to calculate the optical conductivity. The real part of tire conductivity can be expressed as... 

Figure Al.3.30. Theoretical frequency-dependent conductivity for GaAs and CdTe liquids from ab initio molecular dynamics simulations [42].

Yan Y J, Gillilan R E, Whitnell R M, Wilson K R and Mukamel S 1993 Optimal control of molecular dynamics -Liouville space theory J. Chem. Phys. 97 2320... 

Rice S A and Zhao M 2000 Optical Control of Molecular Dynamics (New York Wiley)... 

In equilibrium statistical mechanics, one is concerned with the thennodynamic and other macroscopic properties of matter. The aim is to derive these properties from the laws of molecular dynamics and thus create a link between microscopic molecular motion and thennodynamic behaviour. A typical macroscopic system is composed of a large number A of molecules occupying a volume V which is large compared to that occupied by a molecule ... 

Alavi A 1996 Path integrals and ab initio molecular dynamics Monte Carlo and Molecular Dynamics of Condensed Matter Systems ed K Binder and G Ciccotti (Bologna SIF)... 

The alternative simulation approaches are based on molecular dynamics calculations. This is conceptually simpler that the Monte Carlo method the equations of motion are solved for a system of A molecules, and periodic boundary conditions are again imposed. This method pennits both the equilibrium and transport properties of the system to be evaluated, essentially by numerically solvmg the equations of motion... 

If we wish to know the number of (VpV)-collisions that actually take place in this small time interval, we need to know exactly where each particle is located and then follow the motion of all the particles from time tto time t+ bt. In fact, this is what is done in computer simulated molecular dynamics. We wish to avoid this exact specification of the particle trajectories, and instead carry out a plausible argument for the computation of r To do this, Boltzmann made the following assumption, called the Stosszahlansatz, which we encountered already in the calculation of the mean free path ... 

Radiation probes such as neutrons, x-rays and visible light are used to see the structure of physical systems tlirough elastic scattering experunents. Inelastic scattering experiments measure both the structural and dynamical correlations that exist in a physical system. For a system which is in thennodynamic equilibrium, the molecular dynamics create spatio-temporal correlations which are the manifestation of themial fluctuations around the equilibrium state. For a condensed phase system, dynamical correlations are intimately linked to its structure. For systems in equilibrium, linear response tiieory is an appropriate framework to use to inquire on the spatio-temporal correlations resulting from thennodynamic fluctuations. Appropriate response and correlation functions emerge naturally in this framework, and the role of theory is to understand these correlation fiinctions from first principles. This is the subject of section A3.3.2. 

Progress in the theoretical description of reaction rates in solution of course correlates strongly with that in other theoretical disciplines, in particular those which have profited most from the enonnous advances in computing power such as quantum chemistry and equilibrium as well as non-equilibrium statistical mechanics of liquid solutions where Monte Carlo and molecular dynamics simulations in many cases have taken on the traditional role of experunents, as they allow the detailed investigation of the influence of intra- and intemiolecular potential parameters on the microscopic dynamics not accessible to measurements in the laboratory. No attempt, however, will be made here to address these areas in more than a cursory way, and the interested reader is referred to the corresponding chapters of the encyclopedia. 

Specific solute-solvent interactions involving the first solvation shell only can be treated in detail by discrete solvent models. The various approaches like point charge models, siipennoleciilar calculations, quantum theories of reactions in solution, and their implementations in Monte Carlo methods and molecular dynamics simulations like the Car-Parrinello method are discussed elsewhere in this encyclopedia. Here only some points will be briefly mentioned that seem of relevance for later sections. 

Predicting the solvent or density dependence of rate constants by equation (A3.6.29) or equation (A3.6.31) requires the same ingredients as the calculation of TST rate constants plus an estimate of and a suitable model for the friction coefficient y and its density dependence. While in the framework of molecular dynamics simulations it may be worthwhile to numerically calculate friction coefficients from the average of the relevant time correlation fiinctions, for practical purposes in the analysis of kinetic data it is much more convenient and instructive to use experimentally detemiined macroscopic solvent parameters. 

---