# Open SourceMolecular Dynamics

### Site Tools

Writing /var/www/cp2k.org/www/data/meta/exercises/2015_ethz_mmm/monte_carlo_ice.meta failed

exercises:2015_ethz_mmm:monte_carlo_ice

# Differences

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

 exercises:2015_ethz_mmm:monte_carlo_ice [2015/02/06 17:49]127.0.0.1 external edit exercises:2015_ethz_mmm:monte_carlo_ice [2015/03/18 15:23] (current)sclelia 2015/03/18 15:23 sclelia 2015/03/18 15:18 sclelia [Main Input-File (run this)] 2015/03/18 15:12 sclelia [Initial Coordinates] 2015/03/18 15:11 sclelia 2015/03/18 14:53 sclelia 2015/02/06 17:49 external edit Next revision Previous revision 2015/03/18 15:23 sclelia 2015/03/18 15:18 sclelia [Main Input-File (run this)] 2015/03/18 15:12 sclelia [Initial Coordinates] 2015/03/18 15:11 sclelia 2015/03/18 14:53 sclelia 2015/02/06 17:49 external edit Line 1: Line 1: ====== Properties of Ice from Monte Carlo Simulations ====== ====== Properties of Ice from Monte Carlo Simulations ====== - ​ - - 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. - ​ - 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. \\ - + - + - + - You should run these calculations on 3 cores with ''​bsub -n 3''​. + - ​ + + 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: \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) \ , - 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: \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 ) \ . - 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. + + You should run these calculations on 3 cores with : \\ + ''​bsub -W 01:00 -n 3 mpirun cp2k.popt -i mc_exercise.inp -o mc_exercise.out''​ + ​ + * 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: ​ ​ Line 78: Line 89: 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 ​ Line 87: Line 98: 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 ... You can also share trajectories with your colleges. You can also share trajectories with your colleges. ​ - ===== 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}} ==== Main Input-File (run this) ==== ==== Main Input-File (run this) ==== + This file specifies the type of algorithm do adopt and the simulation parameters. ​ &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 149: &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 Line 157: Line 169: ​ - ==== Auxiliary Input-File (do not run this directly) ==== + ==== Auxiliary Input-File (do NOT run this directly) ==== + This auxiliary input specifies the system properties ​ Line 242: Line 255: ==== Topology File ==== ==== Topology File ==== + This file specifies the topology of the system (which atoms, bonds, angles in the water molecule). It allows the program to distinguish which atoms belong to which water molecule, and therefore discriminate between inter and intramolecular interactions. ​ PSF PSF