Coursehero >> Michigan >> Michigan State University >> PHY 480

Course Hero has millions of student submitted documents similar to the one below including study guides, homework solutions, papers, exam answer keys and textbook solutions.

Dynamics Molecular Lecture 2 As described in the rst lecture, we are using the Verlet algorithm to integrate Newton s equations for particles interacting through the Lennard-Jones interaction. The Lennard-Jones pair potential between two atoms and is given by, V (r ) = 4 [ r 12 r 6 ] (1) For Argon the parameters of the model are m = 6.69 10 26 kg, = 1.654 10 21 J, = 3.405 10 10 m. The minimum of the potential is at rm = 21/6 , and the depth of the minimum is . The time scale of the problem is = (m 2 / ) = 2.17ps. Since we set m = 1, = 1, = 1, one unit of time in our simulations corresponds to 2.17ps. Since the timestep we use is .004, it corresponds to 8.65 f s. This is typical of simulations of atomic motion and in cases where the interactions are strong, the timestep has to be shorter, for example in covalently bonded systems. We want to study the behavior of the model for Argon as a function of temperature. At low temperatures Argon is a fcc crystal and as the temperature is raised, it melts into a liquid. To see this transition, we initialise the system as an fcc crystal and give the atoms a Maxwell-Boltzmann distribution of velocities. We then using the Verlet algorithm you already developed to nd the time evolution of the system. The rst step is thus to set up the fcc lattice, to impose periodic boundary conditions and initialize the velocities. We use the Box-Muller method to x the velocities to a Maxwell-Boltzmann distribution and then scale the velocities to ensure that the center of mass velocity is zero. To set the temperature we can use the relation between the kinetic energy and the velocities given by the Maxwell-Boltzmann distribution, 1 3 ( N 1) k B T = m i | v i |2 2 2 i (2) You will nd that you cannot set the initial temperature accurately, so you need to rescale all of the velocities by a constant factor in order to adjust the temperature. After 4 or 5 temperature rescalings you should reach a value close to the target temperature. Things to calculate 1 From the MD simulation we can calculate a variety of physical properties. Choose two of the following. (i) The total energy as a function of temperature. (ii) The average mean square displacement of each atom as a function of time, for several temperatures. (iii) The pair correlation function, for several temperatures. To nd this function make a histogram of the values of rij as a function of r with a box size r = 0.1 and count the number of atom pairs that lie in each box. Also (iv) The pressure is given by, pV = Nk T B 1 rij Fij 3 i< j (3) where rij is the vector between atoms i and j and Fij is the force between atoms i and j. This equation is due to the virial theorem, and (v) The speci c heat at constant volume, CV , is found from the kinetic energy uctuations using a formula due to Lebowitz, < K2 > 2 3N = (1 ) 2 K 3N 2CV (4) Finally the pair distribution function g(r ) is the probability that two atoms are separated by distance r. To nd this function make a histogram of the values of rij as a function of r with a histogram interval size r = 0.1 and count the number of atom pairs that lie in each interval. Tricks of the trade 1. Setting up the initial fcc lattice. A code for this is posted as is a Mathematica code to visualize the structure. 2. Setting up the initial velocity. For this we use the Box-Muller method. We want to choose velocities from the Maxwell-Boltzmann distribution, m ( v2 + v2 + v2 ) m x y z 3/2 ) exp[ ]. P(v) = ( 2 k B T 2k B T (5) The Box-Muller method uses random numbers drawn from [0,1] to generate random numbers drawn from the standard normal distribution exp[ u2 /2]/(2 )1/2 . See the Wikipedia page on this to see how to do it. Once you have the random numbers from this normalized distribution, you can nd a velocity component using v = (k B T/m)1/2 u. 2 You also need to set the center of mass velocity to zero and to do this, sum up all the velocities (as a vector), divide by N (the number of atoms) and subtract this vector from each vector velocity. During the computation check that the center of mass velocity stays zero. 3. Periodic boundary conditions. These have to be used in two ways for a system which has size box in each direction. (i) First when calculating the distance between two atoms, you have to take the smallest distance taking into consideration the boundary conditions. This can be done for the x,y and z directions separately and is achieved using the following code for each atom pair, dr (:) = r ( j, :) r (i, :) dr (:) = dr (:) N I NT (dr (:)/box ) box (ii) The second way in which the boundary conditions come in is that after the new positions are calculated in the Verlet algorithm, we have to check to see if any particles are outside the box. If they are, they must be moved to their appropriate position on the other side of the box. This is achieved using, x (:, :) = MODULO( x (:, :), box ) 4. Check energy conservation. The total energy (KE + PE) is conserved to six digits and the center of mass momentum stays zero to the same accuracy. These are essential checks on the code 3

Find millions of documents here - Study Guides, Homework Solutions, Papers, Exam Answer Keys and more. Course Hero has millions of course related materials that will enable you to learn better, faster and get an A in all your courses.
Below is a small sample set of documents: