Author Archives: moabdelaziz

Molecular Dynamics Project Conclusions

Molecular Dynamics Project Summary

Sushant Mahat and Mohammed Abdelaziz

Over the last few weeks, we have been working on a computational project on molecular dynamics. This post below summarizes the details of this project, and the overall goals we achieved.


Initial Goals and Motivations:

We wanted to do an experiment in molecular dynamics as it connects to several things we have studied in our previous physics classes. We used our knowledge on Newton’s second law which we learned in classical mechanics; we used Boltzmann’s distribution and the concepts of temperature from thermodynamics and statistical mechanics; and we used the Lennard-Jones potential which has its root in quantum mechanics. All these features made simulating  molecular dynamics an attractive, challenging, and a very educational project.

We had several goals when we started this project. We wanted to be able to simulate a 2-D gas system, and be able to study several of its properties like temperature, energies, speed distributions, and visualization at equilibrium. We then wanted to build up on this system to make it into a 3-D system. We were also thinking about making it possible to increase the temperature of the system and observe how long it takes to come back to equilibrium. This was supposed to be a simulation of periodic heating of gas particles with a laser beam.


Goals We Completed

Due to time constraints, we were unable to fulfill some of our initial goals. However, we were able to finish most of the goals. The goals we met for this project include:

  1. Simulating the motion of a dilute sample of gaseous Argon atoms in two dimensions.
  2. Calculating the mean velocities of each atom and plotting their speed distribution in each simulation.
  3. Determining the temperature of the sample based on initial and equilibrium speeds of each atom, which required a determination of the kinetic energies of all atoms at each time step.
  4. Comparing the calculated speed distribution with the Maxwell-Boltzmann speed distribution curve.
  5. Completing the simulation and the above calculations in three dimensions.
  6. Introduced possibility to simulate temperature change in the system.



Section 9.1 of Giordano and Nakanishi’s Computational Physics was our primary reference for this project. We used MATLAB for the simulations, dividing our code into several functions (subroutines) instead of using one script. The codes and their purposes are outlined below.



This function initializes an evenly spaced grid in 3-D space (2-D in the old version of the program). This grid serves as the basis for the location of each atom, because we want the atoms to start at roughly the same location for each simulation. Two further steps are taken to make sure the simulations are not fully deterministic: each atom was displaced slightly from a point on this evenly spaced grid, and it was given a random velocity (within some limits) in each dimension. This function also allows an option to have different limits of the initial random velocity which can produce an effect similar to starting with gas particles with different temperature.



This function calculates the force on each atom due to every other atom, based on the Lennard-Jones potential, which is ineffective at large distances, attractive at moderate distances, and repulsive at very short distances. This calculation takes exponentially longer for additional particles, because each particle interacts with each other particle.



This function uses the forces from the previous function and each particle’s positions to calculate each particle’s position and velocity at the next time step. It feeds these positions back into TheEnforcer.m and iterates until the positions and velocities of each particle for each time step are calculated. This function also enforces the periodic boundary conditions (that a particle leaving one side of the “box” comes back through the other side, with its velocity unaffected). Particles that are closer to each other through these boundary conditions than they appear must interact as though they were at this close distance. To clarify this, consider a one-dimensional system, ranging from x=0 to x=10. Two particles at x=1 and x=9 should interact as if they are 2 units apart, not 8 units apart. Implementing this system properly was one of the most difficult tasks we faced.



This function takes in each velocity component of each particle at each time step, and first calculates a speed for each particle at each time step. It plots three speed distribution curves by averaging each particle’s speed over each third of the runtime, placing these speeds into bins, and plotting a histogram that is visualized as a curve. The function is configurable to plot any integer D speed distribution curves by averaging each particle’s speed over each 1/D of the runtime instead, but D = 3 effectively showed that the system had reached equilibrium (signified by when two successive speed distribution curves are fairly similar).

This function also plots the kinetic, potential, and total energy averaged for all particles at each time step in the simulation. The kinetic energy was needed to calculate temperature, since the product of Boltzmann’s constant and the system’s temperature is equal to the average kinetic energy of all particles in the system.



This script simply initializes the spatial dimensions of the system and how many atoms will be simulated. This is the code that is directly run by the user to visualize the simulation and look at the speed distributions and energies.

All of these MATLAB files can be found in this Google Drive folder.



We succeeded in producing visualizations of the Argon atoms in two and three dimensions. Even with 100 atoms, the script does not take more than thirty seconds to run, so more particles could feasibly be simulated. Additionally, our speed distributions matched up relatively well with the actual Maxwell-Boltzmann distributions in terms of horizontal location (this shows that our determination of temperature was accurate, since this is what affects the horizontal position of the peak (most common particle speed). However, our generated speed distributions were more narrow than the Maxwell-Boltzmann distribution. It is possible that they could be widened by decreasing the density of the atoms (simulating fewer particles or making the dimensions of the “box” bigger). It is also possible that increasing the maximum initial velocity would increase the width of the speed distribution, since it would give each particle a wider range of initial velocities to have.

We were mostly satisfied with our energy plots; they showed a fairly constant total energy, which we expect (it is not completely constant because the total energy is calculated from averages of kinetic and potential energy of all particles). Kinetic and potential energies became steady as the system reached equilibrium, which takes longer for lower gas densities.

One example of a speed distribution plot (click for full-size):


One example of an energy plot (click for full-size):



Conclusion and Future Potential

Based on the Boltzmann’s distribution plots of our data, we are fairly confident that our simulation is useful at approximating systems of gas molecules. The fact that the total energy of the system remains around a fixed value is also promising. However, the fluctuations in the total energy when the particles start interacting is still not ideal.

One limitation of our problem would be that the system destabilizes when the initial density of our gas molecules are very high. This is a direct result of the steep increase in the Lennard-Jones potential at low distances between particles.

We want our readers to see this code as a work in progress rather than a complete program. The program still has a lot of unrealized potential. For example, the energy graph produced by the simulation for cases with different initial velocity limits could potentially be used to calculate time taken by a system with gases of different temperature to equilibrate. Similarly, with some changes to the initiator and the enforcer function, the program could also be used to simulate a collection of different gases as opposed to a single type of gas molecule. Although we were unable to simulate it, increasing the density and lowering the initial velocity should also make it possible for future users to simulate materials in solid and liquid phases. Since the program already has method to change temperature built into it, some changes could be made to make these temperature changes periodic. This could potentially be used to simulate a laser beam heating the gas particles.

Overall, we are very satisfied with the progress we were able to make throughout the semester. Apart from learning a great deal about different physical concepts, we also developed a deeper appreciation and knowledge of computational approaches towards studying physical systems. We also think that we have developed  a fair understanding of the limitations and the potential of the computational tools available to us, and we are confident that this project has prepared us for many more simulations that we are likely to encounter in our careers as physicists.         


Molecular Dynamics – Results and Analysis

Sushant Mahat

Mohammed Abdelaziz

Since last week, we have managed to greatly fine tune our program, and get more data.

We are using the same basic functions we created last week: TheInitiator and TheUpdater. However, as mentioned in the earlier post, our TheUpdater function was not functioning ideally. The profile of the Lennard-Jones potential is such that when two particles get close together (less than ~1.1 reduced unit in length), they exert a large repulsive force on each other and send each other moving away at large velocities. If you look at the Lennard-Jones potential in Figure 1 below, you can see that the force changes from attractive to slightly repulsive to extremely repulsive very fast (force changes by ~100 reduced units of force over ~0.3 reduced units in length). In real life, two particles repel each other before they get this close, so the force should never become very large. However, our function discretizes time, so there is a small chance that particles get very close to each other before they feel the repulsive force. This creates a repulsive force that catapults one of the particles, which breaks the law of energy conservation. The catapulting effect decreased when we decreased the time step dt (our discretization in time), but it did not go away completely.

To solve this problem, we divided TheUpdater into two parts: TheUpdater and TheEnforcer. TheEnforcer function calculates the force between each pair of particles and the total force on each particle. However, this function is designed such that if two particles accidentally get close to each other during dt, the force they exert on each other will be limited to a certain value (i.e. we cap the max force). We think this change is valid as two particles, with the velocities we use, would never get so close that they would exert more than our max force on each other. Apart from that, we have also made changes to the net force direction. In our last program, we were able to determine repulsive or attractive forces but it turns out we made an error in translating this information to force in the x direction and the y direction. This is especially hard as we have periodic boundary conditions, so it is not always clear whether a particle is “to the left” or “to the right” of another particle, for instance. However, we are confident that we have managed to get the force direction correct this time. We tested it by running a program with just 5 particles in a small box (size of 10* 10), so that it would be easier to monitor their interactions. During our several runs, we saw several different interactions where our particles curve around each other and attract at certain distances and then repel, as expected, which is a significant improvement to our last program. Video 1 showcases this test.

TheUpdater function now calculates positions and updates the velocity of our particle using the output from TheEnforcer. These data are fed into the main body script which calls upon another function: SpeedDist. SpeedDist finds the speed distribution of our particles. Figure 2 below shows that the average speed remains very close to 1 reduced speed unit, which was our initial average speed. In the coming weeks, we will also include another plot that shows the kinetic energy of all particles with respect to time, and the cumulative Kinetic energy of all particles. This information will be interesting because it can directly be used to calculate the temperature of the system, potentially allowing us to test the effects on the system when we manipulate the temperature.

Figure 1. This figure shows the Lennard-Jones potential. The x axis is length in reduced units. The reduced unit of length is based on the size of an argon molecule. The steep increase in potential below a distance of 1 suggests a rapid increase in repulsive force as two particles come close to each other. This essentially creates a “hard sphere” around each particle.



Figure 2. Speed distribution of 100 particles in a 25*25 box. The average speed remains roughly constant over the time interval and closely resembles Maxwell distribution. The total time for the simulation was 10 reduced time units where 1 reduced time units is approximately 2 picoseconds.

Video 1. This is a video of 5 particles in a  10*10 box. You can see the particles interacting in several ways.

Video 2. This is a video of 25 particles in a 10*10 box. It demonstrates that the interactions still work for greater numbers of particles (we have found that even simulations of 200 particles can be calculated in under a minute).


A note on reduced units:

Reduced units are used for convenience and simplicity in the equations used for the simulations because they avoid the need of having to put constants in every calculation. We have been somewhat vague about what each reduced unit represents in terms of common units, but will explicitly define these in a future post.

Edit (4/27/15):

We forgot to upload our most recent working script and functions; they can be found here.


Preliminary Results for Molecular Dynamics Simulations

Sushant Mahat and Mohammed Abdelaziz

We created two functions to simulate the movement of dilute gas molecules in a container.

The first function initiates a grid of equally spaced points on which to initially place each particle. The function gives them random velocities and perturbs their initial position slightly from the equally spaced grid points.

The second function is a routine that updates the position of each particle at each time step. It determines the next position of each particle by determining the force acting on it due to all other particles in the box, based on the Lennard-Jones potential. The Lennard-Jones potential leads to a force that is large and repulsive at short distances, and attractive at larger distances. The repulsive force has been problematic, because two particles that encounter each other in the box get repelled at very large speeds, producing isolated points like those seen in our results below. If we can solve this issue, the results should be much more coherent.

Additionally, the simulation uses periodic boundary conditions, which we successfully implemented. A particle that “exits” one side of the box must return through the other side at the same velocity. The force calculations reflected this; two distances could be used to calculate the force: one that is between two points within the box, and another that takes into account the periodic boundary conditions. The shorter distance of the two was used to calculate each force.

We ran a calculation of 25 particles in a 10×10 box and obtained these results:




Once again, fixing the exponentially increasing forces is our top priority at this point; once we do so, we will be able to accurately calculate speed and velocity distributions of the particles in the box.

The two functions are attached as .pdf files.

Function 1: TheInitiator

Function 2: TheUpdater

A script was used to initialize variables, run the two functions, and plot the results.

EDIT (4/16/15)
The original *.m files have been uploaded and can be found here.


Sushant Mahat and Mohammed Abdelaziz – PHYS 375 Project Proposal

Our project will be aimed at studying the behavior of a sample of N molecules when they are heated by a pulsed laser source. For this project, we will be using Newton’s second law, statistical mechanics, and the information on molecular dynamics simulations in Chapter 9 of Giordano and Nakanishi’s Computational Physics.


In our project we will tackle molecular heating simulations at different levels of difficulty. We will start simply by simulating gas molecules in a box and move on to more complex topics that involve heating and then finally pulsed heating of the gases. From these simulations we will study how long it will take for the gas particles to come to thermal equilibrium and how the speed distributions of the gas particles look when the system is in equilibrium.


After gases, we will try to create similar simulations to study the properties of crystalline solids. If we have time, we will also study the Fermi-Pasta-Ulam problem.
We decided upon this project as this is closely related to what we are studying in our thermal physics class right now and what we studied in classical mechanics in the past. Both branches of science are interesting and very challenging and we hope that undertaking this project will give us a better understanding of the concepts we have learnt so far.