Reviews in Computational Chemistry

204 Spin–Orbit Coupling in Molecules

180.J. Tatchen, Anwendung von Mean-Field Operatoren zur Beschreibung von Spin-Bahn– Effekten in organischen Moleku¨ len, Diploma Thesis, University of Bonn, 1999. URL www.thch.uni-bonn.de/tc/.

181.J. Tatchen, C. M. Marian, M. Waletzke, and S. Grimme, to be published. A Quantum Chemical Study of Spin-Forbidden Transitions in Dithiosuccinimide.

182.B. Minaev, Int. J. Quantum Chem., 17, 367 (1980). Intensities of Spin-Forbidden Transitions in Molecular Oxygen and Selective Heavy-Atom Effects.

183.R. Klotz, C. M. Marian, S. D. Peyerimhoff, B. A. Hess, and R. J. Buenker, Chem. Phys., 89,

223 (1984). Calculation of Spin-Forbidden Radiative Transitions Using Correlated Wave Functions: Lifetimes of b1 þ, a1 States in O2, S2 and SO.

184.H. Ko¨ ppel and W. Domcke, in Encyclopedia of Computational Chemistry, P. v. R. Schleyer, Ed., Wiley, Chichester, 1997, pp. 152–278. Vibronic Dynamics of Polyatomic Molecules.

185.B. Gemein and S. D. Peyerimhoff, Chem. Phys. Lett., 184, 45 (1991). Theoretical Study of the Vibronic Interactions in the Ground and First Excited a3 and a0 3 þ States of the CO Molecule.

186.J. R. Taylor, Scattering Theory: The Quantum Theory on Nonrelativistic Collisions, Wiley, New York, 1972.

¨

187. G. Wentzel, Z. Phys., 43, 524 (1927). Uber strahlungslose Quantenspru¨ nge.

188. G. Wentzel, Phys. Z., 29, 321 (1928). Die unperiodischen Vorga¨nge in der Wellenmechanik.

189. O. K. Rice, Phys. Rev., 33, 748 (1929). Perturbations in Molecules and the Theory of Predissociation and Diffuse Spectra.

190. O. K. Rice, Phys. Rev., 35, 1551 (1930). Perturbations in Molecules and the Theory of Predissociation and Diffuse Spectra. II.

191. J. N. Murrell and J. M. Taylor, Mol. Phys., 16, 609 (1969). Predissociation in Diatomic Spectra with Special Reference to the Schumann–Runge Bands of O2.

192. C. M. Marian, R. Marian, S. D. Peyerimhoff, B. A. Hess, R. J. Buenker, and G. Seger, Mol. Phys., 46, 779 (1982). Ab-Initio Study of Oþ2 Predissociation Phenomena Induced by a Spin– Orbit Coupling Mechanism.

193. D. Edvardsson, S. Lunell, F. Rakowitz, C. M. Marian, and L. Karlsson, Chem. Phys., 229, 203 (1998). Calculation of Predissociation Rates in O22þ by Ab Initio MRD–CI Methods.

194. J. Scho¨ n and H. Ko¨ ppel, J. Chem. Phys., 108, 1503 (1998). Geometric Phases and Quantum Dynamics in Spin–Orbit Coupled Systems.

195. J. Ro¨ melt and R. Runau, Theor. Chim. Acta, 54, 171 (1980). Franck–Condon Matrix Elements for Bound-Continuum Vibrational Transitions Calculated by Numerical Integration and Basis Set Expansion Techniques.

196. D. O. Harris, G. G. Engerholm and W. D. Gwinn, J. Chem. Phys., 43, 1515 (1965). Calculation of Matrix Elements for One-Dimensional Quantum-Mechanical Problems and the Application to Anharmonic Oscillators.

197. Note Added in Proofs: For the evaluation of spin–orbit matrix elements, standard basis sets of at least double-zeta plus polarization quality should be employed. See, e.g., R. Klotz, C. M. Marian, S. D. Peyerimhoff, B. A. Hess, and R. J. Buenker, Chem. Phys., 76, 367 (1983). Study of the Dependence of Spin–Orbit Matrix Elements on AO Basis Set Composition for Inner and Valence Shells: Results for the Multiplet Spiltting of X 3 P and C 3 of SO and X 2 in SOþ. A proper representation of the inner and outer nodes of the valence shell orbitals is essential. The representation of the inner core orbitals was found to be less critical. Scalar relativistic all-electron calculations require a modification of standard Gaussian basis set contractions. For elements H–Kr, scalar relativistic recontractions are available from http://www.emsl.pnl.gov:2080/forms/basisform.html. For the Douglas–Kroll recontracted basis sets, see: W. A. de Jong, R. J. Harrison, and D. A. Dixon, J. Chem. Phys., 114, 48 (2001). Parallel Douglas–Kroll Energy and Gradients in NWChem: Estimating Scalar Relativistic Effects Using Douglas–Kroll Contracted Basis Sets.

206 Cellular Automata Models of Aqueous Solution Systems

are clear limits to this approach, especially when one seeks to explore many of the less detailed, broader aspects of the observations, such as the variations in species populations with time and the statistical and kinetic details of the phenomenon in question.

In the second general approach to this problem, an attempt is made to examine in some manner the overall behavior of the entire ensemble of interacting units. By far, the most common approach here, and the one normally taught from textbooks, is to represent the kinetic behavior of a particular system in terms of an applicable set of coupled differential rate equations. These equations, with their associated rate constants, summarize the bulk behaviors of the ingredients involved in an averaged way. For example, the simple twostep transformation A ! B ! C, can be characterized by the set of rate equations:

where k1 and k2 are the rate constants for the two steps. This formulation can be referred to as the ‘‘traditional’’ approach. Solution of these equations yields the time-dependent populations of the species involved. Although analytical solutions may be impossible to obtain for some very complicated systems of equations, powerful approximate numerical methods make such exact, analytical solutions unnecessary. There are, however, several important—and often unrecognized—limitations of this traditional approach, in addition to its occasional cumbersomeness. The first limitation, which is more a formal than a practical problem, is that despite its widespread employment and obvious utility, such an approach has no solid theoretical justification; its use rests on empiricism.1 A more practical limitation is that a deterministic approach of this nature is truly applicable only at the limit in which an exceedingly large number of interacting particles is involved. Furthermore, this approach gives no information on the nature of the inherent, and often important, fluctuations that can be expected to arise in the behaviors of real, finite systems. Since recent advances in analytical technology have pushed levels of experimental observations to very small samples and even to the single molecule level, the ability to deal with finite samples has become increasingly important.

Stochastic analysis presents an alternative avenue for dealing with the inherently probabilistic and discontinuous microscopic events that underlie macroscopic phenomena. Many processes of chemical and physical interest can be described as random Markov processes.1,2 Unfortunately, solution of a stochastic master equation can present an extremely difficult mathematical challenge for systems of even modest complexity. In response to this difficulty, Gillespie3–5 developed an approach employing numerical Monte Carlo

Introduction 207

techniques that accurately describes the stochastic features of such processes, and several applications have demonstrated the usefulness of this approach.6,7 There is an intermediate level of consideration of chemical phenomena laying somewhere between the intense scrutiny of a single molecule and the averaged, possibly superficial, treatment of a bulk sample with an effectively infinite population. This intermediate level is the level of chemical events for which the phenomena of interest can be modeled or simulated with acceptable accuracy by including a finite number of molecules (greater than one but far less than Avogadro’s number). For example, the behavior of a solution cannot be modeled accurately with just one or two or even a few of molecules. Ideally, however, if we could examine or model a few tens, a hundred, or a thousand ingredients, then some significant insight into the nature of the solution might be achieved. This intermediate level model might, then, be reasonably simpler to create and still give meaningful insights into the phenomena

of interest.

The past three decades have witnessed the development of three broad techniques—molecular dynamics (MD), Monte Carlo (MC), and cellular automata simulations—that approach the study of molecular systems by simulating submicroscopic chemical events at this intermediate level. All three methods focus attention on a modest number of molecules and portray chemical phenomena as being dependent on dynamic, and interactive events (a portrayal consistent with our scientific intuition and a characteristic not intrinsic to either thermodynamics or the traditional deterministic approach based on differential equations). These techniques lend themselves to a visual portrayal of the evolution of the configurations of the systems under study. Because each approach has its own particular advantages and shortcomings, one must take into consideration the pros and cons of each, especially in light of the nature of the problem to be solved.

Both MD and MC methods require either quantum mechanical or force field calculations to determine the energies. Normally, for large systems one must begin with a choice of a force field or potential functions (namely, a set of equations representing the electrostatic, steric, and other interactions that can be expected to apply within the system studied). Regardless of how these energies are computed, in an MD calculation Newton’s laws of motion are applied in a linearized form using a very short (e.g., 0.5 femtosecond) time step. Positions of the atoms are determined at each time step, and the evolution of the system is followed over the time period of the simulation. In an MC calculation, small random changes are made in the positions of the particles in the system. The configurations produced in this way are evaluated and weighted according to their energies as determined by a suitable quantum method or by the potential function chosen for the system. The probability of a change from one configuration to another can thus be determined; some changes can be accepted, and others rejected. Sampling and evaluation of the normally very large number (e.g., 107) of configurations produced in this way allows

208 Cellular Automata Models of Aqueous Solution Systems

estimation of the thermodynamic properties of the system studied. More detailed descriptions of these techniques can be found in recent reviews.8–11 A question arises as to the relationship of cellular automata models to the MC and MD models. Put briefly, cellular automata models are usually

much simpler both in concept and implementation than are MD or MC models. In results and procedures, cellular automata have much in common with the integrated MC technique.11 All three techniques provide visual expressions of the processes involved as well as information on fluctuations in the system’s behavior. The cellular automata models, however, omit the detailed molecular structures and energy calculating methods that are commonly employed in MC and MD simulations, substituting in their place heuristic rules for these features. As a result, cellular automata models are far less demanding computationally and much faster than either of the other techniques, but they do not directly reveal information about the energetics of the process or system. However, the use of rules in place of forces does have the special advantage of providing specific features that can be varied individually, and the consequences of those changes can be queried independent of other variations that are intimately coupled to them. For example, in a cellular automaton model one can independently vary the temperature of a solvent while keeping the temperature of its solute constant in order to assess the effect of environmental temperature, or one can increase one specific reaction or transition rate while holding others constant so as to isolate the influence of a particular process on the overall dynamics of a system.

In the following sections, we describe what cellular automata are, how they work, and then we provide examples of their uses. Because much of chemistry is carried out in condensed phases, especially in water, we limit our applications to this medium but point out that extensions to other types of systems are possible.

CELLULAR AUTOMATA

Historical Background

Cellular automata were first proposed by the mathematician Stanislaw Ulam and the mathematical physicist John von Neumann a half century ago,12–14 although related ideas were put forth even earlier, in the 1940s, by the German engineer Konrad Zuse.15–17 von Neumann’s interest was in the construction of ‘‘self-reproducing automata’’.18 His original idea was to construct a series of mechanical devices or ‘‘automata’’ that would gather and integrate the ingredients that could reproduce themselves. A suggestion by Ulaml3 led von Neumann to consider grids with moving ingredients, operating with rules. The first such system was made up of square cells in a matrix, each with a state, operating with a set of rules in a two-dimensional grid. The system existed with up to 29 different states. With the development of modern

Cellular Automata 209

digital computers, it became increasingly clear that these fairly abstract ideas

could in fact be usefully applied to the examination of real physical systems.19,20 As described by Wolfram,21 cellular automata have five fundamen-

tal defining characteristics:

They consist of a discrete lattice of cells. They evolve in discrete time steps.

Each site takes on a finite number of possible values.

The value of each site evolves according to the same deterministic rules.

The rules for the evolution of a site depend only on a local neighborhood of sites around it.

As we shall see, the fourth characteristic can be modified to include probabilistic rules as well as deterministic rules. An important feature sometimes observed in the evolution of these computational systems was the development of unanticipated patterns of ordered dynamical behavior, or ‘‘emergent properties’’. As Kauffman has expressed it,22 ‘‘Studies of large, randomly assembled cellular automata . . . have now demonstrated that such systems can spontaneously crystallize enormously ordered dynamical behavior. This crystallization hints that hitherto unexpected principles of order may be found, [and] that the order observed may have significant explanatory import in [biology and physics].’’ This proposal has borne considerable fruit, not only in biology and physics, but also in chemistry. Readers are referred to reviews for applications in physics23–26 and biology;27 selected physical and chemical applications have been reviewed by Chopard and co-workers.28

The General Structure

The simulation of a dynamic system using cellular automata requires several parts that make up the process. The cell is the basic model of each ingredient, molecule, or whatever constitutes the system. These cells may have several shapes as part of the matrix or grid of cells. The grid may have boundaries or be part of a topological object that eliminates boundaries. The cells may have rules that apply to all of the edges, or there may be different rules for each edge. This latter plan may impart more detail to the model, as needed for a more detailed study.

The Cells

Cellular automata have been designed for one (1D), two (2D), or threedimensional (3D) arrays. The most commonly used is the 2D grid. The cells may be triangles, squares, hexagons, or other shapes in the 2D grid. The square cell has been the one most widely used over the past 40 years. Each cell in the grid is endowed with a primary state, namely, whether it is empty or occupied with a particle, object, molecule, or whatever the system requires

210 Cellular Automata Models of Aqueous Solution Systems

A

Figure 1 A cellular automata grid showing occupied cells of different states A and B, and unoccupied cells (blank).

to study the dynamic event (see Figure 1). Secondary information is contained in the state description that encodes the differences among cell occupants in a study.

The Cell Shape

The choice of the cell shape is based on the objective of the study. In studies of water and solution phenomena, a square cell is appropriate because the water molecule is quadravalent to hydrogen bonding to other water molecules or solutes. A water molecule donates two hydrogens and two lone-pair electrons in forming the tetrahedral structure that characterizes the liquid state. The four faces of a square cell thus correspond to the bonding opportunities of a water molecule.

Grid Boundary Cells

The ‘‘moving cell’’ may encounter an edge or boundary during its movements. The boundary cell may be treated as any other occupied cell, following rules that permit joining or breaking. A more common practice is to assume that the grid is simulating a small segment of a large dynamic system. In this model, the boundaries do come into play in the results. The grid is then considered to be the surface of a torus (Figure 2); the planar projection of this surface would reveal the movement of a cell off of the edge and reappearing from the opposite edge onto the grid (see Figure 3). In some cases, it is necessary to establish a vertical relationship among occupants, which establishes a gravity

Cellular Automata 211

Figure 2 The grid located on the surface of a torus. Column A of the grid is represented by line A on the torus; row B on the grid is represented by line B on the torus.

effect. In this situation, the grid is on the surface of a cylinder with a boundary condition at the top and bottom, which is an impenetrable boundary.

Variegated Cell Types

Until recently, all of the cellular automata models assumed that each edge of a cell bore the same state and movement rules. Recent work has

Figure 3 Movement of cell occupants at the boundaries of the grid on a torus. The cell occupant A may move right and appear at the opposite side of the grid. The other cell occupant A may move down and appear at the top of the grid.

212 Cellular Automata Models of Aqueous Solution Systems

Figure 4 The possible variegated cells. (a) Possible occupants with two different types of edges allowing for different rules for the edges. (b) Possible occupants with three different types of edges hence possibly different rules for the edges. (c) Possible occupants with two different types of edges hence different rules for the edges.

employed a variegated cell where each edge may have its own state and set of movement rules. Examples of some variegated cells are shown in Figure 4.

Cell Movement

The dynamic character of cellular automata is developed by the simulation of movement of the cells, which may be a simultaneous process or each cell, in turn, may execute a movement. Each cell’s movement is based on rules derived from the states of other nearby cells. These nearby cells constitute a neighborhood. The rules may be deterministic or stochastic, the latter process being driven by probabilities of certain events occurring. These processes are elaborated upon in this section.

Cellular Automata 213

Synchronous–Asynchronous Movement

When we speak of movement of a cell or the movement of cell occupants, we are speaking of the simulation of a movement from one cell to another. Thus a water molecule or some object is postulated to move across space, appearing in a new location at time t þ 1. In the cellular automata models, the actual situation is the exchange of state between two adjacent cells. If we are modeling the movement of a molecule from place A to the adjacent place B, then we must exchange the states of cells A and B. Initially, at time t, cell A has a state corresponding to an occupant molecule, while adjacent cell B is devoid of a molecule (i.e., it is empty). At time t þ 1, the states of the cells A and B have exchanged. Cell A is empty, and cell B has the state of the occupant molecule. This exchange gives the illusion, and the practical consequences, of a movement from A to B. We speak of the movement of cells or of the movement of cell occupants; either way we are describing the process of simulating a movement as stated above.

Cell movement of all occupants in the grid may occur simultaneously (synchronous) or it may occur sequentially (asynchronous). When all cells in the grid have computed their state and have executed their movement (or not) it is called one iteration, a unit of time in the cellular automata simulation. Each cell is identified in the program and is selected randomly for the choice of movement or not. The question of which type of movement to use depends on the system being modeled and the information sought from the model, but it should reflect reality. If the system being studied is a slow process, then synchronous motion may best represent the process. In contrast, if the system is very fast, like proton hopping among water molecules where the cellular automata is using a few thousand cells, then an asynchronous model is desirable. A synchronous execution of the movement rules leads to possible competition for a cell from more than one occupant. A resolution scheme must thus be in place to resolve the competition; otherwise this may interfere with the validity of the model.

Neighborhoods

Cell movement is governed by rules called transition functions. The rules involve the immediate environment of the cell called the neighborhood. The most common neighborhood used in 2D cellular automata is called the von Neumann neighborhood after the pioneer of the method. A cell, i, is in the center of four cells, j, adjoining the four faces of i, [see Figure 5(a) for this pattern]. Another common neighborhood is the Moore neighborhood [Figure 5(b)], where cell i is completely surrounded by j cells. Other neighborhoods include the extended von Neumann neighborhood shown in Figure 5(c), where the k cells beyond j are identified and allowed to participate in movements.

Deterministic/Probabilistic Movement Rules

The rules governing cell movement may be deterministic or probabilistic. Deterministic cellular automata use a fixed set of rules, the values of which are