Skip Navigation | ANU Home | Search ANU | Search FEIT | Feedback

COMP3320/COMP6464 Molecular Dyanmics Project

Worth 20%, due by 18:00 on Friday 18 April 2008

The aim of the project is to develop an elementary molecular dynamics (MD) code and associated visualization software. MD is extremely widely used in computational chemistry, physics, biology, material science etc. It was also one of the first scientific uses of computers going back to work by Alder and Wainwright who, in 1957, simulated the motions of about 150 argon atoms. (See Alder, B. J. and Wainwright, T. E. J. Chem. Phys. 27, 1208 (1957) and J. Chem. Phys. 31, 459 (1959))

In MD each atom is represented as a point in space. The atoms interact with each other via pairwise potentials, in a similar manner to gravity interacting between the planets, sun, moons and stars. In practice this is an approximation since for atoms quantum effects are often important - but we'll ignore this and confine our discussions to "classical" MD.

Gravitational bodies interact via a potential that goes as m₁m₂/r₁₂ where m₁ and m₂ are the masses of the bodies involved and r₁₂ is the distance between their centres of mass. For the purpose of this work we will use a potential of the form

V₁₂ = 1/(r₁₂)¹² - 2/(r₁₂)⁶

To evaluate the interaction for many bodies we must simply perform a sum over all unique two body interactions. Thus for 3 particles we will have interactions 1-2, 1-3 and 2-3.

Consider now a cubic box with a side length of R. Along each dimension of the box we arrange N atoms as indicated by the following figure (not all atoms are shown!), such that there are a total of N*N*N atoms.

For this configuration of atoms we could calculate the total interaction energy arising from the above potential - we term this the "potential" energy (PE) of the system.

---

STEP 1: By Lab in Week 3

---

Your first objective is to produce a basic code that can read as input N (the number of atoms in each dimension) and R (the size of the box in each dimension). Your code should also compute and print out the value for the potential energy.

Develop this code using Visual Python (VPython).

• Visual Python
• Python
Your code should display a visual window showing all the atoms neatly arranged in a cube, and you should printout the total potential energy (in the text window). Note to run VPython in the lab type "vpython". This will give you an idle window. You type code in another window - to obtain this go to the file menu and use the "new window" option. If you select the "open" option it should take you to all the vpython demos that I showed in lectures. (They are in /usr/local/lib/python2.3/site-packages/visual/demos).

I will ask you to demonstrate your code to me in your lab in week 3.

---

STEP 2:

---

Given the pairwise potential there will be a pairwise force. This corresponds to the derivative of the potential with respect to R₁₂. Assuming the location of each atom is given in Cartesian coordinates the Cartesian forces are given by:
F_x = (12/r₁₂¹⁴ - 12/r₁₂⁸)(x₁-x₂)
F_y = (12/r₁₂¹⁴ - 12/r₁₂⁸)(y₁-y₂)
F_z = (12/r₁₂¹⁴ - 12/r₁₂⁸)(z₁-z₂)


The total force on a given atom will be the sum of all the pairwise forces.

Assuming that the forces on at least some of the atoms are non-zero then the atoms will want to move. MD is the process of moving the atoms according to the forces they are subjected to. Moving the atoms equates to integrating the equations of motion over time or solving a differential equation. There are a variety of ways in which we can do this, but a popular method is the so called "velocity verlet" algorithm (W.C Swope et al J. Chem. Physics 76 637 (1982)). This has the following form

coord(t+tstep) = coord(t) + tstep * (velocity(t) + tstep*accel(t)/2)
velocity(t+tstep) = velocity(t) + tstep * (accel(t) + accel(t+tstep))/2

where t is the current time, and tstep is the timestep (the integration step). And by assuming that all particles have unit mass we can replace the acceleration by the force (F = ma).

When the particles begin to move they will gain kinetic energy (KE = 1/2 mv²) with the total KE just a simple sum of all the atomic KEs. If there is no source or sink of energy then the total energy (PE+KE) should remain constant, so the gain in KE is offset by a lose of PE.

Extend your Python code to compute the force on each atom and allow the particles to move. Your code will need to read in a timestep and a total number of timesteps. Your code should also compute and printout the potential, kinetic and total energy at each timestep.

You should have this working by lab in week 5 at the latest.

HERE is some sample output for you to verify that your results are correct

---

STEP 3:

---

Re-write the compute intensive parts of your VPython code in C. The VPython code should read in the various parameters and display the positions of the atoms as they move. The VPhython code should communicate the starting parameters to the C code. This should then run the MD simulation and periodically communicate the coordinates of the atoms back to the VPhython program to be displayed. It is not expected that you communicate the coordinates at every timestep.

In addition to the visual output, you should provide text output that details the input parameters, and provides the total energy, potential energy and kinetic energy periodically during the course of the simulation.

Within the Python portion of your code use time.clock() routine to report the total time taken to perform the simulation. The start time is from the beginning of the first simulation step, while the end time is after the last iteration. This means you do not have to include setup time, but once the clock is started it should not be stopped again until everything is complete.

• COMP3320 People! Your simulation should correspond to an isolated cluster, there should be no distance cutoff used and there should be no use of symmetry.
• COMP6464 and Honours People! As for COMP3320, but your code should also include the ability to do simulations with periodic boundary conditions. ( DON'T PANIC - WE WILL DISCUSS THIS IN CLASS SOON). Whether or not to use periodic boundary conditions should be determined by an input parameter. The size of the periodic box (Sb) should be read as input and should be greater or equal to the size of the original box (R), with the original cube of atoms positioned in the center of the periodic box. A cutoff of Sb/2 should be used to ensure each atom only interacts with one image of any given atom.

  Here is a small code that given two atoms (i and J) determines what the minimum distance between atom i and any of the periodic images of j.

You are required to submit a tar file using the submit tool on one of the student computers:

    submit comp3320 md_project md.tar

or if you are a COMP6464 student

    submit comp6464 md_project md.tar

The tar file should expand to give a directory that contains your source code PLUS a plain text file called README that tells the user (me) exactly how to run your code and either contains or points the user to some sample input/output data so they can verify that your code works (or gives identical results to what you got). The tar file should also contain another plain text file called WRITEUP that details how long it took you to complete this part of the project and details of who you collaborated with (if anyone).
• For all students your WRITEUP file should also give a critical analysis of your code (good things, bad things, why you did it the way you did etc).
• For COMP6464 and honours students your WRITEUP should include some discussion comparing simulations that use periodic boundary conditions with those that don't (e.g. speed, what happens etc).
Your code will be tested on partch - so make sure it works on this system. You can use anything you find generally available (no something you installed in your directory) on partch (including libraries). Marks will be awarded as follows:
• 7/20 marks on the content of your WRITEUP file
• 7/20 marks for a functionally correct code that is documented well (in README and source files).
• 6/20 marks for the measured performance of your code
Your code will be tested on partch - so make sure it works on this system. You can use anything you find generally available (no something you installed in your directory) on partch (including libraries). Marks will be awarded as follows:
• 7/20 marks on the content of your WRITEUP file
• 7/20 marks for a functionally correct code that is documented well (in README and source files).
• 6/20 marks for the measured performance of your code