Molecular dynamics in Python¶

In this TP we will be using the Python language to study one of the major theoretical tools of statistical mechanics and materials science: Molecular dynamics simulation.

Molecular dynamics follows the dynamics of "particles" representing atoms, with an empirical interaction between pairs of particles.

We will be doing this with a highly simplified model of "hard disks" in two dimensions, where it is easy to calculate the trajectories of particles and also easy to visualize the results.

The dynamics of the hard disks will be simulated using the Event-driven strategy in order to make the most of the computer numerical precision.

Class based programming¶

Before starting to look and the simulation, we need to learn a little bit about Classes in python. We will ask that solutions to the problems maintain a clean separation between classes/objects. It makes code easier to read, and easier to reuse or modify.

The philosophy is that similar properties are grouped together in a way that is logical for the physical situation that is being studied, see ClassExample/point2d.html and a explanation from Zorana Zeravcic

Simulation with event-driven methods¶

We will work starting with the two files simul.html and animate.html.

simul.html defines a class that creates a system containing 4 particles in two dimensions __init__() sets the 4 particles in a standard position with the code self.position=. We then choose random velocities, and create sets of indexes self.l, self.m in order to be able to manage in interactions between pairs of particles. Particles are moved forwards in the function md_step(). Read simul.html and understand everything. We will be working with this file for several sessions.

animate.html defines a second class that creates a window and displays the system of particles on the screen. The two classes are written to cooperate, and are used together in run.html. Note how the user of the animate.html does not have to know anything about the window management. Everything is "hidden" in the class. The file run.html uses the two classes to run and animate a simulation. You should read animate.html quickly to see what's there, but you will not have to understand it in detail for the first weeks. animate.html will work with any class which defines TWO properties position and md_step(). You can rewrite the physics entirely in simul.html without modifying the display code.

1. Run the code¶

Try to run the program run.html from a terminal with the command python run.html and understand how the pieces work together. For the first day you should only need to modify the code in simul.html. md_step() is called 200 times, as imposed by the line animate.go(nframes=200) in run.html

The two classes can also be used from a Notebook as in Animation_simulation.html, however while developing code this is less useful since error messages can be difficult to understand from within the Jupyter Notebook. We recommend running from run.html. Classes should only be defined in .html files, then imported into the main program.

2. Write the function wall_time()¶

The code as written just moves the particles at a constant speed. The first task is to modify the file simul.html so that particles bounce off the walls of the enclosing box.

The function wall_time() should find the next collision between a particle and a wall, and return the collision time, particle involved and the direction of the wall : x or y.

Manage the collision with the wall in md_step(), and perform the number of required steps so that the total time treated by the function is simul.time.
md_step() should also keep track of the total momentum transferred to the walls during simul_time, and return the value of the Pressure for the time slice. See Pressure.html for more details. Most of the logic can be treated using the functions np.where() and np.unravel_index(), see Vectors.html

This style of programming is not obvious, it may take some hours to get this part of the code correct.

3. Write the function pair_time()¶

When the particles bounce correctly of the sides of the box you should then implement the "hard disk" collisions between pairs of particles. Use the mathematics in Collisions.html to calculate the collision time between pairs of particles, and then update the velocities in md_step(). numpy.where() and numpy.argmin() can be useful.

3. If you wish, you can write an independent driver routine to replace the graphics display¶

Displaying the motion of the particles on the screen slows down the code (a lot!). In "data acquisition" mode is is normal to have a way of running the code that does not display anything, but rather accumulates data that is written to disk for later analysis. You can call your function md_step() from within a loop:

In [1]:
from simul import Simul
import numpy as np

P_list=[]  # a list with the times series of the pressure
simulation = Simul(simul_time=100, sigma=0.5,L=4)  # New Simul object
for i in range(5):
    p=simulation.md_step()  #  read the result of the pressure from a single time step
    P_list.append(p) 
    
    # save P_list to disk using pickle, and use matplotlib plot to look at the result
print('pressure list', P_list)

Simul::md_step
Simul::md_step
Simul::md_step
Simul::md_step
Simul::md_step
pressure list [-1, -1, -1, -1, -1]

At the end of a simulation write the list to a file using pickle.

To measure the time series of the pressure on the walls look at Pressure.html. pickle and other methods to save and read data are described in HowToSaveObjects.html.

4. Modify for arbitrary number of particles¶

Rather than running with just 4 particles, write a general code for $N$ particles which are placed in a box.

You have used a vector notation for all of your code, so little should need modifying.

Your code should run from the command line: python run.html should produce a correct visual output. Any (non-visual) code producing the pressure time series.

Projects¶

When your code is working a number of projects of varying difficulty can be worked on. The code as described above, should however always continue working. Any modifications that you need for your project should be performed by subclassing the simul.html code see Subclassing.html. If you have a project where this does not seem possible talk to the teaching staff before writing independent/duplicate code. Subclassing allows one to take a class and replace/overwrite functions, while keeping the original code in full working order.

You should produce a Jupyter Notebook to describe which project you decided to run and document your experiences and difficulties. Please include any graphical output/curves that is useful. The project can run as either a command line program with python3, or as a Notebook. Notebook programming can be more complicated so that we do not ask for a fully active notebook.

We are looking for concise and clear descriptions, so there is no reason that these notebooks should be overly long and complicated.

List of Projects¶

The projects are presented in a very shortened form, talk to the teaching staff before choosing one, or propose one yourself. The can be approached with different levels of sophistication on both the physics side and the python side. We will help you choose something which corresponds to your knowledge and interests. The projects as written are starting points, they do not require exhaustive exploration of every question. If you find something that you are curious about explore it!

We do not require projects with dynamic curves, real-time display of histograms....

Some of the projects require changing the graphics, Hints are here. For instance different particles sizes/colors/box shapes

  1. Gas pressure Project/Gas.html, simulation compared to van der Waals

  2. Adiabatic expansion of a gas Project/Adiabatic.html.

  3. Osmotic pressure Osmotic

  4. Growing the particle size for dense disordered packings Project/Growing.html.

  5. Thermal conduction Project/Conduction.html.

  6. Velocity distributions, convergence in time, as a function of N, Project/Distributions.html.

  7. Radial distribution functions, Project/Radial.html.

  8. Bunimovich stadium, Project/Stadium.html.

  9. Time reversal, Project/Reverse.html.

  10. Brownian motion Large particle in bath

  11. Depletion interactions Attraction to wall

  12. Two interacting volumes Moving Wall.`

  13. Spatial mixing Mixing

  14. (harder) Implement periodic boundary conditions

  15. (harder) Implement particles with a more general interaction potential

  16. Have particles which form dimers or polymers