Problems & Solutions :: Molecular Simulations in Mechanics and Physics PDF

Preview

Summary

Problems and Solutions Molecular Simulations Methods in Mechanics and Physics Dibakar Datta Brown University, Providence, USA Sang-Pil Kim Samsung Electronics, Seoul, South Korea Vivek B Shenoy The University of Pennsylvania, Philadelphia, USA Contents 1 Random Numbers and Molecular Simulations …………...

Description

Problems and Solutions

Molecular Simulations Methods in Mechanics and Physics

Dibakar Datta Brown University, Providence, USA

Sang-Pil Kim Samsung Electronics, Seoul, South Korea

Vivek B Shenoy The University of Pennsylvania, Philadelphia, USA

Contents 1 Random Numbers and Molecular Simulations ………………………….... 3 2 Use of Monte Carlo Methods to Do Some Integrals ………………………. 8 3 Use of Monte Carlo Methods for Evaluation of Thermodynamics Properties of One Multistate System ……………………………………....13 4 Use of Monte Carlo Methods to Evaluate the Thermodynamics Properties of the One Dimensional Ising Model ………………………………………18 5 Use of Monte Carlo Methods to Evaluate the Thermodynamics Properties of the Two Dimensional Ising Model …………………………………… 23 6 Kinetic Monte Carlo Method for Langmuir Adsorption-Desorption Problem ……………………………………………………………………. 27 7 Molecular Dynamics for Simple LJ (Two & Three Dimensional Particles) Systems ........................................................................................................... 30 8 Computation of Pressure and Radial Distribution Function Using LAMMPS …………………………………………………………………... 38 9 Computation of Diffusivity and Thermal Conductivity Using LAMMPS 41 10 Indentation Using LAMMPS ………………………………………………43 Appendix A Matlab Code for Chapter 1 ………………………………... 44 Appendix B Matlab Code for Chapter 2 ………………………………... 47 Appendix C Matlab Code and Data File for Chapter 3 ………………... 51 !

1!

Appendix D Matlab Code for Chapter 4 ………………........................... 60 Appendix E Matlab Code for Chapter 5 ………………............................ 65 Appendix F Matlab Code for Chapter 6 ………………............................ 73 Appendix G Lammps Script for Chapter 7 ……………………………... 80 Appendix H Lammps Script and Matlab Code for Chapter 8 ………… 86 Appendix I Lammps Script and Matlab Code for Chapter 9 ………….. 94 Appendix J Lammps Script and Matlab Code for Chapter 10 ………… 98 ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !

2!

!

1 Random Numbers and Molecular Simulations Problem 1.1 Generate 4 million random numbers uniformly distributed between 0 and 1. Compute their mean by distributing them into 100. Estimate the error on the mean and see if it agrees with the deviation of the computed mean from 0.5. Now distribute them into 1000 bins and repeat the procedure.

Procedure of ‘binning’ Take the N events and break them into m sets or ‘bins’ of n = N/m events each (make sure N/m is an integer). Find the mean of the numbers in each of these m bins; let these means be xi , i =1,2,…,m. If the central limit theorem holds, the distribution of xi’s should be Gaussian. Moreover, the uncertainty in the mean of these means, i.e., in

1 m x = ∑ xi m i=1

(1.1)

is expected to be

(x

δx =

2

− x

2

) / ( m − 1)

(1.2)

where

x

2

1 m 2 = ∑ xi m i=1

(1.3)

If m is a large number, we may with very little error replace m-1 by m in this equation. Let’s see how that works for the numbers produced by the generator. We have a lot of latitude in the choice of m. So long as N/m is large enough for the central limit theorem to come into play, all such choices should give much the same value for δx. Also, δx should be comparable, provided the generator is good, to the difference !

3!

between < x > and ½ . What does that mean? It means that roughly two-thirds of the time, δx will be larger than the absolute difference between < x > and ½; the rest of the time the difference will be larger than δx, but not dramatically larger. One to two times as large will be common; two to three times as large will certainly occur but not more than a few percent of the time; and more than three times as large will be rare. No of Bin 100 1000

Deviation With Bin 0.000215 0.000078

Error with Bin 0.000144 0.000145

Deviation Without Bin 0.000215 0.000078

Error Without Bin 0.000144 0.000144

Analytical Error 0.000215 0.000144

Table 1.1 Estimation of error on the mean with and without binning ! !

Problem 1.2 Find a random number generator which claims to produce evenly distributed random numbers on the interval [0,1] or (0,1). Run a test of the generator in which you generate N random numbers with N around ten million and see how they are distributed by dividing the interval into ten subintervals of equal size. Determine the fraction of the numbers in each subinterval and also the uncertainty in this number. Determine the latter by dividing the run up into ten parts (often called “bins"), finding, for each subinterval, the number of events in each bin and then computing, for each subinterval, the mean and the variance of these numbers. Look at the correlations in the sequence of numbers; that is evaluate, for n between 0 and 20 1 N χ ( n ) = ∑ ( xi+n − x N i=1

1 N )( xi − x ) = N ∑ xi+n xi − ( x i=1

)

2

(1.4)

using periodic boundary conditions, meaning x N + j ≡ x j . Notice that

1 N χ ( 0 ) = ∑ ( xi − x N i=1

)( x − x ) = i

x2 − ( x

)

2

is the variance of x. Normalize the correlation function so that χ ( 0 ) = 1 . !

4!

(1.5)

On the basis of your results, does the generator appear to be working, as it should? We have considered 1 million numbers for the speed of computing. !

Figure 1.1 The distribution of random numbers in ten subintervals of equal interval

!

Bin No 1 2 3 4 5 6 7 8 9 10 0.100122 0.099599 0.099753 0.100115 0.099969 0.099820 0.100057 0.100257 0.099806 0.100502 Fraction Table 1.2 Fraction of the numbers in each subinterval !

Bin No 1 2 3 4 5 6 7 8 9 10 100122 99599 99753 100115 99969 99820 100057 100257 99806 100502 No. of Events 0.049953 0.149975 0.249852 0.349914 0.450019 0.550229 0.649917 0.749936 0.849850 0.949997 Mean Variance 0.000833 0.000832 0.000832 0.000831 0.000837 0.000832 0.000831 0.0008309 0.000836 0.000836 ! Table 1.3 Uncertainty quantification of the random numbers

! !

5!

n 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

correlation 1 0.000020259813097612 0.000011675673976863 -0.000060191483037431 0.000021630773839709 -0.000141177154716232 -0.000044922519324975 0.000102002320232830 -0.000013534856363240 0.000073100502788581 0.000071693236619486 0.000029834616016466 -0.000059281554875223 0.000249376029785864 0.000008373749979362 -0.000041418950919470 0.000017963205148352 -0.000063381333307910 -0.000011930398607751 0.000184709093162316 0.000080121154401480 Table 1.4 Estimation of correlation of random numbers

!

We can observe that except

χ (0) ,

all other correlations are close to zero. Hence the

random numbers generated are not correlated. So the generator appears to be working, as it should be.

Problem 1.3 Consider n random numbers x1, x2…. xn that follow the probability distribution

p( x) =

1 1 , −∞≤ x≤∞ π 1+ x 2

(1.6)

If the mean of the distribution computed from these numbers is

1 n µ n = ∑ xi n i=1

(1.7)

derive an expression for the probability distribution of µn . ! !

6!

The given distribution is the Cauchy Distribution with Location Parameter ( x0 = 0 ) and Scale Parameter ( γ = 1 ). The n random numbers x1 , x2 ... xn follow this distribution. If x is a Cauchy distributed random variable, the characteristic function of the Cauchy distribution can be written as:

φx (t ) = E ( e

∞

) = ∫ f ( x )e

ixt

ixt

dx = e− t

(1.8)

−∞

This is the Fourier Transform of the probability density. We can express the original probability density in terms of the characteristic function1 by using the Inverse Fourier Transform:

f ( x) =

1 2π

∫

∞

−∞

φ X ( t ) e−ixt dt

(1.9)

Equation 1.6 gives the mean of the distribution computed from these numbers. The Characteristic Function of the sample mean can be written as: φ µn ( t ) = E ( e

i µnt

)

n

⎛ i1 n ⎞ ⎛ −i ⎞ = E ⎜ e n ∑ xi t ⎟ = ⎜ e n ⎟ = e − t ⎝ i=1 ⎠ ⎝ ⎠

(1.10)

The probability distribution of µn is also Cauchy Distribution. !

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1! The characteristic function is not differential at the origin because the Cauchy distribution does not have an expected value.! !

7!

2 Use of Monte Carlo Methods to Do Some Integrals ! !

Problem 2.1 Use Monte Carlo methods to find the natural logarithm of 2. Also, determine the uncertainty in your result and make a comparison with a precise value of the logarithm. In order to evaluate the value of natural logarithm of 2, we perform the integral: !

I=∫

2

1

!

1 dx x

(2.1)

The precise value of the logarithm: 0.693147181 No of bin

No of points per bin

Result

Uncertainty

Percentage Error

100 100

10000 100000

0.6932258 0.6931653

1.345481e-04 4.988322e-05

0.011342 0.0026140

Table 2.1 Evaluation of natural logarithm of 2

Problem 2.2 Use Monte Carlo methods to do the integrals I=∫

2

1

10 ⎞⎤ ⎡ 10 ⎤ ⎡ ⎛ Πdx sin π nx ∑ n⎟ ⎥ ⎢⎣ i=1 i ⎥⎦ ⎢ ⎜⎝ ⎠⎦ ⎣ n=1

No of bin 100

No of points per bin 1000

(2.2)

Result 0.0003951843

Uncertainty 2.277833e-03

! Table 2.2 Evaluation of 10 dimensional integral using MC Method ! ! !

8!

Problem 2.3 Use importance sampling (not the random walk procedure) to do the integral

I=∫

∞

1

e− x 1+ 2sin 2 x

dx

(2.3)

Use importance sampling according to the random walk procedure to do the same integrals. Find the uncertainty in your results.

Importance and Random Sampling Consider the following example in one dimension: I = f = ∫ dx f ( x ) Q ( x ) = ∫ dx F ( x )

(2.4)

where Q(x) is a probability distribution. We may apply the method just learned to the evaluation of this integral by generating a set of random numbers uniformly on the integral over which the integral is taken and adding up the values of F(x) at these points. But suppose that Q(x) is zero or close to zero on a large part of the interval and that it is much larger on some small part of the interval, while f(x) varies relatively little. Then it would be more efficient to sample the integrand in the range of x where Q is large while largelt ignoring the range of x where Q is very small. A way to do that is to generate a sequence of x’s distributed not uniformly on some interval but rather according to the distribution Q(x). That is, one wants the number of elements in the sequence which fall into an interval dx around x to be proportional to Q(x)dx. Given such a sequence with elements xi , i =1,2,…,N , an approximation to f

is obtained from (assuming Q is normalized to unity) f =

1 ∑ f ( xi ) N i

(2.5)

The approximation improves with increasing N. In practice, there are cases for which this approach doesn’t work well. If f is too wildly varying or becomes very large in places where Q is small, it will not give good answers within a reasonable amount, or perhaps within any amount, of computer time. Of course the same can be said of any !

9!

kind of Monte Carlo evaluation of integrals. Let’s look at a sample example, ∞

1 I = ∫ dx x 2 e− x 20

(2.6)

The probability distribution is

Q ( x ) = exp ( −x ) ,

so we need to generate a sequence of

numbers distributed in this way. There are generators of this kind available in standard software packages, but we can easily make one from a generator, which produces a uniform distribution on the interval [0,1]. Consider two distinct distributions P(u) and Q(x) ans ask what is the relation between x and u if the probability in dx around some x is to be the same as that in du around u. Let that relation be u = u(x). Thus we write

P ( u ) du = Q ( x ) dx

(2.7)

and if we can integrate both sides of this expression we will be able to use the result to find u(x). The particular normalized distributions we have are

P ( u ) = 1 on [ 0,1]

and

Q ( x ) = e− x for all x > 0

(2.8)

Thus we consider

du = e− x dx

(2.9)

which may be integrated to give

u = −e− x + C

(2.10)

where C is a constant. If u = 0 is to coincide with x = 0, then C = 1 and

1− u = e− x

or

x = − log (1− u )

(2.11)

Thus we simply employ the usual generator to give a set of u’s uniformly distributed on [0,1] and turn them into a set of x’s distributed according to exp(-x) by means of the preceding formula. We take that set and evaluate I MC =

1 N 2 ∑ xi 2 N i=1

(2.12)

the result should be close to I. That’s been done for a set of one million x’s (or u’s) broken into 100 bins for purposes of finding the uncertainty in the mean. The result is IMC = 1.0011±0.0024 which is quite reasonable in all respects. !

10!

E = ∑ Ei e− β Ei i

∑e

−βEj

(2.13)

j

the sums being over all states of the system. We can’t possibly sum over all the states of a many-particle system, but we can make an approximate evaluation of

E

using

the MC method. At a finite temperature, there will likely be many contributions to the sums. Consequently, it is much more efficient to proceed by generating a sequence of N states with energies distributed according to the Boltzmann distribution. For that sequence of states,

E =

1 ∑ Ei ≈ E N i

At least one hopes that

(2.14) E

will approximate

E.

The first question is how to generate such a sequence of states. In the example given above, it was easy to generate the sequence of x’s because we could integrate the probability distribution to find the necessary relation between x and u and a simple random number generator gave us a good set of u’s. Usually it isn’t that simple, and one must generate the x’s (or states) some other way. One way is to use a biased random walk. Let’s look at a specific example, which is to evaluate a onedimensional integral, I=

∞

∫ dx f ( x ) P ( x )

(2.15)

−∞

I is the mean value of f over the probability distribution P; for definiteness, let

P( x) =

1 − x 2 /2 e 2π

(2.16)

a Gaussian distribution. Now consider the sequence of x’s produced in a random walk which starts at some plausible value of x, such x = 0, and in which the size of the attempted steps is uniformly distributed between preset limits, say –a < δx < a with a on the order of the distribution’s width i.e., around unity. The basic procedure is as follows: Given that at the ith step the position is xi, one generates a random number δx on the interval [-a,a] and considers a move to the position x′ = xi + δ x . If, P ( x ′ ) > P ( xi ) , the move is accepted and if P ( x′ ) < P ( xi ) , one accepts the move with a !

11!

probability P ( x′ ) / P ( xi ) . The latter is achieved by generating a random number on the interval [0,1] and accepting the move if that random number is smaller than P ( x ′ ) / P ( xi ) . If the move is accepted, xi+1 is set equal to

x ′ ; and if the move is not

accepted, xi+1 is set equal to xi. In this way, one generates a set of x’s which will be distributed according to the probability distribution P(x); that is, the number of points in the sequence in the interval dx around x will be proportional to P(x)dx and so one gets an approximation to I by using these points to evaluate I MC =

1 N ∑ f ( xi ) N i=1

(2.17)

For example, let f ( x ) = x 2 and let N be one million. To get a measure of the uncertainty of the mean, subdivide the simulation into one hundred bins. Also, start from x1 = 0 and let a =1. One such run (a typical one) gives I = 0.9991±0.0048; the exact answer is unity. Now let’s solve the Problem 2.3. !

No of Bin 100

No of points per bin 10000

Result

Uncertainty

0.7884821

1.514067e-04

! Table 2.3 Evaluation of integral using Importance Sampling Method

No of Bin 100

No of points per bin 10000

Result

Uncertainty

0.7813938

4.867925e-04

! ! ! ! ! ! ! ! ! ! ! ! !

Table 2.4 Evaluation of integral using Random Walk Method !

12!

3 Use of Monte Carlo Methods for Evaluation of Thermodynamic Properties of One MultiState System Problem 3.1 Consider an object (a spin, or a segment of a step on a crystal surface with two types of kinks, if you like), which has three internal states with energies εn= (n1)Δ, where n = 1,2,3, and Δ is a constant energy. It is a straightforward matter to evaluate the mean energy and specific heat of this system at a temperature T from the partition function, which is 3

Z (T ) = ∑ e− ε n /kT

(3.1)

n=1

If one does that and tabulates the results, one finds the information given in the file a004.dat given in the appendix. Temperatures are in units of Δ/k; energies, in units of!Δ; and specific heats, in units of k (Boltzmann's constant). (Try to obtain this data on your own) Mean Energy: The Mean Energy of the system is given by

∑εe ε= ∑e

− βε i

i

i

(3.2)

− βε i

i

Specific Heat: The Specific Heat of the system is given by ⎡ ε2 − ε 2⎤ ⎦ C=⎣ 2 kT

!

(3.3)

13!

Figure 3.1 Comparison of Energy and Specific Heat with the reference value

Problem 3.2 Alternatively, one may compute the energy and heat capacity by Monte Carlo methods. Do that for an appropriate set of temperatures using the method of importance sampling. Start from any of the three states and let the system make a biased random walk through the states. Do the latter as follows: At any given point in the simulation, the system is in some state i. Now attempt a move to another state j. First, pick, with equal probability, either of the other two states as the one ( j ) ! to which move is to be attempted. Accept that move with probability unity if ε j < ε i and exp ⎡⎣ − ( ε j − ε i ) / kT ⎤⎦ if ε j > ε i . Accumulate the energy and squared energy of the system after each attempted step; remember, if the attempted move is rejected, the system remains in the state i and the energy and squared energies of this state are added to the corresponding sums. Use a reasonable number of steps (your choice) in this biased random walk and bin the results for purposes of computing uncertainties in the energy and heat !

14!

capacity. Fro...