exercises:2015_ethz_mmm:monte_carlo_ice
Differences
This shows you the differences between two versions of the page.
Next revision | Previous revisionNext revisionBoth sides next revision | ||
exercises:2015_ethz_mmm:monte_carlo_ice [2015/02/06 17:49] – external edit 127.0.0.1 | exercises: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 '' | ||
- | - Add the line '' | ||
- | </ | ||
- | In this exercise we will use Monte Carlo sampling to calculate the [[wp> | + | In this exercise we will use Monte Carlo sampling to calculate the [[wp> |
- | + | ||
- | + | ||
- | <note tip> | + | |
- | You should run these calculations on 3 cores with '' | + | |
- | </ | + | |
+ | In order to speed up the calculation, | ||
+ | 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> | ||
- | ===== 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 | + | 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> | ||
- | Run the input-file '' | + | The dielectric constant of a system describes its response to an external electric field. |
+ | If the dipole moment | ||
- | 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 |
\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 '' | ||
+ | ⇒ It will create a file '' | ||
+ | <note tip> | ||
+ | You should run these calculations on 3 cores with '' | ||
+ | </ | ||
+ | * 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: | ||
< | < | ||
- | you@brutusX | + | you@eulerX |
- | you@brutusX | + | you@eulerX |
- | you@brutusX | + | you@eulerX |
</ | </ | ||
</ | </ | ||
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 '' | You can gather more samples by launching multiple independent runs in parallel with different random number seeds. The seed is given by the '' | ||
< | < | ||
- | you@brutusX | + | you@eulersX |
</ | </ | ||
<note tip> | <note tip> | ||
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}} | ||
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: |
&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 | + | PRINT_COORDS .FALSE. !this avoids the printing of all coordinates and file-size |
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 |
TEMPERATURE 200 | TEMPERATURE 200 | ||
- | & | + | & |
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. |
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