====== Molecular Dynamics simulation of a small molecule====== 

<note warning>
TO USE THE FUNCTION LIBRARY (VERSION UP TO DATE) IN THE INTERACTIVE SHELL:

you@eulerX ~$ module load courses mmm vmd
 
you@eulerX ~$ mmm-init
</note>


<note important> **REMEMBER: this is the command to load the  module for the cp2k program:**

<code bash>
you@eulerX ~$ module load new cp2k
</code>

**and to submit the job:**

<code bash>
you@eulerX ~$ bsub < jobname
</code> 
</note> 




Download the 4.1 exercise into your $HOME folder and unzip it:
<code bash>
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
</code>

<note tip>
All files of this exercise (**all inputs are commented**) can be also downloaded from the wiki: {{exercise_4.1.zip|exercise_4.1.zip}}
</note>

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:


<code - 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
</code>

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

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

<code bash>
you@eulerX exercise_4.1$ bsub cp2k.popt -i inp.nve -o out.nve
</code>
<note tip>
Assignments:
  - We are performing MD at a constant energy, but why we still have to define the temperature?
</note>

  * 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: 
<code bash>
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
</code>
  * Have a look at the corresponding *.ener files (we suggest you to use gnuplot).
<note tip>
Assignments
  - Do you see the energy conservation? Give comments on your observations.
  - Analyse the behavior of potential and kinetic energy, and the temperature.
</note>

<note important>
Hint (plotting with gnuplot).

To plot the Kinetic energy:
<code bash>
gnuplot> plot "nve_md-1.ener" u 1:3 w l t "Kinetic Energy"
</code>
To plot the Potential energy:
<code bash>
gnuplot> plot "nve_md-1.ener" u 1:5 w l t "Potential Energy"
</code>
To plot the Temperature:
<code bash>
gnuplot> plot "nve_md-1.ener" u 1:4 w l t "Temperature"
</code>


</note>

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. 

<note tip>
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 [[doi>10.1063/1.463940|this paper]] by Glenn Martyna (1992).
</note>

In cp2k input files you should again have a look at the following section:
<code - 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
</code>

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: 
<code bash>
you@eulerX exercise_4.1$ bsub cp2k.popt -i inp.100 -o out.100
</code>
  * 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
<code bash>
you@eulerX exercise_4.1$ bsub cp2k.popt -i inp.300 -o out.300
</code>

Now you have the following outputs to study with vmd:
<code>
nve_md-pos-1.pdb
md.100-pos-1.pdb
md.300-pos-1.pdb
</code>

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

<code bash>
vmd nve_md-pos-1.pdb 
</code>

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

<code tcl>
vmd> source "dihedrals.vmd"
</code>

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

<note tip>
Which differences do you notice between the nve, the 100 K and the 300 K case? Can you explain them?
</note>