Simulation method

Simulation refer to models involving a statistical component, i.e. the results are obtained with an uncertainty arising from the finite length of the simulation. A typical application is simulation of a liquid phase, i.e. a detailed description of a solvent or solution. 

A macroscopic quantity X, such as the enthalpy or entropy, can be calculated from the partition function Q, as discussed in Chapter 12. For easy reference we reproduce eqs. (12.6-12.8) here. 

The properties have intentionally been separated in to two groups, those involving derivatives of Q and those which depend directly on Q. 

The partition function Q here describes the whole system consisting of N interacting particles, and the energy states Ei are consequently for all the particles (in Section 12.2 we considered N non-interacting molecules, where the total partition function could be written in terms of the partition function for one molecule, Q — /N[). More correctly 

SIMULATIONS, TIME-DEPENDENT METHODS AND SOLVATION MODELS 

For an isolated molecule in the rigid rotor, harmonic oscillator approximation, the (quantum) energy states are sufficiently regular to allow an explicit construction of the partition function, as discussed in Chapter 12. For a collection of many particles the 

Two major classical simulation techniques, molecular dynamics and Monte Carlo, have been applied to simulation of water-metal interfaces. We first discuss features common to both methodologies and then describe aspects unique to each. The field of computer simulations is an actively evolving one, despite being more than 40 years old. Even for the particular case of water-metal interfaces, many variations exist on the central theme of how best to carry out these calculations. In this chapter, we limit our discussion to the most significant (in our opinion) techniques in use for metal-water interfaces. 

The main goal of simulation methods is to obtain information on the spatial and temporal behavior of a complex system (a material), that is, on its structure and evolution. Simulation methods are subdivided into atomistic and phenomenological methods. Atomistic methods directly consider the evolution of the system of interest at the atomic level with regard to the microscopic structure of the substance. These methods include classical and quantum MD and various modifications of the MC technique. Phenomenological methods are based on macroscopic equations in which the atomistic nature of the material is not directly taken into account. Within the multiscale approach, both groups of methods mutually complement each other, which permits the physicochemical system under study to be described most comprehensively. 

The MD technique permits thin films to be investigated under both equilibrium and nonequilibrium conditions. Therefore, in addition to equilibrium properties, such as structure, density, etc., MD predicts the nonequilibrium behavior of thin films, for example, energy relaxation for an ion interacting 

As noted before, simulation methods may be considered to be intermediate between theory and experiment. In fact, various kinds of simulations have been devised for the purpose of computing g R) as well as other properties of the liquid. We now make a distinction between three of these methods. 

Morrell and Hildebrand (1934) described a simple and interesting method of computing g R) for a system of balls of macroscopic size. They used gelatin balls suspended in liquid gelatin that were mechanically shaken in their vessel. The coordinates of some of the balls could be measured by taking photographs of the system at various times. From these data, they were able to compute the distribution of distances between the balls, and from this, the radial distribution function (see Hildebrand et al., 1970). 

A similar method was employed by Bernal and his collaborators (Bernal and King, 1968), who used steel balls in a random configuration to compute g(R) (as well as other distribution functions of interest in understanding the random character of the liquid state). The form of g(R) obtained from 

A closely related experiment can be carried out on a computer. Instead of shaking a system of steel balls by mechanical means, one can generate random configurations of particles by computer. The latter method can be used to compute the RDF as well as other thermodynamic quantities of a system. 

The specific method devised by Metropolis et aL (1953) to compute the properties of liquids is now known as the Monte Carlo (MC) method. In fact, this is a special procedure to compute multidimensional integrals numerically. 

The quantity F can be approximated by a sum over discrete events (configurations) 

A glance at (5.5.13) reveals that two very severe difficulties arise in any attempt to approximate the average F by a sum of the form (5.5.14). First, a configuration means a specification of 3N coordinates (for spherical particles), and, clearly, such a number of coordinates cannot be handled in a computer if N is of the order of 10 . This limitation forces us to choose N of the order of a few hundred. The question of the suitability of such a small sample of molecules to represent a macroscopic system immediately arises. The second difficulty concerns the convergence of the sum in (5.5.14). Suppose we have already chosen N. The question is How many configurations n should we select in order to ensure the validity of the approximation in (5.5.14) Again, time limitations impose restrictions on the number n that we can afford in an actual experiment.  

The simulation techniques used for polyelectrolytes in solution are extensions of the standard methods used for neutral polymers. The polymer chain is modeled as a set of connected beads. The beads are charged depending on the charge fraction, but otherwise the details of the monomer structure are neglected. Various means of connecting the bonded monomers are used. In lattice Monte Carlo the bonds are of course fixed. Two sets of simulations have used the rotational isomeric state model. Other simulations have used Hookean springs or the finite-extendable-nonlinear-elastic (FENE) potential. No important dependence on the nature of the bonds is expected at this level of modeling the polymer chain. 

Treating the solvent is more complex for charged polymers in solution than for neutral polymers in solution. In neutral polymer simulations the solvent quality can be treated by effectively altering the monomer-monomer interaction. For polyelectrolytes in solution other considerations are also 

The long-range Coulomb interaction requires special treatment. 

The MI image has been chosen over the Ewald method, because for the system sizes studied the MI method is faster than the Ewald method. Even so, the MI method is slow compared to neutral system simulations. The main obstacle with charged systems is poor scaling with the number of particles. Since all particle pairs interact the computation time scales as N for the MI method. (For large enough N the Ewald method changes from scaling as N to Thus for even rather small system sizes the 

Following common practice, the expectation values in the derivation of the equations for free energy methods are written as ensemble averages. 

In this section, we will briefly describe the static and dynamic simulation methods and then give some results for the liquid of greatest geochemical importance water. Applications of such methods to mineral structures will be discussed in later chapters. 

Static simulations are generally based on an energy minimization procedure that is, the energy of the system is written as a function of structural variables that include atomic coordinates and cell dimensions. The 

Combined local-density-functional molecular dynamics approach 

An important new development within solid-state theory is the combination of self-consistent band structure, structure determination, and molecular dynamics within the local-density approximation as developed by Car and Parrinello (1985). Our discussion follows that of Srivastava and Weaire (1987). 

In the Car and Parrinello (1985) scheme, ion dynamics is combined with a fictitious classical electron dynamics, with nuclei assigned real masses and the electron wave functions arbitrary fictitious masses. One starts the molecular-dynamics simulation at high temperature and cools progressively to zero temperature to find the ground state of both electrons and ions simultaneously. Although this approach at first seems strange and unphysical, it has yielded excellent results for amorphous Si (Car and Parrinello, 1988) and recently for SiOj (Allan and Teter, 1987) and S clusters (Hohl et al., 1988) and will probably play an important role in the future development of the field. 

Most molecular simulation techniques can be categorized as being among three main types (1) quantum mechanics, (2) molecular dynamics (MD) and (3) kinetic Monte Carlo (KMC) simulation. Quantum mechanics methods, which include ah initio, semi-empirical and density functional techniques, are useful for understanding chemical mechanisms and estimating chemical kinetic parameters for gas-phase 

The probability distribution for each configuration is described by the Master equation [96]  

The chapter contains three main sections. In the first, we present a detailed but not exhaustive survey of simulation methods for the calculation of free energy, entropy, and heat capacity. This is followed by a section that discusses results from several hydrophobic hydration and hydrophobic interaction simulations. These applications are then used in the last section of this chapter to illustrate the simulation methods and to highlight several important conceptual developments in the theory of hydrophobicity. 

Molecular dynamics (MD) and Monte Carlo (MC) simulation techniques have been used now for decades to characterize aqueous solutions. The most basic elements that underlie these techniques, such as numerical integration algorithms and the Metropolis method, are discussed thoroughly else-where, so they are not included here. Our intention here is to survey methods used for determining the thermodynamic and structural quantities most closely tied to hydrophobicity. 

Most of the applications discussed below involve atomic solutes with either hard-sphere (HS) or Lennard-Jones (LJ) solute-water interaction potentials. The potential energy of a HS is infinite if the solute-oxygen distance is less than the sum of the solute and water HS radii and zero otherwise. An effective HS radius is assigned to water in some reasonable fashion. The solute-water LJ interaction has the form 

As was mentioned earlier the master equation (5.37) generally cannot be solved. To get some experience of the behaviour of chemical systems we might do stochastic simulation experiments using Monte-Carlo techniques (Introductions to Monte-Carlo methods are given in Hammersby Hand-scomb 1964, and Srejder 1965. Their applications in chemical physics are discussed in Binder (1979.) 

Based on a theorem of Doob (1953, pp. 244) we can assume that the waiting time of the system in state j can be considered as an exponentially 

Simulation experiments are performed by determining the duration while the system does not change its state, and then the reaction what occurs at the end of the interval. 

Set the time variable / = 0. Prescribe the initial numbers of molecules (a deterministic set of initial conditions is assumed). Specify the elementary reactions and their rate constants. Specify the stopping time . 

Generate two (pseudo)random numbers to select a time interval A/, and a reaction occurring at time / H- A/. 

The bond fluctuation model [72] is used to simulate the motion of the polymer chains on the lattice. In this model, each segment occupies eight lattice sites of a unit cell, and each site can be a part of only one segment (self-avoiding walk condition). This condition is necessary to account for the excluded volume of the polymer chains. For a given chain, the bond length between two successive seg- 

To monitor the time evolution of the long-range ordering, the collective structure factor at the constant time interval of phase separation is calculated by Eq. (20)  

The first moment q t) of the structure factor is also computed in order to observe the coarsening processes in the later stage of phase separation more clearly. This quantity q t), the inverse of which is a measure of the average domain size, is defines as 

For typical measure of domain size, one usually considers either the location qmax(t) of the peak of the spherically averaged structure factor or some moment of S(q,t). Here only the first moment ql will be considered as a measure of average domain size, since qmax cannot be precisely determined due to the discretization of the scattering vector q in a finite lattice. Five independent runs are performed for each case, and all the results are reported by averaging the data from five independent runs. 

1 The intermolecutar atom-atom model for organic crystals 

In the early 1970s, the availability of an ever increasing amount of crystal structure data for organic compounds prompted the development of empirical schemes for the 

The physical interpretation of the necessary parameters is difficult or impossible, if one considers how difficult it must be to reconcile the simple expressions 4.38-4.40 with the intricate electronic effects that have been sketched in the preceding sections. A recommendable attitude [46] is then to consider the atom-atom model as a useful numerical machinery for the calculation of lattice energies, without attaching too much physical significance to it. Thus, terms are attractive terms but cannot be considered to represent dispersion or polarization, or exponential terms are repulsive terms but do not necessarily represent Pauli repulsion, and even the terms are of a coulomb-law type but, as discussed in former sections of this chapter, are but a poor representation of the true coulombic potential energy. 

The calibration of the parameters in an atom-atom interaction energy curve is better carried out in terms of the position of the minimum, 7 °, the well depth, s, and the curvature parameter, X [47]. An atom-atom potential can be written as  

For a 12-6 potential curve, (dE /dz) = 72 (a constant) by comparison with equation 4.44, a 12-6 curve corresponds to a 6-exp curve with X = 13.8. 

---

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. 

This chapter concentrates on describing molecular simulation methods which have a counectiou with the statistical mechanical description of condensed matter, and hence relate to theoretical approaches to understanding phenomena such as phase equilibria, rare events, and quantum mechanical effects. 

One of the flexibilities of eomputer simulation is that it is possible to define the themiodynamie eonditions eorresponding to one of many statistieal ensembles, eaeh of whieh may be most suitable for the purpose of the study. A knowledge of the underlying statistieal meehanies is essential in the design of eorreet simulation methods, and in the analysis of simulation results. Flere we deseribe two of the most eommoir statistieal ensembles, but examples of the use of other ensembles will appear later in the ehapter. 

Although this section has concentrated on MD, it should not be forgotten that lattice-based MC codes may be parallelized very efficiently for more infomration on parallel simulation methods see [220. 221. 222 and 223] and references therein. 

With the rapid development of computer power, and the continual hmovation of simulation methods, it is impossible to predict what may be achieved over the next few years, except to say that the outlook is very promising. The areas of rare events, phase equilibria, and quantum simulation continue to be active. 

Kremer K 1996 Computer simulation methods for polymer physics Monte Carlo and Molecular Dynamics of Condensed Matter Systems vol 49, ed K Binder and G Ciccotti (Bologna Italian Physical Society) pp 669-723... 

Swope W C and Andersen H C 1995 A computer simulation method for the calculation of chemical potentials of liquids and solids using the bicanonical ensemble J. Chem. Phys. f02 2851-63... 

Consta S, Wilding N B, Frenkel D and Alexandrowicz Z 1999 Recoil growth an efficient simulation method for multi-... 

Kotelyanskii M 1997 Simulation methods for modelling amorphous polymers Trends Polym. Sol. 5 192-8... 

Handy, N.C. Density functional theory. In Quantum mechanical simulation methods for studying biological systems, D. Bicout and M. Field, eds. Springer, Berlin (1996) 1-35. 

P. Bala, P. Grochowski, B. Lesyng, and J. A. McCammon Quantum-classical molecular dynamics. Models and applications. In Quantum Mechanical Simulation Methods for Studying Biological Systems (M. Fields, ed.). Les Houches, France (1995)... 

Other methods which are applied to conformational analysis and to generating multiple conformations and which can be regarded as random or stochastic techniques, since they explore the conformational space in a non-deterministic fashion, arc genetic algorithms (GA) [137, 1381 simulation methods, such as molecular dynamics (MD) and Monte Carlo (MC) simulations 1139], as well as simulated annealing [140], All of those approaches and their application to generate ensembles of conformations arc discussed in Chapter II, Section 7.2 in the Handbook. 

D. W Heermann, Computer Simulation Methods in Theoretical Physics, Springer, Berlin, 1986. 

Monte Carlo sim u lat ion s pro vide an altern ate approach to the generation of stable con form ation s. As with HyperCh ern s o th er simulation methods, the effects of temperature changes and solvation arc easily incorporated into th c ealcii lation s. 

In this chapter we shall discuss some of the general principles involved in the two most common simulation techniques used in molecular modelling the molecular dynamics and the Monte Carlo methods. We shall also discuss several concepts that are common to both of these methods. A more detailed discussion of the two simulation methods can be found in Chapters 7 and 8. 

A ll ide variety of thermodynamic properties can be calculated from compufer simulations a Comparison of experimental and calculated values for such properties is an important way in which the accuracy of the simulation and the underlying energy model can be quantified. Simulation methods also enable predictions to be made of the thermodynamic properties of V.stems for which there is no experimental data, or for which experimental data is difficult or impossible to obtain. Simulations can also provide structural information about the... 

The correct treatment of boundaries and boundary effects is crucial to simulation methods because it enables macroscopic properties to be calculated from simulations using relatively small numbers of particles. The importance of boundary effects can be illustrated by considering the following simple example. Suppose we have a cube of volume 1 litre which is filled with water at room temperature. The cube contains approximately 3.3 X 10 molecules. Interactions with the walls can extend up to 10 molecular diameters into the fluid. The diameter of the water molecule is approximately 2.8 A and so the number of water molecules that are interacting with the boundary is about 2 x 10. So only about one in 1.5 million water molecules is influenced by interactions with the walls of the container. The number of particles in a Monte Carlo or molecular dynamics simulation is far fewer than 10 -10 and is frequently less than 1000. In a system of 1000 water molecules most, if not all of them, would be within the influence of the walls of the boundary. Clecirly, a simulation of 1000 water molecules in a vessel would not be an appropriate way to derive bulk properties. The alternative is to dispense with the container altogether. Now, approximately three-quarters of the molecules would be at the surface of the sample rather than being in the bulk. Such a situation would be relevcUit to studies of liquid drops, but not to studies of bulk phenomena. 

A number of simulation methods based on Equation (7.115) have been described. Thess differ in the assumptions that are made about the nature of frictional and random forces A common simplifying assumption is that the collision frequency 7 is independent o time and position. The random force R(f) is often assumed to be uncorrelated with th particle velocities, positions and the forces acting on them, and to obey a Gaussiar distribution with zero mean. The force F, is assumed to be constant over the time step o the integration. 

---