.. -*- encoding: utf-8 -*- ==================================== 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. :ref:`install Gromacs ` 2. :ref:`run a MD simulation ` 3. :ref:`analyze the output ` You should *record your results in your own Wiki page*. Therefore, go to the `LJ Fluid in Gromacs Wiki`_ on Blackboard (:menuselection:`Course Documents --> Practicals --> Practical 11 --> LJ Fluid in Gromacs Wiki`) and create your own results page in the wiki. .. _Gromacs: http://www.gromacs.org .. _LJ Fluid in Gromacs Wiki: https://myasucourses.asu.edu/webapps/Bb-wiki-bb_bb60/viewWikiLink?course_id=_250406_1&content_id=_7779965_1 Installation ============ All files can be found on the page for `Practical 11`_. .. _Practical 11: http://becksteinlab.physics.asu.edu/pages/courses/2013/SimBioNano/11/ .. _install-gromacs: 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 :envvar:`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 .. _MD-simulation: 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. .. _grompp: http://manual.gromacs.org/online/grompp.html .. _mdrun: http://manual.gromacs.org/online/mdrun.html 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 :command:`mdrun` you can also write a PDB file.) .. _analysis: Analysis ======== *Record your results in your own Wiki page* in the `LJ Fluid in Gromacs Wiki`_ on Blackboard (:menuselection:`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. .. _`g_energy`: http://manual.gromacs.org/current/online/g_energy.html .. _xmgrace: http://plasma-gate.weizmann.ac.il/Grace/ 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. .. _`g_rdf`: http://manual.gromacs.org/current/online/g_rdf.html 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 .. _`g_msd`: http://manual.gromacs.org/current/online/g_msd.html 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 .. |cm2/s| replace:: cm\ :sup:`2`/s .. |g/cm3| replace:: g/cm\ :sup:`3` Compare to an experimental literature value of 2.83 x 10\ :sup:`-5` |cm2/s| @ 94.4 K and 1.374 |g/cm3| (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\ :sup:`-5` |cm2/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/cm3| D = 2.43x10\ :sup:`-5` |cm2/s| Simulation: * Rahman 1964: D = 2.43x10\ :sup:`-5` |cm2/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?