Differences

This shows you the differences between two versions of the page.

--- exercises:2021_uzh_acpc2:ex02 [2021/04/13 11:02] – external edit 127.0.0.1
+++ exercises:2021_uzh_acpc2:ex02 [2021/05/19 14:30] (current) – jglan
@@ Line 2: / Line 2: @@
 ===== Water =====
-[[http://www1.lsbu.ac.uk/water/water_models.html|Water molecular models]] are computational techniques that have been developed in order to help discover the structure of water. In this section, you will be asked to calculate some physical properties based on classical molecular dynamics simulation. The TIP3/Fw model will be usded in the simulations.
+[[http://www.idc-online.com/technical_references/pdfs/chemical_engineering/Water_models.pdf|Water molecular models]] are computational techniques that have been developed in order to help discover the structure of water. In this section, you will be asked to calculate some physical properties based on classical molecular dynamics simulation. The TIP3/Fw model will be used in the simulations.
+We have prepared a CP2K input file ''water.inp'' for running a MD simulation of liquid water using the force field from the first exercise (parameterized by [[https://aip.scitation.org/doi/pdf/10.1063/1.1884609|Praprotnik et al.]]).
+<note important>Download {{ :exercises:2021_uzh_acpc2:water.zip |}} and extract it.</note>
+<code>wget https://www.cp2k.org/_media/exercises:2021_uzh_acpc2:water.zip
+unzip water.zip
+</code>
-We have prepared a CP2K input file ''water.inp'' for running a MD simulation of liquid water using the force field from the first exercise (parametrized by [[https://aip.scitation.org/doi/pdf/10.1063/1.1884609|Praprotnik et al.]]).
-<note important>Download {{ :exercises:2019_uzh_acpc2:water.zip |water.zip}} and extract it.</note>
 <note>**TASK 1**
-  * Check that the MD is energy conserving and //well-behaved//.
+  * Check that the MD simulation is energy conserving and //well-behaved//.
   * What are the final average temperatures of the simulation?
   * The initial atomic configuration stems from an equilibration run. At which temperature was the system (approximately) equilibrated?
 </note>
-Next we are going to analyze the trajectories in order to calculate the [[http://en.wikipedia.org/wiki/Radial_distribution_function|radial distribution function]] (rdf,  $g(r)$ ) as a function of distance  $r$ .
+Next, we are going to analyze the trajectories in order to calculate the [[http://en.wikipedia.org/wiki/Radial_distribution_function|radial distribution function]] (rdf,  $g(r)$ ) as a function of the distance  $r$ .
-VMD comes with an extension for exactly this purpose: In the VMD Main window open "Extensions -> Analysis" click on "Radial Pair Distribution function  $g(r)$ ". In the appearing window use "Utilities -> Set unit cell dimensions" to let VMD know the simulation box you used. After that use Selection 1 and 2 to define the atomic types that you want to calculate the rdf for, for example "element H".
+VMD comes with an extension for exactly this purpose: In the VMD Main window open "Extensions -> Analysis" click on "Radial Pair Distribution function  $g(r)$ ". In the appearing window use "Utilities -> Set unit cell dimensions" to let VMD know the size of the periodic simulation box which you used. After that use Selection 1 and 2 to define the atomic types that you want to calculate the rdf for, for example "element H".
 <note>**TASK 2**
-  * Plot  $g_{O-O}(r)$  at 300K and experimental value ''goo.ALS'' taken at 300 K into same graph.
+  * Plot  $g_{O-O}(r)$  at 300K and the experimental value provided in ''goo.ALS'' taken at 300 K into the same graph.
 </note>
-Then we will calculate diffusion coefficient. It is a proportionality constant between the molar flux due to molecular diffusion and the gradient in the concentration of the species (or the driving force for diffusion), which is defined by:
+Now, we will calculate the diffusion coefficient which is a proportionality constant between the molar flux due to molecular diffusion and the gradient in the concentration of the species (or the driving force for diffusion) and is defined by:
 ** $6D=\lim_{t\to\infty}  \ \frac{\delta <r^2(t)>}{\delta t}$ **
-To evaluate this expression, all that is needed is to evaluate at each point in time in the calculation the average of the square of the distance that each atom has traveled since the start of the production phase of the dynamics, and examining the slope of this function at a long time. By storing the initial coordinates, it is straightforward to evaluate the square of the distance. However, some care is needed due to the use of periodic boundary conditions: the program stores x, the coordinates, but in many programs, during the dynamics, if any atom has its x, y, or z coordinate become larger than box size or smaller than zero, then it is moved back into the other side of the box. This has the effect of making the raw distance traveled meaningless. The value of D is obtained from the slope, at a long time, of the right-hand side of the above equation (you need to divide by six to obtain D, take the slope, and also be careful with units).
+To evaluate this expression, all that is needed is to evaluate the average of the square of the distance that each atom has traveled since the start of the production phase of the dynamics for each point in time during the simulation, and examining the slope of this function in the long time limit. By storing the initial coordinates, it is straightforward to evaluate the square of the traveled distance. However, some care is needed due to the use of periodic boundary conditions: the program stores x, the coordinates, but in many programs, during the dynamics, if any atom has its x, y, or z coordinate become larger than the box size or smaller than zero, it is moved back to the other side of the box. This has the effect of making the raw distance traveled meaningless.
+Finally, if you take care of the above, the value of D is obtained from the slope, at a long time, of the right-hand side of the above equation (also be careful with the units).
-VMD comes with an extension for exactly this purpose: In the VMD Main window open “Extensions → Analysis” click on “RMSD Trajectory Tool”. In the appearing window use “all” to let VMD know the molecule you want to track. Tick "Plot", and press "RMSD", you will have the RMSD plot for the water system.
+Once again, VMD comes with an extension for exactly this purpose: In the VMD
+Main window open “Extensions → Analysis” click on “RMSD Trajectory Tool”. In the
+appearing window change “protein” to “all” to let VMD know the molecule you want
+to track. Press "RMSD" to run the analysis (Note: this might take a while!).
+Finally, use "File -> Plot data" to plot the RMSD for the water system.
 <note>**TASK 3**
-  * Plot RMSD and MSD for the water at 300K and calculate corresponding diffusion coefficient from the slope of MSD, are they expected?
+  * Plot the RMSD and MSD for the water at 300K and calculate the corresponding diffusion coefficient from the slope of the MSD. Does the value match your expectation?
 </note>
-<note important>The diffusion coefficient is calculated using MSD but NOT RMSD.</note>
+<note important>The diffusion coefficient is calculated using the MSD, **NOT** the RMSD.</note>
+<note tip>RMSD stands for [[https://en.wikipedia.org/wiki/Root-mean-square_deviation|Root-mean-square-deviation]]. Accordingly, MSD simply does NOT take the //square root//.</note>
-We will compute the vibrational spectrum, and dielectric constant of water based on molecular dynamics. The spectra for water are available in this paper [[https://aip.scitation.org/doi/pdf/10.1063/1.1884609|https://doi.org/10.1063/1.1884609]]. The provided program computes the correlation function of the (derivative of) the dipole moment and performs the Fourier transform.
+Now, you will compute the vibrational spectrum and dielectric constant of water based on the previous MD simulation. Reference spectra for water are available from [[https://aip.scitation.org/doi/pdf/10.1063/1.1884609|Praprotnik et al.]]. We provide you with a small FORTRAN program which computes the correlation function of the (derivative of) the dipole moment and performs the Fourier transform:
 \begin{equation}
-A(\omega)\propto{\int\langle{\dot{\mu}}({\tau}){\dot{\mu}}(t+{\tau})\rangle_{\tau}e^{-i{\omega}t}d{t}},
+A(\omega)\propto{\int\langle{\dot{\mu}}({\tau}){\dot{\mu}}(t+{\tau})\rangle_{\tau}e^{-i{\omega}t}d{t}} \ .
 \label{eq:auto}
 \end{equation}
@@ Line 55: / Line 65: @@
 \end{equation}
-Compile the FORTRAN code, and execute the program
+To perform this calculation, compile the FORTRAN code, and execute the program:
 <code>
 gfortran cpt_ir_diele.f90 -o cpt_ir_diele.o
@@ Line 62: / Line 72: @@
 <note>**TASK 4**
-  * Compute the IR spectrum and plot it, match the frequencies with vibrational mode.
+  * Compute the IR spectrum and plot it. Match the frequencies to their vibrational modes.
   * Compute the dielectric constant of water at 300K.
-  * Does IR or dielectric constant match the experimenal value? If not, why?
+  * Does the IR or dielectric constant match the experimental value? If not, why?
 </note>
@@ Line 72: / Line 82: @@
 In particular, it has more than one long-lived conformation, which we will identify in this exercise by mapping out its //potential energy surface//.
-The conformations of glyala dipeptide are characterized by the dihedral angles of the backbone.
+The conformations of the glyala dipeptide are characterized by the dihedral angles of the backbone.
-Below, we color carbons in green, hydrogens in white, oxygen in red and nitrogen in blue, i.e.
+Below, we color carbons in green, hydrogens in white, oxygen in red and nitrogen in blue, showing that the torsional angle  $\phi$  is N-C-C-N , while  $\psi$  is C-N-C-C along the backbone.
-the torsional angle  $\phi$  is N-C-C-N , while  $\psi$  is C-N-C-C along the backbone.
 <note important>Please download {{ exercises:2017_uzh_acpc2:glyala-epot.tar.gz |glyala-epot.tar.gz}} and extract this file using <code> tar -xvf glyala-epot.tar.gz </code></note>
@@ Line 81: / Line 90: @@
 {{  :exercises:2017_uzh_acpc2:vmdscene.png ?direct&400 |}}
-<note>**TASK 1**
+<note>**TASK 5**
 Visualize the structure ''glyala.pdb'' with VMD and determine the atomic indices of the atoms defining the dihedral angles.
 </note>
-<note important>//Note:// While VMD starts counting atoms from 0, CP2K starts counting from 1, i.e. the VMD indices need to be increased by 1.</note>
+<note important>//Note:// While VMD starts counting atoms from 0, CP2K starts counting from 1. Thus, the VMD indices need to be increased by 1 when inserted into your CP2K input file.</note>
-With this knowledge at hand,
+With this knowledge at hand, we will fix the dihedral angles and perform geometry optimization for all remaining degrees of freedom.
-we will fix the dihedral angles and perform geometry optimization for all remaining degrees of freedom.
-<note>**TASK 2**
+<note>**TASK 6**
-  - The atomic indices defining the dihedral indices in the input file ''geo.in'' are missing. Replace ''I1'' to ''I4'' by the atomic indices determined previously.
+  - The atomic indices defining the dihedral indices in the input file ''geo.in'' are missing. Replace ''I1'' to ''I4'' by the atomic indices determined in //Task 5//.
-  - Use ''perform-gopt.sh'' to perform the grid of geometry optimizations.
+  - Use the provided bash script, ''perform-gopt.sh'', to perform the grid of geometry optimizations.
-  - Use gnuplot to plot the potential energy surface (we have provided a script ''epot.gp''). Which are the two most favoured conformations? <code> $ gnuplot</code><code> gnuplot > load "epot.gp"</code>
+  - Use gnuplot to plot the potential energy surface (we have provided a script ''epot.gp'') <code> gnuplot > load "epot.gp"</code>What are the two most favored conformations?
 </note>
 ===== Glyala in water =====
-Now, we will move to a more realistic system - Glyala in water. We will preformed a MD of glyala in water and save the trajectory.
-The initial geometry provided in the PDB file is a glyala molecule solvated by 73 water molecules. The geometry is not equilibrated. You need first to equilibrate the system at 300K. When the system is equilibrated, you need to analysis the result.
+Now, we will move to a more realistic system - Glyala in water. We will preform a MD simulation of glyala in water and save the trajectory.
-<note important>Download the {{ exercises:2017_uzh_acpc2:glyala_water.tar.gz |glyala_water.tar.gz}} and extract it. </note>
+The initial geometry provided in the PDB file is a glyala molecule solvated by 73 water molecules. The geometry is not equilibrated. Thus, you first need to equilibrate the system at 300K. When the system is equilibrated, you need to analyze the result.
-<note>**TASK 6**
+<note important>Download the {{ exercises:2017_uzh_acpc2:glyala_water.tar.gz |glyala_water.tar.gz}} and extract it. See above on how to extract a ''tar''-archive.</note>
-   - Perform the molecular dynamics simulation using NVT ensemble at 300K. Change TIMECON (i.e.500, 2000 fs) in the &THERMOSTAT section.
+<note>**TASK 7**
+   - Perform the MD simulation using an NVT ensemble at 300K. Change TIMECON (i.e.500, 2000 fs) in the &THERMOSTAT section.
    - Determine from which step the system is equilibrated, plot the calculated properties and explain why.
    - Compute the O-O radial distribution function for water with acceptable statistics using 20 ps (after equilibration) of simulated time.
-   - Determine the solvation shell by calculating RDF of g $_{CO}$  (carbon atoms from glyala and oxygen atoms from water)
+   - Determine the solvation shell by calculating the RDF of g $_{CO}$  (carbon atoms from glyala and oxygen atoms from water)
 </note>
 <note tip>**Tip for O-O RDF for water**
-In last exercise, one already knew how to calculate the RDF for the Argon system. In TASK3, you need to calculate the RDF only for water instead of whole system. Since the glyala contain two oxygen atoms, it is not reasonable to include the oxygen atoms in glyala molecule if we are only interested in O-O RDF for water.
+From the last exercise, you already know how to calculate the RDF for the Argon system. However, in TASK 7 you need to calculate the RDF only for water instead of the whole system. Since the glyala molecule contains two oxygen atoms itself, it is not reasonable to include these oxygen atoms of glyala if we are only interested in the O-O RDF for water.
-Using VMD, the O-O RDF for the water can be easily calculated. In the <code> Selection 1, Selection 2 </code>, one need to specify <code> element O and not same residue as element C </code> The frames should start from the beginning of production run.
+However, using VMD, the O-O RDF for the water can still be easily calculated. In the <code> Selection 1, Selection 2 </code>, you need to specify <code> element O and not same residue as element C </code> in order to exclude the oxygen atoms present in glyala. The frames should start from the beginning of production run.
 </note>