4.2. MD simulation from scratch#

The ToyMD code provided is a minimalistic python program to run MD simulations. By clicking this link you can open the files for this exercise on noto.

The main program is toy_md.py takes care of the MD propagation and runs different sub-programs for differrent calculations needed for a MD simulations. For example, toy_md_params.py takes care of reading the parameters from the params.txt file, or toy_md_forces.py calculates bonded and non-bonded forces.

For this exercise, we will focus of the main program toy_md.py, implementing a propagation algorithm, and toy_md_forces.py, implementing PBCs. Take some time to familarize with the program structure in the ToyMD folder. In that folder, you will also find a carbon-dioxide folder, containing parameters for a toy system: a box of 216 carbon dioxide molecules. At the end of the exercise, you will be able to run a MD simulation of the system.

Exercise 3

How does an MD program work? Describe schematically (either describe steps in words thoroughly and/or provide an annotated sketch of the steps) what are the main steps to perform a molecular dynamics simulation from start to finish. Point out the main differences between your scheme and Monte Carlo methods.

In this exercise, you will implement a propagation algorithm and PBCs. Unfortunately, you will be able to verify your implementations are correct only when you will have modified both pieces of code needed (Ex4 and Ex6). However, you will find tips and suggestions as comments in the python codes to help you.

import sys
sys.path.append("..")
import helpers
from helpers import *
helpers.set_style()

Exercise 4

Implement the velocity verlet algorithm in the Toy MD program (toy_md.py).

Exercise 5

Why are most MD and MC simulations based on periodic systems? Explain the main purpose of periodic boundary conditions in these schemes.

Exercise 6

Describe a possible implementation of periodic boundary conditions and provide an implementation in toy_md_forces.py to compute the distance of two particles taking into account the perodic boundary conditions in the distance_PBC function (i.e a particle at the edge of the primary cell has a very short distance to the particle in the next cell).

With a box size of [2,2,2] and three points the elementwise distances between point [0.5,0.5,0.5] should be the same to the points [-0.5,-0.5,-0.5] and [1.5,1.5,1.5]. The distance should be [1,1,1].

To test whether your implementation works you can run a short MD run of a box of a toy system consisting in CO2 molecules.

You can run the Toy MD code in two ways:

  1. via teminal:

    • open a new terminal

    • navigate to the ToyMD/carbon-dioxide folder

    • run the toy_md.py script via the following bash command

      python3 ../toy_md.py -c co2.pdb -p params.txt -ff force_field.txt -o traj.pdb -w co2-output.pdb

  2. via jupyter notebook:

    • you can run bash scripts in jupyter notebook cells if you add a ! at the beginning of the command

    • in this case, you will not move to the carbon-dioxide folder, so you will need to correctly specify all the paths relative to the current one, where the notebook is located. You can know what is the current folder by running the next cell, using the pwd (print working directory) bash command

Note: if you run the code without any modification, you will get a ZeroDivisionError: float division by zero error.

!pwd

/home/mdmc/mdmc/Ex4

# If you choose option 2, add a bash command to run toy_md.py in this cell
# start cell with a ! symbol

An example output should look like this.

$ ../toy_md.py -c co2.pdb -p params.txt -f force_field.txt -o traj.pdb -w co2-output.pdb
There are 432 bonds
There are 216 angles
There are 648 conect
There are 648 exclusions
Step:     0 Epot   -588.570 Ekin     72.833 Etot   -515.736 T    9.01 lambda 1.00
Step:     1 Epot   -794.971 Ekin    251.656 Etot   -543.316 T   31.14 lambda 1.00
Step:     2 Epot  -1004.530 Ekin    440.865 Etot   -563.665 T   54.55 lambda 1.00
Step:     3 Epot  -1189.030 Ekin    543.860 Etot   -645.170 T   67.30 lambda 1.00
Step:     4 Epot  -1170.759 Ekin    514.693 Etot   -656.066 T   63.69 lambda 1.00
Step:     5 Epot   -949.895 Ekin    376.968 Etot   -572.926 T   46.65 lambda 1.00
Step:     6 Epot   -820.466 Ekin    204.490 Etot   -615.976 T   25.30 lambda 1.00
Step:     7 Epot   -695.645 Ekin     80.842 Etot   -614.804 T   10.00 lambda 1.00
Step:     8 Epot   -710.875 Ekin     61.903 Etot   -648.973 T    7.66 lambda 1.00
Step:     9 Epot   -795.046 Ekin    155.754 Etot   -639.292 T   19.27 lambda 1.00

Exercise 7 - Bonus

What ensemble does the code in the ToyMD program sample in? Which of the quantities linear momentum \(p\) and angular momentum \(L\) are conserved when running with and without periodic boundary conditions?

4.3. Visualize trajectory#

Once your code is ready, running it will produce some output files. In particular, traj.pdb will contain the trajectory of your toy system.

Note that the params.txt file contains the parameters for the MD run. Among them, number-of-steps will indicate the number of integration time steps to be performed and output-frequency will indicate how ofter to store coordinates in the output file. Of course, what you will see in the viewer depends on how long the MD simulation is and how often you saved the coordinates.

Using the cells below you can visualize the trajectory generated. You first need to finish and run the code before you can look at the trajectory.

import MDAnalysis

u = MDAnalysis.Universe("./ToyMD/carbon-dioxide/traj.pdb")

show_trajectory(u)

Exercise 8

What do you observe during the trajectory? Try to relate what you see with the parameters set for the simulation (i.e. simulation length, time step, temperature). Explain the changes, if any, you observe in the simulation if you change one or two of those parameters (note: you will need to delete your traj.pdb