CP2K Open Source Molecular Dynamics

User Tools

Site Tools

Trace:
exercises:2015_ethz_mmm:monte_carlo_ice

Differences

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

Link to this comparison view

Next revision		Previous revision
Next revisionBoth sides next revision
exercises:2015_ethz_mmm:monte_carlo_ice [2015/02/06 17:49] – external edit 127.0.0.1exercises:2015_ethz_mmm:monte_carlo_ice [2015/03/18 15:18] – [Main Input-File (run this)] sclelia
Line 1: Line 1:
 ====== Properties of Ice from Monte Carlo Simulations ====== ====== Properties of Ice from Monte Carlo Simulations ======
-<note important> 
-  - Add the line ''PRINT_COORDS .FALSE.'' in the ''TMC'' section. It disables the output of the coordinate trajectory and thus avoids the files-size problem. 
-  - Add the line ''RND_DETERMINISTIC 42'' to the ''TMC'' section to choose a random number seed. You have to replace ''42'' with different values to obtain different samplings. 
-</note> 
  
-In this exercise we will use Monte Carlo sampling to calculate the [[wp>Relative_permittivity|dielectric constant]] and the [[wp>Thermal_expansion#Coefficient_of_thermal_expansion| thermal expansion coefficient]] of the disordered hexagonal phase (Ih) of water ice. In order to speed up the calculation we are going to use the simple SPC/E water model. For more information see: [[doi>10.1063/1.1568337]]. A DFT-version of the calculation can be found here: [[doi>10.1021/jp4103355]]. +In this exercise we will use Monte Carlo sampling to calculate the [[wp>Relative_permittivity|dielectric constant]] of the disordered hexagonal phase (Ih) of water ice. \\ 
- +
- +
-<note tip> +
-You should run these calculations on 3 cores with ''bsub -n 3''. +
-</note>+
  
 +In order to speed up the calculation, we are going to use the SPC/E water model.
 +The model is a non polarizable forcefield model, with parameters for:
 +  * Charge of the atomic species
 +  * Harmonic O-H bond elongation $ U_{bond}= k_{bond}⋅(r-r_{eq})^2$ 
 +  * Harmonic H-O-H angle bending $ U_{angle}= k_{angle}⋅(θ-θ_{eq})^2$ 
 +  * Lennard Jones interaction between non bonded species (interaction between atoms belonging to different water molecules)
 + (A DFT-version of the calculation can be found here: [[doi>10.1021/jp4103355]]).
  
-===== Task 1: Calculate the dielectric constant ===== +===== Introduction =====
-The dielectric constant describes the response of a system to an external electric field. Within a linear response approximation the //Kubo Formula// is applicable. The Kubo Formula is an equation which expresses the linear response of an observable quantity due to a time-dependent perturbation. Hence, one can calculate the dielectric constant based on a sampling of the dipole moment.+
  
-The advantage of Monte Carlo is that we can employ special hydrogen reordering moves, which lead to an effective sampling of the dipole. +Water molecules are arranged according to the so called //"ice rules"//:  
 +  - Water molecules are present as neutral H2O  
 +  - Each molecule makes four hydrogen bonds with its four nearest neighbors, two as a hydrogen bond donor and two as an acceptor.  
 +All the possible proton arrangements in the ice strucrure are not energetically equivalent and we sample and weight them with a Monte Carlo approach. For this purpose, a specific algorithm was created (see: [[doi>10.1063/1.1568337]]), that allows to employ special proton reordering moves, which lead to an effective sampling of the ice dipole (N.B: different proton arrangement = different dipole). \\
  
-Run the input-file ''mc_exercise.inp''. It will create a file ''tmc_trajectory_T200.00.dip_cl'', which contains the dipole moment for each accepted step. Plot the histogram of the z-component of the dipole moment. The distribution should be symmetric, because the simulation cell is also symmetric in the z-direction.+The dielectric constant of a system describes its response to an external electric field.  
 +If the dipole moment is properly sampled, one can compute the dielectric constant of ice, by applying the //Kubo Formula//. This is valid in the approximation that the response of the system to the time-dependent perturbation (the field) is linear. \\
  
-The dielectric constant can then calculated from the dipole moments via:+The dielectric constant can then be calculated from the dipole moments via:
 \begin{equation} \begin{equation}
 \epsilon = 1 + \left(\frac{4 \pi}{3 \epsilon_0 V  k_B T } \right ) \operatorname{Var}(M) \ , \epsilon = 1 + \left(\frac{4 \pi}{3 \epsilon_0 V  k_B T } \right ) \operatorname{Var}(M) \ ,
 \end{equation} \end{equation}
  
-where $M$ denotes the dipol moment of the entire simulation cell and $\operatorname{Var}(M)$ denotes the variance of the dipol moment of the sampling:+where $M$ denotes the dipole moment of the entire simulation cell and $\operatorname{Var}(M)$ denotes the variance of the dipole moment of the sampling:
 \begin{equation} \begin{equation}
 \operatorname{Var}(M) = (\langle M \cdot M\rangle - \langle M\rangle\langle M\rangle ) \ . \operatorname{Var}(M) = (\langle M \cdot M\rangle - \langle M\rangle\langle M\rangle ) \ .
 \end{equation} \end{equation}
  
-To simplify this task you can use the following python-script:+===== Task 1: Calculate the dielectric constant from one simulation ===== 
 +      * As usual, download the 4 required files (at the end of this page: input, auxiliary input, ice coordinates and topology) in the same directory. \\ 
 +      * Run the input-file ''mc_exercise.inp'' 
 +⇒ It will create a file ''tmc_trajectory_T200.00.dip_cl'', which contains the dipole moment for each accepted step.  
 +<note tip> 
 +You should run these calculations on 3 cores with ''bsub -n 3''. 
 +</note> 
 +    * Plot the histogram of the z-component of the dipole moment. \\ 
 +⇒ The distribution should be symmetric, because the simulation cell is also symmetric in the z-direction. 
 +  * Compute the dielectric constant of ice. To simplify this task you can use the following python-script: 
  
 <code python calc_dielectric_constant.py> <code python calc_dielectric_constant.py>
Line 78: Line 88:
 Before you can run the python-script you have to load a newer python-module and set the executable-bit: Before you can run the python-script you have to load a newer python-module and set the executable-bit:
 <code> <code>
-you@brutusX ~$ module load new python/2.7.6 +you@eulerX ~$ module load python/2.7.6 
-you@brutusX ~$ chmod +x calc_dielectric_constant.py +you@eulerX ~$ chmod +x calc_dielectric_constant.py 
-you@brutusX ~$ ./calc_dielectric_constant.py 200 tmc_trajectory_T200.00.dip_cl+you@eulerX ~$ ./calc_dielectric_constant.py 200 tmc_trajectory_T200.00.dip_cl
 </code> </code>
 </note> </note>
Line 87: Line 97:
 You can gather more samples by launching multiple independent runs in parallel with different random number seeds. The seed is given by the ''RND_DETERMINISTIC'' keyword. The gathered trajectories can then be analyzed collectively with the python-script: You can gather more samples by launching multiple independent runs in parallel with different random number seeds. The seed is given by the ''RND_DETERMINISTIC'' keyword. The gathered trajectories can then be analyzed collectively with the python-script:
 <code> <code>
-you@brutusX ~$ ./calc_dielectric_constant.py 200 trj_1_T200.00.dip_cl trj_1_T200.00.dip_cl ...+you@eulersX ~$ ./calc_dielectric_constant.py 200 trj_1_T200.00.dip_cl trj_1_T200.00.dip_cl ...
 </code> </code>
 <note tip> <note tip>
 You can also share trajectories with your colleges. You can also share trajectories with your colleges.
 </note> </note>
-===== Task 3: Calculate the thermal expansion coefficient ===== +
-Run the MC sampling again at temperature of 150K . Plot the cell volume of each sampling. Read off the converged cell volume and plot it vs. the temperature. Determine the expansion coefficient via a linear fit. Do you trust your result? Why (not)?+
  
 ===== Required Files ===== ===== Required Files =====
 ==== Initial Coordinates ==== ==== Initial Coordinates ====
 +This file contains the initial ice coordinates. 
 {{ice_ih_96.xyz.gz| Download here}} {{ice_ih_96.xyz.gz| Download here}}
  
Line 104: Line 114:
 &GLOBAL &GLOBAL
   PROJECT H2O_MC   PROJECT H2O_MC
-  PROGRAM TMC+  PROGRAM TMC         ! Tree Monte Carlo algorithm 
   RUN_TYPE TMC   RUN_TYPE TMC
-  PRINT_LEVEL LOW +  PRINT_LEVEL LOW     ! Low amount of information written in the output 
-  WALLTIME 1:00:00+  WALLTIME 1:00:00    ! The simulation will last one hour and then will stop 
 &END GLOBAL &END GLOBAL
 &MOTION &MOTION
   &TMC   &TMC
       RND_DETERMINISTIC 42 !<=== Change this number to obtain different samplings       RND_DETERMINISTIC 42 !<=== Change this number to obtain different samplings
-      PRINT_COORDS .FALSE. !this avoids the file-size problem+      PRINT_COORDS .FALSE. !this avoids the printing of all coordinates and file-size problems
       GROUP_CC_SIZE 0       GROUP_CC_SIZE 0
       NUM_MC_ELEM 100000       NUM_MC_ELEM 100000
-      ENERGY_FILE_NAME H2O_ice.inp+      ENERGY_FILE_NAME H2O_ice.inp ! refers to the auxiliary input for the force specification
       TEMPERATURE 200       TEMPERATURE 200
-      &MOVE_TYPE      MOL_TRANS+      &MOVE_TYPE      MOL_TRANS       ! specifies "proton moves" for better sampling
         SIZE          0.05         SIZE          0.05
         PROB          3         PROB          3
Line 137: Line 147:
       &END       &END
       PRESSURE 0.01       PRESSURE 0.01
-      ESIMATE_ACC_PROB .TRUE.+      ESIMATE_ACC_PROB .TRUE.        !accuracy parameters, do not change
       RESTART_OUT 0       RESTART_OUT 0
       NUM_MV_ELEM_IN_CELL 5       NUM_MV_ELEM_IN_CELL 5
exercises/2015_ethz_mmm/monte_carlo_ice.txt · Last modified: 2020/08/21 10:15 by 127.0.0.1