MD MSM Protocol

Parametrization of organic molecules

The parametrization is not required for peptides, for cyclic peptides a workaround is needed to overcome the cyclization problem.

The parametrization is carried out on the minimum energy conformer (MEC) resulting from Omega conformational sampling after the minimization with Szybki.

WARNING: Omega changes the atoms order, thus it is MANDATORY to use the files obtained by the Szybki minimization of Omega output to avoid problems when the MD simulation protocol is used to study all conformations. In principle, the output of Omega is always the same and it does not depend on the epsilon used to simulate the environment.

The parametrization is carried out with acpype tool.

acpype on site

Set up of the environment

source /home/germondi/gromacs_setup.source

In the setup file the commands are the following

# GROMACS environment export PATH=/opt/gromacs/bin:$PATH

# VMD environment (from BiKi)
export PATH=/opt/biki/BiKiLifeSciences/vmd-1.9.2/vmd_bin/:$PATH

# Babel setup (installation in user directory)
export PATH=/opt/openbabel/openbabel-openbabel-2-4-0/bin:$PATH
export BABEL_LIBDIR=/opt/openbabel/openbabel-openbabel-2-4-0/lib
export BABEL_DATADIR=/opt/openbabel/openbabel-openbabel-2-4-0/data

# acpype (and Antechamber) environment
export PATH=/home/germondi/acpype/scripts:$PATH
export PATH=/home/germondi/acpype/amber19-0_linux/bin:$PATH
export PATH=/home/germondi/acpype/amber19-0_linux/dat:$PATH
export PATH=/home/germondi/acpype/amber19-0_linux/lib:$PATH

Copy the file with the MEC coordinates (mol_solv_final_mec.mol2) in the working directory and then submit it to acpype

acpype.py -i mol_solv_final_mec.mol2 -b ach -n molcharge

Where molcharge is the charge of the molecule.

The output is saved in a directory mol_solv_final_mec.acpype.

Files that we are looking for are the files with the GAFF parametrization:

mol_solv_final_GMX.itp
mol_solf_final_GMX.top

acpype on-line

Open an Ineternet browser and then go to the Bio2Bite web site (a registration is required)

Click on the submit link, login with your credentials and then load mol_solv_final_mec.mol2

Select the correct charge and submit.

You will receive a mail when the parametrization is finished.

Download and unzip the file with the results, a directory acpype is created with all the needed files, in particular:

mol_solv_final_GMX.itp
mol_solf_final_GMX.top

Topology test

Firstly, it is possible to test the topology files using the em.mdp file included in the acpype output following the instruction reported in em.mdp.

A second test is performed in the directory where the simulations will be launched.

The new tests will be launched directly in the topology directory top. Thus, copy the following files in top:

mol_solv_final_GMX.itp
mol_solf_final_GMX.top
if the simulation is run in chloroform

CHL_Box.pdb
CHL.itp

the file mol2 that will be used in the next simulations (in principle the same used by acpype or a file obtained by the conformational analysis

Create a file em.mdp to test the topology with a minimization. The file has to contain the following instructions:

; Parameters describing what to do, when to stop and what to save
integrator = steep         ; Algorithm (steep = steepest descent minimization)
emtol       = 1000.0        ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep      = 0.01          ; Minimization 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
cutoff-scheme   = Verlet    ; Buffered neighbor searching
ns_type         = grid      ; Method to determine neighbor list (simple, grid)
coulombtype     = PME       ; Treatment of long range electrostatic interactions
rcoulomb        = 1.0       ; Short-range electrostatic cut-off
rvdw            = 1.0       ; Short-range Van der Waals cut-off
pbc             = xyz       ; Periodic Boundary Conditions in all 3 dimensions

Now we test the topology using the original mol2 file to be sure that there will be no problems when other conformations will be used

obabel -imol2 mol_solv_final_mec.mol2 -opdb -O mol_solv_final_mec.pdb
gmx editconf -f mol_solv_final_mec.pdb -o mol_test.gro -bt cubic -d 0.5

Now run the minimization with Gromacs

gmx grompp -f em.mdp -c mol_test.gro -p ach_wat_fin_mec_GMX.top -o em.tpr -v
gmx mdrun -ntmpi 1 -v -deffnm em

Problems that could be encountered during grompp run:

1) Non-matching atom names

Message:

Fatal error:
Too many warnings (1), gmx terminated.
WARNING 1 [file arv_GMX.top, line 15]:
1 non-matching atom name
atom names from arv_GMX.top will be used
atom names from arv852_chl_01.gro will be ignored

Solution:

Check the itp file and correct atom names. In particular, some atoms, i.e. chlorine, are reported in capital letters in the pdb file, i.e. CL, whereas in the itp file not, i.e. Cl. In these cases, the CL has to be changed in Cl, check the definition in the pdb file. This solution should guarantee compatibility with all the conformers.

2) If grompp exits with an error about a cut-off, the box dimensions has to be increased until the procedure finishes without errors (i.e. change -d 0.5 to -d 1.0 or larger values)

Check em.gro with VMD to be sure that all is ok

vmd em.gro

From the VMD console:

pbc box_draw

Check the dimensions of the box to be sure that the molecule is well positioned.

If all is ok, move to the next step.

MD setup in a box of water

Create a new directory mol_wat_setup and copy herein the following files from the test directory

mol_solv_final_GMX.itp
mol_solf_final_GMX.top
mol_test.gro

the gro file could be obtained using a file resulting from a conformational sampling following the instruction in Topology test section

Solvate the molecule with water

gmx solvate -cp ach_test.gro -o ach_test_wat.gro -p ach_wat_fin_mec_GMX.top

Check the result with VMD

The topology (.top) file is updated with the addition of a new row in [ molecules ] section

Edit the topology file

delete [ defaults ] section
insert itp information about deafult ff, water model (tip3p) and ions (even if it is necessary only for charged molecules)

; ach_wat_fin_mec_GMX.top created by acpype (v: 2019-07-10T18:04:00UTC) on Fri Apr 10 12:05:13 2020

; Include default ff
#include "amber99sb.ff/forcefield.itp"

; Include mol_solv_fin_mec_GMX.itp topology
#include "mol_solv_fin_mec_GMX.itp"

; Include the water model
#include "amber99sb.ff/tip3p.itp"

; Include ions (not alway necessary)
#include "amber99sb.ff/ions.itp"

The order is mandatory because forcefield.itp contains the [ defaults ] section that must be unique whereas the topology file of our molecule has a section [ atomtypes ] to define new atoms that must be located before the definitions of other molecules.

Prepare the file tpr that will be directly used in the minimization for neutral molecules whereas it is used to add ion(s) when the molecule is charged

gmx grompp -f em.mdp -c mol_test_wat.gro -p mol_wat_fin_mec_GMX.top -o em.tpr

Now if the molecule is ionized, it is mandatory to balance the number of charges by the addition of a counter ion, To add ions, the file tpr is required

gmx genion -s em.tpr -o mol_test_wat_fin.gro -p ach_GMX_1.top -pname NA -nname CL -neutral

Select SOL option

Check the resulting gro file to be sure that ion(s) was(were) added and check that topology was modified according to the added ion(s)

If ions were added, rebuild the tpr file

gmx grompp -f em.mdp -c mol_test_wat.gro -p mol_wat_fin_mec_GMX.top -o em.tpr

and run the test minimization

gmx mdrun -ntmpi 1 -v -deffnm em

Check the output em.gro with VMD

Check the result with VMD

MD setup in chloroform

Check that in the top file the following files are present:

mol_GMX.itp
mol_GMX.top
CHL_Box.pdb
CHL.itp
mol.gro

gmx solvate -cp mol.gro -cs CHL_Box.pdb -o mv1_model1_solv.gro -p topol_cyc.top

Check the result with VMD

The topology (.top) file is updated with the addition of a new row in [ molecules ] section

Edit the topology file

delete [ defaults ] section
insert itp information about deafult ff, water model (tip3p) and ions (even if it is necessary only for charged molecules)

; arv_GMX.top created by acpype (v: 2020-03-15T13:27:00CET) on Fri Sep 25 10:01:49 2020

; Include default ff
#include "amber99sb.ff/forcefield.itp"

; Include mol_GMX.itp topology
#include "mol_GMX.itp"

; Include chloroform
#include "CHL.itp"

; Include ions (not alway necessary)
#include "amber99sb.ff/ions.itp"

The order is mandatory because forcefield.itp contains the [ defaults ] section that must be unique whereas the topology file of our molecule has a section [ atomtypes ] to define new atoms that must be located before the definitions of other molecules.
h3 atom type of chloroform is not present in the [Atom types] section of the mol.itp file but this section cannot be called more than one times. To resolve the error it is mandatory to add in the [Atom types] section of the mol.itp a row with h3 definition

[ atomtypes ]
;name bond_type mass charge ptype sigma epsilon Amb
... ...
h3 h3 0.00000 0.00000 A 2.11499e-01 6.56888e-02 ; 1.19 0.0157

Prepare the file tpr that will be directly used in the minimization for neutral molecules whereas it is used to add ion(s) when the molecule is charged

gmx grompp -f em.mdp -c mol_test_wat.gro -p mol_wat_fin_mec_GMX.top -o em.tpr

Now if the molecule is ionized, it is mandatory to balance the number of charges by the addition of a counter ion, To add ions, the file tpr is required

gmx genion -s em.tpr -o mol_test_wat_fin.gro -p ach_GMX_1.top -pname NA -nname CL -neutral

Select the option related to chloroform option

Check the resulting gro file to be sure that ion(s) was(were) added and check that topology was modified according to the added ion(s)

If ions were added, rebuild the tpr file

gmx grompp -f em.mdp -c mol_test_wat.gro -p mol_wat_fin_mec_GMX.top -o em.tpr

and run the test minimization

gmx mdrun -ntmpi 1 -v -deffnm em

Check the output em.gro with VMD

Check the result with VMD

Setup MD runs for MSM

Now, the setup of MD simulations has to be carried out to all the conformations selected from the Conformational Sampling / Cluster step.

In principle, every MD simulation will run for 100 ns, thus 100 conformations will be used to arrive to 100 x 100 ns = 10000 ns = 10 μs (suggested total run time from literature).

Probably, a first multi run of about 50 conformations could be a good starting point for a tentative MSM analysis. The MSM analysis could suggest the new conformations for the remaining simulations.

Firstly, select mol_01-50.mol2 conformations.

Then, check the consistency between the selected conformations and the MEC used in the parametrization

Check the name of the residue, that it must be the same
Check the order of the atoms, that must be the same

All steps are collected in bash scripts that allow the automation of the process.

Step 0. Preparation of directory tree

Create the working directory and the run the script md_step0_setupdir.sh

# Usage: ./md_step0.sh
# setup of the directories

mkdir diagn
mkdir md_out
mkdir mol2

mkdir traj
mkdir traj_msm

mkdir scripts

# File list
#   em.mdp
#   eq-1.mdp
#   eq-2.mdp
#   eq-3.mdp
#   eq-4.mdp
#   production.mdp
#   production_BriefForDebug.mdp

mkdir top

# File list
#   mol.itp
#   mol.top
#   CHL_Box.pdb (if the simulation is run in chloroform)
#   CHL.itp (the should be move in the root simulation directory)

The script makes all the directories.

Copy the mol2 files in the directory mol2.

If you only need to run simulations on some mol2 files, you can create subdirectories to move unnecessary files to. The simulation scripts run on the mol2 file located in the mol2 root directory

Step 1. Setup and minimization

Run the script: ./md_step1.sh basename

# Usage: ./md_step1.sh root_name
# root_name: is the root name of the generated topology file

export GMX_MAXBACKUP=-1

for f in ./mol2/*.mol2; do

    name=$(basename $f .mol2)

    echo $name
    cp mol2/$name.mol2 .
    obabel -imol2 $name.mol2 -opdb -O $name.pdb
    gmx editconf -f $name.pdb -o $name.gro -bt cubic -box 3.5 3.5 3.5
    cp top/$1_GMX.top $name.top
    cp top/$1_GMX.itp .
    gmx solvate -cp $name.gro -o ${name}_sol.gro -p $name.top -cs top/CHL_Box.pdb
    gmx grompp -f scripts/em.mdp -c ${name}_sol.gro -p $name.top -o em.tpr
    echo 4 | gmx genion -s em.tpr -o ${name}_sol.gro -p $name.top -pname NA -nname CL -neutral
    gmx grompp -f scripts/em.mdp -c ${name}_sol.gro -p $name.top -o em.tpr
    mpirun -np 1 mdrun -deffnm em -cpt 1 -cpo em_restart.cpt
    cp em.gro md_out/${name}_sol_em.gro
    mv ${name}_sol.gro md_out

done

mv *.gro md_out
rm *.mol2
rm *.pdb
rm em*

Some remarks:

the basename must be the name used for the conformational sampling/minimization procedure
in the editconf the cubic box dimensions are fixed to 3.5 but the dimensions have to be checked to be sure that the molecule is well positioned inside the box
finally, the number of solvent molecules could be different in the gro file producing an error because in the root topology file the number of chloroform molecules is fixed, in this case it is necessary to modify manually the number of solvent molecules

To check if the run has finished without error, it is possible to count top files, their number must be the same of the number of mol2 files.

ls *.top | wc -l
ls mol2/*.mol2 | wc -l