Molecular dynamics of liquid $\text{H}_2\text{O}$

Radial distribution function

We have prepared a CP2K input file md.in for running a MD simulation of liquid water using the force field from the first exercise (parametrized by Praprotnik et al.).

TASK 1
  • Check that the MD is energy conserving and well-behaved.

Repeat the MD using initial temperatures 200 and 400 K. In order not to overwrite any of your previous files, it is advisable to run the new simulations in different folders.

TASK 2
  • What are the final average temperatures in each of the simulations?
  • Why are they different from the initial ones?
  • The initial atomic configuration stems from an equilibration run. At which temperature was the system (approximately) equilibrated?

Next we are going to analyze the trajectories in order to calculate the radial distribution function (rdf, $g(r)$) as a function of temperature.

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

TASK 3
  • Plot $g_{O-O}(r)$ at 200, 300 and 400 K into the same graph.
  • What are the differences in the height of the first peak?
  • What does this say about the structure of the liquid and is this expected? (2P)
  • Compare to experimental data goo.ALS taken at 300 K.

Infrared spectrum

Due to the partial charges on the oxygen and hydrogen atoms, both the stretching and the bending motion of the $\text{H}_2\text{O}$ molecule give rise to oscillations in its dipole moment. In MD simulations, these frequencies can be extracted from the autocorrelation function of the total dipole moment of all charges in the simulation box.

Repeat the MD at 300 K, but now uncomment the &DIPOLE section in the input file in order to write the total dipole moment to dipole.traj.

We have provided a short Fortran program dipole_correlation.f90 to calculate the autocorrelation function. Use the gfortran compiler to generate the executable and, once the MD simulation is finished, use it to calculate the dipole-autocorrelation function (adjusting its input file dipole.in as needed): 

gfortran dipole_correlation.f90 -o dipole_correlation.x  # compile fortran program
./dipole_correlation.x < dipole.in                       # calculate dipole autocorrelation
TASK 4
  • What are the approximate frequencies of the stretching and bending modes?
  • How do they compare to the normal mode frequencies of the isolated molecule calculated previously? Use a longer simulation time (e.g. 40 ps) to obtain a clearer spectrum.
  • Perform a simulation with a larger time step (e.g. 1.5 fs). What is the effect on the spectrum? Provide graphs of both spectra in the report.
  • Compare with Figure 3 in the paper of Praprotnik et al.