Molecular Dynamics project
Date: Spring semester 2017
In this project we will implement a model called molecular dynamics (MD). Molecular dynamics allows us to study the dynamics of atoms and explore the phase space. For those of you who are taking Thermal and Statistical Physics PHY410, this project could be an interesting choice, as well as the one on the Ising model. The atoms interact through a force given by the negative gradient of a potential energy function. With this force, we can integrate Newton's laws. While exploring the phase space, we will sample statistical properties of the system like energy, temperature, pressure, diffusion constant and so on.
We will also learn how to do write object oriented code. You are given a code skeleton that contains the basic structure of an MD code, but you will have to implement many of the functions yourself. We think that learn by example is a great way to learn how to proper object orient your code, but you can of course start from scratch if you want. You can find the code here. You can clone or fork this repository into your own. The class structure which you developed for project 3 can also be used again here, as well as the Velocity Verlet algorithm of project 3.
After you have downloaded the repository, open the project in Qt Creator and see if it runs. The program should create 100 argon atoms and place them uniformly inside a cubic box with sides 10 Angstroms. Each atom is given a velocity according to the Maxwell-Boltzmann distribution so that
where $m$ is the mass of the atom, $k_B$ is Boltzmann's constant and $T$ is the temperature. We recognize this as a normal distribution with zero mean and standard deviation $\sigma = \sqrt{k_B T/m}$. We recognize this also as the Boltzmann distribution $P(E) \propto \exp{(-\beta E)}$, with $E = 0.5mv^2$. The program will evolve the system in time with no forces so that all atoms move in a straight line. It will create a file called movie.xyz containing all the timesteps ready to visualize with for example Ovito. You can download Ovito from this site. You should start with the fun part - to look at the simulation in Ovito. We then strongly recommend that, before you start, you look at the code structure and try to understand how the different classes are connected to each other and how a typical timestep is computed. You can skip understanding the contents of unitconverter.cpp and vec3.cpp.
Before you start, we would like that you pay attention to the following two steps:
First, spend some 30 minutes figuring out what molecular dynamics is and which problems it can be applied to. What areas of physics can it be used in? What about chemistry and biology? What you find here should be part of the introductory section in your report.
Then run the program, visualize the output with for example Ovito (as mentioned above). Now, spend some 30-60 minutes looking at the code and understand how the output is produced and try to understand every step in the code. The best way to do this is to start at the top of the main() function and follow every line and every function call, step by step. As a part of your report, you should explain the code structure and draw a flow chart (a diagram showing how the program is executed). This part is extremely important and will save you a lot of time. And of course, please ask us teachers if you have any questions!
When working with MD, the system size is usually limited by available computer resources. A typical large MD simulation (with a parallelized code) contains a few million atoms corresponding to a system much smaller than a cubic micron. In order to get rid of boundary effects, we usually apply periodic boundary conditions so that we simulate a system of infinite size. You need to implement the applyPeriodicBoundaryConditions function in the System class. After doing so, run the program again and notice how the atoms now are contained inside a box. Discuss the benefits of this strategy.
The atoms are usually given velocities according to the Maxwell-Boltzmann distribution (you should discuss why this is the case) which will result in a nonzero net momentum in the system. Implement the removeMomentum function in the System class so that the net momentum is zero.
The atoms are now uniformly distributed in space. This is not very physical, so we should place the atoms in a crystal structure. When we later implement the Lennard-Jones potential, we will see that the face-centered cubic (FCC) lattice is a stable structure for the potential (it is actually the crystalline structure of argon). A lattice is built up by unit cells - a group of atoms - so that a larger system can be created by repeating these cells in space. An FCC lattice unit cell of size $b$ in Angstroms (this is the lattice constant) consists of four atoms with local coordinates
You can now create $N_x \times N_y \times N_z$ such unit cells next to each other to form a larger system so that the origin of unit cell $(i,j,k)$ is
where $i=0,1,..., N_x-1, j=0,1,..., N_y-1, k=0,1,..., N_z-1$. The unit vectors of the unit cells are scaled with the lattice constant $b$ so that
Implement the function createFCCLattice(int numberOfUnitCellsEachDimension, double latticeConstant) in the System class. Remember to remove the 100 atoms that are created as the code is right now. Use lattice constant $b=5.26$ $\mathrm{Angstrom}$. Remember to update the system size so the applyPeriodicBoundaryConditions function works properly. If you now visualize the result of the program, you should see the nice crystalline structure in the beginning of the simulation. What is the density $\rho$ in your system?
While this looks nice, we need to implement the computation of forces in order to see interesting physics. In this project, we will use the Lennard-Jones potential which calculates the energy between two atoms $i$ and $j$ as
where $r_{ij} = \vert\mathbf{r}_i - \mathbf{r}_j\vert$ is the distance from atom $i$ to atom $j$, $\epsilon$ is the depth of the potential well (with dimension energy) and $\sigma$ is the distance at which the potential is zero. For argon, optimal values of the parameters are:
In our code, we use the so called MD units (see the units definition at the end). This potential reproduces thermodynamic equilibrium properties that are in good agreement with experimental values for argon. The equation of state for this system is the van der Waals equation of state. Phases such as solid, liquid and gas are reproduced from this simple model with phase transitions and other properties as well with a remarkable agreement with data.
The total potential energy $V$ in the system is computed by summing over all pairs of atoms (counting each pair only once)
The force is calculated by taking the negative gradient of the potential
giving the $x$-component (the other components are calculated the same way)
where $x_{ij}$ is the $x$-component of $\mathbf{r}_{ij}$.
Calculate (analytically) the force and implement the calculateForces function in the LennardJones class, using the Velocity Verlet algorithm that you developed for project 3. You will have to take care of the periodic boundary conditions. You could use the minimum image convention as described here. Now run the simulation with 5 unit cells in each dimension (500 atoms) for different initial temperatures. Can you tell if the system is in a solid state just by looking at the system in Ovito? At what initial temperature is the crystal melting?
We now have a working Molecular dynamics code! The next step is to calculate some physical properties of the system. The easiest properties to measure are the kinetic and potential energy (calculated in the LennardJones class). The kinetic energy is defined as
where $m_i$ and $v_i$ is the mass and the speed of atom $i$. You can use the function atom->velocity.lengthSquared() to calculate the dot product of the velocities.
You can also calculate an estimate of the temperature through the equipartition theorem, see for example D. Schroeder's An Introduction to Thermal Physics for details,
We can use this to define an instantaneous temperature
All of these quantities should be calculated in the StatisticsSampler class. All the statistical quantities should be saved to a file which needs to be implemented in the saveToFile function in the same class.
As you will see, since the temperature is proportional to the kinetic energy, its initial value will drop in the beginning when the system is initialized in the crystal structure. Why does this happen? What happens to the total energy? In order to have a system with final temperature $T$, what initial temperature $T_i$ is needed? What is the ratio $T/T_i$? Do multiple simulations with different random seeds and plot this ratio as a function of time. Again visualize with for example Ovito and try to determine at what temperature the system actually melts at. After the temperature drops, it will reach a more or less stable value. The system is then said to be in an equilibrium state.
The last exercise is to compute the diffusion constant and use this to measure the melting temperature. We use the Einstein relation that relates the so called mean square displacement (MSD) to the diffusion constant
where $D$ is the diffusion constant, $t$ is time. The mean square displacement for atom $i$ is calculated as
where $\mathbf{r}_i(t)$ is the position of atom $i$ at time $t$. When you implement this you need to add the initial position as a new property in the Atom class. You also need to update the applyPeriodicBoundaryConditions so that $r_i^2(t)$ is the actual displacement as if no periodic boundary condition has occured. Why is this important?
We can use the diffusion constant to find the melting temperature. Atoms in a solid will not diffuse much, meaning that a solid will have a diffusion constant close to zero. Now solve the Einstein relation for the diffusion constant $D$ and add a function in StatisticsSampler thats measures this quantity. Plot the diffusion constant as a function of temperature $T$ and use this to find an estimate for the melting temperature. Use what you found earlier to make sure you simulate temperatures both above and below the melting temperature. Also remember not to use the initial temperature, but the measured temperature after the system has reached equilibrium. You can use the builtin UnitConverter to convert the temperature to SI units as illustrated in the top of the main function.
The idea of an MD simulation is to sample microstates from a statistical ensemble. In our simulation, we have a constant number of atoms $N$, a constant volume $V$ and a (more or less, depending on our integrator) constant energy $E$. This corresponds to the microcanonical ensemble (NVE).
Although we produce the dynamics of atoms, we should not see MD as a model that will give the true trajectories of the atoms, but rather a model that allows us to explore the phase space according to the probabilities we can calculate with statistical mechanics. A fundamental assumption in any MD simulation is the ergodic hypothesis. It states that over long periods of time, the time spent by a system in some region of the phase space of microstates with the same energy is proportional to the volume of this region, i.e., that all accessible microstates are equiprobable over a long period of time. This means that by evolving the atoms through time will visit microstates of the ensemble according to the true probabilities given by the ensemble.
The program uses by default a set of units. We have chosen the following four units
Boltzmann's constant is then equal to one (convince yourself about this), and other units can be derived by using known relations. We can for example find the unit of time by using $E=mc^2$. We have that
which can be solved for time
By inserting the values above, we get
This way, we can continue working out the values of the other units. What are the units of force? And pressure? In the project, a class that calculates these values is given.
Here follows a brief recipe and recommendation on how to write a report for each project.
Give a short description of the nature of the problem and the eventual numerical methods you have used.
Describe the algorithm you have used and/or developed. Here you may find it convenient to use pseudocoding. In many cases you can describe the algorithm in the program itself.
Include the source code of your program. Comment your program properly.
If possible, try to find analytic solutions, or known limits in order to test your program when developing the code.
Include your results either in figure form or in a table. Remember to label your results. All tables and figures should have relevant captions and labels on the axes.
Try to evaluate the reliabilty and numerical stability/precision of your results. If possible, include a qualitative and/or quantitative discussion of the numerical stability, eventual loss of precision etc.
Try to give an interpretation of you results in your answers to the problems.
Critique: if possible include your comments and reflections about the exercise, whether you felt you learnt something, ideas for improvements and other thoughts you've made when solving the exercise. We wish to keep this course at the interactive level and your comments can help us improve it.
Try to establish a practice where you log your work at the computerlab. You may find such a logbook very handy at later stages in your work, especially when you don't properly remember what a previous test version of your program did. Here you could also record the time spent on solving the exercise, various algorithms you may have tested or other topics which you feel worthy of mentioning.
The preferred format for the report is a PDF file. You can also use DOC or postscript formats or as an ipython notebook file. As programming language we prefer that you choose between C/C++, Fortran2008 or Python. The following prescription should be followed when preparing the report:
Use your github repository to upload your report. Indicate where the report is by creating for example a Report folder. Please send us as soon as possible your github username.
Place your programs in a folder called for example Programs or src, in order to indicate where your programs are. You can use a README file to tell us how your github folders are organized.
In your git repository, please include a folder which contains selected results. These can be in the form of output from your code for a selected set of runs and input parameters.
In this and all later projects, you should include tests (for example unit tests) of your code(s).
Comments from us on your projects, with score and detailed feedback will be emailed to you.
Finally, we encourage you to work two and two together. Optimal working groups consist of 2-3 students. You can then hand in a common report.