Practical: GROMACS + CP2K Part I

Overview

Teaching: 25 min
Exercises: 50 min

Questions

• What is GROMACS-CP2K QMMM Interface?

• How it could be used?

Objectives

• Getting started with GROMACS-CP2K Interface

• Learning how to prepare your system for a simple QM calculation

• Umbrella sampling using GROMACS-CP2K Interface

The video recording for this practical session can be found here. Note that this is the full recording with solutions to the exercises shown. We strongly advise that you pause the video and try the exercises when they come up.

The slides from the recording can be found here.

Preparing for the tutorial

Everything, which is written inside the gray box are a commands, that should be executed in the terminal window, string-by-string, each following with the ENTER button.
Please note that <...> in the commands means, that everything, including <> symbols, must be replaced with your own specific information. Be careful!

Helpful utilities and commands

Some exercises will require usage of less Linux tool for looking up into the content of files. In case you are not familiar with it, here is a short list of hotkeys, which could be used inside LESS editor:

• q – exit
• / – search for a pattern which will take you to the next occurrence
• n – for next match in forward
• N (SHIFT+n) – for previous match in backward
• g – go to the start of file
• G (SHIFT+g) – go to the end of file

All exercises will require you to submit job for computing using sbatch run.sh command. To check status of your job following commands would be useful.

• squeue -u <your login name> - checks status of all your jobs. Output will look like that:
• scancel <JOBID> will remove the job, if you occasionally submitted it.

Setting up tutorial environment

Let’s start the tutorial with the following steps

1. Execute commands in the terminal (making sure to use the budget/project ID you were assigned when signing up to this course where appropriate):

     module use /work/y07/shared/archer2-lmod/training/
     module load gromacs/2021.1+cp2k
     cd /work/<your project code>/<your project code>/<your login name>
     svn checkout https://github.com/EPCCed/qmmm_gromacs_cp2k_sc/trunk/tutorial
   

2. And go to the tutorial directory
   cd tutorial

Exercise 1: Setting up simple QM system

1) Go to nma directory:
cd nma

In the directory located forcefield and nma.pdb file with a geometry of simple compuond N-Methylacetamide (NMA)
You can download it and inspect structure with VMD or PyMOL.

2) Make topology for the system using the following command:
gmx_cp2k pdb2gmx -f nma.pdb
choose the following forcefield and water model:

Select the Force Field:
From current directory:
1: CHARMM36 all-atom force field (March 2019)
....
....
Select the Water Model:
1: TIP3P       TIP 3-point, recommended, by default uses CHARMM TIP3 with LJ on H

Files topol.top, conf.gro and posre.itp should appear in the directory

3) Look into Gromacs input file em.mdp:
less em.mdp
The following lines contain QMMM MdModule options:

; CP2K QMMM parameters
qmmm-active              = true   ; Activate QMMM MdModule
qmmm-qmgroup             = System ; Index group of QM atoms
qmmm-qmmethod            = PBE    ; Method to use
qmmm-qmcharge            = 0      ; Charge of QM system
qmmm-qmmultiplicity      = 1      ; Multiplicity of QM system

4) Lets perform energy minimization first for that molecule using QMMM interface
Generate Gromacs-CP2K simulation file:
gmx_cp2k grompp -f em.mdp -p topol.top -c conf.gro -o nma-em.tpr
files nma-em.tpr, nma-em.inp and nma-em.pdb should appear in the directory

5) Edit run-em.sh and change [budget code] to your budget code (e.g. ta036-username) and run QMMM simulation: sbatch run-em.sh

6) While job is running you can check the content of nma-em.inp
less nma-em.inp

7) At the end of the job use the following command to extract potential energy:
gmx_cp2k energy
and choose 5 Potential
File energy.xvg should appear in the directory. It contains data with Potential energy (kJ/mol) against optimization step. You can open it in Grace or copy data into any other software (i.e. Excel).

8) Next we will perform short (100 fs) MD simulation with QM. At first look into the md-nvt.mdp file:
less md-nvt.mdp
It contains parameters for performing dynamics with QM forces in the NVT ensemble at 300K

9) Generate Gromacs-CP2K simulation file:
gmx_cp2k grompp -f md-nvt.mdp -p topol.top -c conf.gro -o nma-nvt.tpr
files nma-nvt.tpr, nma-nvt.inp and nma-nvt.pdb should appear in the directory

10) Edit run-nvt.sh and change [budget code] to your budget code (e.g. ta036-username) and run QMMM simulation:
sbatch run-nvt.sh

11) At the end of the simulation you can download trajectory file traj.trr and render it using your favorite software (e.g. VMD, PyMOL).
Also you could check temperature as a function of time with the following command:
gmx_cp2k energy
and choose 9 Temperature
File energy.xvg will contain data about Temperature (K) against simulation time (ps).
Notice, how temperature fluctuates around 300K.

Exercise 2: Umbrella sampling simulation with QMMM MdModule

1) Go to stilbene_vacuum directory:
cd ../stilbene_vacuum

2) Look up in the table and pick-up starting structure and dihedral angle value:

3) Copy chosen starting structure:
cp eq_gro/<your starting gro> ./conf.gro

4) Modify Gromacs input file qmmm_md_umbrella.mdp with value of your chosen Dihedral angle:
sed -i "s/@umbr@/<your dihedral angle>/" qmmm_md_umbrella.mdp
You can also modify pull-coord1-init option in the qmmm_md_umbrella.mdp file with vim or any other editor.

5) Add group of atoms which will be treated with QM to the index file (in that case all atoms are QM):
gmx_cp2k make_ndx -f conf.gro -n index.ndx

> 0
> name 7 QMatoms
> q

6) Generate Gromacs-CP2K simulation file:
gmx_cp2k grompp -f qmmm_md_umbrella.mdp -p topol.top -c conf.gro -n index.ndx -o stilbene.tpr
files stilbene.tpr, stilbene.inp and stilbene.pdb should appear in the directory

7) Edit run.sh and change [budget code] to your budget code (e.g. ta036-username) and run QMMM simulation:
sbatch run.sh
(Note that this might take up to one hour to run).

8) While job is running you can check the content of stilbene.inp
less stilbene.inp
and of qmmm_md_umbrella.mdp
less qmmm_md_umbrella.mdp

9) Check the free energy profiles generated from 100 steps (100 fs) and 10000 steps (10 ps) of QMMM MD simulation for each frame from the eq_gro directory: profile-100fs.xvg and profile-10ps.xvg. You can open them with Grace or copy data into any other software (i.e. Excel).

Key Points

• QM simulation could be activated by adding several parameters into the mdp file

• Most of the simulation techniques from GROMACS are available also within QMMM

• When doing advanced sampling with QMMM one should be aware of the distribution and final profile quality