Practical 11: Liquid Argon (again)

In this practical we learn to use the Gromacs MD package. As the first example we study (again) liquid argon.

Today you should

1. install Gromacs
2. run a MD simulation
3. analyze the output

You should record your results in your own Wiki page. Therefore, go to the LJ Fluid in Gromacs Wiki on Blackboard (Course Documents ‣ Practicals ‣ Practical 11 ‣ LJ Fluid in Gromacs Wiki) and create your own results page in the wiki.

Installation

All files can be found on the page for Practical 11.

Gromacs 4.5.5 Mac OS X 10.6.8

Install the pre-compiled Gromacs 4.5.5 distribution on your iMac. Get files and unpack in your home directory:

curl -O http://becksteinlab.physics.asu.edu/pages/courses/2013/SimBioNano/11/gromacs455_static_MacOSX_10.6.8.tgz
cd ~
tar -zxvf gromacs455_static_MacOSX_10.6.8.tgz

Add the $HOME/opt/bin directory to your PATH in ~/.profile:

echo 'export PATH=$HOME/opt/bin:$PATH' >> ~/.profile
. ~/.profile

so that you can automatically find the Gromacs commands in ~/opt/bin. (This only needs to be done once!)

Input files for the Argon simulations

Go to your work directory for this practical and get the input files:

mkdir -p NAME/p11 && cd NAME/p11
curl -O http://becksteinlab.physics.asu.edu/pages/courses/2013/SimBioNano/11/Argon_GMX.tar.gz
tar -zxvf Argon_GMX.tar.gz

You can also put the contents of scripts into ~/opt/bin:

cp scripts/* ~/opt/bin

NVE MD simulation of liquid Argon

Make a directory for the simulation, e.g.

mkdir argon_94K
cd argon_94K

Input structure

Generate cubic system of Ar, 512 atoms with scripts/generate_lattice.py script:

generate_lattice.py -l cubic 512

(You might have to give the full path to scripts/generate_lattice.py, depending on where you installed it.)

MD simulation

Run MD at

• 94 K
• dt = 0.01 ps
• 250 ps (or longer)

Start with the parameters in LJ.mdp (suitable for a NVE simulation of a liquid of neutral molecules) and the topology file argon.top (defines the LJ parameters for argon and describes the system to Gromacs). Copy the input files to your simulation directory:

cp ../Argon/{LJ.mdp,argon.top} .

Prepare the Gromacs environment with GMXRC, create the “binary run input file” md.tpr with the Gromacs preprocessor grompp, and run the simulation using the mdrun command:

. ~/opt/bin/GMXRC
grompp -p argon.top -c Ar_rho=0.0207125_N=512_cubic.pdb -f LJ.mdp -o md.tpr
mdrun -v -deffnm md

This should not take longer than 2 minutes.

Gromacs basics

Input files:

• conformation = coordinates (pdb, gro)
• topology (top)
• run parameters (mdp)

From those files you create a single, portable run input file (“tpr”). You can copy this file to another machine (e.g. a supercomputer) and run the simulation with only this file as an input.

Output files:

• trajectory, “XTC file” (md.xtc; one can also write a “trr” file that can contain velocities and forces). The XTC is a very efficient storage format due to a lossy compression algorithm, that only stores coordinates to a given precision, by default 1000th of a nm. TRR files are full precision but take up a lot of space. In most cases we only need XTCs. Both XTC and TRR are binary formats that are not human readable.
• energy file (md.edr); must be read with the g_energy tool
• log file (md.log); human readable.
• last frame (md.gro, but using the -c md.pdb option to mdrun you can also write a PDB file.)

Analysis

Record your results in your own Wiki page in the LJ Fluid in Gromacs Wiki on Blackboard (Course Documents ‣ Practicals ‣ Practical 11 ‣ LJ Fluid in Gromacs Wiki).

Visualization

Visualize in VMD by using the md.gro file as the topology (load it first) and the load the md.xtc trajectory into the same molecule (number 0).

Show a representative snapshot of your system.

Tip: You can show the unitcell with the tcl commands (type at the VMD> prompt):

pbc box

Energy and Temperature

Use the g_energy tool to analyze the md.edr file, which contains the energy terms collected over the simulation.

g_energy -f md.edr -o energy  -xvg none
4 5 6 7 9

Select the terms you want from the following list by selecting either (part of) the name or the number or a combination.:

End your selection with an empty line or a zero.
-------------------------------------------------------------------
  1  LJ-(SR)          2  Disper.-corr.    3  Coulomb-(SR)     4  Potential
  5  Kinetic-En.      6  Total-Energy     7  Temperature      8  Pres.-DC
  9  Pressure        10  Vir-XX          11  Vir-XY          12  Vir-XZ

Example output:

Statistics over 25001 steps [ 0.0000 through 250.0000 ps ], 5 data sets
All statistics are over 25001 points

Energy                      Average   Err.Est.       RMSD  Tot-Drift
-------------------------------------------------------------------------------
Potential                  -2947.63        0.3    17.8312   -1.73157  (kJ/mol)
Kinetic En.                 643.805       0.31    17.8422    1.72684  (kJ/mol)
Total Energy               -2303.83     0.0075   0.187319 -0.00471748  (kJ/mol)
Temperature                  101.02      0.048    2.79962   0.270959  (K)
Pressure                    167.781        1.3    60.0152   -7.39932  (bar)

One needs to write Python code to read the xvg files (or use the xmgrace program (not installed on the iMacs!)).

The following code plots energies in temperatures in a 2x2 layout:

import numpy as np
energy = np.loadtxt("energy.xvg", unpack=True)

t = energy[0]
U, KE, Etot, T = energy[1:5]

import matplotlib.pyplot as plt

# layout with 3 graphs in 2x2
fig = plt.figure(figsize=(6,5))

axE = fig.add_subplot(221)
axT = fig.add_subplot(223)
axhT = fig.add_subplot(224)

# energies

axE.plot(t, Etot, 'k-', lw=2, label="total")
axE.plot(t, U, 'b-', lw=2, label="potential")
axE.plot(t, KE, 'r-', lw=2, label="kinetic")
axE.set_ylabel("energy (kJ/mol)")
axE.legend(loc="best")

# plot T

# timeseries
axT.plot(t, T, 'r-', lw=2, label="temperature")
axT.set_ylabel("temperature (K)")
axT.set_xlabel("time (ps)")

# histogram

h,edges = np.histogram(T, bins=20, density=True)
midpoints = 0.5*(edges[1:] + edges[:-1])

# plot the histogram as "marginal" over the time series
axhT.plot(h, midpoints, 'r-')
axhT.set_ylim(axT.get_ylim())
axhT.get_xaxis().set_visible(False)    # hide whole axis (counts/density)

# save figure as high resolution PNG
fig.savefig("energies.png", dpi=300)

# mean and fluctuations
print Etot.mean(), Etot.std()

Collect the averages and standard deviations and the plots. Write a short discussion of your results. Comment on the values that you obtained.

Radial distribution function

The radial distribution function g(r) can be calculated in VMD or use g_rdf:

g_rdf -s md.tpr -f md.xtc -o rdf -xvg none

Plot the resulting rdf.xvg file using matplotlib as above.

Self diffusion coefficient

To analyze the self diffusion of Ar, use g_msd:

g_msd -s md.tpr -f md.xtc -o msd -xvg none

Plot:

diffusion = np.loadtxt("msd.xvg", unpack=True)
plot(diffusion[0], diffusion[1], lw=3)

Fitting from 25 to 225 ps

D[        Ar] 2.7806 (+/- 0.2003) 1e-5 cm^2/s

Compare to an experimental literature value of 2.83 x 10^-5 cm²/s @ 94.4 K and 1.374 g/cm³ (Cook GA, Argon, Helium and the Rare Gases. Intersciences: NY, 1961) (via Lee et al, Bull. Korean Chem Soc 24 (2003), 178)

Older experimental values:

• G. Cini-Castagnoli and F. P. Ricci, Self-Diffusion in Liquid Argon, J. Chem. Phys. 32, 19 (1960); http://dx.doi.org/10.1063/1.1700899

  The self-diffusion coefficient of liquid argon has been measured by the capillary technique at 84.56±0.04 °K and at about 70 cm Hg. The average value of D(A37-A40) is

  D = (1.53±0.03)x10^-5 cm²/s.

  70 cm Hg == 93325.4 Pa = 0.9211 atm = 0.9333 bar

  http://www.unitarium.com/pressure

• Naghizade and Rice, J Chem Phys 36 (1962), 2710 at T 90K and 1.374 g/cm³

  D = 2.43x10^-5 cm²/s

Simulation:

• Rahman 1964:

  D = 2.43x10^-5 cm²/s

Extensions

• How do your results change when simulating at T = 120 K or 200 K or 60 K?
• How do the results differ when starting from a fcc system, especially below the melting temperature?