MacKinnon - Computational Physics, страница 10

For example,, where is the numberthe error in a trapezium rule calculation of a –dimensional integral falls as of different values of the integrand used. In a Monte–Carlo calculation the error falls as independently of the dimension. Hence forMonte-Carlo integration will usually converge faster.We consider the expression for the average of a statistic,when is a random number distributedaccording to a distribution \, then e"<Ky&¾¾Ž ×jk"]KÑ&jk"]KÑ&gf;> \ "<Ky&„jk"]KÑ&%1K :eKKK10eTR’ 0(4.6)hfkjwhich is just a generalisation of the well known results for (e.g.

)or i, where we are using thenotationto denote averaging. Now consider an integral of the sort which might arise while using`Laplace transforms.—Ž> J!²"„$qKÑ&Rjk"<KÑ&„K["]KÑ&;hJ!²"„$qKÑ&kj "]KÑ& … jk"<K &\[TThis integral can be evaluated by generating a set of random numbers, ,distribution, \, and calculating the mean ofas]K²T : [[ [ : K (4.7)- , from a Poisson(4.8)eThe error in this mean is evaluated as usual by consideringthe corresponding standard error of themeanGjlf(4.9)iAs mentioned earlier, Monte–Carlo integration can be particularly efficient in the case of multi–dimensionalintegrals. However this case is particularly susceptible to the flaws in random number generators. It is acommon feature that when a random set of coordinates in a –dimensional space is generated (i.e.

a setof random numbers), the resulting distribution contains (hyper)planes on which the probability is eithersignificantly higher or lower than expected.¾0$ ° j 0 "<Ky& m$ jk"<Ky& 0 ± [¾Monte Carlo Methods and Simulation404.3 The Metropolis Algorithm e ¿ Þiäe(4.10)ä ¿ Þº is the energy of the system in state i and ã 7 Ëon Q .

Such problems can be solved using thewhereMetropolis et al. (1953) algorithm.ÒLet us suppose the system is initially in a particular state and we change it to another state â . Thedetailed balancecondition demands that in equilibrium q p'€ € the/€ p9flow from Ò to â must be balanced by the flowÒfrom â to .

This can be expressed as \ Q  p€\ Q(4.11)\where is the probability of finding the system in state i and Qis the probability (or rate) that a system p€€in state i will make a transition to state j. (4.11) can be rearranged to read/€ p9 \\  5 m mQ(4.12)Qeeƒ ¿ ‚Ï Þ[(4.13)Generally the right–hand–side of (4.3) is known and we want to generate a set of states which obey the€ rates such€ that distribution \ . This can be achieved by choosing the transition€ p€c ` \  e 5 m e m ƒ if \ € Ž— \  or ºº € —lºº  [‚ Ï Þ if \ \ or Ž(4.14)€  Q\ or ¿d 5m m— \ € a randome number,e Þ ƒ r , is chosen between and and the system is moved to state âIn practice if \7‚Ï\\only if r is less thanor ¿.5mIn statistical mechanics we commonly wante to evaluate thermodynamicaverages of the form5mTf(4.14) is not the only way in which the condition (4.3) can be fulfilled, but it is by far the most commonlyused.An important feature of the procedure is that it is never necessary to evaluate the partition function, thedenominator in (4.10) but only the relative probabilities of the different states.

This is usually much easierto achieve as it only requires the calculation of the change of energy from one state to another.Note that, although we have derived the algorithm in the context of thermodynamics, its use is by nomeans confined to that case. See for example section 4.4.4.3.1 The Ising model €€¯$ … ß ßAs a simple example of the Metropolis Method we consider the Ising model of a ferromagnets9ß} ( 0T(4.15)€€,º L$ … ßwhere J is a positive energy,, and i and j are nearest neighbours on a lattice.

In this case we changefrom one state to another by flipping a single spin and the change in energy is simplys(4.16)where the sum is only over the nearest neighbours of the flipped spin.The simulation proceeds by choosing a spin (usually at random) and testing whether the energy wouldbe increased or decreased by flipping the spin.5utvm If it is decreased the rules (4.14) say that the spin shoulddefinitely be flipped. If, on the other hand, the energy is increased, a uniform random number, r , betweenand is generated and compared with.

If it is smaller the spin is flipped, otherwise the spin isunchanged.Further information can be found in section 4.7.¿eMonte Carlo Methods and Simulation414.3.2 Thermodynamic AveragesTo average over a thermodynamic quantity it suffices to average over the values for the sequence of statesgenerated by the Metropolis algorithm. However it is usually wise to carry out a number of Monte–Carlo steps before starting to do any averaging.

This is to guarantee that the system is in thermodynamicequilibrium while the averaging is carried out.The sequence of random changes is often considered as a sort of time axis. In practice we think (e.g.) about the time required to reach equilibrium. Sometimes, however, the transition rate becomes very lowand the system effectively gets stuck in a non–equilibrium state. This is often the case at low temperatureswhen almost every change causes an increase in energy.4.4 Quantum Monte–CarloThe term Quantum Monte–Carlo does not refer to a particular method but rather to any method usingMonte–Carlo type methods to solve quantum (usually many–body) problems. As an example consider theevaluation of the energy of a trial wave function w r where r is a 3N dimensional position coordinate ofN particles,yxrw ! r Hw r(4.17)rw ! r w rxwhere H is the Hamiltonian operator.

This can be turned into an appropriate form for Monte–Carlo integration (see section 4.2) by rewriting as"&º û " " & & " " & &ûºÃ> r Ó z w " r& " z 0 & 0 Ô · w " r& eT Hx w " r&C¸û rzw r z(4.18)such that the quantity in braces ( ,o- ) has the form of a probability distribution. This integral can now easilybe evaluated by Monte–Carlo integration.

Typically a sequence of r’s is generated using the Metropolis(section 4.3) algorithm, so that it is not even necessary to normalise the trial wave function, and the quantityin square brackets (4.18) is averaged over the r’s.This method, variational quantum Monte–Carlo, presupposes we have a good guess for the wave function and want to evaluate an integral over it.

It is only one of several different techniques which are referredto as quantum Monte Carlo. Others include, Diffusion Monte–Carlo, Green’s function Monte–Carlo andWorld Line Monte–Carlo.uv4.5 Molecular DynamicsAn alternative approach to studying the behaviour of large systems is simply to solve the equations ofmotion for a reasonably large number of molecules.

Usually this is done by using one of the methodsdescribed for Ordinary Differential Equations (chapter 1) to solve the coupled equations of motion forthe atoms or molecules to be described. Molecular solids in particular can be studied by consideringthe molecules as rigid objects and using phenomenological classical equations of motion to describe theinteraction between them. As a simple example the Lennard–Jones potential" r &; | Œ ° œ ± T 0 $ f ° œ a± { rrz(4.19)has been successfully used to describe the thermodynamics of noble gases.4.5.1 General Principles Consider a set of particles interacting through a 2–body potential,written in the formrvz" r & . The equations of motion can be(4.20)Monte Carlo Methods and Simulation42 v ¯$   …€}|  z "  $ € &y[ ~ / r r ] z r r z(4.21)r"U& z 7 r to avoid unnecessary numerical differentiation.In practice it is better to use the scalar force rr rA common feature of such problems is that the time derivative of one variable only involves the othervariable as with r and v in the above equations of motion (4.5.1).

In such circumstances a leap–frog likeQ appropriate.QeQeT we writemethod suggests itself as the mostHence r  0 f O v ~  € Q €Qr(4.22)QS#T QegT f O …}€ |  Q $ Q È rQ $ r€Q [ / ] Ç _r r _ r $ rv(4.23)v_ _____This method has the advantage of simplicity as well as the merit of being properly conservative (section 2.6).The temperature is defined from the kinetic energy of N particles viau~T‡0 Ë n Q … 0 / v0—Ž(4.24)where the averagesare taken with respect to time.€As an example of another thermodynamic quantity consider the specific heat at constant volume.This can be calculated by changing the total energy by multiplying all the velocities by a constant amount() and running for some time to determine the temperature.

Note that, as the temperature is definedin terms of a time average, multiplying all the velocities by does not necessarily imply a simple changeQ . At a 1st order phase transition, for example, the system might equilibrate to theof temperature, Qsame temperature as before with the additional energy contributing to the latent heat.When the specific heat is known for a range of temperatures it becomes possible, at least in principle,to calculate the entropy from the relationship1(4.25)Q€€QQ‚1The pressure is rather more tricky as it is defined in terms of the free energy,Q , using F› 0ßr  r r~r\ $ r ½NSƒr~ $ ß(4.26)and hence has a contribution from the entropy as well as the internal energy. Nevertheless methods existfor calculating this.It is also possible to define modified equations of motion which, rather than conserving energy andvolume, conserve temperature or pressure.