Quick and dirty GROMACS how-to
This is all executed on sierra.
First. the PDB file must be "cleaned" for use in Gromacs:
pdb2gmx -ignh -f 1UXC.pdb -o uxc.pdb -p uxc.top
Here, the RCSB pdb is 1UXC.pdb, and the output file I named uxc.pdb. The parameter file created is specified with "-p". You may also specify a force-field with the "-ff" switch, but for me, it's easier to pick from the list of available force-fields. In this case, I selected 13 / GROMOS96 53a6. pdb2gmx will also ask you about a water type. Gromacs utilizes SPC water, so I selected 1 / SPC water.
In order to create a solvent box around the protein, two steps must be made. First, the box is set up with "editconf":
editconf -bt cubic -f uxc.pdb -o uxc.pdb -c -d 1.5
Then the actual box is generated. Notice that a new parameter file is also created/overwritten:
genbox -cp uxc.pdb -cs spc216.gro -o uxc_b4em.pdb -p uxc.top
"-cp uxc.pdb" specifies the input file, "-cs spc216.gro" specifies the library file containing the SPC water, and "-o uxc_b4em.pdb" specifies the output file.
Now you will pre-process the files using grompp. You must create a input file "em.mdp" for the simulation/minimization settings. Here is an example:
title = UXC
cpp = /usr/bin/cpp
define = -DFLEXIBLE
constraints = none
integrator = steep
dt = 0.002 ; ps
nsteps = 500
nstlist = 5
ns_type = grid
rlist = 0.9
coulombtype = PME
rcoulomb = 0.9
rvdw = 1.4
fourierspacing = 0.12
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
pme_order = 4
ewald_rtol = 1e-5
optimize_fft = yes
;
; Energy minimizing
;
emtol = 1000
emstep = 0.01
Run grompp ("-c file" specifies the input PDB file):
grompp -f em.mdp -c uxc_b4em.pdb -p uxc.top -o uxc_em.tpr
Pay attention to the output, as it will give you useful information. For example, if the system has net charge, it must be neutralized by adding chloride and/or sodium ions. This involves adding the ions with genion, and modifying the topology file to adjust counts of solvent molecules removed, and ions added. Here the net charge is +4, so we are going to add four chlorides:
genion –s uxc_em.tpr –o uxc_ion.pdb –nname CL- –nn 4 –g uxc_ion.log
Make sure to use the SOL group for replacements!
Examine the output, and make note of how many ions were added/how many solvent molecules were replaced. Edit your topology file (uxc.top) and find the [molecules] section. Make the necessary adjustments in the counts, and add a CL entry (here we see it after the adjustment):
[ molecules ]
; Compound #mols
Protein_chain_A 1
SOL 7829
CL 4
Now the pre-processor must be run again to reflect these changes. (Also remember to re-run grompp if you change things in xx.mdp!)
grompp -f em.mdp -c uxc_ion.pdb -p uxc.top -o uxc_em.tpr
The minimization is then performed. If it doesn't converge, add steps in em.mdp, re-pre-process and run again.
mdrun –s uxc_em.tpr –o uxc_em.trr –c uxc_b4pr.pdb –g em.log –e em.edr &
"–c uxc_b4pr.pdb" specifies the output from this step.
Now we may run a fixed-coordinates simulation, relaxing the water (this is also known as "soaking" the water into the protein). The protein coordinates are fixed in this instance. First, a control file is written, very similar to the energy-minimzing file, above. Here is an example for a restrained soaking (pr.mdp): Comments are preceded by a semicolon.
title = UXC
cpp = /usr/bin/cpp
define = -DPOSRES
constraints = all-bonds
integrator = md
dt = 0.002;ps!
nsteps = 10000 ; total 20.0 ps.
nstcomm = 1
nstxout = 250 ; collect data every 0.5 ps
nstvout = 1000
nstfout = 0
nstlog = 10
nstenergy = 10
nstlist = 5
ns_type = grid
rlist = 0.9
coulombtype = PME
rcoulomb = 0.9
rvdw = 1.4
fourierspacing = 0.12
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
pme_order = 4
ewald_rtol = 1e-5
optimize_fft = yes
; Berendsen temperature coupling is on in two groups
Tcoupl = berendsen
tau_t = 0.1 0.1 0.1 ; last entry if you have Cl ions
tc-grps = protein sol CL
ref_t = 300 300 300
; Pressure coupling is on
Pcoupl = berendsen
tau_p = 0.5
compressibility = 4.5e-5
ref_p = 1.0
; Generate velocites is on at 300 K.
gen_vel = yes
gen_temp = 300.0
gen_seed = 173529
In grompp, "-r uxc_b4pr.pdb" specifies the restraint file. In this case, since the whole PDB file is specified, and because of the "constraints = all-bonds" entry in pr.mdp, the protein is frozen, while the waters are allowed to relax.
grompp -f pr.mdp -c uxc_b4pr.pdb -r uxc_b4pr.pdb -p uxc.top -o uxc_pr.tpr
mdrun –s uxc_pr.tpr –o uxc_pr.trr –c uxc_b4md.pdb –g pr.log –e pr.edr &
Here is an md.mdp for an explictily-solvated system. Note that all bonds are still constrained. Is this realistic? You may wish to experiment with what constraints (if any) are necessary. These links to mdp and mdrun/bonds might be useful.
title = UXC
cpp = /usr/bin/cpp
constraints = all-bonds
integrator =md
dt = 0.002 ;ps!
nsteps = 25000 ; total 50 ps.
nstcomm =1
nstxout = 500 ; collect data every 1 ps nstvout =0
nstfout =0
nstlist =5
ns_type = grid
rlist = 0.9
coulombtype = PME
rcoulomb = 0.9
rvdw = 1.4
fourierspacing = 0.12
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
pme_order = 4
ewald_rtol = 1e-5
optimize_fft = yes
; Berendsen temperature coupling is on in two groups Tcoupl = berendsen
tau_t =0.1 0.1 0.1
tc-grps = protein sol CL-
ref_t = 300 300 300
; Pressure coupling is on
Pcoupl = berendsen
tau_p = 0.5
compressibility = 4.5e-5
ref_p = 1.0
; Generate velocites is on at 300 K.
gen_vel = yes
gen_temp = 300.0
gen_seed = 173529
grompp –f md.mdp –c uxc_b4md.pdb –r uxc_b4md.pdb –p uxc.top –o uxc_md.tpr
mdrun –s uxc_md.tpr –o uxc_md.trr –c uxc_pmd.pdb –g md.log –e md.edr &
Once you have run the actual MD, you may compress the trajectory using trjconv to save on disk space.
trjconv –f filename.trr –o filename.xtc
Once you have made the *.xtc file, you may delete the *.trr file.
This is all executed on sierra.
First. the PDB file must be "cleaned" for use in Gromacs:
pdb2gmx -ignh -f 1UXC.pdb -o uxc.pdb -p uxc.top
Here, the RCSB pdb is 1UXC.pdb, and the output file I named uxc.pdb. The parameter file created is specified with "-p". You may also specify a force-field with the "-ff" switch, but for me, it's easier to pick from the list of available force-fields. In this case, I selected 13 / GROMOS96 53a6. pdb2gmx will also ask you about a water type. Gromacs utilizes SPC water, so I selected 1 / SPC water.
In order to create a solvent box around the protein, two steps must be made. First, the box is set up with "editconf":
editconf -bt cubic -f uxc.pdb -o uxc.pdb -c -d 1.5
Then the actual box is generated. Notice that a new parameter file is also created/overwritten:
genbox -cp uxc.pdb -cs spc216.gro -o uxc_b4em.pdb -p uxc.top
"-cp uxc.pdb" specifies the input file, "-cs spc216.gro" specifies the library file containing the SPC water, and "-o uxc_b4em.pdb" specifies the output file.
Now you will pre-process the files using grompp. You must create a input file "em.mdp" for the simulation/minimization settings. Here is an example:
title = UXC
cpp = /usr/bin/cpp
define = -DFLEXIBLE
constraints = none
integrator = steep
dt = 0.002 ; ps
nsteps = 500
nstlist = 5
ns_type = grid
rlist = 0.9
coulombtype = PME
rcoulomb = 0.9
rvdw = 1.4
fourierspacing = 0.12
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
pme_order = 4
ewald_rtol = 1e-5
optimize_fft = yes
;
; Energy minimizing
;
emtol = 1000
emstep = 0.01
Run grompp ("-c file" specifies the input PDB file):
grompp -f em.mdp -c uxc_b4em.pdb -p uxc.top -o uxc_em.tpr
Pay attention to the output, as it will give you useful information. For example, if the system has net charge, it must be neutralized by adding chloride and/or sodium ions. This involves adding the ions with genion, and modifying the topology file to adjust counts of solvent molecules removed, and ions added. Here the net charge is +4, so we are going to add four chlorides:
genion –s uxc_em.tpr –o uxc_ion.pdb –nname CL- –nn 4 –g uxc_ion.log
Make sure to use the SOL group for replacements!
Examine the output, and make note of how many ions were added/how many solvent molecules were replaced. Edit your topology file (uxc.top) and find the [molecules] section. Make the necessary adjustments in the counts, and add a CL entry (here we see it after the adjustment):
[ molecules ]
; Compound #mols
Protein_chain_A 1
SOL 7829
CL 4
Now the pre-processor must be run again to reflect these changes. (Also remember to re-run grompp if you change things in xx.mdp!)
grompp -f em.mdp -c uxc_ion.pdb -p uxc.top -o uxc_em.tpr
The minimization is then performed. If it doesn't converge, add steps in em.mdp, re-pre-process and run again.
mdrun –s uxc_em.tpr –o uxc_em.trr –c uxc_b4pr.pdb –g em.log –e em.edr &
"–c uxc_b4pr.pdb" specifies the output from this step.
Now we may run a fixed-coordinates simulation, relaxing the water (this is also known as "soaking" the water into the protein). The protein coordinates are fixed in this instance. First, a control file is written, very similar to the energy-minimzing file, above. Here is an example for a restrained soaking (pr.mdp): Comments are preceded by a semicolon.
title = UXC
cpp = /usr/bin/cpp
define = -DPOSRES
constraints = all-bonds
integrator = md
dt = 0.002;ps!
nsteps = 10000 ; total 20.0 ps.
nstcomm = 1
nstxout = 250 ; collect data every 0.5 ps
nstvout = 1000
nstfout = 0
nstlog = 10
nstenergy = 10
nstlist = 5
ns_type = grid
rlist = 0.9
coulombtype = PME
rcoulomb = 0.9
rvdw = 1.4
fourierspacing = 0.12
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
pme_order = 4
ewald_rtol = 1e-5
optimize_fft = yes
; Berendsen temperature coupling is on in two groups
Tcoupl = berendsen
tau_t = 0.1 0.1 0.1 ; last entry if you have Cl ions
tc-grps = protein sol CL
ref_t = 300 300 300
; Pressure coupling is on
Pcoupl = berendsen
tau_p = 0.5
compressibility = 4.5e-5
ref_p = 1.0
; Generate velocites is on at 300 K.
gen_vel = yes
gen_temp = 300.0
gen_seed = 173529
In grompp, "-r uxc_b4pr.pdb" specifies the restraint file. In this case, since the whole PDB file is specified, and because of the "constraints = all-bonds" entry in pr.mdp, the protein is frozen, while the waters are allowed to relax.
grompp -f pr.mdp -c uxc_b4pr.pdb -r uxc_b4pr.pdb -p uxc.top -o uxc_pr.tpr
mdrun –s uxc_pr.tpr –o uxc_pr.trr –c uxc_b4md.pdb –g pr.log –e pr.edr &
Here is an md.mdp for an explictily-solvated system. Note that all bonds are still constrained. Is this realistic? You may wish to experiment with what constraints (if any) are necessary. These links to mdp and mdrun/bonds might be useful.
title = UXC
cpp = /usr/bin/cpp
constraints = all-bonds
integrator =md
dt = 0.002 ;ps!
nsteps = 25000 ; total 50 ps.
nstcomm =1
nstxout = 500 ; collect data every 1 ps nstvout =0
nstfout =0
nstlist =5
ns_type = grid
rlist = 0.9
coulombtype = PME
rcoulomb = 0.9
rvdw = 1.4
fourierspacing = 0.12
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
pme_order = 4
ewald_rtol = 1e-5
optimize_fft = yes
; Berendsen temperature coupling is on in two groups Tcoupl = berendsen
tau_t =0.1 0.1 0.1
tc-grps = protein sol CL-
ref_t = 300 300 300
; Pressure coupling is on
Pcoupl = berendsen
tau_p = 0.5
compressibility = 4.5e-5
ref_p = 1.0
; Generate velocites is on at 300 K.
gen_vel = yes
gen_temp = 300.0
gen_seed = 173529
grompp –f md.mdp –c uxc_b4md.pdb –r uxc_b4md.pdb –p uxc.top –o uxc_md.tpr
mdrun –s uxc_md.tpr –o uxc_md.trr –c uxc_pmd.pdb –g md.log –e md.edr &
Once you have run the actual MD, you may compress the trajectory using trjconv to save on disk space.
trjconv –f filename.trr –o filename.xtc
Once you have made the *.xtc file, you may delete the *.trr file.