Monte Carlo Simulation Methods Lecture Notes

1

MONTE CARLO SIMULATION METHODS Monte Carlo simulation employs random sampling. That is configurations are generated by making random changes in positions of the atoms / molecules. The orientations and conformations are also changed randomly when required. Monte Carlo simulations use a technique called ‘importance sampling’. Importance sampling generates low energy configurations and hence the properties can be calculated accurately. Calculating properties Properties are calculated by integration. For example, the average potential energy can be determined by evaluating the following integral:

= ∫ drN V(rN)ρ(rN)

-------------------- (1)

ρ(rN) is the probability of getting the configuration rN. ρ(rN) = exp(-V(rN) / kBT) / Z ----------------------- (2) where ‘Z’ is the configurational integral ZNVT = ∫drN exp(-V(rN) / kBT)

----------------------- (3)

These integrals cannot be evaluated analytically. Hence, numerical methods are used to evaluate these integrals, for e.g, by using trapezium rules. The area under the curve is approximated as the sum of the areas under the trapeziums. The trapezium rule is a way of estimating the area under a curve. We know that the area under a curve is given by integration, so the trapezium rule gives a method of estimating integrals. This is useful when we come across integrals that we don't know how to evaluate. The trapezium rule works by splitting the area under a curve into a number of trapeziums, which we know the area of. If we want to find the area under a curve between the points x0 and xn, we divide this interval up into smaller intervals, each of which has length h (see diagram above). Then we find that:

2

where y0 = f(x0) and y1 = f(x1) etc If the original interval was split up into n smaller intervals, then h is given by: h = (xn - x0)/n

3

Random method This is an alternative method.

We start the familiar example of finding the area of a circle. Figure 1 below shows a circle with radius inscribed within a square. The area of the circle is , and the area of the square is . The ratio of the area of the circle to the area of the square is

Figure 1. If we could compute ratio, then we could multiple it by four to obtain the value . One particularly simple way to do this is to pick lattice points in the square and count how many of them lie inside the circle, see Figure 2. Suppose for example that the points are selected, then there are 812 points inside the circle and 212

4

points outside the circle and the percentage of points inside the circle is following calculation

. Then the area of the circle is approximated with the

Monte Carlo Method for Monte Carlo methods can be thought of as statistical simulation methods that utilize a sequences of random numbers to perform the simulation. The name "Monte Carlo'' was coined by Nicholas Constantine Metropolis (1915-1999) and inspired by Stanslaw Ulam (1909-1986), because of the similarity of statistical simulation to games of chance, and because Monte Carlo is a center for gambling and games of chance. In a typical process one compute the number of points in a set A that lies inside box R. The ratio of the number of points that fall inside A to the total number of points tried is equal to the ratio of the two areas (or volume in 3 dimensions). The accuracy of the ratio depends on the number of points used, with more points leading to a more accurate value. A simple Monte Carlo simulation to approximate the value of selecting points number of points that satisfy were points satisfying and

could involve randomly

in the unit square and determining the ratio , where is . In a typical simulation of sample size there , shown in Figure 3. Using this data, we obtain

5

Figure 3.

The idea is to use random points for the numerical evaluation of an integral, that is, using random points to determine the area under the function, see the picture below

6

The integral of the function f(x) is approximately the total area times the fraction of points that fall under the curve of f(x). Naturally this method for evaluation of an integral (using "random" points) is competitive only for the multi-dimensional case and/or complicated functions. Area under the curve = Total area(A) x [No. of points under the curve / total number of points Generated] Metropolis Monte Carlo Method Consider the simulation of an atomic fluid. At each step of the simulation, a new configuration is generated by making a random change in the Cartesian coordinates of a single randomly chosen particle using a random number generator. A unique random number is generated individually for the three directions. Then the energy of the new configuration is calculated. The entire system need not be considered for the energy calculation but only those contributions involving the particle which has been just moved are considered. [The Metropolis algorithm is as follows. The first conformation is randomly generated. At each point in the construction of the chain of conformations a move is attempted to the current conformation. The move is rejected immediately if the local chain conformation is not compatible with the attempted move or if it violates the excluded volume condition. If these two conditions are satisfied then the so called Metropolis criterion is applied. If the difference between the energy of the resulting conformation and the energy of the current conformation, , is negative (i.e. the energy of the resulting conformation is smaller than the energy of the current conformation), then the resulting conformation is accepted and it becomes the new conformation in the chain. If is positive, however, a (pseudo)random number between 0 and

7

1, 0
. If

then the resulting conformation is refused. Whenever the conformation resulting from the attempted move is refused for any of the three possible reasons, then the new conformation of the chain is the same current conformation. For sequence selection the same algorithm is used but the ``moves'' correspond to switching the position of two monomers while the conformation is kept fixed ] If the new configuration has lower energy than the previous one it is accepted and it is the starting point for the next iteration. If the new configuration has a higher energy than the previous one, then the Boltzman factor exp(-ΔV / kBT), is compared to a random number between 0 and 1. If the Boltzman factor is greater than the random number, the new configuration is accepted. If it is not, the new move is rejected and the original configuration is retained for the next move. Monte Carlo simulation of rigid molecule For rigid and non-spherical molecules, their orientations should also be changed in addition to their position. During each Monte Carlo step, one molecule is subjected to translation and rotation. Positions are expressed in terms of the centre of the masses of the particles and the orientation is changed by rotating around a chosen axis by a random angle. For example, if rotated about x-axis by an angle δw, then the new positions will be:

x’ y’

1 =

z’

0

0

x

0

cosδw -sinδw

y

0

sinδw

z

cosδw

Euler’s angles A better way of representing the rotations of molecules is to use Euler’s angles, φ, θ and ψ. φ is the rotation about the z-axis; x- and y-axes are moved. θ is the rotation about the new x-axis; y- and z-axes are moved.

8

ψ is the rotation about new z-axis. If the Euler angles are randomly changed by δφ, δθ and δψ, then a vector Vold is moved according to the following matrix equation: Vnew = AVold Where matrix A is complex matrix in sin and cos of the Euler’s angles Disadvantages of the Euler’s angles The rotation matrix contains a total of six trigonometric functions for each of the three Euler’s angles. The trigonometric functions are computationally expensive to calculate. The alternative is to use quaternions. Quaternions It is a four-dimensional vector such that the sum of the components = 1. q02 + q12 + q22 + q32 = 1 The quaternion components are related to the Euler angles as follows: q0 = cos(1/2)θ cos(1/2)(φ + ψ) q1 = sin(1/2)θ cos(1/2) (φ + ψ) q2 = sin(1/2)θ sin(1/2) (φ + ψ) q3 = cos(1/2)θ sin(1/2) (φ + ψ) Then the rotational matrix A can be written in terms of the quaternion. The quaternion vector is rotated to new (random) orientation to generate a new orientation. Since it is a four-dimensional vector, the orientation must be performed in fourdimensional space.

Monte Carlo Simulation of Flexible Molecules 1. System must be small

9

2. Some of the internal degrees of freedom are frozen or 3. Specials models or methods are employed. •

Random changes are performed to the Cartesian coordinates of individual atoms, translations and rotations of the entire molecule to generate a new configuration of the flexible molecule.

•

Very small atomic displacements are required to achieve an acceptable acceptance ratio. Hence, the phase space is covered slowly.

•

Even small movements away from an equilibrium bond length will cause a large increase in the energy. Hence, bond lengths, bond angles etc. are frozen. This works well for small molecules such as butane.

•

But for large molecules even relatively small bond rotations may cause large movements of atoms along the chain. This will lead to high energy configurations as shown below:

ω

•

A bond rotation in the middle of a molecule may lead to a large movement at the end.

Models used in Monte Carlo simulation of polymers •

Polymer is a macro molecule constructed by chemically linking together a sequence of molecular fragments.

•

In polyethylene or polystyrene, all of the molecular fragments comprise of the same basic unit or monomer.

•

Other polymers may contain mixtures of monomers.

•

Proteins are polypeptide chains in which each unit is one of the twenty amino acids.

10

•

Cross-linking between different chains leads to further variations in the constitution and structure of a polymer.

•

All these features may affect the overall properties of the molecule.

•

The properties of the polymer under different conditions such as in solution, in a polymer melt or in the crystalline state are of interest.

•

We can develop theories based on Molecular Modeling to understand these properties of polymers and we can also predict their properties from Molecular Modeling.

Lattice Models of polymers (a) Growing a polymer onto the lattice 1. This model is simple because many states can be generated and studied very

rapidly. 2. 2D and 3D lattices are used. 3. The simplest models are cubic and tetrahedral lattices. 4. Successive monomers occupy adjacent lattice points

Cubic lattice

Tetrahedral lattice

Complex Models