Molecular Dynamics simulation of a small molecule

TO USE THE FUNCTION LIBRARY (VERSION UP TO DATE) IN THE INTERACTIVE SHELL:

you@eulerX ~$ module load courses mmm vmd

you@eulerX ~$ mmm-init

REMEMBER: this is the command to load the module for the cp2k program:

you@eulerX ~$ module load new cp2k

and to submit the job:

you@eulerX ~$ bsub < jobname

Download the 4.1 exercise into your $HOME folder and unzip it:

you@eulerX ~$ wget http://www.cp2k.org/_media/exercises:2015_ethz_mmm:exercise_4.1.zip
you@eulerX ~$ unzip exercises:2015_ethz_mmm:exercise_4.1.zip
you@eulerX ~$ cd exercise_4.1

All files of this exercise (all inputs are commented) can be also downloaded from the wiki: exercise_4.1.zip

You will start from a configuration already computed in the second lecture (inp.a.pdb) which is included in the repository of this exercise as well. Update the following part of the file inp.nve for the first simulation:

md_part.inp.nve

&MD                                           ! This section defines the whole set of parameters needed perform an MD run.
  ???????? ???                                ! Please specify the appropriete ensemble for you MD simulation
  ????? ??????                                ! Please specify the number of MD steps to perform
  ???????? ???? ???                           ! Please specify the length of an integration step
  ??????????? ?????                           ! Please specify the initial temperature
&END MD

To get more information, please visit cp2k reference manual, section Molecular Dynamics: http://manual.cp2k.org/trunk/CP2K_INPUT/MOTION/MD.html

Perform a constant energy simulation, 100000 time steps, with a time step of 1 fs. Use 100 K as an initial temperature!

you@eulerX exercise_4.1$ bsub cp2k.popt -i inp.nve -o out.nve

Assignments:

We are performing MD at a constant energy, but why we still have to define the temperature?

Make four copies of the previous input file (say inp.nve_0.1, inp.nve_2.0, inp.nve_3.0, inp.nve_4.0), in each input file modify the time step (use 0.1, 2, 3, 4 fs respectively) and the name of the project.
Perform the simulations with all these input files:

you@eulerX exercise_4.1$ bsub cp2k.popt -i inp.nve_0.1 -o out.nve_0.1
you@eulerX exercise_4.1$ bsub cp2k.popt -i inp.nve_2.0 -o out.nve_2.0
you@eulerX exercise_4.1$ bsub cp2k.popt -i inp.nve_3.0 -o out.nve_3.0
you@eulerX exercise_4.1$ bsub cp2k.popt -i inp.nve_4.0 -o out.nve_4.0

Have a look at the corresponding *.ener files (we suggest you to use gnuplot).

Assignments

Do you see the energy conservation? Give comments on your observations.
Analyse the behavior of potential and kinetic energy, and the temperature.

Hint (plotting with gnuplot).

To plot the Kinetic energy:

gnuplot> plot "nve_md-1.ener" u 1:3 w l t "Kinetic Energy"

To plot the Potential energy:

gnuplot> plot "nve_md-1.ener" u 1:5 w l t "Potential Energy"

To plot the Temperature:

gnuplot> plot "nve_md-1.ener" u 1:4 w l t "Temperature"

Now you will perform a constant Temperature simulations, where the system is in contact with a thermostat, and the conserved quantity includes the thermostat degrees of freedom.

Concerning temperature control, in these exercises we will use the NOSE-HOOVER chains method. This has been briefly described in the lecture, and is presented in this paper by Glenn Martyna (1992).

In cp2k input files you should again have a look at the following section:

md_part.inp.300

  &MD                                           ! This section defines the whole set of parameters needed perform an MD run.
    ???????? ???                                ! Please specify the appropriete ensemble for you MD simulation
    ????? ??????                                ! Please specify the number of MD steps to perform
    ???????? ???                                ! Please specify the length of an integration step
    ??????????? ???                             ! Please specify the temperature of the simulation
    &??????????                                 ! Please specify a thermostat section here
      &????                                     ! Please put here a section which specfies Nose-Hoover thermostat chain
        TIMECON 50                              ! Timeconstant of the thermostat chain
        LENGTH 3                                ! Length of the Nose-Hoover chain 
        YOSHIDA 3                               ! Order of the yoshida integretor used for the thermostat
      &??? 
    &???
  &END MD

Edit the inp.100 file (Put there: NVT ensemble, 100000 steps of simulation, 100 K, Nose-Hoover thermostat and 1.0 fs of timestep). The first simulation is done at 100 K:

you@eulerX exercise_4.1$ bsub cp2k.popt -i inp.100 -o out.100

Then, perform another simulation, using the restart file from the previous simulation: inp.300. But first you have to edit it exactly like in the previous case, but put 300 K instead of 100 K

you@eulerX exercise_4.1$ bsub cp2k.popt -i inp.300 -o out.300

Now you have the following outputs to study with vmd:

nve_md-pos-1.pdb
md.100-pos-1.pdb
md.300-pos-1.pdb

Open (for example) nve_md-pos-1.pdb with VMD:

vmd nve_md-pos-1.pdb

Open Tk Console (Extensions menu > Tk console). And to define the two dihedrals (PHI and PSI) from there, you can enter:

vmd> source "dihedrals.vmd"

You can also pick from the extensions the “RMSD trajectory tool” and use it to align the molecule along the trajectory (Extensions>Analysis>RMSC Trajectory Tool). Replace the word “protein” with “all” in the selection, and then use “align”. You will see that now the molecule is well aligned along the path.

Now using “Labels” menu, plot the graph of two dihedral angles:

Go to Graphics > Labels
In the drop-down list chose Dihedrals
Chose both dihedrals in the list
Go to the “Graph” section
Press on the “Graph…” button
(Optional) save these graps in a text file (File > Export to ASCII matrix…)

Which differences do you notice between the nve, the 100 K and the 300 K case? Can you explain them?