This is an old revision of the document!

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 ; 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

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).

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 (input and scripts are all 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. Use the file inp.nve for the first simulation, which is a constant energy simulation.

md_part.inp.nve

&MD                                           ! This section defines the whole set of parameters needed perform an MD run.
  ENSEMBLE NVE                                ! The ensemble/integrator that you want to use for MD propagation
  STEPS 100000                                ! The number of MD steps to perform
  TIMESTEP [fs] 1.0                           ! The length of an integration step
  TEMPERATURE 100.0                           ! The temperature in K used to initialize the velocities with init and pos restart velocities
&END MD

Perform a constant energy simulation, 100000 time steps, with a time step of 1 fs.

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

Using a different input file, modify the time step and the name of the project. Do it for 0.1, 2, 3, 4 fs.
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.
Describe 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. In cp2k input file it can be specified as follows:

md_part.inp.300

  &MD                                           ! This section defines the whole set of parameters needed perform an MD run.
    ENSEMBLE NVT                                ! The ensemble/integrator that you want to use for MD propagation
    STEPS 100000                                ! The number of MD steps to perform
    TIMESTEP 1.0                                ! The length of an integration step
    TEMPERATURE 300.0                           ! The temperature in K used to initialize the velocities with init and pos restart velocities
    &THERMOSTAT                                 ! This section specifies thermostat type and parameters controlling the thermostat
      &NOSE                                     ! This section specifies paramameters of the Nose Hoover thermostat chain
        TIMECON 50                              ! Timeconstant of the thermostat chain
      &END NOSE
    &END
  &END MD

The first simulation is done at 100 K:

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

Then, perform a simulation at 300 K, using the restart file from the previous simulation: inp.300.

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:

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.

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