COMP3320/COMP6464 Molecular Dyanmics Project

Worth 20%, due by 18:00 on Wed 21 April 2010

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:

Develop this code using Visual Python (VPython).

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 that your python code, call it md1.py should accept the following arguments in order, i) N - number of atoms in each dimension of the box ii) R - the size of the box in each dimension, iii) tstep - the time step, and iv) iter - the number of time steps in your simulation. It should therefore run as follows: python md1.py N R tstep iter.

STEP 2:

₁₂

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:

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.

STEP 4:

Step4_Report.pdf

Assuming that the total energy before the first iteration is exact, present an analysis of numerical errors due to the following factors: i) increasing number of atoms, ii) increasing number of iterations, and iii) increasing timestep size. For each of these factors, identify and discuss the most likely source of numerical errors. (4 marks)
What is the per-iteration (i.e. each timestep) execution time of your code as a function of different problem dimensions? (2 marks)
Furthermore, how does your code scale as a function of the problem size? (e.g. is it O(NlogN),O(N)?) (2 marks)
What is the most computationally intensive component of your code in terms of time? (2 marks)
During each iteration, what proportion of the total execution time is spent calculating the new i) coordinates, ii) forces, and iii) velocities? How does this proportion change as problem dimension increases? (2 marks)
Does the most computationally intensive component correspond to the lines of code that is executed most often? (2 marks)
Show some hardware performance counter data of your C code, these may include metrics such as Cycle counts, L1/L2 cache miss rates, Instr/Cycle, and branch mispredictions, whatever appears worthwhile to you. Explain the significance of each metric that you have chosen, and analyse how they vary with different problem dimensions and how they correlate with observed execution times (if applicable). [You can do this on either Wallaman or Partch] (4 marks) (Download perfctr.tar.gz for perf counter code)
Develop a modification of your code (e.g. change the array layout or access pattern) and compare the performance with your main version using any of the tools you have been introduced to during the labs and lectures. (2 marks)
For COMP6464 and Honours Students Only: Suppose now you want to convert your molecular dynamics code into a benchmark similar to the LINPACK benchmmarks.
- Briefly describe the implementation of the LINPACK benchmark, and how the TOP500 list is compiled using LINPACK. (1 marks)
- Using how LINPACK calculates the number of floating point operations (FLOP) as a model, devise a method for measuring the MFLOP/s in relation to your MD code. Specifically, outline how you will calculate the number of FLOPs for a given problem dimension. (2 marks)
- Using your method for calculating MFLOP/s, determine the values of N_max, and R_max for your MD code on Partch and Wallaman. Discuss why the particular dimension N_max yielded the best performance. (4 marks)
- What is the theoretical peak MFLOP/s (R_peak) of Partch and Wallaman? Discuss how this compares with the actual maximum MFLOP/s (R_max)? (1 marks)
- What aspects of the computational architecture does your benchmark assess? What aspects does it not assess? (For example the STREAM benchmark is designed to assess memory bandwidth, but not cache performance) (2 marks)

To complete this step, we envisage that you will need to modify your C code and its Makefile substantially (ie. add calls to timers and performance counters, change compiler flags to enable gprof, etc) you must ensure that you provide sufficient details in your report so that we can easily reproduce the performance results that you reported WITHOUT having to modify your code or Makefile. Marks will be deducted if this requirement is not met.

Formatting of your program

All your submitted code should include this header at the top:

Name:
Student Number:
Course:
Assignment Number:
Name of this file:
Lab Group:

I declare that the material I am submitting in this file is entirely
my own work.  I have not collaborated with anyone to produce it, nor
have I copied it, in part or in full, from work produced by someone else.
I have not given access to this material to any other student in this course.

By Stage 3, you should have at least two python files, one completely written in python called md1.py and another called md2.py which calls C code.
I must be able to run your code from the command line by typing python md1.py N R tstep iter as explained above.
Your code should output to stdout the total energy, potential energy, and kinetic energy at every timestep.

The output for md1.py and md2.py should be in the same format as the sample output: HERE, ie.

      Dimension of Cube 2.000000
      Total atoms       27
      Time step         0.001000
      Number of steps   10

      TIME         Total Energy            PE                 KE
      -------------------------------------------------------------------
      0.000000      -7.577011e+01      -7.577011e+01       0.000000e+00 
      0.001000      -7.577011e+01      -7.577026e+01       1.504815e-04 
      0.002000      -7.577011e+01      -7.577071e+01       6.018923e-04 
      0.003000      -7.577011e+01      -7.577147e+01       1.354131e-03 
      0.004000      -7.577011e+01      -7.577252e+01       2.407030e-03 
      0.005000      -7.577011e+01      -7.577387e+01       3.760354e-03 
      0.006000      -7.577011e+01      -7.577552e+01       5.413798e-03 
      0.007000      -7.577011e+01      -7.577748e+01       7.366993e-03 
      0.008000      -7.577011e+01      -7.577973e+01       9.619501e-03 
      0.009000      -7.577011e+01      -7.578228e+01       1.217082e-02

Submission and Assessment:

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). This need not be as detailed as the report for Step 4. 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).

Your code will be tested on partch (and on wallaman if you use performance counters) - 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 and wallaman (including libraries).

Your assignment will be marked out of 40 for COMP3320 students, and out of 50 for COMP6464 students. This is worth 20% of your total overall course score. The marks will be distributed as follows:

Appearance of python simulation (4 marks)
Correctness of your code in calculating energy values, forces, velocities (8 marks)
Coding style, commenting, interfacing in C and Python code (8 marks)
The distribution of marks for the report is given above (20 marks COMP3320, 30 marks COMP6464)

	Skip Navigation CECS Home \| ANU Home \| Search ANU
	ANU College of Engineering and Computer Science School of Computer Science