* I'm not an expert at GROMACS -- I'm writing this up as I am learning myself. Use what I show here as a starting point for your own explorations. Don't treat it like a scientifically complete piece of work.
1. Building the .gro
Draw a methanol molecule in e.g. Avogadro (available in the debian repos) and save as .xyz.
6You should edit the file to give the atoms unique names -- keep them short though, as each string field is limited to 5 chars in the .gro file. You will refer to these names later in your metahnol.itp file.
C -1.60284 -0.81294 0.00066
O -0.18719 -0.81849 -0.00206
H -1.97028 -1.66300 0.58111
H -1.97029 -0.87387 -1.02686
H -1.95386 0.11552 0.45631
H 0.09473 -1.65393 -0.41207
Here's my new methanol.xyz:
6
meth
COH -1.60284 -0.81294 0.00066
OHC -0.18719 -0.81849 -0.00206
HC1 -1.97028 -1.66300 0.58111
HC2 -1.97029 -0.87387 -1.02686
HC3 -1.95386 0.11552 0.45631
HO 0.09473 -1.65393 -0.41207
Next we convert it to a .gro file and create a 2.5 nm x 2.5 nm x 2.5 nm box around it.
editconf -f methanol.xyz -o methanol.gro -box 2.5 2.5 2.5
methanol.gro:
meth
6
1 ??? COH 1 1.216 1.264 1.257
1 ??? OHC 2 1.358 1.263 1.257
1 ??? HC1 3 1.179 1.179 1.315
1 ??? HC2 4 1.179 1.258 1.154
1 ??? HC3 5 1.181 1.357 1.302
1 ??? HO 6 1.386 1.180 1.216
2.50000 2.50000 2.50000
2. Building the .top
We're going to make it a bit easier and more modular for ourselves by using the #include directive and by (re)using .itp files -- this way we can pull instructions and data from stock files.
step1.top:
[defaults]
1 1 yes 1.0 1.0
#include "atomtypes.itp"
#include "methanol.itp"
#include "water.itp"
[system]
wood liquor
[molecules]
meth 1
SOL 1
It now remains for us to define the molecule meth in methanol.itp, SOL in water.itp and define the proper atomtypes in atomtypes.itp.
First create a .gzmat file from the .xyz using openbabel
babel methanol.xyz methanol.gzmat
methanol.gzmat:
#Put Keywords Here, check Charge and Multiplicity.
meth
0 1
Xx
Xx 1 r2
Xx 1 r3 2 a3
Xx 1 r4 2 a4 3 d4
Xx 1 r5 2 a5 3 d5
Ho 2 r6 1 a6 3 d6
Variables:
r2= 1.4157
r3= 1.0929
a3= 109.52
r4= 1.0929
a4= 109.52
d4= 239.22
r5= 1.0922
a5= 109.00
d5= 119.61
r6= 0.9724
a6= 107.10
d6= 60.39
The zmat format gives an easy overview over bond-lengths and angles, which is what I use it for. For something as small as methanol you don't really need a zmat file, but it can come in handy for large molecules. We'll also include #ifndef, #else and #endif -- this way we have an easy way of switching between using constraints and not doing so. ifndef checks to see whether FLEXIBLE is defined or not (ifndef=If Not Defined).
Remember that the atom order in methanol.itp should be the same as in the .gro file.
Remember that the atom order in methanol.itp should be the same as in the .gro file.
The partial charges are addressed under atomtypes. Be aware that while gzmat outputs bond lengths in Ångström you should use nanometres in your .itp file or things will get weird.
methanol.itp:
;;;;;;;;;;;;;;;;;;
; Methanol
;;;;;;;;;;;;;;;;;;
[ moleculetype ]
meth 3
[atoms]
;nr type resnr res atom cgnr chrg mass
1 COH 1 meth COH 1 0.1450 12.011
2 OHC 1 meth OHC 1 -0.683 15.9994
3 HC 1 meth HC1 1 0.0400 1.00079
4 HC 1 meth HC2 1 0.0400 1.00079
5 HC 1 meth HC3 1 0.0400 1.00079
6 HO 1 meth HO 1 0.4180 1.00079
#ifndef FLEXIBLE
[constraints]
;ai aj f nm
3 1 1 0.110929 ; H-C
4 1 1 0.110929 ; H-C
5 1 1 0.110929 ; H-C
1 2 1 0.14157 ; C-O
6 2 1 0.09724 ; H-O
[exclusions]
1 2 3 4 5 6
2 1 3 4 5 6
3 1 2 4 5 6
4 1 2 3 5 6
5 1 2 3 4 6
6 1 2 3 4 5
#else
[bonds]
;ai aj f nm kb
3 1 1 0.110929 10 ; H-C, r3
4 1 1 0.110929 10 ; H-C, r4
5 1 1 0.110929 10 ; H-C, r5
2 1 1 0.14157 100 ; C-O, r2
6 2 1 0.09724 1 ; H-O, r6
;[ pairs ]
; i j funct
;2 6
;3 6
;4 6
[angles]
;ai aj ak f angle k_a
3 1 2 1 109.52 500 ; H-C-O a3
4 1 2 1 109.52 500 ; H-C-O a4
5 1 2 1 109.52 500 ; H-C-O a5
6 2 1 1 107.10 500 ; H-O-C a6
3 1 4 1 120.00 500
4 1 5 1 120.00 500
3 1 5 1 120.00 500
[dihedrals];proper - f=1
6 2 1 3 1 60 300 2 ;H(O)-O-C-H(C) d6
;6 2 1 4 1
;6 2 1 5 1
[dihedrals];improper - f=2
4 1 2 3 2 239.22 300 ;H-C-O-H(C) d4
5 1 2 3 2 119.61 300 ;H-C-O-H(C) d5
#endif
3. Building the water.itp
The partial charges are addressed under atomtypes, but this is essentially a copy of the (/usr/local/gromacs/share/gromacs/)top/oplsaa.ff/spce.itp file.
water.itp:
4. Building the atomtypes.itp
The partial charges are addressed under atomtypes, but this is essentially a copy of the (/usr/local/gromacs/share/gromacs/)top/oplsaa.ff/spce.itp file.
water.itp:
[moleculetype]
sol 2
[atoms]
;nr type resnr res atom cgnr charge mass
1 OW 1 SOL OW 1 -0.8476 15.9994
2 HW 1 SOL HW1 1 0.4238 1.00079
3 HW 1 SOL HW2 1 0.4238 1.00079
#ifndef FLEXIBLE
[settles]
;OW fun doh dhh
1 1 0.1 0.16330
[exclusions]
1 2 3
2 1 3
3 1 2
#else
[bonds]
;i j f l(nm) kb
1 2 1 0.1 34.5e3
1 3 1 0.1 34.5e3
[angles]
;ai aj ak f angle k_a
2 1 3 1 109.47 383
#endif
4. Building the atomtypes.itp
This is the tricky bit -- you need to pick the correct parameters for your atoms, and as they are empirical it is difficult to know what is 'right'.
A quick way of getting decent starting values is to
cat /usr/local/gromacs/share/gromacs/top/oplsaa.ff/ffnonbonded.itp
and look through the OPLS-AA for different atoms.
e.g.
cat /usr/local/gromacs/share/gromacs/top/oplsaa.ff/ffnonbonded.itp | grep OH
5. Our em.mdp
This is from one of Justin Lemkul's tutorials
6. Energy minimisation
I'll keep it brief:
Add water
genbox_dd -cp methanol.gro -cs /usr/local/gromacs/share/gromacs/top/spc216.gro -o step1.gro -p step1.top
Prepare mdrun input
grompp_dd -f em.mdp -po step2.mdp -c step1.gro -p step1.top -pp step2.top -o step2.tpr
Do energy minimsation
mdrun_dd -v -s step2.tpr -o step3.trr -x step3.xtc -c step3.gro -e step3.edr -g step3.log
Using a run.sh
I get tired of typing the same thing over and over - especially when troubleshooting
Create a run.sh file and call it with sh run.sh every time
7. Equilibration and production runs
See previous tutorial for mdp files.
I know it's a bit of a cop-out, but no time to make a more detailed tutorial right now.
A quick way of getting decent starting values is to
cat /usr/local/gromacs/share/gromacs/top/oplsaa.ff/ffnonbonded.itp
and look through the OPLS-AA for different atoms.
e.g.
cat /usr/local/gromacs/share/gromacs/top/oplsaa.ff/ffnonbonded.itp | grep OH
gives
; name bond_type mass charge ptype sigma epsilon
opls_023 OH 8 15.99940 -0.700 A 3.07000e-01 7.11280e-01 ; SIG
opls_078 OH 8 15.99940 -0.700 A 3.07000e-01 7.11280e-01 ; SIG
opls_154 OH 8 15.99940 -0.683 A 3.12000e-01 7.11280e-01
opls_162 OH 8 15.99940 -0.635 A 3.07000e-01 7.11280e-01
opls_167 OH 8 15.99940 -0.585 A 3.07000e-01 7.11280e-01
opls_169 OH 8 15.99940 -0.700 A 3.07000e-01 7.11280e-01
opls_171 OH 8 15.99940 -0.730 A 3.07000e-01 7.11280e-01
opls_187 OH 8 15.99940 -0.700 A 3.07000e-01 7.11280e-01
opls_268 OH 8 15.99940 -0.530 A 3.00000e-01 7.11280e-01
opls_420 OH 8 15.99940 -0.980 A 3.15000e-01 1.04600e+00
opls_434 OH 8 15.99940 -1.300 A 3.20000e-01 1.04600e+00
For most definitions sigma=3.07000e-01 and epsilon=7.11280e-01, so we'll use that. However, we want to use c6 and c12 notation, where c6=4*epsilon*sigma^6 and c12=4*epsilon*sigma^12.
An alternative is to look through the literature for suitable values.
An advantage with making an atomtypes.itp file yourself is that you can slowly expand it between experiments.
Here's my atomtypes.itp:
[atomtypes]
;atom no mass charge ptype c6 c12
;water
OW 8 15.9994 -0.8476 A 0.261710e-02 0.26331e-5
HW 1 1.00079 0.42380 A 0.000000e-02 0.00000e-5
;methanol
; Methanol, Jorgensen et al. JACS 118 pp. 11225 (1996)
;opls-aa/ffnonbonded.itp
COH 6 12.01100 0.145 A 0.20305e-2 0.37326e-5 ;s/e 3.50000e-01 2.76144e-01
OHC 8 15.99940 -0.683 A 0.26244e-2 0.24208e-5 ;s/e 3.12000e-01 7.11280e-01
HC 1 1.00800 0.040 A 0.012258e-2 0.0029926e-5 ;s/e 2.50000e-01 1.25520e-01
HO 1 1.00800 0.418 A 0.00000000 0.000000000 ;s/e 0.00000e+00 0.00000e+00
;CO2
OC 8 15.9994 -0.35 A 2.5438e-1 2.0478e-4
CO 6 12.0107 0.700 A 5.2044e-2 2.5080e-5
;CO3
OCO2 8 15.9994 -1.0450 A 0.261710e-2 0.26331e-05
CO3 6 12.0107 1.13500 A 0.234020e-2 0.33740e-05
5. Our em.mdp
This is from one of Justin Lemkul's tutorials
; Parameters describing what to do, when to stop and what to save
define = -DFLEXIBLE ; relaxes constraints in the .itp files
integrator = steep ; Algorithm (steep = steepest descent minimization)
emtol = 1000.0 ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep = 0.01 ; Energy step size
nsteps = 50000 ; Maximum number of (minimization) steps to perform
; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist = 1 ; Frequency to update the neighbor list and long range forces
ns_type = grid ; Method to determine neighbor list (simple, grid)
rlist = 1.2 ; Cut-off for making neighbor list (short range forces)
coulombtype = PME ; Treatment of long range electrostatic interactions
rcoulomb = 1.2 ; Short-range electrostatic cut-off
rvdw = 1.2 ; Short-range Van der Waals cut-off
pbc = xyz ; Periodic Boundary Conditions (yes/no)
6. Energy minimisation
I'll keep it brief:
Add water
genbox_dd -cp methanol.gro -cs /usr/local/gromacs/share/gromacs/top/spc216.gro -o step1.gro -p step1.top
Prepare mdrun input
grompp_dd -f em.mdp -po step2.mdp -c step1.gro -p step1.top -pp step2.top -o step2.tpr
Do energy minimsation
mdrun_dd -v -s step2.tpr -o step3.trr -x step3.xtc -c step3.gro -e step3.edr -g step3.log
Using a run.sh
I get tired of typing the same thing over and over - especially when troubleshooting
Create a run.sh file and call it with sh run.sh every time
editconf_dd -f methanol.xyz -o methanol.gro -box 2.5 2.5 2.5
genbox_dd -cp methanol.gro -cs /usr/local/gromacs/share/gromacs/top/spc216.gro -o step1.gro -p step1.top
grompp_dd -f em.mdp -po step2.mdp -c step1.gro -p step1.top -pp step2.top -o step2.tpr
mdrun_dd -v -s step2.tpr -o step3.trr -x step3.xtc -c step3.gro -e step3.edr -g step3.log
7. Equilibration and production runs
See previous tutorial for mdp files.
I know it's a bit of a cop-out, but no time to make a more detailed tutorial right now.
