MD simulation from scratch

MD simulation from scratch#

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

The main program toy_md.py takes care of the MD propagation and runs different sub-programs for different calculations needed for an MD simulation 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 on the main program toy_md.py, implementing a propagation algorithm, and on toy_md_forces.py, implementing PBCs. Take some time to familarize yourself with the program structure in the ToyMD folder. In that folder, you will also find a carbon-dioxide folder, containing parameters for our 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 this 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 and main parameters/inputs to perform a molecular dynamics simulation from start to finish? Point out the main difference(s) between your scheme and Monte Carlo methods.

In this exercise, you will implement a propagation algorithm and PBCs. Unfortunately, you will only be able to verify if your implementations are correct after modifying 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] 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] and this distance should be [1,1,1]. Furthermore, the distance to the identical point in an adjacent box (e.g [-1.5,-1.5,-1.5]) should be [0, 0, 0].

You can test if your implementation of the distance_PBC function in toy_md_forces.py is running correctly by replacing the expected values for the elementwise distances in the following cell.

import numpy as np
sys.path.append("ToyMD")
import toy_md_forces
from toy_md_forces import distance_pbc
# Your cubic box
box = np.array([2.0, 2.0, 2.0])

# Test 1
dx = distance_pbc(box, np.array([0.5, 0.5, 0.5]), np.array([-0.5, -0.5, -0.5]))
print(f"Test 1: dx = {dx}")
assert np.allclose(dx, ### Expected distance [x,y,z] ###), f"Failed: {dx}"

# Test 2
dx = distance_pbc(box, np.array([0.5, 0.5, 0.5]), np.array([1.5, 1.5, 1.5]))
print(f"Test 2: dx = {dx}")
assert np.allclose(dx, ### Expected distance [x,y,z] ###), f"Failed: {dx}"

# Test 2
dx = distance_pbc(box, np.array([0.5, 0.5, 0.5]), np.array([-1.5, -1.5, -1.5]))
print(f"Test 2: dx = {dx}")
assert np.allclose(dx, ### Expected distance [x,y,z] ###), f"Failed: {dx}"

print("\nAll tests passed!")

To test whether your whole implementation works you can run a short MD run 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

# 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   -589.493 Ekin     73.741 Etot   -515.751 T    9.12 lambda 1.00
Step:     1 Epot   -806.967 Ekin    254.801 Etot   -552.167 T   31.53 lambda 1.00
Step:     2 Epot  -1010.214 Ekin    446.395 Etot   -563.818 T   55.24 lambda 1.00
Step:     3 Epot  -1191.846 Ekin    550.711 Etot   -641.135 T   68.15 lambda 1.00
Step:     4 Epot  -1147.327 Ekin    521.203 Etot   -626.124 T   64.49 lambda 1.00
Step:     5 Epot   -934.489 Ekin    381.774 Etot   -552.715 T   47.24 lambda 1.00
Step:     6 Epot   -835.432 Ekin    207.229 Etot   -628.203 T   25.64 lambda 1.00
Step:     7 Epot   -663.648 Ekin     82.305 Etot   -581.344 T   10.18 lambda 1.00
Step:     8 Epot   -696.184 Ekin     63.622 Etot   -632.562 T    7.87 lambda 1.00
Step:     9 Epot   -782.979 Ekin    159.262 Etot   -623.717 T   19.71 lambda 1.00

Exercise 7 - Bonus

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

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 often 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: rename your traj.pdb file if you want to keep your previous run).

Contents