GROMACS How-tos

author

2019-09-11 18:40 (updated 2022-07-19 05:02)

Tutorials

Introduction to molecular modeling (slides SharePoint folder, video)] - theoretical framework, modeling levels, limitations and possibilities, systems and methods (Erik Lindahl).
Molecular dynamics simulations and GROMACS introduction (slides SharePoint folder, video) - force fields, integrators, control of temperature and pressure (Berk Hess).
Step by step GROMACS Tutorials (Justin A. Lemkul).
MD tutorials in SharePoint folder.
Papers in the group Zotero library, inside folders labeled classics.
Classic texts listed under Computational modeling.

Example simulations

Zeolite Adsorption

Generate the structure file conf.gro.
The system configuration from MC (e.g., topmon output file box1config.xyz) provides a good starting point to initialize GROMACS MD simulations, especially when you need to load a large number of adsorbate molecules into a porous material. You should modify the xyz2gro script to suit your setup. At a minimum, double check variable assignment for nbead, confG, target_atom, residuenameG, atomnameG. You may also need to adapt confH, residuenameH, atomnameH if you simulate a different type of host framework. Then run the script and the last parameter specifies how many adsorbate molecules you want in the conf.gro file to be generated.
```
/path/to/xyz2gro.py /path/to/zeolite.pdb /path/to/box1config.xyz conf.gro NUM_ADSORBATE
```
Set up the topology file topol.top, using #include mechanism to source the parameters from the forcefield folder. TraPPE Forcefield files can be found in /home/common/sample-input/gromacs/trappe.ff. topol.top can be shared among simulations where the same numbers of molecules for each species are used.

Set up the MD parameter file grompp.mdp. For simulating a rigid framework, we do not remove center-of-mass motion, write out only adsorbate coordinates, and freeze framework motion:

comm-mode = None
compressed-x-grps = Water
freezegrps = ZEO
freezedim = Y Y Y

We can then exclude all host-host interactions with group-based pair lists

NOTE: As of GROMACS 2020, the group cutoff scheme support has been removed. If you are using GROMACS 2020 or newer, use Verlet scheme

cutoff-scheme = group
nstlist = 10
nstcalclr = 1
ns_type = grid
rlist = 1.6
energygrps = ZEO
energygrp_excl = ZEO ZEO

However, for GPU simulations, Verlet pair lists have to be used, which does not yet support energy exclusion.

cutoff-scheme = Verlet
verlet-buffer-tolerance = 0.005
coulomb-modifier = potential-shift-Verlet
vdw-modifier = potential-shift-Verlet

Example settings are given below, but you should choose appropriate values for other parameters: dt, nsteps, nstvout, nstlog, nstenergy, nstxout-compressed, Tcoupl, ref_t, Pcoupl, ref_p, gen_temp. grompp.mdp can be shared for the same molecule at same state points (e.g., temperature and pressure).

;dt = 0.0005
;nsteps = 1000000
nstxout = 0
;nstvout = 40000
nstfout = 0
;nstlog = 40000
;nstenergy = 40000
;nstxout-compressed = 4000
compressed-x-precision = 1000
coulombtype = PME
coulomb-modifier = potential-shift
rcoulomb-switch = 1.3
rcoulomb = 1.4
vdw-type = cut-off
vdw-modifier = potential-shift
rvdw-switch = 1.3
rvdw = 1.4
DispCorr = EnerPres
fourierspacing = 0.12
pme_order = 8
ewald_rtol = 1e-06
ewald-rtol-lj = 1e-03
Tcoupl = nose-hoover ;berendsen no v-rescale
nsttcouple = 1
nh-chain-length = 1
tc-grps = System
tau_t = 0.1
;ref_t = 308
Pcoupl = no ;MTTK
Pcoupltype = isotropic
nstpcouple = 1
tau_p = 0.2
compressibility = 4.5e-5
;ref_p = 1.013
gen_vel = yes ;no
;gen_temp = 308
gen_seed = 1993
constraint-algorithm = LINCS
continuation = no ;yes
shake-tol = 1e-06
lincs-order = 8
lincs-iter = 2
lincs-warnangle = 30

If you simulate non-charged adsorbates, then you can use a larger time step and an otherwise unacceptable Coulomb method.

dt = 0.001
coulombtype = cut-off

Generate the binary run input file equil1.tpr
```
/path/to/gmx_mpi_d grompp -f /path/to/grompp.mdp -c conf.gro -p /path/to/topol.top -o equil1.tpr -po equil1.mdp -maxwarn 1
```
-maxwarn 1 is needed because you will get a warning “Can not exclude the lattice Coulomb energy between energy groups”.
Prepare submission scripts. The best practice is to write a generic submission script and put it in the local asset directory (e.g., /path/to/projects/ABC/bin/). Then create symbolic links:
```
ln -sf ../../bin/sub-gromacs.script .  # Run in current directory /path/to/projects/ABC/moleculeX/123K/RUN1
```
NOTE: these symlinks must use relative paths instead of absolute paths, because the local asset directory will be synced to the remote target under your home directory, which may have a different absolute path.

Submit the job using shotgun

shotgun submit -n 1 -c 4 -m 5200 -t 24 -r ghpcc06 -q long -j ../../bin/sub-gromacs-ghpcc.script  # submitting the job

If the job completes successfully, you would want to modify grompp.mdp for production runs.

;nsteps = 20000000
Tcoupl = nose-hoover ;if you have used berendsen during equilibration
gen_vel = no ;yes
continuation = yes ;no

Finally, generate the binary input file prod1.tpr and make sure that state.cpt is removed if you have previously used checkpoint.

rm state.cpt
/path/to/gmx_mpi_d grompp -f /path/to/grompp.mdp -c equil1.tpr -p /path/to/topol.top -o prod1.tpr -po prod1.mdp -t equil1.cpt -maxwarn 1;
shotgun submit -n 1 -c 4 -m 5200 -t 24 -r ghpcc06 -q long -j ../../bin/sub-gromacs-ghpcc.script  # submitting the job

If you job runs out of time or is otherwise killed, you can use the checkpoint file to continue.
```
ln -sf equil1.cpt state.cpt
```
Using the checkpoint will cause new data to be appended to the same set of files equil1.log, equil1.xtc, etc., as if the simulation had not been interrupted.

In this case, we do not need to generate a new .tpr file. Simply submit the job using shotgun:
```
shotgun submit -n 1 -c 4 -m 5200 -t 24 -r ghpcc06 -q long -j ../../bin/sub-gromacs-ghpcc.script  # submitting the job
```

On the other hand, if you wish to continue from where the simulation was stopped, but wrote to a different set of files (e.g., if you want to change parameters), you can do so by generating a new .tpr file.

rm state.cpt  # make sure that there is no leftover checkpoint file
/path/to/gmx_mpi_d grompp -f /path/to/grompp.mdp -c equil1.tpr -p /path/to/topol.top -o equil2.tpr -po equil2.mdp -t equil1.cpt -maxwarn 1;

And submit the job using shotgun:

shotgun submit -n 1 -c 4 -m 5200 -t 24 -r ghpcc06 -q long -j ../../bin/sub-gromacs-ghpcc.script  # submitting the job

Visualize trajectories where some atoms are not written

In the procedure above, we use compressed-x-grps to output coordinates of only sorbate molecules (group Water), the .tpr file will not be able to work with the resulting .xtc files, as it contains references to group ZEO that is not present in the .xtc files. To use the various GROMACS post-processing tools, you need to prepare alternative .mdp, .top, and .gro files. As an example, we can prepare grompp-analysis.mdp by removing references to group ZEO in energygrps, energygrp_excl, freezegrps, and freezedim. Then prepare analysis.top by setting molecules that are not in the .xtc file to zero:

OSiZ 0
SiZ  0

You can do this using sed

sed -E -e "s/SiZ +[0-9]+/SiZ 0/" -e "s/SiOH +[0-9]+/SiOH 0/" topol.top > analysis.top

After that, generate analysis.gro by removing the same set of atoms:

gawk '$2!="ZEO"{print}' conf.gro > analysis.gro
nguest=$(wc -l analysis.gro | cut -d " " -f 1)
((nguest-=3))
sed -i -E "2s/ *[0-9]+/$nguest/" analysis.gro;

Then, generate analysis.tpr to use with Gromacs tools for analysis or trajectory format conversion:

/path/to/gmx_mpi_d grompp -f /path/to/grompp-analysis.mdp -c analysis.gro -p /path/to/analysis.top -o analysis.tpr -po analysis.mdp -maxwarn 1

Tutorials

Example simulations

Zeolite Adsorption

Visualize trajectories where some atoms are not written

Comments