Contents

5.2. How to run the exercises

Access the NVE version of the code using this link and activate the virtual environment in the terminal again.

source ~/my_venvs/mdmc/bin/activate

Now you can proceed to modify and run the code files.

5.3. Exercises

5.3.1. MD Initialization and Temperature

  1. Implement the initialization described in the Theory section in our ToyMD code in toy_md.py in the section ##INITIALIZATION HERE##. For this you should pick random velocities using np.random.randn() which will produce random numbers according to a gaussian distribution \(\mathcal{N}(0,1)\). \(\sigma\) of \(\mathcal{N}\) needs to be adjusted such that the width of the velocity distribution matches the kinetic energy at the target temperature. Use the variables masses[i] for the mass of particle i and Boltzmanns constant 0.00831415. Remember that for each degree of freedom (e.g velocity in x direction) (5.22) holds.

  2. Implement the Berendsen thermostat in the toy_md.py and toy_md_integrate.py (change the compute_lambda_T function) files.

  3. A better thermostat is the Andersen thermostat. It can be implemented as follows below. Describe what problems this thermostat will present to us e.g for sampling of diffusion coefficients. What advantage does this thermostat have compared to the Berendsen thermostat?

#Andersen Thermostat
for i in range(N_part): # Loop over all particles
    if (r.random()<float(md_params["nu-T"])*float(md_params['time-step'])): # pick random particles according to nu
        sd=np.sqrt(0.00831415*float(md_params['temperature'])/masses[i])  # Velocity distribution at target temperature
        velocities[i]=0+np.random.randn(3)*sd

5.3.2. Force field

Now we can run a MD simulation in our desired NVT ensemble using our code on simple systems. Let’s simulate a box of CO2 molecules.

Investigate the force_field.txt file in the carbon-dioxide folder and then run a short molecular dynamics simulation using the following parameters:

number-of-steps  2000 # Number of integration time steps
time-step        0.001  # Integration time step (picoseconds)
temperature      50  # Simulation temperature
tau-T            0.05 # Temperature coupling time (picoseconds)
output-frequency  10   # Store coordinates every N steps

The simulation is run using

path/to/toy_md.py -c co2.pdb -p params.txt -f force_field.txt -o co2-traj.pdb -w co2-final.pdb

where co2-traj.pdb is a trajectory written each output-frequency steps and co2-final.pdb is the final geometry after the simulation has run the specified number of steps.

Open the resulting trajectory in VMD and plot the bond distance of a C-O bond. What to you observe?

  1. Plot the distribution of the CO and OO bond lengths during the simulation. How do the sampled values correspond to the values set in the force field?

  2. Visualize the radial distribution function of the CO2 trajectory. What do you observe? Bonus: What would you observe for a heterogenous system (i.e a polar molecule solvated in liquid CO2)?

How to open trajectory in VMD and plot the distance

Linux

To open a trajectory or structure in VMD, navigate your terminal to the file of interest (using the cd command) and type the following command:

vmd name_of_file.pdb

You can do this e.g on vdi.epfl.ch

MacOS, Windows

Open VMD and select New Molecule

How to navigate

You can explore the system by holding your left mouse button and moving. You can translate the system by pressing T or selecting Mouse \rightarrow Translate Mode in the VMD Main panel. Rotation and scaling can be performed from this menu as well, or you can just type r and s respectively. To change the centre of the viewport (i.e the position about which geometry operations occur), press c and select a particle to move the centre to.

To view the trajectory as a movie, press the play button from the VMD Main panel, or manually move the counter bar to scroll through the trajectory. Experiment with the buttons on this panel to become more comfortable - you will need this program for the following exercises.

Changing the Drawing Mode

Click on Graphics in the VMD Main panel, and go to Representations.... On the Draw style tab, select ColourID from the Colouring Method drop-down box, and select a colour of your choosing. Next, select CPK from the Drawing Method drop-down box, and change Sphere Scale to 0.5 and Sphere Resolution to 20. Click Apply and close the Graphical Representations window.

Changing the Viewport Background

Click on the Graphics item in the VMD Main panel, and select the Colours... option. Next, click on the Display item within the Categories list, and select Background. Change the background of the viewport to 8 white. Close the Colour Controls window.

Rendering and Saving

To save a snapshot of the current frame, go to File in the VMD Main panel, and click the Render... option. Change the filename to something meaningful, and click Start Rendering. The mac program Preview will immediately open, and from there you can export the snapshot as a .png image.

Plotting the bond distance Press the 2 button on your keyboard and click the atoms of a single CO2 molecule. After clicking two atoms the bond distance will be plotted for the current frame.

Click on the Graphics item in the VMD Main panel, click on Labels and select Bonds in the dropdown instead of Atoms. Click on the first bond shown in the list e.g ( SOL1:C   SOL1:O1) and select the Graph tab. Click on the Graph button and a plot of the distance will be shown.

Visualizing the Radial Distribution Function In the VMD main panel go to Extensions, Analysis, Radial Pair Distribution Function g(r). In Use Molecule select your trajectory file (e.g traj.pdb). Enter as Selection1 and Selection2 name C and leave the checkmarks below as is. Click on Compute g(r).

5.3.3. Timestep and coupling

  1. Test the influence of the coupling parameter \(\tau\) and timestep \(dt\). You set both of those parameters in params.txt Test 3 different settings of \(\tau\) and \(dt\) and describe what you observe.

5.1. Theory 6. Properties from Molecular Dynamics