Analysis of molecular dynamics simulation of water

 

Objective: to demonstrate the analysis tools and to test the accuracy of models of water.

 

First you will want to examine the output file to see that the run was successful.

Once the molecular dynamics DISCOVER3 job is complete you can analyze the results using the DECIPHER module. Go to the left corner to find the pull-down menu and select DECIPHER. In the following I will assume that your file has the name water_1.*. Use the actual name of the *.car, *.mdf and other files wherever you see [water_1].

 

Read in the water_1.car file using the Molecule/Get command in INSIGHT. After running dynamics you will see an apparently unstructured argon solvent as shown below.

 

Use the side pull-down menu to run DECIPHER. In decipher use the Configurations/Get command to generate a menu. Use the Files list on the right-hand side to get the water_1.his file.

Configurations File [water_1.his]

Assembly/molecule [WATER_1]

Load_Cell [on] (check the box to turn it yellow)

Load_Pressure [off] (do not check the box)

Load_Velocities [on] (check the box to turn it yellow)

Frame Spec [*] (reads all of the frames)

Select execute.


It is interesting to observe the trajectory visually. You can do this by using the Configurations/Animate command. When you execute this command the frames will be loaded into memory and played step by step on the monitor so that you can observe the motion of the molecules.

 

Radial Distribution Functions

 

To calculate the radial distribution function averaged over the entire trajectory (i.e. averaging all of the frames in the history file) use the

Spectra/RadialDistFunc command. For argon there was only one radial distribution function and that depended on the Ar-Ar distances. For water, there are three radial distribution functions, the O-O, O-H, and H-H functions.

For example, to generate the O – H pair distribution in the RadialDistFunc menu

Distribution_type [pair]

Molecule_Spec [::o*]

Molecule_2_Spec [::h*]

Coord_Number [on] (check the box to turn it yellow)

R_Min [0.0]

R_Max [10.0]

R_Step [0.1]

Bulk_Density [on] (check the box to turn it yellow)

Data_Sources [trajectory] (check the box to turn it yellow)

Background [on] (check the box to turn it yellow)

Create_TBL_File [rdf] (check the box to turn it yellow)

 

In the Parameters menu box (on the right) the Assembly/Molecule level should be set to Molecule.

 

The procedure for the other radial distribution functions is similar except that you need to replace ::o* by ::h* to obtain the H-H distribution function and you need to replace ::h* by ::o* to obtain the O-O distribution function.

 

Mean Square Displacement

 

To calculate the mean squared displacement use the Spectra/Mean_Sq_Disp command.

Molecule_Spec [WATER_1]

Data_Source [trajectory]

Frame_Spec [*]

Output

MSD [on] (check the box to turn it yellow)

Dist_Traveled [on] (check the box to turn it yellow)

Background (check the box to turn it yellow)

Create_TBL_File [msd] (check the box to turn it yellow)

 

Velocity Autocorrelation Function

 

To calculate the velocity correlation function and associated spectral density first use the Functions/Velocity menu.

Function/Name [velocity1]

Definition levels

Per_Atom (check the diamond to turn it yellow)

Molecule_Spec [WATER_1]

 

Then use the Spectra/TimeCorrFunc pulldown

Function Spec [velocity 1] (you can find this under Parameters/Function Levels Vector on the right side or you can simply type it in).

Function 2 Spec [same]

Output

Corr_Func [on] (check the box to turn it yellow)

Spectral_Density [on] (check the box to turn it yellow)

Zero padding [on] (check the box to turn it yellow)

Calculation method

Direct [on] (check the box to turn it yellow)

 

t0_step [1]

t_length [200]

Frame Spec [*]

Background (check the box to turn it yellow)

Create_TBL_File [vacf] (check the box to turn it yellow)

 

Dipole Autocorrelation Function

 

To calculate the dipole correlation function and associated spectral density first use the Functions/Dipole menu.

Function/Name [dipole1] (note the 1 at the end is arbitrary in the definition)

Definition levels

Per_Molecule [on] (check the diamond to turn it yellow)

Molecule_Spec [WATER_1]

 

Then use the Spectra/TimeCorrFunc pulldown

Function Spec [dipole1] (you can find this under Parameters/Function Levels Vector on the right side or you can simply type it in).

Function 2 Spec [same]

Output

Corr_Func [on] (check the box to turn it yellow)

Spectral_Density [on] (check the box to turn it yellow)

Zero padding [on] (check the box to turn it yellow)

Calculation method

Direct [on] (check the box to turn it yellow)

t0_step [1]

t_length [200]

Frame Spec [*]

Background (check the box to turn it yellow)

Create_TBL_File [vacf] (check the box to turn it yellow)

 

Average Quantities

 

The average internal energy, pressure, are temperature are found at the bottom of the *.out output file. In addition, to the average values please note the root-mean-square deviations of these values. These will be interesting when comparing different simulation conditions, e.g. different time steps, potential cutoff choice, number of atoms (size of system) etc.

 

File Handling

 

Once you have obtained the results they should be labeled using your personal code at the # sign as follows

h2o_#_1.oordf - O-O radial distribution function

h2o_#_1.ohrdf - O-H radial distribution function

h2o_#_1.hhrdf - H-H radial distribution function

h2o_#_1.msd - mean square displacement

h2o_#_1.vacf - velocity autocorrelation function

h2o_#_1.vsd - velocity spectral density

h2o_#_1.dacf - dipole autocorrelation function

h2o_#_1.dsd - dipole spectral density

h2o_#_1.ave - average energy, pressure, and temperature from output file

 

If you have more than one condition as part of your assignment the ending number may need to be changed accordingly.

 

These data will be used for comparisons in class.