Difference between revisions of "Code:gmxliquid"

From Predictive Chemistry
Jump to: navigation, search
(analyze dynamic data)
 
(3 intermediate revisions by the same user not shown)
Line 1: Line 1:
In this tutorial, I'm assuming you already have the gromacs commands (grompp, mdrun, etc.) loaded. On circe at USF, this is done with:
+
In this tutorial, I'm assuming you already have the gromacs commands (grompp, mdrun, etc.) loaded.
+
I'm also assuming you have gromacs 5, which made significant changes to the way
module load apps/gromacs/4.5.5
+
command-line tools are used.
   
 
You'll also need the starting files [[Media:liq_files.tgz|here]].
 
You'll also need the starting files [[Media:liq_files.tgz|here]].
  +
Or, on circe, just run
  +
  +
cp -R /shares/mri_workshop/Liquids .
  +
cd Liquids
   
 
= Create an input system =
 
= Create an input system =
Line 21: Line 25:
 
Next, we need to generate an 8 nm<math>^3</math> box with 33 EtOH molecules
 
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
+
gmx insert-molecules -box 2 2 2 -ci etoh.pdb -o box.pdb -nmol 33
 
nano topol.top # check for line "ETH 33" (no changes required)
 
nano topol.top # check for line "ETH 33" (no changes required)
   
Line 27: Line 31:
 
these molecules. We'll do this just to check that grompp works
 
these molecules. We'll do this just to check that grompp works
   
grompp -c box.pdb -f equ.mdp -o equ.tpr
+
gmx grompp -c box.pdb -f equ.mdp -o equ.tpr
   
  +
If this completes without errors, we know that we ''could'' simulate the ETH molecules we have now.
  +
That would give a high-density gas, rather than the liquid we want.
 
Next, the complete water+EtOH mixture can be made by filling the voids with water. We'll use TIP4P here.
 
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
+
gmx solvate -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)
+
Since this last command added waters, you have to update the topology file, adding "SOL 130" (or your number of waters)
   
 
nano topol.top
 
nano topol.top
Line 43: Line 49:
 
First, we have to create a file describing how gromacs should do the minimization.
 
First, we have to create a file describing how gromacs should do the minimization.
   
grompp -f min.mdp -c sys.pdb -o min
+
gmx grompp -p topol.top -f min.mdp -c sys.pdb -o min
   
 
A lot of backup files (starting with '#') accumulate, and we remove them like so:
 
A lot of backup files (starting with '#') accumulate, and we remove them like so:
Line 51: Line 57:
 
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.
 
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
+
gmx mdrun -deffnm min
  +
  +
The latest gromacs has a bug of some sort or can't work with OPLS-AA for EtOH, since
  +
I get a segfault here. Hopefully this will be fixed soon. Otherwise, you
  +
can manually delete the EtOH from sys.pdb and topol.top and
  +
then re-run the grompp and mdrun steps.
   
 
Now we're ready to take the minimized structure and run dynamics. This first round is called equilibration,
 
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.
 
since we intend to let the system settle into a thermal equilibrium state.
   
grompp -c min.gro -f equ.mdp -o equ.tpr
+
gmx 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.
 
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
 
nano equ.sh
qsub -q devel equ.sh
+
sbatch equ.sh
qstat -u `whoami`
+
squeue -u `whoami`
   
 
As the job is running you can read through the log file.
 
As the job is running you can read through the log file.
Line 77: Line 83:
   
 
nano run.mdp # no changes required
 
nano run.mdp # no changes required
grompp -c equ.gro -f run.mdp -o run.tpr
+
gmx grompp -c equ.gro -f run.mdp -o run.tpr
   
   
Line 83: Line 89:
   
 
nano run.sh
 
nano run.sh
qsub -q devel run.sh
+
sbatch run.sh
qstat -u `whoami`
+
squeue -u `whoami`
   
 
Read through run.log
 
Read through run.log
Line 98: Line 104:
 
Let's first check what happened during equilibration. Since the volume was allowed to change, we should be able to plot it vs. time:
 
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
+
gmx energy -f equ.edr -o equ.en.xvg
 
# Select: Potential Kinetic-En. Total-Energy Volume
 
# Select: Potential Kinetic-En. Total-Energy Volume
   
Line 121: Line 127:
 
calculate here shows the density of molecules within spherical shells around the central molecule.
 
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
+
gmx rdf -rdf mol_com -ng 2 -s run -f run.trr -o run.w-rdf.xvg
 
# select Water as "reference group", select "Water", then "ETH"
 
# select Water as "reference group", select "Water", then "ETH"
 
# re-run for ETH-ETH rdf
 
# re-run for ETH-ETH rdf
g_rdf -rdf mol_com -ng 1 -f run.trr -o run.e-rdf.xvg
+
gmx rdf -rdf mol_com -ng 1 -s run -f run.trr -o run.e-rdf.xvg
   
 
gnuplot can plot these together with:
 
gnuplot can plot these together with:
Line 134: Line 140:
 
The dielectric of a liquid measures its ability to be polarized by an applied field. High dielectric materials make good capacitors.
 
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
+
gmx dipoles -f run.trr -s run.tpr -c dipcorr.xvg -corr total -P 1
   
 
<pre><nowiki>
 
<pre><nowiki>
Line 171: Line 177:
 
but it usually holds over longer times.
 
but it usually holds over longer times.
   
g_msd -f run.trr -ngroup 2 -o run.msd.xvg
+
gmx msd -f run.trr -ngroup 2 -o run.msd.xvg
   
 
<pre><nowiki>
 
<pre><nowiki>
Line 184: Line 190:
 
More detail on velocity correlation:
 
More detail on velocity correlation:
   
g_velacc -mol -f run.trr -o run.w-acf.xvg # Choose group 3 (Water)
+
gmx velacc -mol -s run -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)
+
gmx velacc -mol -s run -f run.trr -o run.e-acf.xvg # Choose group 2 (ETH)
   
 
you can compare these to some of the above references.
 
you can compare these to some of the above references.
  +
  +
  +
== Appendix ==
  +
  +
Here are cut and paste versions of the 3 most important files for running this tutorial.
  +
  +
TITLE An empty 20 Ang. cube simulation box (box.pdb)
  +
CRYST1 20.000 20.000 20.000 90.00 90.00 90.00 P 1 1
  +
MODEL 1
  +
ENDMDL
  +
  +
A topology using the opls-aa forcefield (topol.top).
  +
  +
#include "oplsaa.ff/forcefield.itp"
  +
#include "oplsaa.ff/ethanol.itp"
  +
#include "oplsaa.ff/tip4p.itp"
  +
  +
[ system ]
  +
80 Proof Water
  +
  +
[ molecules ]
  +
ETH 33
  +
  +
An ethanol molecule from the gromacs shares directory (etoh.pdb):
  +
  +
HETATM 1 C ETH 1 -0.000 -0.000 0.000 0.00 0.00 C
  +
HETATM 2 H ETH 1 -0.785 0.244 -0.653 0.00 0.00 H
  +
HETATM 3 H ETH 1 0.322 -0.981 -0.190 0.00 0.00 H
  +
HETATM 4 H ETH 1 -0.335 0.073 0.993 0.00 0.00 H
  +
HETATM 5 C ETH 1 1.171 0.975 -0.220 0.00 0.00 C
  +
HETATM 6 H ETH 1 2.014 0.675 0.402 0.00 0.00 H
  +
HETATM 7 H ETH 1 1.468 0.957 -1.268 0.00 0.00 H
  +
HETATM 8 OA ETH 1 0.763 2.299 0.138 0.00 0.00 O
  +
HETATM 9 HO ETH 1 0.490 2.313 1.058 0.00 0.00 H
  +
  +
An mdp file for a short equilibration run (equ.mdp):
  +
  +
constraints = none
  +
continuation = no
  +
gen_vel = yes
  +
gen_temp = 300
  +
gen_seed = 9875945
  +
  +
  +
integrator = md
  +
  +
tinit = 0
  +
dt = 0.001
  +
nsteps = 100000 ; 100 ps
  +
nstcomm = 500
  +
  +
nstlog = 5000
  +
nstenergy = 500
  +
nstcalcenergy = 500
  +
nstxout = 5000
  +
nstvout = 5000
  +
nstfout = 5000
  +
  +
comm-mode = Linear
  +
  +
emtol = 100
  +
emstep = 0.001
  +
  +
nstlist = 10
  +
ns_type = grid
  +
pbc = xyz
  +
rlist = 0.9
  +
  +
coulombtype = pme
  +
rcoulomb = 0.9
  +
vdw-type = Cut-Off
  +
rvdw = 0.9
  +
  +
DispCorr = EnerPres
  +
fourierspacing = 0.1
  +
  +
pme_order = 4
  +
ewald_rtol = 1e-5
  +
ewald_geometry = 3d
  +
optimize_fft = yes
  +
  +
Tcoupl = v-rescale
  +
tc-grps = System
  +
tau_t = 1.0
  +
ref_t = 300
  +
  +
Pcoupl = Berendsen
  +
; Pcoupl = Parrinello-Rahman
  +
Pcoupltype = isotropic
  +
tau_p = 1.8
  +
compressibility = 4.5e-5
  +
ref_p = 1.01325
  +
  +
  +
* Install instructions for [[Code:libxdrfile]]

Latest revision as of 13:51, 17 April 2017

In this tutorial, I'm assuming you already have the gromacs commands (grompp, mdrun, etc.) loaded. I'm also assuming you have gromacs 5, which made significant changes to the way command-line tools are used.

You'll also need the starting files here. Or, on circe, just run

cp -R /shares/mri_workshop/Liquids .
cd Liquids

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

gmx insert-molecules -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

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

If this completes without errors, we know that we could simulate the ETH molecules we have now. That would give a high-density gas, rather than the liquid we want. Next, the complete water+EtOH mixture can be made by filling the voids with water. We'll use TIP4P here.

gmx solvate -cs tip4p -cp box.pdb -o sys.pdb

Since this last command added waters, you have to update the topology file, adding "SOL 130" (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.

gmx grompp -p topol.top -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.

gmx mdrun -deffnm min

The latest gromacs has a bug of some sort or can't work with OPLS-AA for EtOH, since I get a segfault here. Hopefully this will be fixed soon. Otherwise, you can manually delete the EtOH from sys.pdb and topol.top and then re-run the grompp and mdrun steps.

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.

gmx 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
sbatch equ.sh
squeue -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
gmx grompp -c equ.gro -f run.mdp -o run.tpr


Check and start the job (no changes required).

nano run.sh
sbatch run.sh
squeue -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:

gmx 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
    • Gas-phase TIP4P energy = 0 (no self interactions in the model)
    • Gas-phase OPLSAA-EtOH energy = 25 kJ/mol (<math>\pm</math>1)
    • (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.

gmx rdf -rdf mol_com -ng 2 -s run -f run.trr -o run.w-rdf.xvg
# select Water as "reference group", select "Water", then "ETH"
# re-run for ETH-ETH rdf
gmx rdf -rdf mol_com -ng 1 -s run -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.

 gmx dipoles -f run.trr -s run.tpr -c dipcorr.xvg -corr total -P 1
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

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.

gmx msd -f run.trr -ngroup 2 -o run.msd.xvg
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

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

More detail on velocity correlation:

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

you can compare these to some of the above references.


Appendix

Here are cut and paste versions of the 3 most important files for running this tutorial.

TITLE     An empty 20 Ang. cube simulation box (box.pdb)
CRYST1   20.000   20.000   20.000  90.00  90.00  90.00 P 1           1
MODEL        1
ENDMDL

A topology using the opls-aa forcefield (topol.top).

#include "oplsaa.ff/forcefield.itp"
#include "oplsaa.ff/ethanol.itp"
#include "oplsaa.ff/tip4p.itp"

[ system ]
80 Proof Water

[ molecules ]
ETH             33

An ethanol molecule from the gromacs shares directory (etoh.pdb):

HETATM    1  C   ETH     1      -0.000  -0.000   0.000  0.00  0.00           C
HETATM    2  H   ETH     1      -0.785   0.244  -0.653  0.00  0.00           H
HETATM    3  H   ETH     1       0.322  -0.981  -0.190  0.00  0.00           H
HETATM    4  H   ETH     1      -0.335   0.073   0.993  0.00  0.00           H
HETATM    5  C   ETH     1       1.171   0.975  -0.220  0.00  0.00           C
HETATM    6  H   ETH     1       2.014   0.675   0.402  0.00  0.00           H
HETATM    7  H   ETH     1       1.468   0.957  -1.268  0.00  0.00           H
HETATM    8  OA  ETH     1       0.763   2.299   0.138  0.00  0.00           O
HETATM    9  HO  ETH     1       0.490   2.313   1.058  0.00  0.00           H

An mdp file for a short equilibration run (equ.mdp):

constraints              = none
continuation             = no
gen_vel                  = yes
gen_temp                 = 300
gen_seed                 = 9875945


integrator               = md

tinit                    = 0
dt                       = 0.001
nsteps                   = 100000 ; 100 ps
nstcomm                  = 500

nstlog                   = 5000
nstenergy                = 500
nstcalcenergy		 = 500
nstxout                  = 5000
nstvout                  = 5000
nstfout                  = 5000

comm-mode                = Linear

emtol                    = 100
emstep                   = 0.001

nstlist                  = 10
ns_type                  = grid
pbc                      = xyz
rlist                    = 0.9

coulombtype              = pme
rcoulomb                 = 0.9
vdw-type                 = Cut-Off
rvdw                     = 0.9

DispCorr                 = EnerPres
fourierspacing           = 0.1

pme_order                = 4
ewald_rtol               = 1e-5
ewald_geometry           = 3d
optimize_fft             = yes

Tcoupl                   = v-rescale
tc-grps                  = System
tau_t                    = 1.0
ref_t                    = 300

Pcoupl                   = Berendsen
; Pcoupl                   = Parrinello-Rahman
Pcoupltype               = isotropic
tau_p                    = 1.8
compressibility          = 4.5e-5
ref_p                    = 1.01325