exercises:2016_ethz_mmm:monte_carlo_ice
Differences
This shows you the differences between two versions of the page.
— | exercises:2016_ethz_mmm:monte_carlo_ice [2020/08/21 10:15] (current) – created - external edit 127.0.0.1 | ||
---|---|---|---|
Line 1: | Line 1: | ||
+ | ====== Properties of Ice from Monte Carlo Simulations ====== | ||
+ | |||
+ | In this exercise we will use Monte Carlo sampling to calculate the [[wp> | ||
+ | |||
+ | 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> | ||
+ | |||
+ | ===== Introduction ===== | ||
+ | |||
+ | 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> | ||
+ | |||
+ | 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 be calculated from the dipole moments via: | ||
+ | \begin{equation} | ||
+ | \epsilon = 1 + \left(\frac{4 \pi}{3 \epsilon_0 V k_B T } \right ) \operatorname{Var}(M) \ , | ||
+ | \end{equation} | ||
+ | |||
+ | 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} | ||
+ | \operatorname{Var}(M) = (\langle M \cdot M\rangle - \langle M\rangle\langle M\rangle ) \ . | ||
+ | \end{equation} | ||
+ | |||
+ | ===== 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> | ||
+ | # | ||
+ | |||
+ | import sys | ||
+ | import numpy as np | ||
+ | |||
+ | if(len(sys.argv) < 3): | ||
+ | print " | ||
+ | sys.exit() | ||
+ | |||
+ | T = float(sys.argv[1]) | ||
+ | print " | ||
+ | |||
+ | weights = []; M_vec = []; cell_vol = [] | ||
+ | for fn in sys.argv[2: | ||
+ | print " | ||
+ | assert(fn.endswith(" | ||
+ | raw_data = np.loadtxt(fn) # [C*Angstrom] | ||
+ | weights.append(raw_data[1:, | ||
+ | M_vec.append(raw_data[: | ||
+ | cell_vol.append(np.loadtxt(fn.replace(" | ||
+ | |||
+ | weights = np.concatenate(weights) | ||
+ | M_vec = np.concatenate(M_vec) | ||
+ | cell_vol = np.concatenate(cell_vol) | ||
+ | |||
+ | V = np.mean(cell_vol) # [Angstrom^3] | ||
+ | print "Total number of samples: %d" | ||
+ | print "Mean cell volume: %.2f Angstrom^3" | ||
+ | |||
+ | M = np.sqrt(np.sum(np.square(M_vec), | ||
+ | var = np.average(np.square(M), | ||
+ | |||
+ | epsilon_0 = 8.8541878176e-12 #[F/m] vacuum permittivity | ||
+ | e = 1.602176565e-19 # [C] elementary charge | ||
+ | kb = 1.3806488e-23 #[J/K] Boltzmann constant | ||
+ | angstrom2meter = 1e-10 | ||
+ | s = (e*e*4*np.pi)/ | ||
+ | |||
+ | epsilon = 1 + s*var | ||
+ | print " | ||
+ | </ | ||
+ | |||
+ | <note tip> | ||
+ | Before you can run the python-script you have to load a newer python-module and set the executable-bit: | ||
+ | < | ||
+ | you@eulerX ~$ module load python/ | ||
+ | you@eulerX ~$ chmod +x calc_dielectric_constant.py | ||
+ | you@eulerX ~$ ./ | ||
+ | </ | ||
+ | </ | ||
+ | |||
+ | ===== Task 2: Gather more samplings ===== | ||
+ | You can gather more samples by launching multiple independent runs in parallel with different random number seeds. The seed is given by the '' | ||
+ | < | ||
+ | you@eulersX ~$ ./ | ||
+ | </ | ||
+ | <note tip> | ||
+ | You can also share trajectories with your colleges. | ||
+ | </ | ||
+ | |||
+ | |||
+ | ===== Required Files ===== | ||
+ | ==== Initial Coordinates ==== | ||
+ | This file contains the initial ice coordinates. | ||
+ | {{ice_ih_96.xyz.gz| Download here}} | ||
+ | |||
+ | |||
+ | ==== Main Input-File (run this) ==== | ||
+ | This file specifies the type of algorithm do adopt and the simulation parameters. | ||
+ | <code - mc_exercise.inp> | ||
+ | &GLOBAL | ||
+ | PROJECT H2O_MC | ||
+ | PROGRAM TMC ! Tree Monte Carlo algorithm | ||
+ | RUN_TYPE TMC | ||
+ | PRINT_LEVEL LOW ! Low amount of information written in the output | ||
+ | WALLTIME 1: | ||
+ | &END GLOBAL | ||
+ | &MOTION | ||
+ | &TMC | ||
+ | RND_DETERMINISTIC 42 !<=== Change this number to obtain different samplings | ||
+ | PRINT_COORDS .FALSE. !this avoids the printing of all coordinates and file-size problems | ||
+ | GROUP_CC_SIZE 0 | ||
+ | NUM_MC_ELEM 100000 | ||
+ | ENERGY_FILE_NAME H2O_ice.inp ! refers to the auxiliary input for the force specification | ||
+ | TEMPERATURE 200 | ||
+ | & | ||
+ | SIZE 0.05 | ||
+ | PROB 3 | ||
+ | &END | ||
+ | & | ||
+ | SIZE 0.01 | ||
+ | PROB 3 | ||
+ | &END | ||
+ | & | ||
+ | SIZE 5 | ||
+ | PROB 3 | ||
+ | &END | ||
+ | & | ||
+ | PROB 5 | ||
+ | &END | ||
+ | & | ||
+ | SIZE 0.1 | ||
+ | PROB 1 | ||
+ | &END | ||
+ | PRESSURE 0.01 | ||
+ | ESIMATE_ACC_PROB .TRUE. | ||
+ | RESTART_OUT 0 | ||
+ | NUM_MV_ELEM_IN_CELL 5 | ||
+ | PRINT_ONLY_ACC .FALSE. | ||
+ | & | ||
+ | CLASSICAL_DIPOLE_MOMENTS | ||
+ | RESTART .FALSE. | ||
+ | &CHARGE | ||
+ | ATOM O | ||
+ | | ||
+ | &END CHARGE | ||
+ | &CHARGE | ||
+ | ATOM H | ||
+ | | ||
+ | &END CHARGE | ||
+ | &END TMC_ANALYSIS | ||
+ | &END TMC | ||
+ | &END MOTION | ||
+ | </ | ||
+ | |||
+ | ==== Auxiliary Input-File (do NOT run this directly) ==== | ||
+ | This auxiliary input specifies the system properties | ||
+ | <code - H2O_ice.inp> | ||
+ | |||
+ | & | ||
+ | METHOD FIST | ||
+ | &MM | ||
+ | & | ||
+ | GEO_CHECK OFF | ||
+ | VERLET_SKIN 4.0 | ||
+ | &END NEIGHBOR_LISTS | ||
+ | & | ||
+ | &SPLINE | ||
+ | EMAX_SPLINE 10.0 | ||
+ | &END | ||
+ | &CHARGE | ||
+ | ATOM O | ||
+ | CHARGE -0.8476 | ||
+ | &END CHARGE | ||
+ | &CHARGE | ||
+ | ATOM H | ||
+ | CHARGE 0.4238 | ||
+ | &END CHARGE | ||
+ | &BOND | ||
+ | ATOMS O H | ||
+ | K [nm^-2kjmol] 502080.0 | ||
+ | R0 [nm] 0.09572 | ||
+ | KIND G87 | ||
+ | &END BOND | ||
+ | &BEND | ||
+ | ATOMS H | ||
+ | THETA0 | ||
+ | K [rad^2kjmol] 627.600 | ||
+ | KIND G87 | ||
+ | &END BEND | ||
+ | & | ||
+ | & | ||
+ | ATOMS O O | ||
+ | EPSILON [kjmol] 0.650 | ||
+ | SIGMA [angstrom] 3.166 | ||
+ | RCUT 8.0 | ||
+ | &END LENNARD-JONES | ||
+ | & | ||
+ | ATOMS O H | ||
+ | EPSILON 0.0 | ||
+ | SIGMA 3.1655 | ||
+ | RCUT 8.0 | ||
+ | &END LENNARD-JONES | ||
+ | & | ||
+ | ATOMS H H | ||
+ | EPSILON 0.0 | ||
+ | SIGMA 3.1655 | ||
+ | RCUT 8.0 | ||
+ | &END LENNARD-JONES | ||
+ | &END NONBONDED | ||
+ | &END FORCEFIELD | ||
+ | & | ||
+ | &EWALD | ||
+ | EWALD_TYPE ewald | ||
+ | ALPHA .45 | ||
+ | GMAX 19 | ||
+ | &END EWALD | ||
+ | &END POISSON | ||
+ | &END MM | ||
+ | &SUBSYS | ||
+ | &CELL | ||
+ | ABC 13.5084309945 15.6001709945 14.7072509945 | ||
+ | &END CELL | ||
+ | & | ||
+ | COORD_FILE_NAME ./ | ||
+ | COORD_FILE_FORMAT xyz | ||
+ | CONNECTIVITY MOL_SET | ||
+ | & | ||
+ | & | ||
+ | NMOL 96 | ||
+ | CONN_FILE_NAME topology_H2O.psf | ||
+ | &END | ||
+ | &END | ||
+ | &END TOPOLOGY | ||
+ | &END SUBSYS | ||
+ | &END FORCE_EVAL | ||
+ | &GLOBAL | ||
+ | &END GLOBAL | ||
+ | </ | ||
+ | |||
+ | ==== 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. | ||
+ | <code - topology_H2O.psf> | ||
+ | PSF | ||
+ | |||
+ | 1 !NTITLE | ||
+ | | ||
+ | |||
+ | | ||
+ | 1 WAT 1 WAT O O 0.000000 | ||
+ | 2 WAT 1 WAT H H 0.000000 | ||
+ | 3 WAT 1 WAT H H 0.000000 | ||
+ | |||
+ | | ||
+ | | ||
+ | |||
+ | | ||
+ | | ||
+ | |||
+ | | ||
+ | |||
+ | | ||
+ | |||
+ | | ||
+ | |||
+ | | ||
+ | |||
+ | | ||
+ | </ | ||