View source for Code:gmxliquid

In this tutorial, I'm assuming you already have the gromacs commands (grompp, mdrun, etc.) loaded.  On circe at USF, this is done with:

 module load apps/gromacs/4.5.5

You'll also need the starting files [[Media:liq_files.tgz|here]].

= Create an input system =

First, we do some size calculations in the python prompt.

 >>> # Compute number of EtOH molecules
 >>> # 8 nm^3 = 8e-21 cc
 >>> # 0.4 % By vol.
 >>> # rho_EtOH = 0.78924 g / cc
 >>> # FW = 46 g/mol
 >>>
 >>> 8e-21*0.4*0.78924/46.0*6.022e23
 33


Next, we need to generate an 8 nm<math>^3</math> box with 33 EtOH molecules

 genbox -box 2 2 2 -ci etoh.pdb -o box.pdb -nmol 33
 nano topol.top # check for line "ETH  33" (no changes required)

If the topology file and forcefields are 'there', we can create a test simulation using just
these molecules.  We'll do this just to check that grompp works

 grompp -c box.pdb -f equ.mdp -o equ.tpr

Next, the complete water+EtOH mixture can be made by filling the voids with water.  We'll use TIP4P here.

 genbox -cs tip4p -cp box.pdb -o sys.pdb

Since this last command added waters, you have to update the topology file, adding "SOL 133" (or your number of waters)

 nano topol.top

= run minimization =

Molecular structures can be touchy, since contacting atoms causes large forces and 'blows up' a system.  The risk of this can be reduced by minimizing the energy of initial starting systems.

First, we have to create a file describing how gromacs should do the minimization.

 grompp -f min.mdp -c sys.pdb -o min

A lot of backup files (starting with '#') accumulate, and we remove them like so:

 rm -f \#* mdout.mdp

Now we're ready to run minimization.  Since the system is small, minimization is cheap and fast - so we do it on the head node.  For larger systems, this command would be put in a job-script.

 mdrun -s min -o min -g min -e min -c min

Now we're ready to take the minimized structure and run dynamics.  This first round is called equilibration,
since we intend to let the system settle into a thermal equilibrium state.

 grompp -c min.gro -f equ.mdp -o equ.tpr

To run this one, we'll use the cluster by writing down the command in a script file.  That script gets sent to the cluster.  No changes to equ.sh should be required.

 nano equ.sh
 qsub -q devel equ.sh
 qstat -u `whoami`

As the job is running you can read through the log file.

 less equ.log
 tail -f equ.log
 # enter Ctl-C to stop tailing

Note that the energy is output at every time-step.  Why is the total energy changing with time?  What about the total volume?

= run dynamics =

Running dynamics uses the same procedure as before, but now we need to worry about how much and which output to produce.

 nano run.mdp # no changes required
 grompp -c equ.gro -f run.mdp -o run.tpr


Check and start the job (no changes required).

 nano run.sh
 qsub -q devel run.sh
 qstat -u `whoami`

Read through run.log

 less run.log
 tail -f run.log
 # enter Ctl-C to stop tailing

Notice that the total energy remains relatively constant now, but the individual energies are still very noisy.  What does this say about the shape of the molecules?

= analyze dynamic data =

Let's first check what happened during equilibration.  Since the volume was allowed to change, we should be able to plot it vs. time:

 g_energy -f equ.edr -o equ.en.xvg
 # Select: Potential Kinetic-En. Total-Energy Volume

* At what time does the potential energy stabilize?
* What is the density?
** answer: (33*46 + 18*133)/Avogadro / 7.0333834916722191e-21 = 0.9236 g/mL
** These two quantities inform on the heat and volume change of mixing.
** Pure H2O (TIP4P): rho = 1.001, H ~ 11.6 kcal/mol, D = 3.9e-5 cm^2/s, eps = 52
** Pure H2O (expt): rho = 0.9971, H = 10.5 kcal/mol, D = 2.3e-5 cm^2/s, eps = 78.4
** (J. Chem. Phys. 123, 234505, 2005)
** Pure EtOH: rho = 0.7873, H = 10.0 kcal/mol, D = 1.0e-5 cm^2/s, eps(expt) = 23
** (J. Phys. Chem. B 1997, 101, 78-86, J. Phys. Chem. B, 2014, 118 (34), 10156-66)

Gnuplot can make this plot with

 plot 'equ.en.xvg' u 1:2 w l, 'equ.en.xvg' u 1:5 axis x1y2 w l


The radial distribution function is a classical measure of the structuring of liquid water.  The one we
calculate here shows the density of molecules within spherical shells around the central molecule.

 g_rdf -rdf mol_com -ng 2 -f run.trr -o run.w-rdf.xvg
 # select Water as "reference group", select "Water", then "ETH"
 # re-run for ETH-ETH rdf
 g_rdf -rdf mol_com -ng 1 -f run.trr -o run.e-rdf.xvg

gnuplot can plot these together with:

 plot 'run.w-rdf.xvg' u 1:2 w l, 'run.w-rdf.xvg' u 1:3 w l, 'run.e-rdf.xvg' u 1:2 w l

Notice how waters stack much more closely together than EtOH?  How is this related to their liquid densities?

The dielectric of a liquid measures its ability to be polarized by an applied field.  High dielectric materials make good capacitors.

  g_dipoles -f run.trr -s run.tpr -c dipcorr.xvg -corr total -P 1

<pre><nowiki>
Dipole moment (Debye)
---------------------
Average  =   2.2234  Std. Dev. =   0.1095  Error =   0.0001

The following averages for the complete trajectory have been calculated:

 Total < M_x > = 7.79882 Debye
 Total < M_y > = -1.11977 Debye
 Total < M_z > = 9.49511 Debye

 Total < M_x^2 > = 882.516 Debye^2
 Total < M_y^2 > = 1239.78 Debye^2
 Total < M_z^2 > = 925.457 Debye^2

 Total < |M|^2 > = 3047.75 Debye^2
 Total |< M >|^2 = 152.233 Debye^2

 < |M|^2 > - |< M >|^2 = 2895.52 Debye^2

Finite system Kirkwood g factor G_k = 3.52836
Infinite system Kirkwood g factor g_k = 2.37948

Epsilon = 43.1739
</nowiki></pre>

Water and EtOH are constantly undergoing random collisions in solution.  The net effect
of this is diffusion of the two molecules.  Einstein showed that this can be tracked
by watching the mean squared displacement of the molecules over time.  Since
the distance traveled over a random walk with diffusion constant <math>D</math>
has a Gaussian distribution with variance <math>2D\Delta t</math>, where <math>\Delta t</math>
is the elapsed time, the mean squared displacement vs time should be a straight line with
slope <math>2D</math>.  Of course, over short times, this picture is only approximate,
but it usually holds over longer times.

 g_msd -f run.trr -ngroup 2 -o run.msd.xvg

<pre><nowiki>
Fitting from 100 to 900 ps

D[     Water] 1.4336 (+/- 0.0850) 1e-5 cm^2/s
D[       ETH] 0.7312 (+/- 0.1279) 1e-5 cm^2/s
</nowiki></pre>

How do these compare with the literature values (above) for the pure liquids?

More detail on velocity correlation:

 g_velacc -mol -f run.trr -o run.w-acf.xvg # Choose group 3 (Water)
 g_velacc -mol -f run.trr -o run.e-acf.xvg # Choose group 2 (ETH)

you can compare these to some of the above references.