22.53 Problem Set I solution

TA: Jingli Liu

November 6, 2000

1 Problem 1 Solution

1.1 Simulate solid at low T without T scaling

1.1.1 Energy equipartition

Statement The average kinetic energy of molecules in the system stands for the temperature. In our system( Lennard-Jones potential interaction), there are three degrees of freedom for each molecular: x,y and z direction. Thus, when having reached equilibrum, we should expect that the average kinetic energy should be divided equally by these three degrees. That means each componet of the average kinetic energy makes the same contribution to the temperature. In mathematics, we can express as:

< v_{x}^{2} > + < v_{y}^{2} > + < v_{z}^{2} >= \frac{kT}{m}  (1)

where k is Boltzman constant (1.38x10^{−23}J/K), T is the dimensional temperature and m is the mass of atom ( for argon, the atomic weight is 39.948g/mol, so the mass of the atom is 6.634x10^{−26}kg). We use the Haile code to calculate these three averages with some small modifications. Our parameters are as follows:

• NP (Number of Particles): 108

Figure 1: Energy Equipartition

• NEQ (Number of Equilibration Timesteps): 1000 • MAXKB (Number of Timesteps): 3000

• TR (Reduced Temperature): 0.1

• DR (Reduced density): 1.1

Results Here the unit of horizon axis is δt=δT imeUnits From our calculateions obtain (reduced)

< v_{x}^{2} >= 0.0983  (2)
< v_{y}^{2} >= 0.0952  (3)
< v_{z}^{2} >= 0.0900  (4)

Temperature = 0.10. So,< v_{x}^{2} >≈< v_{y}^{2} >≈< v_{z}^{2} >≈ T.

Discussion Actually this is a special case of a more general equipartion theorem of statistical mechanics. If there are more dgrees of freedom, such as rotation, vibration, etc., each degree wil have the average energy which equals to \frac{1}{2}KT, and total average kinetic energy will be a product of degrees of freedom and \frac{1}{2}KT. Thus, the total kinetic energy in the system can be evaluated. And it is a little different from the potential energy, which can only be concerned with the configuration of the molecules.

1.1.2 Reversibility

Statements Let’s investigate such a solid which internal molecules aare assumed to obey Newton’s law of motion. In this case, we will think of their reversibility. Having chosen 32 particles, we will simulate a solid; after running 10 time steps, we suddenly change the signs of the velocities. then, we can expect all the particles go right back. The same time later, whole system will be back to its original position except all the velocities are inverted.

• NP: 32

• NEQ: 0

• MAXKB:22

• TR: 0.6

• DR: 1.1 Result We record all these 32 particles original positions and velocities. We suddenly change the signs of velocities and also X3,X5,Y3,Y5,Z3,Z5 at the 11th step. We then will see at the 20th step, the velocities of the particles are reversed. The following two tables are the positions and velocities of particles of the 1st step and the 20th separately.

\begin{matrix}X & Y & Z & V_{x} & V_{y} & V_{z} \\ 3.073499 & 3.074764 & 0.006527 & -0.403826 & -0.151734 & 1.301584 \\ 0.001265 & 0.770241 & 0.775623 & 0.252414 & 0.273412 & 1.345026 \\ 0.770666 & 0.000646 & 0.775632 & 0.356531 & 0.127730 & 1.346892 \\ 0.766012 & 0.762923 & 3.075463 & -0.573789 & -1.188352 & -0.011627 \\ 1.533708 & 0.002945 & 0.004475 & -0.807977 & 0.585949 & 0.893388 \\ 1.536895 & 0.763510 & 0.765245 & -0.171153 & -1.070617 & -0.725430 \\ 2.311994 & 3.072210 & 0.766303 & 1.065738 & -0.660597 & -0.515103 \\ 2.311720 & 0.766110 & 3.072026 & 1.013313 & -0.553555 & -0.698098 \\ 0.001275 & 1.533760 & 3.070478 & 0.254681 & -0.798802 & -1.005693 \\ 3.071372 & 2.311412 & 0.770797 & -0.825548 & 0.949798 & 0.381618 \\ 0.771730 & 1.534116 & 0.774141 & 0.567580 & -0.727539 & 1.050257 \\ 0.764010 & 2.305584 & 3.071556 & -0.972209 & -0.211097 & -0.789458 \\ 1.533540 & 1.536959 & 3.070822 & -0.843675 & -0.161316 & -0.938241 \\ 1.532865 & 2.304212 & 0.772720 & -0.978470 & -0.484450 & 0.765391 \\ 2.310007 & 1.543589 & 0.768089 & 0.670678 & 1.161128 & -0.159013 \\ 2.302862 & 2.304667 & 0.005257 & -0.755682 & -0.394555 & 1.047448 \\ 0.001221 & 0.001858 & 1.531650 & 0.243843 & 0.370387 & -1.217968 \\ 0.002794 & 0.771460 & 2.301258 & 0.557214 & 0.513466 & -1.075796 \\ 0.771803 & 0.001659 & 2.312815 & 0.581739 & 0.331035 & 1.231975 \\ 0.766615 & 0.762668 & 1.537440 & -0.450725 & -1.238472 & -0.062982 \\ 1.542602 & 0.002791 & 1.542013 & 0.966641 & 0.556860 & 0.846990 \\ 1.536196 & 0.763930 & 2.302721 & -0.312573 & -0.987248 & -0.782716 \\ 2.307537 & 3.070804 & 2.302152 & 0.180669 & -0.942747 & -0.894878 \\ 2.307256 & 0.764556 & 1.532883 & 0.121800 & -0.861581 & -0.974415 \\ 3.073480 & 1.543000 & 1.541513 & -0.408008 & 1.046169 & 0.748560 \\ 0.005977 & 2.305783 & 2.303527 & 1.190499 & -0.169798 & -0.621606 \\ 0.766213 & 1.543752 & 2.307568 & -0.530750 & 1.194072 & 0.182839 \\ 0.770337 & 2.312963 & 1.539869 & 0.290295 & 1.262377 & 0.420401 \\ 1.539862 & 1.544159 & 1.538647 & 0.419477 & 1.274989 & 0.177056 \\ 1.537355 & 2.311009 & 2.301856 & -0.080381 & 0.872038 & -0.954273 \\ 2.301967 & 1.535125 & 2.310641 & -0.933737 & -0.525780 & 0.797142 \\ 2.308222 & 2.309693 & 1.532199 & 0.315391 & 0.608834 & -1.109270 \end{matrix}

5

\begin{matrix}X & Y & Z & V_{x} & V_{y} & V_{z} \\ 3.073500 & 3.074765 & 0.006528 & 0.403601 & 0.151681 & -1.300589 \\ 0.001265 & 0.770241 & 0.775624 & -0.252260 & -0.273748 & -1.344167 \\ 0.770666 & 0.000646 & 0.775633 & -0.356392 & -0.127338 & -1.346074 \\ 0.766012 & 0.762922 & 3.075464 & 0.573776 & 1.187509 & 0.011519 \\ 1.533707 & 0.002945 & 0.004475 & 0.807243 & -0.585194 & -0.892975 \\ 1.536895 & 0.763509 & 0.765245 & 0.170554 & 1.069627 & 0.724990 \\ 2.311996 & 3.072210 & 0.766303 & -1.064643 & 0.660107 & 0.515004 \\ 2.311720 & 0.766110 & 3.072027 & -1.012723 & 0.553364 & 0.697702 \\ 0.001275 & 1.533760 & 3.070478 & -0.254625 & 0.798342 & 1.004754 \\ 3.071370 & 2.311414 & 0.770797 & 0.824290 & -0.948726 & -0.381190 \\ 0.771731 & 1.534116 & 0.774142 & -0.566935 & 0.727058 & -1.049735 \\ 0.764010 & 2.305584 & 3.071554 & 0.971704 & 0.210939 & 0.788378 \\ 1.533541 & 1.536960 & 3.070822 & 0.843478 & 0.161545 & 0.937661 \\ 1.532866 & 2.304211 & 0.772721 & 0.978253 & 0.484004 & -0.764759 \\ 2.310007 & 1.543590 & 0.768089 & -0.670075 & -1.160001 & 0.159231 \\ 2.302861 & 2.304667 & 0.005258 & 0.755473 & 0.394369 & -1.046398 \\ 0.001221 & 0.001858 & 1.531649 & -0.243828 & -0.370081 & 1.216802 \\ 0.002794 & 0.771461 & 2.301258 & -0.556816 & -0.512780 & 1.075487 \\ 0.771804 & 0.001659 & 2.312816 & -0.581005 & -0.330863 & -1.231342 \\ 0.766614 & 0.762666 & 1.537439 & 0.450076 & 1.237325 & 0.062600 \\ 1.542603 & 0.002791 & 1.542014 & -0.966306 & -0.556547 & -0.846172 \\ 1.536196 & 0.763929 & 2.302722 & 0.312440 & 0.986471 & 0.782320 \\ 2.307537 & 3.070803 & 2.302152 & -0.181136 & 0.942434 & 0.894012 \\ 2.307255 & 0.764555 & 1.532883 & -0.121606 & 0.860678 & 0.974029 \\ 3.073481 & 1.543001 & 1.541513 & 0.407803 & -1.045766 & -0.748166 \\ 0.005978 & 2.305782 & 2.303527 & -1.189182 & 0.169249 & 0.621143 \\ 0.766212 & 1.543753 & 2.307568 & 0.530013 & -1.193003 & -0.182293 \\ 0.770337 & 2.312964 & 1.539868 & -0.290062 & -1.261960 & -0.420130 \\ 1.539863 & 1.544160 & 1.538647 & -0.419287 & -1.273736 & -0.177030 \\ 1.537355 & 2.311010 & 2.301855 & 0.080124 & -0.871770 & 0.953437 \\ 2.301967 & 1.535124 & 2.310642 & 0.933348 & 0.525304 & -0.796448 \\ 2.308222 & 2.309693 & 1.532199 & -0.315295 & -0.608494 & 1.108398 \end{matrix}

6 Discussion Since our equations of motion are symmetric in time, that will inevitably cause the reversibility. For a simple example, let’s look at some particle moving under Newton’s law:

\frac{d^{2}r}{dt^{2}} = F  (5)

And after solving that equation, one would obtain a trajectory of this particle as r = r(t). Suppose we can let time go back, which allows us to live in a time of t′ = −t. Then we are very surprised to find out that the form of Newton’s law doesn’t change. And that’s the origin of the reversibility . It’s thus very naturally to check our code to see if it can perform the rever sibility. this is a very simple but convenient way. Also , I have more to speak is about the signs of velocities and higher derivat ives. Our transformation is t′ = −t, and \frac{d^{n}r}{dt^{′n}} = (−1)^{n} \frac{d^{n}r}{dt^{n}}, therefore, we should not forget that we have to change 3rd and 5th order derivatives, not only the 1st derivatives( which are velocities).

1.1.3 Checking Energy Conservation by varying δt

Statement

• NP: 108

• NEQ: 1000

• MAXKB: 1000

• TR: 0.6

• DR: 1.1

In order to be compared, four different values of δt( 0.001, 0.005, 0.007, 0.01) are chosen.

From these four different values of δt(0.001, 0.005, 0.007, 0.01), we find that 0.005 and 0.001 are better. At δt =0.01, we can see that the total energy increase with time and it is very abnormal. It shows that the time step is too large that the displacement is too far. At Figure 2: Energy Conservation

8 this time, the conservation of energy is increasingly damaged. However, if we would like to extend our time scale, we have to improve the δt as large as possible. Here, we will have to deal with this dilemma. And I think the most important thing is to have a direct image what this simulation is, then we will have to get some practices from some tests to know which value is the best.

1.2 Study a solid at several different Temperatures

Our simulation are carried out at the following conditions:

• NP: 108

• NEQ: 1000

• MAXKB: 4000

• DR: 1.2

And we choose four different temperatures as 0.5, 1.0, 2.0 and 1.2(reduced temperature). Running results are the follows. From figure 3 we can see that the configuration energy increases as the temperature goes up. However, at high temperatures, for example T R = 5.0, the configuration energy no longer conserves. This can hint us that something has changed inthe system. At least, our system is no longer stable at high temperatures. Recalling our knowledge of phase transition, we can inger from the configuration energy’s changes that there exists some similar changes in the system, because at this time temperature is constant while the energy of the system increases. The fig 4 shows the variation of pressure. From Figure 5, we can see that the higher the temperature is, the biger the Mean Square Displacement is. this is reasonable since the temperature indicates the average kinetic energy of the molecules. when the molecule has bigger kinetic energy, it will have moretendency to break out the control otf the bound potential. At low temperature, it might only vibrate around its balance position. But as temperature goes sufficiently high, it has the chance to Figure 3: < U_{c} >∼ t

10 Figure 4: < P >∼ t

11 Figure 5: < r^{2} >∼ t

12

Figure 6: Distribution of G(R)

go out to move freely. And that is a basic feature of liquids. so, in some sense, we can say that the solid is melting. Next, we will compare the g(r) of under these four temperatures with the ideal FCC structure, which means T=0 and the atom is strictly bound in its position. As the figure 6 has shown, the function g(r) becomes more and more flat as temperature goes up. The ideal FCC structure is a sum of δ functions. So its figure is feally very sharp and high. when temperatuer increases, the particles moves back and through. Therefore, g(r) becomes more and more flat just as the figures have shown. The number of peaks decreases and it is more difficult to find which position is peak.

1.3 Determine the thermal expansion coefficient at approximately zero pressure

Thermal expansion coefficient is defined as:

ß = \frac{1}{V} (\frac{∂V}{∂T} )_{N,P}  (6)

According to this definition, we should simulate our Lennard-Jones system at constant pressure if we want to calculate from the above definition directly. However, our simulationis molecular dynamics simulation at constant temperature and constant temperature or energy. And, we cannot keep the pressure constant in Haile code. Thus we need to transform this definitionof thermal expansion coefficient to another form. First since we use the density as our parameter as our input parameter in Haile code, it is more convenient to use the density (ρ) instead of volume (V). Substituting V = \frac{N}{ρ} into (6) yields that

ß = −\frac{1}{ρ}(\frac{∂ρ}{∂T})_{N,P}  (7)

Second, we would like to change the derivative at constant pressure. In thermodynamics, we know that there is a relation among the state quantities (pressure, temperature, density) which is called the state equation, i.e.,

f(P, T, ρ) = 0  (8)

Thus,

(\frac{∂ρ}{∂T})_{N,P} = − \frac{\frac{∂f}{∂T}}{\frac{∂f}{∂ρ}} = (− \frac{\frac{∂f}{∂T}}{\frac{∂f}{∂P}} )(− \frac{\frac{∂f}{∂P}}{\frac{∂f}{∂ρ}} ) = −\frac{(\frac{∂P}{∂T})_{N,ρ}}{(\frac{∂P}{∂ρ})_{N,T}}  (9)

Thermal expansion coefficient can then be written as:

ß = \frac{1}{ρ} \frac{(\frac{∂P}{∂T})_{N,P}}{(\frac{∂P}{∂ρ})_{N,T}}  (10)

In Simulation, we will determine the quantities in (10)to get the final results of thermal expansion coefficient. In order to be at approximately zero pressure, we would choose high temperature and low density as follows:

• T R = T^{∗} = 2.0

• DR = ρ^{∗} = 0.01

In simulation, we have done as folowing procedures ( a star on the variable means it is in reduced units.):

1. Determine the current pressure at given temperature and density. The input parameters of Haile code are set as

NP=864, NEQ=4000, MAXKB=4000, TR=2.0, DR=0.01.

Then we got tge result of pressure as

P^{∗} = 0.020

2. Determine the term (\frac{∂P}{∂T})_{N,T} At this step we should vary the density instead. Let Haile code run at following

NP=864, NEQ=4000, MAXKB=4000, TR=2.0, DR=0.015

The resultant pressure is 0.030. So,

(\frac{∂P}{∂ρ})^{∗} = \frac{0.030 − 0.020}{0.015 − 0.010} = 0.5  (11)

And ß = 0.5/120.0 = 0.004167(K^{−1}).

This is one case calculation of thermal expansion coefficient at approximately zero pressure. From the table handed out in class, we can find that ß for argon at zero pressure is: ß = 0.3 The two data matches each other.