Properties of Ice from Monte Carlo Simulations

Properties of Ice from Monte Carlo Simulations

In this exercise we will use Monte Carlo sampling to calculate the dielectric constant of the disordered hexagonal phase (Ih) of water ice.

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: 10.1021/jp4103355).

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: 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).

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 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:

calc_dielectric_constant.py

#!/usr/bin/env python
 
import sys
import numpy as np
 
if(len(sys.argv) < 3):
    print "Usage: calc_dielectric_constant.py <temperature> <traj_1.dip_cl> ... <traj_N.dip_cl>"
    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@eulerX ~$ module load python/2.7.6
you@eulerX ~$ chmod +x calc_dielectric_constant.py
you@eulerX ~$ ./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@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.

Required Files

Initial Coordinates

This file contains the initial ice coordinates. Download here

Main Input-File (run this)

This file specifies the type of algorithm do adopt and the simulation parameters.

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:00:00    ! The simulation will last one hour and then will stop 
&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
      &MOVE_TYPE      MOL_TRANS       ! specifies "proton moves" for better sampling
        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.        !accuracy parameters, do not change
      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)

This auxiliary input specifies the system properties

H2O_ice.inp

&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

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.

topology_H2O.psf

 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

Table of Contents