====== 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]]. You should run these calculations on 3 cores with ''bsub -n 3''. ===== Task 1: Calculate the dielectric constant ===== 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. 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 can then 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 dipol moment of the entire simulation cell and $\operatorname{Var}(M)$ denotes the variance of the dipol moment of the sampling: \begin{equation} \operatorname{Var}(M) = (\langle M \cdot M\rangle - \langle M\rangle\langle M\rangle ) \ . \end{equation} To simplify this task you can use the following python-script: #!/usr/bin/env python import sys import numpy as np if(len(sys.argv) < 3): print "Usage: calc_dielectric_constant.py ... " sys.exit() T = float(sys.argv[1]) print "Temperature: %.2f K"%T weights = []; M_vec = []; cell_vol = [] for fn in sys.argv[2:]: print "Reading: "+fn assert(fn.endswith(".dip_cl")) raw_data = np.loadtxt(fn) # [C*Angstrom] weights.append(raw_data[1:,0]-raw_data[:-1,0]) M_vec.append(raw_data[:-1,1:4]) cell_vol.append(np.loadtxt(fn.replace(".dip_cl", ".cell"), usecols=(10,))) 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"%np.sum(weights) print "Mean cell volume: %.2f Angstrom^3"%V M = np.sqrt(np.sum(np.square(M_vec), axis=1)) # take norm of dipol moment var = np.average(np.square(M), weights=weights) - np.square(np.average(M, weights=weights)) 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)/(3*V*kb*T*angstrom2meter*epsilon_0) epsilon = 1 + s*var print "Dielectric constant: %.2f"%epsilon Before you can run the python-script you have to load a newer python-module and set the executable-bit: you@brutusX ~$ module load new python/2.7.6 you@brutusX ~$ chmod +x calc_dielectric_constant.py you@brutusX ~$ ./calc_dielectric_constant.py 200 tmc_trajectory_T200.00.dip_cl ===== 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 ''RND_DETERMINISTIC'' keyword. The gathered trajectories can then be analyzed collectively with the python-script: you@brutusX ~$ ./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. ===== 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 ===== ==== Initial Coordinates ==== {{ice_ih_96.xyz.gz| Download here}} ==== Main Input-File (run this) ==== &GLOBAL PROJECT H2O_MC PROGRAM TMC RUN_TYPE TMC PRINT_LEVEL LOW WALLTIME 1:00:00 &END GLOBAL &MOTION &TMC RND_DETERMINISTIC 42 !<=== Change this number to obtain different samplings PRINT_COORDS .FALSE. !this avoids the file-size problem GROUP_CC_SIZE 0 NUM_MC_ELEM 100000 ENERGY_FILE_NAME H2O_ice.inp TEMPERATURE 200 &MOVE_TYPE MOL_TRANS SIZE 0.05 PROB 3 &END &MOVE_TYPE ATOM_TRANS SIZE 0.01 PROB 3 &END &MOVE_TYPE MOL_ROT SIZE 5 PROB 3 &END &MOVE_TYPE PROT_REORDER PROB 5 &END &MOVE_TYPE VOL_MOVE 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. &TMC_ANALYSIS CLASSICAL_DIPOLE_MOMENTS RESTART .FALSE. &CHARGE ATOM O CHARGE -0.8476 &END CHARGE &CHARGE ATOM H CHARGE 0.4238 &END CHARGE &END TMC_ANALYSIS &END TMC &END MOTION ==== Auxiliary Input-File (do not run this directly) ==== &FORCE_EVAL METHOD FIST &MM &NEIGHBOR_LISTS GEO_CHECK OFF VERLET_SKIN 4.0 &END NEIGHBOR_LISTS &FORCEFIELD &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 O H THETA0 [deg] 104.500 K [rad^2kjmol] 627.600 KIND G87 &END BEND &NONBONDED &LENNARD-JONES ATOMS O O EPSILON [kjmol] 0.650 SIGMA [angstrom] 3.166 RCUT 8.0 &END LENNARD-JONES &LENNARD-JONES ATOMS O H EPSILON 0.0 SIGMA 3.1655 RCUT 8.0 &END LENNARD-JONES &LENNARD-JONES ATOMS H H EPSILON 0.0 SIGMA 3.1655 RCUT 8.0 &END LENNARD-JONES &END NONBONDED &END FORCEFIELD &POISSON &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 &TOPOLOGY COORD_FILE_NAME ./ice_ih_96.xyz COORD_FILE_FORMAT xyz CONNECTIVITY MOL_SET &MOL_SET &MOLECULE NMOL 96 CONN_FILE_NAME topology_H2O.psf &END &END &END TOPOLOGY &END SUBSYS &END FORCE_EVAL &GLOBAL &END GLOBAL ==== Topology File ==== PSF 1 !NTITLE Topology file for water 3 !NATOM 1 WAT 1 WAT O O 0.000000 15.999400 0 2 WAT 1 WAT H H 0.000000 1.008000 0 3 WAT 1 WAT H H 0.000000 1.008000 0 2 !NBOND 1 2 1 3 1 !NTHETA 2 1 3 0 !NPHI 0 !NIMPHI 0 !NDON 0 !NACC 0 !NNB