Monte Carlo Simulation Methods for Nanosystems

"Random Numbers. How are the various decisions made? To start with, the computer must have a source of uniformly distributed psuedo-random numbers. A much used algorithm for generating such numbers is the so-called von Neumann "middle-square digits." Here, an arbitrary n-digit integer is squared, creating a 2n-digit product. A new integer is formed by extracting the middle n-digits from the product. This process is iterated over and over, forming a chain of integers whose properties have been extensively studied. Clearly this chain of numbers repeats after some point. H. Lehmer has suggested a scheme based on the Kronecker-Weyl theorem that generates all possible numbers of n digits before it repeats. (See "Random-Number Generators" for a discussion of various approaches to the generation of random numbers.)

Once one has an algorithm for generating a uniformly distributed set of random numbers, these numbers must be to transformed into the nonuniform distribution g desired for the property of interest. It can be shown that the function f needed to achieve this transformation is just the inverse of the nonuniform distribution function, that is f=g-1." Nicholas C. Metropolis


The Metropolis Monte Carlo (MC) simulation methods can be used in nanoscience to simulate various complex physical phenomena including property calculations, prediction of phase transitions, self-assembly, thermally averaged structures and charge distributions, just to name a few [1]. There exist a variety of MC simulations, which are used depending on the nano system under consideration and the kind of computational results in mind. They include, but not limited to, Classical MC, Quantum MC and Volumetric MC.

In the Classical MC the classical Boltzmann distribution is used as the starting point to perform various property calculations. Through

Quantum MC simulations one can compute quantum mechanical energies, wave functions and electronic structure using the Schrodinger equation. The Volumetric MC is used to calculate molecular volumes and sample molecular phase-space surfaces.

The basis of MC simulations is random numbers. The most important problem for a successful MC simulation is generating random numbers [1].

Generating Random Numbers

One typically needs to generate random numbers distributed according to some given distribution function. Let us first see how uniformly distributed random numbers can be generated. For the sake of completeness, we will briefly mention how one can generate a set of pseudo random numbers. The prefix "pseudo" is used because the generated numbers are never truly random, since they can be regenerated exactly if the same initial "seed" and computer program are chosen. More importantly, the generated set is always periodic, but one tries to make its period as large as possible to avoid any correlations among the generated numbers [1-3].

In general, however, we recommend the reader to use the compiler's own random, or rather "pseudo random" numbers package, or subroutines written by professionals.

Generating Uniformly Distributed Random Numbers in [0,1)

If there is no stochastic process occurring in the computer chip, it will not be able to generate truly random numbers. All it could produce are "pseudo" random numbers generated according to some specific algorithm. The most widely used algorithm is the "Linear Congruential Random Number Generator" presented by D. J. Lehmer in 1949 [4]. The set of pseudo random numbers rk may be generated according to the following rule:

so = seed; Sk = (g * Sk _i + i) mod(m), for k > 0 0 < Sk < m rk = Sk/m (1)

The numbers rk thus defined are "pseudo-randomly" distributed in [0,1). The starting point s0 is called the seed, so that the same set of numbers can be reproduced if one uses the same seed. The symbols g, i, and m are called generator (or multiplier), increment and modulus, respectively. The "mod" function, which C++ implements with %, returns the remainder. For example, 90 mod 13, written as 90 % 13 in C++, is 12, the remainder in the division of 90 into 13.

As an example of the application of Eq. (1) let us choose the initial value in the sequence of random numbers, s0 = 201 and g = 9, i = 2, m = 13, then Sj=(9*201+2)mod13= 1811 mod13=4 and r;=4/13, (4 is the remainder in the division of 1811 into 13). After we calculate one random number "4/13", to evaluate the next random number, we substitute sj into the expression on the right hand side. Consequently we calculate, s2 = (9 * 4+2) mod 13 = 38 mod 13 = 12 and the next random number is r2=12/13. This can be continued repeatedly and the result will be the following set of periodic random numbers,


The sequence produced in this way is periodic of 3, but one can make the period very large by choosing appropriate values for the above parameters. The seed, s0, should preferably be a large odd number, otherwise, one needs to drop the first few generated points. The choice of m = 2e simplifies the mod operation: in the binary representation of a given number, one only needs to keep the modulus m with most digits. Exponent E should be chosen as large as possible: E>35 would be a good choice to make the period as large as possible. Usually, it is taken as the number of bits in the largest positive integer. A large value for g will reduce serial correlations. The number i must be prime to m; g-1 should be multiple of every prime factor of m. In this case the period of the sequence would be 2E. For more details see [1,3].

Generating Random Numbers in Ta,b) According to a Given Distribution Function P(x)

Importance Sampling: Also known as "biased sampling" can be used for generating random numbers according to a given distribution function P(x). This is because the efficiency of random sampling methods can be appreciably enhanced by a suitable choice of sample points. The procedure described is called "importance sampling".

The purpose behind importance sampling is to concentrate on the sample points where the value of the function is appreciable and avoid taking sample points where the value of the function is quite negligible. It is necessary to allow for this bias in the sample points by weighting the sample values appropriately. The general approach to importance sampling consists of two steps:

(1) The first step is to carry out a preliminary evaluation of the function. The region of interest is divided into fairly large cells, the value of the function is determined at the center of each cell, and this value and the cumulative sum are recorded.

(2) The next step is to take as many sample points as is desired, but to adapt the sampling procedure so that the number of sample points in each cell is proportional to the function value at the center of that cell, as determined in step (1). However, each sample is given a weight inversely proportional to the central value of its cell.

There are a few ways to generate random numbers through importance sampling. We will mention here four such methods.

(i) For the sake of simplicity, we can mention one simple algorithm, which can easily be implemented and used in one dimension. For a given interval [a,b) and desired number of random points N, one will first divide the interval into VN pieces (this number is somehow arbitrary), and then generate a uniform set of points in each piece, the number of generated points being np( x,. y x P( x), where xi is a point in the middle of each piece. Clearly, the generated points are distributed according to P and they are almost N in number. It must be noted, however, that the generated points are strongly correlated in position and need a reshuffling before being used.

(ii) It is possible to analytically map the uniform random numbers in [0,1) to [a,b) if the distribution function P(x) is simple enough. If r is uniform in [0,1), then we have x dr = P( x)dx ^ r = J P(t )dt = g (x, a); g (b, a) = 1 (2)

This relation needs to be inverted if possible to obtain x as a function of r. If r is uniform in [0,1), then x will be distributed in [a,b) according to P.

(iii) For an arbitrarily complicated function P, a third way is to use the well-known Metropolis algorithm [5]. In this algorithm, one can start from an arbitrary point in [a,b), say x0=(a+b)/2, and add a random number to it: x1=x0+d(2r-1) where r is the random number in [0,1) and d is the magnitude of the step. If P(x1)/P(x0)>1, then the move is accepted; otherwise the ratio is compared to a random number in [0,1). If the ratio is greater, then again, the move is accepted, otherwise it is rejected and the old point x0 is kept, and one starts another trial from this point. There is also the possibility that the generated points go out of [a,b) range in which case, they will be rejected. Compared to the first method, this is by far less efficient because there will be many rejections. The efficiency of this method depends strongly on the choice of step size d. With a small d most trials will be accepted, but a good sampling will require many trials. Whereas if d is large, one can sample the [a,b) interval quicker, but there might be many rejections.

(iv) A more efficient method due to Von Neumann proposed in 1951 consists of considering a function f {for x in [a,b); f(x) > P(x)} for which one can use method (ii) presented above, and generate points distributed according to f (see Figure 1). The simplest choice would be to take f=constant; but larger than the maximum of P. Once such point x is generated in [a,b), one compares a random number r in [0,1) to the ratio P(x)/f(x). If r is less than the ratio, the number x is accepted, otherwise it is rejected and one should repeat the process with another point x distributed according to f. For constant f, the probability of accepting x is clearly proportional to P(x).

Figure 1. Von Neumann algorithm involves selecting 2 random numbers in [0,1], the first one, distributed according to f, is the x-coordinate; and the second one, uniform, the y-coordinate between 0 and f(x). If it falls below P(x) the point is accepted [1].

Monte Carlo Integration Method

Standard deterministic methods of numerical integration of a function are well known (see for example [6]). The simplest method is the trapezoidal rule, which adds up the function values at equally separated points and multiplies the result by the width of the intervals. A more accurate method is the Simpson's rule, which, instead of approximating the function by linear segments on each interval, replaces it by a quadratic function. Even a more sophisticated method is the Gaussian quadrature, where the integral is the weighted sum of the function evaluated at certain special points. These methods, even though quite accurate for one-dimensional (single) integrals, soon become very prohibitive in multi-dimensional (multiple) integrals as the required number of function evaluations grow quite rapidly as Nd. N is the number of function evaluations required for a one-dimensional integral and d is the number of dimensions in a multi-dimensional integral. In this situation, Monte Carlo integration using random numbers can become very handy.

Regardless of the dimension, one can approximate the integral by a sum in which the function is evaluated at random points:

The points x{ are randomly and uniformly distributed on the [a,b) interval. As the number of points N grows, the central limit theorem (CLT) guarantees the convergence of the sum to the real value of the integral. The error being of the order of . More precisely, the CLT

states that if we repeat this process of generating N random numbers and computing the above sum S, we generate a set of random variables S. These variables have a Gaussian distribution of mean I (the real value of the integral) and variance o = v

/4n where v is the mean deviation of

(b - a) f( x) from the mean value I . The value of v can be estimated, however, by considering the deviation from the mean in a given sample

(the usual standard deviation). Of course v defined as above is sample-dependent; one therefore can take its average over all the taken samples. The CLT implies that as the number of samples is increased, the mean obtained from the samples converges to the real mean I ~ 1/ y[m , where m is the number of samples, and the uncertainty on the accuracy of the result is reduced. The error may be estimated from the CLT:

0 0

Post a comment