24 January 2012

56. Gromacs -- setting up and running energy minimisation of methanol in water. Using itp files

* First browse through this tutorial.
* 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.
6

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
You 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.

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.

3. Building the methanol.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.

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:
[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

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.

No comments:

Post a Comment