GROMACS

Setting up SIRAH

Required Software

GROMACS

  • GROMACS 4.5.5 or later version properly installed in your computer.

Tip

Download the lastest version of GROMACS and check its Installation Guide.

VMD

The molecular visualization program VMD, version 1.9.3 or later (freely available download). Check VMD’s Installation guide for more instructions.

Prior knowledge

How to perform a standard atomistic molecular dynamics simulations with GROMACS and basic usage of VMD.

Tip

GROMACS

If you are not familiar with atomistic molecular dynamics simulations, we strongly recommend you to perform the basic usage tutorials of GROMACS: Introduction to Molecular Dynamics, Introduction to Membrane-Protein Simulation and others. For more details and documentation, check GROMACS Manual.

VMD

If you are not familiar with VMD, we strongly recommend you to perform the basic usage tutorial of VMD (using VMD). For more details and documentation, check VMD User Guide and the official VMD documentation.

Download and setting SIRAH

Note

Please report bugs, errors or enhancement requests through Issue Tracker or if you have a question about SIRAH open a New Discussion.

Download SIRAH for GROMACS and uncompress it into your working directory.

tar -xzvf sirah_x2.2_20-07.ff.tar.gz

In case of using another version of sirah:

tar -xzvf ``sirah_[version].tgz``

You will get a folder sirah_[version].ff/ containing the force field definition, the SIRAH Tools in sirah_[version].ff/tools/, molecular structures to build up systems in sirah_[version].ff/PDB, the required material to perform the tutorial in sirah_[version].ff/tutorial/X.X/.

Caution

Remember to change X.X to the number that corresponds to the tutorial you are going to do.

Make a new folder for this tutorial in your working directory:

mkdir tutorialX.X; cd tutorialX.X

Create the following symbolic links in the folder tutorialX.X:

ln -s ../sirah_[version].ff sirah.ff
ln -s sirah.ff/residuetypes.dat
ln -s sirah.ff/specbond.dat

Warning

Files residuetypes.dat and specbond.dat are essential for the correct definition of molecular groups and auto-detection of disulfide bonds and cyclic DNA polymers.

1. DNA molecule in explicit solvent

Note

Please report bugs, errors or enhancement requests through Issue Tracker or if you have a question about SIRAH open a New Discussion.

This tutorial shows how to use the SIRAH force field to perform a coarse grained (CG) simulation of a double stranded DNA in explicit solvent (called WatFour, WT4). The main references for this tutorial are: Dans et al. (latest parameters are those reported here), Darré et al., and Machado & Pantano. We strongly advise you to read these articles before starting the tutorial.

Important

Check the Setting up SIRAH section for download and set up details before starting this tutorial. Since this is Tutorial 1, remember to replace X.X, the files corresponding to this tutorial can be found in: sirah.ff/tutorial/1/

1.1. Build CG representations

Map the atomistic structure of a 20-mer DNA to its CG representation:

./sirah.ff/tools/CGCONV/cgconv.pl -i sirah.ff/tutorial/1/dna.pdb -o dna_cg.pdb

The input file -i dna.pdb contains all the heavy atoms composing the DNA molecule, while the output -o dna_cg.pdb preserves a few of them.

Please check both PDB structures using VMD:

vmd -m ./sirah.ff/tutorial/1/dna.pdb dna_cg.pdb

Tip

This is the basic usage of the script cgconv.pl, you can learn other capabilities from its help by typing:

./sirah.ff/tools/CGCONV/cgconv.pl -h

From now on it is just normal GROMACS stuff!

Caution

In GROMACS versions prior to 5.x, the “gmx” command should not be used.

1.2. PDB to GROMACS format

Use pdb2gmx to convert your PDB file into GROMACS format:

gmx pdb2gmx -f dna_cg.pdb -o dna_cg.gro -merge all

When prompted, choose SIRAH force field and then SIRAH solvent models.

Note

Warning messages about long, triangular or square bonds are fine and expected due to the CG topology of some residues.

Caution

However, missing atom messages are errors which probably trace back to the mapping step. In that case, check your atomistic and mapped structures and do not carry on the simulation until the problem is solved.

Note

Merging both DNA chains is convenient when planning to apply restraints between them.

During long simulations of DNA, capping residues may eventually separate. If you want to avoid this effect, which is called helix fraying, add Watson-Crick (WC) restraints at terminal base pairs. To add these restraints edit topol.top to include the file WC_RST.itp at the end of the [ moleculetype ] section:

Topology without WC restraints

Topology with WC restraints

; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif



; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif

; Watson-Crick restraints
#include "./sirah.ff/tutorial/1/WC_RST.itp"

1.3. Solvate the system

Define the simulation box of the system

gmx editconf -f dna_cg.gro -o dna_cg_box.gro -bt octahedron -d 2 -c

Add WT4 molecules:

gmx solvate -cp dna_cg_box.gro -cs sirah.ff/wt416.gro -o dna_cg_sol.gro

Note

In GROMACS versions prior to 5.x, the command gmx solvate was called genbox.

Edit the [ molecules ] section in topol.top to include the number of added WT4 molecules:

Topology before editing

Topology after editing

[ molecules ]
; Compound #mols
DNA_chain_A 1

[ molecules ]
; Compound #mols
DNA_chain_A 1
WT4 3179

Hint

If you forget to read the number of added WT4 molecules from the output of solvate, then use the following command line to get it

grep -c WP1 dna_cg_sol.gro

Caution

The number of added WT4 molecules, 3179, may change according to the software version.

Add CG counterions and 0.15M NaCl:

gmx grompp -f sirah.ff/tutorial/1/GPU/em_CGDNA.mdp -p topol.top -c dna_cg_sol.gro -o dna_cg_sol.tpr
gmx genion -s dna_cg_sol.tpr -o dna_cg_ion.gro -np 113 -pname NaW -nn 75 -nname ClW

When prompted, choose to substitute WT4 molecules by ions.

Note

The available electrolyte species in SIRAH force field are: Na⁺ (NaW), K⁺ (KW) and Cl⁻ (ClW) which represent solvated ions in solution. One ion pair (e.g., NaW-ClW) each 34 WT4 molecules results in a salt concentration of ~0.15M (see Appendix for details). Counterions were added according to Machado et al..

Edit the [ molecules ] section in topol.top to include the CG ions and the correct number of WT4.

Before running the simulation it may be a good idea to visualize your molecular system. CG molecules are not recognized by molecular visualizers and will not display correctly. To fix this problem you may generate a PSF file of the system using the script g_top2psf.pl:

./sirah.ff/tools/g_top2psf.pl -i topol.top -o dna_cg_ion.psf

Note

This is the basic usage of the script g_top2psf.pl, you can learn other capabilities from its help:

./sirah.ff/tools/g_top2psf.pl -h

Use VMD to check how the CG system looks like:

vmd dna_cg_ion.psf dna_cg_ion.gro -e sirah.ff/tools/sirah_vmdtk.tcl

Tip

VMD assigns default radius to unknown atom types, the script sirah_vmdtk.tcl sets the right ones, according to the CG representation. It also provides a kit of useful selection macros, coloring methods and backmapping utilities. Use the command sirah_help in the Tcl/Tk console of VMD to access the manual pages. To learn about SIRAH Tools’ capabilities, you can also go to the SIRAH Tools tutorial.

1.4. Run the simulation

Important

By default in this tutorial we will use input files for GROMACS on GPU (sirah.ff/tutorial/1/GPU). Example input files for using GROMACS on CPU can be found at: sirah.ff/tutorial/1/CPU.

The folder sirah.ff/tutorial/1/GPU/ contains typical input files for energy minimization (em_CGDNA.mdp), equilibration (eq_CGDNA.mdp) and production (md_CGDNA.mdp) runs. Please check carefully the input flags therein.

Create an index file:

echo "q" | gmx make_ndx -f dna_cg_ion.gro -o dna_cg_ion.ndx

Note

WT4 and CG ions (NaW and ClW) are automatically set to the group SIRAH-Solvent.

Make a new folder for the run:

mkdir run; cd run

Energy Minimization:

gmx grompp -f ../sirah.ff/tutorial/1/GPU/em_CGDNA.mdp -p ../topol.top -po em.mdp -n ../dna_cg_ion.ndx -c ../dna_cg_ion.gro -o dna_cg_em.tpr
mdrun -deffnm dna_cg_em &> EM.log &

Equilibration:

gmx grompp -f ../sirah.ff/tutorial/1/GPU/eq_CGDNA.mdp -p ../topol.top -po eq.mdp -n ../dna_cg_ion.ndx -c dna_cg_em.gro -o dna_cg_eq.tpr
mdrun -deffnm dna_cg_eq &> EQ.log &

Production (100ns):

gmx grompp -f ../sirah.ff/tutorial/1/GPU/md_CGDNA.mdp -p ../topol.top -po md.mdp -n ../dna_cg_ion.ndx -c dna_cg_eq.gro -o dna_cg_md.tpr
mdrun -deffnm dna_cg_md &> MD.log &

Note

GPU flags have been set for GROMACS 4.6.7; however, different versions may object to certain specifications.

1.5. Visualizing the simulation

That’s it! Now you can analyze the trajectory.

Process the output trajectory at folder run/ to account for the Periodic Boundary Conditions (PBC):

gmx trjconv -s dna_cg_em.tpr -f dna_cg_md.xtc -o dna_cg_md_pbc.xtc -n ../dna_cg_ion.ndx -ur compact -center -pbc mol

When prompted, choose DNA for centering and System for output.

Now you can check the simulation using VMD:

vmd ../dna_cg_ion.psf ../dna_cg_ion.gro dna_cg_md_pbc.xtc -e ../sirah.ff/tools/sirah_vmdtk.tcl

Note

The file sirah_vmdtk.tcl is a Tcl script that is part of SIRAH Tools and contains the macros to properly visualize the coarse-grained structures in VMD. Use the command sirah-help in the Tcl/Tk console of VMD to access the manual pages. To learn about SIRAH Tools’ capabilities, you can also go to the SIRAH Tools tutorial.

2. Hybrid solvation

Note

Please report bugs, errors or enhancement requests through Issue Tracker or if you have a question about SIRAH open a New Discussion.

This tutorial shows how to apply the hybrid solvation approach of SIRAH force field to speed up the simulation of an atomistic solute. The example system contains a DNA fiber surrounded by a shell of atomistic waters, which are embedded in coarse-grained (CG) molecules, called WT4, to represent bulk water. The general procedure is extensible to any other solute. The hybrid solvation methodology is well tested for SPC, SPC/E, TIP3P atomistic water models. The main references for this tutorial are: Darré et al., Machado et al., Gonzalez et al., and Machado & Pantano. We strongly advise you to read these articles before starting the tutorial.

Important

Check the Setting up SIRAH section for download and set up details before starting this tutorial. Since this is Tutorial 2, remember to replace X.X, the files corresponding to this tutorial can be found in: sirah.ff/tutorial/2/

2.1. PDB to GROMACS format

Use pdb2gmx to convert the PDB file of the solute dna.pdb into GROMACS format:

gmx pdb2gmx -f sirah.ff/tutorial/2/dna.pdb -o dna.gro

Caution

In GROMACS versions prior to 5.x, the “gmx” command should not be used.

When prompted, choose AMBER99SB and then TIP3P as water model.

For GROMACS to recognize SIRAH, edit your topology file topol.top adding the following lines after the force field definition:

Topology before editing

Topology after editing

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


; Include forcefield parameters
#include "amber99sb.ff/forcefield.itp"
#include "./sirah.ff/hybsol_comb2.itp"
#include "./sirah.ff/solv.itp"

Important

hybsol_comb2.itp is a parameter file, while solv.itp links to the topologies of WT4 and SIRAH ions. The choice of the right parameter file depends on the chosen atomistic force field (see Table). The plug-in works smoothly with the implemented force fields of the GROMACS distribution. When using a customized force field (e.g. for lipids) choose the parameter file according to the combination rule and check that the atom type you are using for the atomistic water match the DEFAULT in the SIRAH parameter file, which is OW. If they don’t match then rename OW accordingly to your definition.

Table. SIRAH parameter file to include in the topology according to the chosen atomistic force field.

Atomistic force field

Parameter file

Combination rule

GMX

GROMOS

AMBER

CHARMM

OPLS

hybsol_comb1.itp

1

X

X

hybsol_comb2.itp

2

X

X

hybsol_comb3.itp

3

X

2.2. Solvate the system

Define the simulation box of the system:

gmx editconf -f dna.gro -o dna_box.gro -bt octahedron -d 2 -c

Then add an atomistic water shell:

gmx solvate -cp dna_box.gro -cs spc216.gro -o dna_shell.gro -shell 1

Note

In GROMACS versions prior to 5.x, the command gmx solvate was called genbox.

Finally add the CG solvent:

gmx solvate -cp dna_shell.gro -cs ./sirah.ff/wt4tip3p.gro -o dna_sol.gro

Edit the [ molecules ] section in topol.top to include the number of SOL and WT4 molecules:

Topology before editing

Topology after editing

[ molecules ]
; Compound #mols
DNA_chain_A 1
DNA_chain_B 1


[ molecules ]
; Compound #mols
DNA_chain_A 1
DNA_chain_B 1
SOL 3580
WT4 2659

Hint

If you forget to look at the result of solvate to see how many SOL and WT4 molecules were added, you can use the following command line to get these numbers:

grep -c OW dna_sol.gro; grep -c WP1 dna_sol.gro

Caution

The number of added SOL (3580) and WT4 (2659) molecules, may change according to the software version.

Add CG counterions:

gmx grompp -f sirah.ff/tutorial/2/GPU/em_HYBSOL.mdp -p topol.top -c dna_sol.gro -o dna_sol.tpr
gmx genion -s dna_sol.tpr -o dna_ion.gro -np 38 -pname NaW

When prompted, choose to substitute WT4 molecules by NaW ions.

Note

The available electrolyte species in SIRAH force field are: Na⁺ (NaW), K⁺ (KW) and Cl⁻ (ClW) which represent solvated ions in solution. One ion pair (e.g., NaW-ClW) each 34 WT4 molecules results in a salt concentration of ~0.15M (see Appendix for details). Be aware that SIRAH ions remain within the CG phase. So, if the presence of atomistic electrolytes in close contact with the solute is important to describe the physics of the system you will have to add them.

Edit the [ molecules ] section in topol.top to include the 38 NaW ions and the correct number of WT4 molecules.

Before running the simulation it may be a good idea to visualize your molecular system. CG molecules are not recognized by molecular visualizers and will not display correctly. To fix this problem you may generate a PSF file of the system using the script g_top2psf.pl:

./sirah.ff/tools/g_top2psf.pl -i topol.top -o dna_ion.psf

Note

This is the basic usage of the script g_top2psf.pl, you can learn other capabilities from its help:

./sirah.ff/tools/g_top2psf.pl -h

Use VMD to check how the hybrid system looks like:

vmd dna_ion.psf dna_ion.gro -e sirah.ff/tools/sirah_vmdtk.tcl

Note

VMD assigns default radius to unknown atom types, the script sirah_vmdtk.tcl sets the right ones, according to the CG representation. It also provides a kit of useful selection macros, coloring methods and backmapping utilities. Use the command sirah_help in the Tcl/Tk console of VMD to access the manual pages. To learn about SIRAH Tools’ capabilities, you can also go to the SIRAH Tools tutorial.

2.3. Run the simulation

Important

By default in this tutorial we will use input files for GROMACS on GPU (sirah.ff/tutorial/2/GPU). Example input files for using GROMACS on CPU can be found at: sirah.ff/tutorial/2/CPU.

The folder sirah.ff/tutorial/2/GPU/ contains typical input files for energy minimization (em_HYBSOL.mdp), equilibration (eq_HYBSOL.mdp) and production (md_HYBSOL.mdp) runs. Please check carefully the input flags therein.

Create an index files adding a group for WT4 and NaW:

echo -e "r WT4 | r NaW\nq\n" | make_ndx -f dna_ion.gro -o dna_ion.ndx

Note

WT4 and CG ions (NaW and ClW) are automatically set to the group SIRAH-Solvent.

Make a new folder for the run:

mkdir run; cd run

Energy Minimization:

gmx grompp -f ../sirah.ff/tutorial/2/GPU/em_HYBSOL.mdp -p ../topol.top -po em.mdp -n ../dna_ion.ndx -c ../dna_ion.gro -o dna_em.tpr
mdrun -deffnm dna_em &> EM.log &

Equilibration:

gmx grompp -f ../sirah.ff/tutorial/2/GPU/eq_HYBSOL.mdp -p ../topol.top -po eq.mdp -n ../dna_ion.ndx -c dna_em.gro -o dna_eq.tpr
mdrun -deffnm dna_eq &> EQ.log &

Production (100ns):

gmx grompp -f ../sirah.ff/tutorial/2/GPU/md_HYBSOL.mdp -p ../topol.top -po md.mdp -n ../dna_ion.ndx -c dna_eq.gro -o dna_md.tpr
mdrun -deffnm dna_md &> MD.log &

Note

GPU flags have been set for GROMACS 4.6.7; however, different versions may object to certain specifications.

2.4. Visualizing the simulation

That’s it! Now you can analyze the trajectory.

Process the output trajectory at folder run/ to account for the Periodic Boundary Conditions (PBC):

gmx trjconv -s dna_em.tpr -f dna_md.xtc -o dna_md_pbc.xtc -n ../dna_ion.ndx -ur compact -center -pbc mol

When prompted, choose DNA for centering and System for output.

Now you can check the simulation using VMD:

vmd ../dna_ion.psf ../dna_ion.gro dna_md_pbc.xtc -e ../sirah.ff/tools/sirah_vmdtk.tcl

Note

The file sirah_vmdtk.tcl is a Tcl script that is part of SIRAH Tools and contains the macros to properly visualize the coarse-grained structures in VMD. Use the command sirah-help in the Tcl/Tk console of VMD to access the manual pages. To learn about SIRAH Tools’ capabilities, you can also go to the SIRAH Tools tutorial.

3. Proteins in explicit solvent

Note

Please report bugs, errors or enhancement requests through Issue Tracker or if you have a question about SIRAH open a New Discussion.

This tutorial shows how to use the SIRAH force field to perform a coarse grained (CG) simulation of a protein in explicit solvent (called WatFour, WT4). The main references for this tutorial are: Dans et al., Darré et al., Machado et al., and Machado & Pantano. We strongly advise you to read these articles before starting the tutorial.

Important

Check the Setting up SIRAH section for download and set up details before starting this tutorial. Since this is Tutorial 3, remember to replace X.X, the files corresponding to this tutorial can be found in: sirah.ff/tutorial/3/

3.1. Build CG representations

Caution

The mapping to CG requires the correct protonation state of each residue at a given pH. We recommend using the CHARMM-GUI server and use the PDB Reader & Manipulator to prepare your system. An account is required to access any of the CHARMM-GUI Input Generator modules, and it can take up to 24 hours to obtain one.

Other option is the PDB2PQR server and choosing the output naming scheme of AMBER for best compatibility. This server was utilized to generate the PQR file featured in this tutorial. Be aware that modified residues lacking parameters such as: MSE (seleno MET), TPO (phosphorylated THY), SEP (phosphorylated SER) or others are deleted from the PQR file by the server. In that case, mutate the residues to their unmodified form before submitting the structure to the server.

Map the protonated atomistic structure of protein 1CRN to its CG representation:

./sirah.ff/tools/CGCONV/cgconv.pl -i sirah.ff/tutorial/3/1CRN.pqr -o 1CRN_cg.pdb

The input file -i 1CRN.pqr contains the atomistic representation of 1CRN structure at pH 7.0, while the output -o 1CRN_cg.pdb is its SIRAH CG representation.

Tip

This is the basic usage of the script cgconv.pl, you can learn other capabilities from its help by typing:

./sirah.ff/tools/CGCONV/cgconv.pl -h

Note

Pay attention to residue names when mapping structures from other atomistic force fields or experimental structures. Although we provide compatibility for naming schemes in PDB, GMX, GROMOS, CHARMM and OPLS, there might always be some ambiguity in the residue naming, specially regarding protonation states, that may lead to a wrong mapping. For example, SIRAH Tools always maps the residue name “HIS” to a Histidine protonated at the epsilon nitrogen (\(N_{\epsilon}\)) regardless the actual proton placement. Similarly, protonated Glutamic and Aspartic acid residues must be named “GLH” and “ASH”, otherwise they will be treated as negative charged residues. In addition, protonated and disulfide bonded Cysteines must be named “CYS” and “CYX” respectively. These kind of situations need to be carefully checked by the users. In all cases the residues preserve their identity when mapping and back-mapping the structures. Hence, the total charge of the protein should be the same at atomistic and SIRAH levels. You can check the following mapping file to be sure of the compatibility: sirah.ff/tools/CGCONV/maps/sirah_prot.map.

Please check both PDB and PQR structures using VMD:

vmd -m sirah.ff/tutorial/3/1CRN.pqr 1CRN_cg.pdb

From now on it is just normal GROMACS stuff!

Caution

In GROMACS versions prior to 5.x, the “gmx” command should not be used.

3.2. PDB to GROMACS format

Use pdb2gmx to convert your PDB file into GROMACS format:

gmx pdb2gmx -f 1CRN_cg.pdb -o 1CRN_cg.gro

When prompted, choose SIRAH force field and then SIRAH solvent models.

Note

By default charged terminal are used, but it is possible to set them neutral with option -ter.

Note

Warning messages about long, triangular or square bonds are fine and expected due to the CG topology of some residues.

Caution

However, missing atom messages are errors which probably trace back to the mapping step. In that case, check your atomistic and mapped structures and do not carry on the simulation until the problem is solved.

3.3. Solvate the system

Define the simulation box of the system

gmx editconf -f 1CRN_cg.gro -o 1CRN_cg_box.gro -bt octahedron -d 2.0 -c

Add WT4 molecules:

gmx solvate -cp 1CRN_cg_box.gro -cs sirah.ff/wt416.gro -o 1CRN_cg_sol1.gro

Note

In GROMACS versions prior to 5.x, the command gmx solvate was called genbox.

Edit the [ molecules ] section in topol.top to include the number of added WT4 molecules:

Topology before editing

Topology after editing

[ molecules ]
; Compound #mols
Protein_chain_A 1

[ molecules ]
; Compound #mols
Protein_chain_A 1
WT4 756

Hint

If you forget to read the number of added WT4 molecules from the output of solvate, then use the following command line to get it

grep -c WP1 1CRN_cg_sol1.gro

Caution

The number of added WT4 molecules, 756, may change according to the software version.

Remove WT4 molecules within 0.3 nm of protein:

echo q | gmx make_ndx -f 1CRN_cg_sol1.gro -o 1CRN_cg_sol1.ndx
gmx grompp -f sirah.ff/tutorial/3/GPU/em1_CGPROT.mdp -p topol.top -po delete1.mdp -c 1CRN_cg_sol1.gro -o 1CRN_cg_sol1.tpr
gmx select -f 1CRN_cg_sol1.gro -s 1CRN_cg_sol1.tpr -n 1CRN_cg_sol1.ndx -on rm_close_wt4.ndx -select 'not (same residue as (resname WT4 and within 0.3 of group Protein))'
gmx editconf -f 1CRN_cg_sol1.gro -o 1CRN_cg_sol2.gro -n rm_close_wt4.ndx

Note

In GROMACS versions prior to 5.x, the command gmx select was called g_select.

Edit the [ molecules ] section in topol.top to include the correct number of WT4 molecules:

grep -c WP1 1CRN_cg_sol2.gro

Add CG counterions and 0.15M NaCl:

gmx grompp -f sirah.ff/tutorial/3/GPU/em1_CGPROT.mdp -p topol.top -po delete2.mdp -c 1CRN_cg_sol2.gro -o 1CRN_cg_sol2.tpr
gmx genion -s 1CRN_cg_sol2.tpr -o 1CRN_cg_ion.gro -np 22 -pname NaW -nn 22 -nname ClW

When prompted, choose to substitute WT4 molecules by ions.

Note

The available electrolyte species in SIRAH force field are: Na⁺ (NaW), K⁺ (KW) and Cl⁻ (ClW) which represent solvated ions in solution. One ion pair (e.g., NaW-ClW) each 34 WT4 molecules results in a salt concentration of ~0.15M (see Appendix for details). Counterions were added according to Machado et al..

Edit the [ molecules ] section in topol.top to include the CG ions and the correct number of WT4.

Before running the simulation it may be a good idea to visualize your molecular system. CG molecules are not recognized by molecular visualizers and will not display correctly. To fix this problem you may generate a PSF file of the system using the script g_top2psf.pl:

./sirah.ff/tools/g_top2psf.pl -i topol.top -o 1CRN_cg_ion.psf

Note

This is the basic usage of the script g_top2psf.pl, you can learn other capabilities from its help:

./sirah.ff/tools/g_top2psf.pl -h

Use VMD to check how the CG system looks like:

vmd 1CRN_cg_ion.psf 1CRN_cg_ion.gro -e sirah.ff/tools/sirah_vmdtk.tcl

Note

VMD assigns default radius to unknown atom types, the script sirah_vmdtk.tcl sets the right ones, according to the CG representation. It also provides a kit of useful selection macros, coloring methods and backmapping utilities. Use the command sirah_help in the Tcl/Tk console of VMD to access the manual pages. To learn about SIRAH Tools’ capabilities, you can also go to the SIRAH Tools tutorial.

Create an index file including a group for the backbone GN and GO beads:

echo -e "a GN GO\n\nq" | gmx make_ndx -f 1CRN_cg_ion.gro -o 1CRN_cg_ion.ndx

Note

WT4 and CG ions (NaW, ClW) are automatically set to the group SIRAH-Solvent.

Generate restraint files for the backbone GN and GO beads:

gmx genrestr -f 1CRN_cg.gro -n 1CRN_cg_ion.ndx -o bkbres.itp
gmx genrestr -f 1CRN_cg.gro -n 1CRN_cg_ion.ndx -o bkbres_soft.itp -fc 100 100 100

When prompted, choose the group GN_GO

Add the restraints to topol.top:

Topology before editing

Topology after editing

; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif








; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif
; Backbone restraints
#ifdef GN_GO
#include "bkbres.itp"
#endif
; Backbone soft restrains
#ifdef GN_GO_SOFT
#include "bkbres_soft.itp"
#endif

3.4. Run the simulation

Important

By default in this tutorial we will use input files for GROMACS on GPU (sirah.ff/tutorial/3/GPU). Example input files for using GROMACS on CPU can be found at: sirah.ff/tutorial/3/CPU.

The folder sirah.ff/tutorial/3/GPU/ contains typical input files for energy minimization (em1_CGPROT.mdp, em2_CGPROT.mdp), equilibration (eq1_CGPROT.mdp, eq2_CGPROT.mdp) and production (md_CGPROT.mdp) runs. Please check carefully the input flags therein.

Make a new folder for the run:

mkdir -p run; cd run

Energy Minimization of side chains by restraining the backbone:

gmx grompp -f ../sirah.ff/tutorial/3/GPU/em1_CGPROT.mdp -p ../topol.top -po em1.mdp -n ../1CRN_cg_ion.ndx -c ../1CRN_cg_ion.gro -r ../1CRN_cg_ion.gro -o 1CRN_cg_em1.tpr
gmx mdrun -deffnm 1CRN_cg_em1 &> EM1.log &

Energy Minimization of whole system:

gmx grompp -f ../sirah.ff/tutorial/3/GPU/em2_CGPROT.mdp -p ../topol.top -po em2.mdp -n ../1CRN_cg_ion.ndx -c 1CRN_cg_em1.gro -o 1CRN_cg_em2.tpr
gmx mdrun -deffnm 1CRN_cg_em2 &> EM2.log &

Solvent equilibration:

gmx grompp -f ../sirah.ff/tutorial/3/GPU/eq1_CGPROT.mdp -p ../topol.top -po eq1.mdp -n ../1CRN_cg_ion.ndx -c 1CRN_cg_em2.gro -r 1CRN_cg_em2.gro -o 1CRN_cg_eq1.tpr
gmx mdrun -deffnm 1CRN_cg_eq1 &> EQ1.log &

Soft equilibration to improve side chain solvation:

gmx grompp -f ../sirah.ff/tutorial/3/GPU/eq2_CGPROT.mdp -p ../topol.top -po eq2.mdp -n ../1CRN_cg_ion.ndx -c 1CRN_cg_eq1.gro -r 1CRN_cg_eq1.gro -o 1CRN_cg_eq2.tpr
gmx mdrun -deffnm 1CRN_cg_eq2 &> EQ2.log &

Production (1000 ns):

gmx grompp -f ../sirah.ff/tutorial/3/GPU/md_CGPROT.mdp -p ../topol.top -po md.mdp -n ../1CRN_cg_ion.ndx -c 1CRN_cg_eq2.gro -o 1CRN_cg_md.tpr
gmx mdrun -deffnm 1CRN_cg_md &> MD.log &

Note

GPU flags have been set for GROMACS 4.6.7; however, different versions may object to certain specifications.

3.5. Visualizing the simulation

That’s it! Now you can analyze the trajectory.

Process the output trajectory at folder run/ to account for the Periodic Boundary Conditions (PBC):

gmx trjconv -s 1CRN_cg_em1.tpr -f 1CRN_cg_md.xtc -o 1CRN_cg_md_pbc.xtc -n ../1CRN_cg_ion.ndx -ur compact -center -pbc mol

When prompted, choose Protein for centering and System for output.

Now you can check the simulation using VMD:

vmd ../1CRN_cg_ion.psf ../1CRN_cg_ion.gro 1CRN_cg_md_pbc.xtc -e ../sirah.ff/tools/sirah_vmdtk.tcl

Note

The file sirah_vmdtk.tcl is a Tcl script that is part of SIRAH Tools and contains the macros to properly visualize the coarse-grained structures in VMD. Use the command sirah-help in the Tcl/Tk console of VMD to access the manual pages. To learn about SIRAH Tools’ capabilities, you can also go to the SIRAH Tools tutorial.

4. Closed circular DNA in explicit solvent

Note

Please report bugs, errors or enhancement requests through Issue Tracker or if you have a question about SIRAH open a New Discussion.

This tutorial shows how to use the SIRAH force field to perform a coarse grained (CG) simulation of a closed circular DNA in explicit solvent (called WatFour, WT4). The main references for this tutorial are:Dans et al. (latest parameters are those reported here), and Machado & Pantano. We strongly advise you to read these articles before starting the tutorial.

Important

Check the Setting up SIRAH section for download and set up details before starting this tutorial. Since this is Tutorial 4, remember to replace X.X, the files corresponding to this tutorial can be found in: sirah.ff/tutorial/4/

4.1. Build CG representations

Map the atomistic structure of the closed circular DNA to its CG representation:

./sirah.ff/tools/CGCONV/cgconv.pl -i sirah.ff/tutorial/4/ccdna.pdb | sed -e 's/DCX A 1/CW5 A 1/' -e 's/DCX B 1/CW5 B 1/' > ccdna_cg.pdb

The input file -i ccdna.pdb contains all the heavy atoms composing the DNA molecule, while the output -o ccdna_cg.pdb preserves a few of them.

Important

5’ end residues mast be renamed to AW5, TW5, GW5 or CW5 to represent the corresponding Adenine, Thymine, Guanine or Cytosine extremes in a closed circular DNA.

Please check both PDB structures using VMD:

vmd -m ./sirah.ff/tutorial/4/ccdna.pdb ccdna_cg.pdb

Tip

This is the basic usage of the script cgconv.pl, you can learn other capabilities from its help by typing:

./sirah.ff/tools/CGCONV/cgconv.pl -h

From now on it is just normal GROMACS stuff!

Caution

In GROMACS versions prior to 5.x, the “gmx” command should not be used.

4.2. PDB to GROMACS format

Use pdb2gmx to convert your PDB file into GROMACS format:

gmx pdb2gmx -f ccdna_cg.pdb -o ccdna_cg.gro

When prompted, choose SIRAH force field and then SIRAH solvent models.

Note

Warning messages about long, triangular or square bonds are fine and expected due to the CG topology of some residues.

Caution

However, missing atom messages are errors which probably trace back to the mapping step. In that case, check your atomistic and mapped structures and do not carry on the simulation until the problem is solved.

Note

The distance cut-off to close the DNA molecule is read from specbond.dat. This parameter may need to be modified according to your system’s topology.

4.3. Solvate the system

Define the simulation box of the system

gmx editconf -f ccdna_cg.gro -o ccdna_cg_box.gro -bt octahedron -d 2 -c

Add WT4 molecules:

gmx solvate -cp ccdna_cg_box.gro -cs sirah.ff/wt416.gro -o ccdna_cg_sol.gro

Note

In GROMACS versions prior to 5.x, the command gmx solvate was called genbox.

Edit the [ molecules ] section in topol.top to include the number of added WT4 molecules:

Topology before editing

Topology after editing

[ molecules ]
; Compound #mols
DNA_chain_A 1
DNA_chain_B 1

[ molecules ]
; Compound #mols
DNA_chain_A 1
DNA_chain_B 1
WT4 10240

Hint

If you forget to read the number of added WT4 molecules from the output of solvate, then use the following command line to get it

grep -c WP1 ccdna_cg_sol.gro

Caution

The number of added WT4 molecules, 10240, may change according to the software version.

Add CG counterions and 0.15M NaCl:

gmx grompp -f sirah.ff/tutorial/4/GPU/em_CGDNA.mdp -p topol.top -c ccdna_cg_sol.gro -o ccdna_cg_sol.tpr
gmx genion -s ccdna_cg_sol.tpr -o ccdna_cg_ion.gro -np 401 -pname NaW -nn 201 -nname ClW

When prompted, choose to substitute WT4 molecules by ions.

Note

The available electrolyte species in SIRAH force field are: Na⁺ (NaW), K⁺ (KW) and Cl⁻ (ClW) which represent solvated ions in solution. One ion pair (e.g., NaW-ClW) each 34 WT4 molecules results in a salt concentration of ~0.15M (see Appendix for details). Counterions were added according to Machado et al..

Edit the [ molecules ] section in topol.top to include the CG ions and the correct number of WT4.

Before running the simulation it may be a good idea to visualize your molecular system and check if the DNA circle is closed. CG molecules are not recognized by molecular visualizers and will not display correctly. To fix this problem you may generate a PSF file of the system using the script g_top2psf.pl:

./sirah.ff/tools/g_top2psf.pl -i topol.top -o ccdna_cg_ion.psf

Note

This is the basic usage of the script g_top2psf.pl, you can learn other capabilities from its help:

./sirah.ff/tools/g_top2psf.pl -h

Use VMD to check how the CG system looks like:

vmd ccdna_cg_ion.psf ccdna_cg_ion.gro -e sirah.ff/tools/sirah_vmdtk.tcl

Note

VMD assigns default radius to unknown atom types, the script sirah_vmdtk.tcl sets the right ones, according to the CG representation. It also provides a kit of useful selection macros, coloring methods and backmapping utilities. Use the command sirah_help in the Tcl/Tk console of VMD to access the manual pages. To learn about SIRAH Tools’ capabilities, you can also go to the SIRAH Tools tutorial.

4.4. Run the simulation

Important

By default in this tutorial we will use input files for GROMACS on GPU (sirah.ff/tutorial/4/GPU). Example input files for using GROMACS on CPU can be found at: sirah.ff/tutorial/4/CPU.

The folder sirah.ff/tutorial/4/GPU/ contains typical input files for energy minimization (em_CGDNA.mdp), equilibration (eq_CGDNA.mdp) and production (md_CGDNA.mdp) runs. Please check carefully the input flags therein.

Create an index file:

echo "q" | gmx make_ndx -f ccdna_cg_ion.gro -o ccdna_cg_ion.ndx

Note

WT4 and CG ions (NaW and ClW) are automatically set to the group SIRAH-Solvent.

Make a new folder for the run:

mkdir run; cd run

Energy Minimization:

gmx grompp -f ../sirah.ff/tutorial/4/GPU/em_CGDNA.mdp -p ../topol.top -po em.mdp -n ../ccdna_cg_ion.ndx -c ../ccdna_cg_ion.gro -o ccdna_cg_em.tpr
mdrun -deffnm ccdna_cg_em &> EM.log &

Equilibration:

gmx grompp -f ../sirah.ff/tutorial/4/GPU/eq_CGDNA.mdp -p ../topol.top -po eq.mdp -n ../ccdna_cg_ion.ndx -c ccdna_cg_em.gro -o ccdna_cg_eq.tpr
mdrun -deffnm ccdna_cg_eq &> EQ.log &

Production (100ns):

gmx grompp -f ../sirah.ff/tutorial/4/GPU/md_CGDNA.mdp -p ../topol.top -po md.mdp -n ../ccdna_cg_ion.ndx -c ccdna_cg_eq.gro -o ccdna_cg_md.tpr mdrun
mdrun -deffnm ccdna_cg_md &> MD.log &

Note

GPU flags have been set for GROMACS 4.6.7; however, different versions may object to certain specifications.

4.5. Visualizing the simulation

That’s it! Now you can analyze the trajectory.

Process the output trajectory at folder run/ to account for the Periodic Boundary Conditions (PBC):

gmx trjconv -s ccdna_cg_em.tpr -f ccdna_cg_md.xtc -o ccdna_cg_md_pbc.xtc -n ../ccdna_cg_ion.ndx -ur compact -center -pbc mol

When prompted, choose DNA for centering and System for output.

Now you can check the simulation using VMD:

vmd ../ccdna_cg_ion.psf ../ccdna_cg_ion.gro ccdna_cg_md_pbc.xtc -e ../sirah.ff/tools/sirah_vmdtk.tcl

Note

The file sirah_vmdtk.tcl is a Tcl script that is part of SIRAH Tools and contains the macros to properly visualize the coarse-grained structures in VMD. Use the command sirah-help in the Tcl/Tk console of VMD to access the manual pages. To learn about SIRAH Tools’ capabilities, you can also go to the SIRAH Tools tutorial.

5. Lipid bilayers in explicit solvent

Note

Please report bugs, errors or enhancement requests through Issue Tracker or if you have a question about SIRAH open a New Discussion.

This tutorial shows how to use the SIRAH force field to perform a coarse grained (CG) simulation of a DMPC bilayer in explicit solvent (called WatFour, WT4). The main references for this tutorial are: Barrera et al. and Machado & Pantano. We strongly advise you to read these articles before starting the tutorial. You may also find interesting this book chapter.

Important

Check the Setting up SIRAH section for download and set up details before starting this tutorial. Since this is Tutorial 5, remember to replace X.X, the files corresponding to this tutorial can be found in: sirah.ff/tutorial/5/

5.1. Build CG representations

Map the atomistic structure of the preassembled DMPC bilayer to its CG representation:

./sirah.ff/tools/CGCONV/cgconv.pl -i sirah.ff/tutorial/5/DMPC64.pdb -o DMPC64_cg.pdb -a sirah.ff/tools/CGCONV/maps/tieleman_lipid.map

The input file -i DMPC64.pdb contains the atomistic representation of the DMPC bilayer, while the output -o DMPC64_cg.pdb is its SIRAH CG representation. The flag -a is a mapping scheme for lipids.

Important

By default, no mapping is applied to lipids, as there is no standard naming convention for them. So users are requested to append a MAP file from the list in Table 1, by setting the flag -a in cgconv.pl. We recommend using PACKMOL for building the system. Reference building-block structures are provided at folder sirah.ff/PDB/, which agree with the mapping scheme in sirah.ff/tools/CGCONV/maps/tieleman_lipid.map. The provided DMPC bilayer contains 64 lipid molecules per leaflet distributed in a 6.4 * 6.4 nm surface, taking into account an approximate area per lipid of 0.64 nm2 at 333 K. The starting configuration was created with the input file sirah.ff/tutorial/5/DMPC_bilayer.pkm. See FAQs for cautions on mapping lipids to SIRAH and tips on using fragment-based topologies. If you do not find your issue please start a discussion in our github discussion page F&Q.

Tip

This an advanced usage of the script cgconv.pl, you can learn other capabilities from its help by typing:

./sirah.ff/tools/CGCONV/cgconv.pl -h

The input file DMPC64.pdb contains all the heavy atoms composing the lipids, while the output DMPC64_cg.pdb preserves a few of them. Please check both PDB structures using VMD:

vmd -m sirah.ff/tutorial/5/DMPC64.pdb DMPC64_cg.pdb

From now on it is just normal GROMACS stuff!

Caution

In GROMACS versions prior to 5.x, the “gmx” command should not be used.

5.2. PDB to GROMACS format

Use pdb2gmx to convert your PDB file into GROMACS format:

gmx pdb2gmx -f DMPC64_cg.pdb -o DMPC64_cg.gro

When prompted, choose SIRAH force field and then SIRAH solvent models.

Note

Warning messages about long, triangular or square bonds are fine and expected due to the CG topology of some residues.

Caution

However, missing atom messages are errors which probably trace back to the mapping step. In that case, check your atomistic and mapped structures and do not carry on the simulation until the problem is solved.

5.3. Solvate the system

Define the simulation box of the system

gmx editconf -f DMPC64_cg.gro -o DMPC64_cg_box.gro -box 6.6 6.6 10 -c

Note

As PACKMOL does not consider periodicity while building up the system, increasing the box size a few Angstroms may be required to avoid bad contacts between images.

Add WT4 molecules:

gmx solvate -cp DMPC64_cg_box.gro -cs sirah.ff/wt416.gro -o DMPC64_cg_sol1.gro

Note

In GROMACS versions prior to 5.x, the command gmx solvate was called genbox.

Edit the [ molecules ] section in topol.top to include the number of added WT4 molecules:

Topology before editing

Topology after editing

[ molecules ]
; Compound #mols
Lipid_chain_A 1
Lipid_chain_B 1

[ molecules ]
; Compound #mols
Lipid_chain_A 1
Lipid_chain_B 1
WT4 850

Hint

If you forget to read the number of added WT4 molecules from the output of solvate, then use the following command line to get it

grep -c WP1 DMPC64_cg_sol1.gro

Caution

The number of added WT4 molecules, 850, may change according to the software version.

Remove misplaced WT4 molecules inside the bilayer:

gmx grompp -f sirah.ff/tutorial/5/CPU/em_CGLIP.mdp -p topol.top -c DMPC64_cg_sol1.gro -o DMPC64_cg_sol1.tpr
echo -e "a BE* BC1* BC2* BCT*\n name 7 tail\n\nq" | gmx make_ndx -f DMPC64_cg_sol1.gro -o DMPC64_cg_sol1.ndx
gmx select -f DMPC64_cg_sol1.gro -s DMPC64_cg_sol1.tpr -n DMPC64_cg_sol1.ndx -on rm_close_wt4.ndx -select "not (same residue as (resname WT4 and within 0.4 of group tail))"
gmx editconf -f DMPC64_cg_sol1.gro -o DMPC64_cg_sol2.gro -n rm_close_wt4.ndx

Note

Consult sirah.ff/0ISSUES and FAQs for information on known solvation issues.

Edit the [ molecules ] section in topol.top to correct the number of WT4 molecules:

Hint

If you forget to read the number of added WT4 molecules from the output of solvate, then use the following command line to get it

grep -c WP1 DMPC64_cg_sol2.gro

Add CG counterions and 0.15M NaCl:

gmx grompp -f sirah.ff/tutorial/5/GPU/em_CGLIP.mdp -p topol.top -c DMPC64_cg_sol2.gro -o DMPC64_cg_sol2.tpr
gmx genion -s DMPC64_cg_sol2.tpr -o DMPC64_cg_ion.gro -np 21 -pname NaW -nn 21 -nname ClW

When prompted, choose to substitute WT4 molecules by ions.

Note

The available electrolyte species in SIRAH force field are: Na⁺ (NaW), K⁺ (KW) and Cl⁻ (ClW) which represent solvated ions in solution. One ion pair (e.g., NaW-ClW) each 34 WT4 molecules results in a salt concentration of ~0.15M (see Appendix for details). Counterions were added according to Machado et al..

Edit the [ molecules ] section in topol.top to include the CG ions and the correct number of WT4.

Before running the simulation it may be a good idea to visualize your molecular system. CG molecules are not recognized by molecular visualizers and will not display correctly. To fix this problem you may generate a PSF file of the system using the script g_top2psf.pl:

./sirah.ff/tools/g_top2psf.pl -i topol.top -o DMPC64_cg_ion.psf

Note

This is the basic usage of the script g_top2psf.pl, you can learn other capabilities from its help:

./sirah.ff/tools/g_top2psf.pl -h

Use VMD to check how the CG system looks like:

vmd DMPC64_cg_ion.psf DMPC64_cg_ion.gro -e sirah.ff/tools/sirah_vmdtk.tcl

Tip

VMD assigns default radius to unknown atom types, the script sirah_vmdtk.tcl sets the right ones, according to the CG representation. It also provides a kit of useful selection macros, coloring methods and backmapping utilities. Use the command sirah_help in the Tcl/Tk console of VMD to access the manual pages. To learn about SIRAH Tools’ capabilities, you can also go to the SIRAH Tools tutorial.

5.4. Run the simulation

Important

By default in this tutorial we will use input files for GROMACS on GPU (sirah.ff/tutorial/5/GPU). Example input files for using GROMACS on CPU can be found at: sirah.ff/tutorial/5/CPU.

The folder sirah.ff/tutorial/5/GPU/ contains typical input files for energy minimization (em_CGLIP.mdp), equilibration (eq_CGLIP.mdp) and production (md_CGLIP.mdp) runs. Please check carefully the input flags therein.

Create an index file:

echo "q" | gmx make_ndx -f DMPC64_cg_ion.gro -o DMPC64_cg_ion.ndx

Note

WT4 and CG ions (NaW, ClW) are automatically set to the group SIRAH-Solvent while DMPC (named CMM at CG level) is assigned to group Lipid.

Make a new folder for the run:

mkdir run; cd run

Energy Minimization:

gmx grompp -f ../sirah.ff/tutorial/5/GPU/em_CGLIP.mdp -p ../topol.top -po em.mdp -n ../DMPC64_cg_ion.ndx -c ../DMPC64_cg_ion.gro -o DMPC64_cg_em.tpr
gmx mdrun -deffnm DMPC64_cg_em &> EM.log &

Equilibration:

gmx grompp -f ../sirah.ff/tutorial/5/GPU/eq_CGLIP.mdp -p ../topol.top -po eq.mdp -n ../DMPC64_cg_ion.ndx -c DMPC64_cg_em.gro -o DMPC64_cg_eq.tpr
gmx mdrun -deffnm DMPC_cg_eq &> EQ.log &

Production (100ns):

gmx grompp -f ../sirah.ff/tutorial/5/GPU/md_CGLIP.mdp -p ../topol.top -po md.mdp -n ../DMPC64_cg_ion.ndx -c DMPC64_cg_eq.gro -o DMPC64_cg_md.tpr
gmx mdrun -deffnm DMPC64_cg_md &> MD.log &

Note

GPU flags have been set for GROMACS 4.6.7; however, different versions may object to certain specifications.

5.5. Visualizing the simulation

That’s it! Now you can analyze the trajectory.

Process the output trajectory at folder run/ to account for the Periodic Boundary Conditions (PBC):

gmx trjconv -s DMPC64_cg_em.tpr -f DMPC64_cg_md.xtc -o DMPC64_cg_md_pbc.xtc -n ../DMPC64_cg_ion.ndx -pbc mol

When prompted, choose System for output.

Now you can check the simulation using VMD:

vmd ../DMPC64_cg_ion.psf ../DMPC64_cg_ion.gro DMPC64_cg_md_pbc.xtc -e ../sirah.ff/tools/sirah_vmdtk.tcl

Note

The file sirah_vmdtk.tcl is a Tcl script that is part of SIRAH Tools and contains the macros to properly visualize the coarse-grained structures in VMD. Use the command sirah-help in the Tcl/Tk console of VMD to access the manual pages. To learn about SIRAH Tools’ capabilities, you can also go to the SIRAH Tools tutorial.

Outside VMD, calculate the area per lipid:

gmx energy -f DMPC64_cg_md.edr -o box_XY.xvg

When prompted, choose Box-X and Box-Y for output. End your selection with zero.

Note

To calculate the area per lipid, divide the membrane’s area by the DMPC molecules per leaflet:

\[\frac{Area}{Lipid} = \frac{Box(x) * Box(y)}{64}\]

To calculate density profiles and bilayer thickness, include a new group for phosphate beads in the index file:

echo -e "a BFO\nq\n" | gmx make_ndx -f DMPC64_cg_em.gro -n ../DMPC64_cg_ion.ndx -o DMPC64_cg_ion.ndx
gmx density -sl 1000 -ng 5 -f DMPC64_cg_md_pbc.xtc -s DMPC64_cg_em.tpr -n DMPC64_cg_ion.ndx -o density_profile.xvg

When prompted, choose Lipid, WT4, NaW, ClW, and BFO.

Use Grace to plot the results:

xmgrace -nxy density_profile.xvg

Note

The thickness of the bilayer is the distance between the two peaks corresponding to the position of phosphate beads (BFO) along the z-axis.

6. Membrane proteins in explicit solvent

Note

Please report bugs, errors or enhancement requests through Issue Tracker or if you have a question about SIRAH open a New Discussion.

This tutorial shows how to use the SIRAH force field to perform a coarse grained (CG) simulation of a protein embedded in a lipid bilayer using explicit solvent (called WatFour, WT4). The main references for this tutorial are: Machado et al., Barrera et al., and Machado & Pantano. We strongly advise you to read these articles before starting the tutorial.

Note

We strongly advise you to read and complete Tutorial 3 and Tutorial 5 before starting.

Important

Check the Setting up SIRAH section for download and set up details before starting this tutorial. Since this is Tutorial 6, remember to replace X.X, the files corresponding to this tutorial can be found in: sirah.ff/tutorial/6/

6.1. Build CG representations

Caution

The mapping to CG requires the correct protonation state of each residue at a given pH. We recommend using the CHARMM-GUI server and use the PDB Reader & Manipulator to prepare your system. An account is required to access any of the CHARMM-GUI Input Generator modules, and it can take up to 24 hours to obtain one.

Other option is the PDB2PQR server and choosing the output naming scheme of AMBER for best compatibility. This server was utilized to generate the PQR file featured in this tutorial. Be aware that modified residues lacking parameters such as: MSE (seleno MET), TPO (phosphorylated THY), SEP (phosphorylated SER) or others are deleted from the PQR file by the server. In that case, mutate the residues to their unmodified form before submitting the structure to the server.

See Tutorial 3 for cautions while preparing and mapping atomistic proteins to SIRAH.

Map the atomistic structure of protein 2KYV to its CG representation:

./sirah.ff/tools/CGCONV/cgconv.pl -i sirah.ff/tutorial/6/2kyv.pqr -o 2kyv_cg.pdb

The input file -i 2kyv.pqr contains the atomistic representation of 2KYV structure at pH 7.0, while the output -o 2kyv_cg.pdb is its SIRAH CG representation.

Important

If you already have an atomistic protein within a membrane, then you can simply map the entire system to SIRAH (this is highly recommended) and skip the step of embedding the protein into a lipid bilayer, however clipping the membrane patch may be required to set a correct solvation box (see bellow). By default no mapping is applied to lipids, as there is no standard naming convention for them. So users are requested to append a MAP file from the list in Table 1, by setting the flag -a in cgconv.pl. We recommend using PACKMOL for building the system. Reference building-block structures are provided at folder sirah.ff/PDB/, which agree with the mapping scheme sirah.ff/tools/CGCONV/maps/tieleman_lipid.map. See FAQs for cautions on mapping lipids to SIRAH and tips on using fragment-based topologies.

Tip

This an advanced usage of the script cgconv.pl, you can learn other capabilities from its help by typing:

./sirah.ff/tools/CGCONV/cgconv.pl -h

6.1.1. Embed the protein in a lipid bilayer

We will show one possible way to do it by starting from a pre-equilibrated CG membrane patch.

In sirah_[version].ff/tutorial/6/ you will find a pre-stabilized CG DMPC bilayer patch, concatenate it with the previously generated CG representation of the phospholamban (PLN) pentamer (PDB code: 2KYV):

head -qn -1 2kyv_cg.pdb ./sirah.ff/tutorial/6/DMPC_cg.pdb > 2kyv_DMPC_cg_init.pdb

Luckily, we already oriented the protein inside the membrane. For setting up your own system you can go to Orientations of Proteins in Membranes (OPM) database and (if your structure is available) use the dummy atoms provided there to make them match with your membrane model (see Figure 1).

_images/Tuto7.png

Figure 1. Protein oriented inside the membrane from the OPM database with dummy atoms represented as orange spheres.

6.1.2. Delete close contact lipid molecules

We need to use VMD to to delete lipid molecules in close contact with the protein. For a proper treatment and visualization of the system in VMD you must first generate the molecular topology and initial coordinate files.

Use pdb2gmx to convert your PDB file into GROMACS format:

gmx pdb2gmx -f 2kyv_DMPC_cg_init.pdb -o 2kyv_DMPC_cg_init.gro -p init_topol -i init_posre

When prompted, choose SIRAH force field and then SIRAH solvent models.

Caution

In GROMACS versions prior to 5.x, the “gmx” command should not be used.

CG molecules are not recognized by molecular visualizers and will not display correctly. To fix this problem you may generate a PSF file of the system using the script g_top2psf.pl:

./sirah.ff/tools/g_top2psf.pl -i init_topol.top -o 2kyv_DMPC_cg_init.psf

Note

This is the basic usage of the script g_top2psf.pl, you can learn other capabilities from its help:

./sirah.ff/tools/g_top2psf.pl -h

Use VMD to check how the CG system looks like:

vmd 2kyv_DMPC_cg_init.psf 2kyv_DMPC_cg_init.gro -e sirah.ff/tools/sirah_vmdtk.tcl

Tip

VMD assigns default radius to unknown atom types, the script sirah_vmdtk.tcl sets the right ones, according to the CG representation. It also provides a kit of useful selection macros, coloring methods and backmapping utilities. Use the command sirah_help in the Tcl/Tk console of VMD to access the manual pages. To learn about SIRAH Tools’ capabilities, you can also go to the SIRAH Tools tutorial.

In the VMD main window, select Graphics > Representations. In the Selected Atoms box, type:

not (same residue as (sirah_membrane within 3.5 of sirah_protein))

To save the refined protein-membrane system, in the VMD main window click on 2kyv_DMPC_cg_init.psf, then select File > Save Coordinates. In the Selected atoms option choose the selection you have just created and Save as 2kyv_DMPC_cg.pdb.

From now on it is just normal GROMACS stuff!

6.2. PDB to GROMACS format

Use pdb2gmx to convert your PDB file into GROMACS format:

gmx pdb2gmx -f 2kyv_DMPC_cg.pdb -o 2kyv_DMPC_cg.gro

When prompted, choose SIRAH force field and then SIRAH solvent models.

Note

By default charged terminal are used but it is possible to set them neutral with option -ter

Note

Warning messages about long, triangular or square bonds are fine and expected due to the CG topology of some residues.

Caution

However, missing atom messages are errors which probably trace back to the mapping step. In that case, check your atomistic and mapped structures and do not carry on the simulation until the problem is solved.

6.3. Solvate the system

Define the simulation box of the system

gmx editconf -f 2kyv_DMPC_cg.gro -o 2kyv_DMPC_cg_box.gro -box 12 12 12 -c

Add WT4 molecules:

gmx solvate -cp 2kyv_DMPC_cg_box.gro -cs sirah.ff/wt416.gro -o 2kyv_DMPC_cg_solv1.gro

Note

Before GROMACS version 5.x, the command gmx solvate was called genbox.

Edit the [ molecules ] section in topol.top to include the number of added WT4 molecules:

Topology before editing

Topology after editing

[ molecules ]
; Compound #mols
Protein_chain_A 1
Protein_chain_B 1
Protein_chain_C 1
Protein_chain_D 1
Protein_chain_E 1
Lipid_chain_F 1

[ molecules ]
; Compound #mols
Protein_chain_A 1
Protein_chain_B 1
Protein_chain_C 1
Protein_chain_D 1
Protein_chain_E 1
Lipid_chain_F 1
WT4 3875

Hint

If you forget to read the number of added WT4 molecules from the output of solvate, then use the following command line to get it

grep -c WP1 2kyv_DMPC_cg_solv1.gro

Caution

The number of added WT4 molecules, 3875, may change according to the software version.

Remove misplaced WT4 molecules inside the bilayer:

gmx grompp -f sirah.ff/tutorial/6/GPU/em1_CGLIPROT.mdp -p topol.top -c 2kyv_DMPC_cg_solv1.gro -o 2kyv_DMPC_cg_solv1.tpr
echo "q" | gmx make_ndx -f 2kyv_DMPC_cg_solv1.gro -o 2kyv_DMPC_cg_solv1.ndx
gmx select -f 2kyv_DMPC_cg_solv1.gro -s 2kyv_DMPC_cg_solv1.tpr -n 2kyv_DMPC_cg_solv1.ndx -on rm_close_wt4.ndx -sf sirah.ff/tutorial/6/rm_close_wt4.dat
gmx editconf -f 2kyv_DMPC_cg_solv1.gro -o 2kyv_DMPC_cg_solv2.gro -n rm_close_wt4.ndx

Caution

New GROMACS versions may complain about the non-neutral charge of the system, aborting the generation of the TPR file by command grompp. We will neutralize the system later, so to overcame this issue, just allow a warning message by adding the following keyword to the grompp command line: -maxwarn 1

Note

Consult sirah.ff/0ISSUES and FAQs for information on known solvation issues.

Edit the [ molecules ] section in topol.top to correct the number of WT4 molecules:

Hint

If you forget to read the number of added WT4 molecules from the output of solvate, then use the following command line to get it

grep -c WP1 2kyv_DMPC_cg_solv2.gro

Add CG counterions and 0.15M NaCl:

gmx grompp -f sirah.ff/tutorial/6/GPU/em1_CGLIPROT.mdp -p topol.top -c 2kyv_DMPC_cg_solv2.gro -o 2kyv_DMPC_cg_solv2.tpr
gmx genion -s 2kyv_DMPC_cg_solv2.tpr -o 2kyv_DMPC_cg_ion.gro -np 95 -pname NaW -nn 110 -nname ClW

When prompted, choose to substitute WT4 molecules by ions.

Note

The available electrolyte species in SIRAH force field are: Na⁺ (NaW), K⁺ (KW) and Cl⁻ (ClW) which represent solvated ions in solution. One ion pair (e.g., NaW-ClW) each 34 WT4 molecules results in a salt concentration of ~0.15M (see Appendix for details). Counterions were added according to Machado et al..

Edit the [ molecules ] section in topol.top to include the CG ions and the correct number of WT4.

Before running the simulation it may be a good idea to visualize your molecular system. CG molecules are not recognized by molecular visualizers and will not display correctly. To fix this problem you may generate a PSF file of the system using the script g_top2psf.pl:

./sirah.ff/tools/g_top2psf.pl -i topol.top -o 2kyv_DMPC_cg_ion.psf

Note

This is the basic usage of the script g_top2psf.pl, you can learn other capabilities from its help:

./sirah.ff/tools/g_top2psf.pl -h

Use VMD to check how the CG system looks like:

vmd 2kyv_DMPC_cg_ion.psf 2kyv_DMPC_cg_ion.gro -e sirah.ff/tools/sirah_vmdtk.tcl

Tip

VMD assigns default radius to unknown atom types, the script sirah_vmdtk.tcl sets the right ones, according to the CG representation. It also provides a kit of useful selection macros, coloring methods and backmapping utilities. Use the command sirah_help in the Tcl/Tk console of VMD to access the manual pages. To learn about SIRAH Tools’ capabilities, you can also go to the SIRAH Tools tutorial.

6.4. Generate position restraint files

To achive a proper interaction between the protein and bilayer, we will perform a equilibration step applying restraints over the protein backbone and lipids’ phosphate groups.

Important

GROMACS enumerates atoms starting from 1 in each topology file (e.i. topol_*.itp). In order to avoid atom index discordance between the .gro file of the system and each topology, we will need independent .gro files for the PLN monomer and the DMPC bilayer.

Crate an index file of the system with a group for PLN monomer, then generate the .gro file of it.

Note

WT4 and CG ions (NaW, ClW) are automatically set to the group SIRAH-Solvent while DMPC (named CMM at CG level) is assigned to group Lipid.

echo -e "ri 1-52\nq" | gmx make_ndx -f 2kyv_DMPC_cg_ion.gro -o 2kyv_DMPC_cg_ion.ndx
gmx editconf -f 2kyv_DMPC_cg_ion.gro -n 2kyv_DMPC_cg_ion.ndx -o 2kyv_DMPC_cg_monomer.gro

When prompted, choose r_1-52.

Create an index file for the monomer and add a group for the backbone GO and GN beads.

echo -e "a GO GN \nq" | gmx make_ndx -f 2kyv_DMPC_cg_monomer.gro -o 2kyv_DMPC_cg_monomer.ndx

Create a position restraint file for the monomer backbone.

gmx genrestr -f 2kyv_DMPC_cg_monomer.gro -n 2kyv_DMPC_cg_monomer.ndx -fc 100 100 100 -o posre_BB.itp

When prompted, choose GO_GN.

Edit each topol_Protein_chain_*.itp (A to E) to include the new position restraints:

Topology before editing

Topology after editing

; Include Position restraint file
#ifdef POSRES
#include "posre_Protein.itp"
#endif




; Include Position restraint file
#ifdef POSRES
#include "posre_Protein.itp"
#endif

#ifdef POSREBB
#include "posre_BB.itp"
#endif

Use a similar procedure to set the positional restraints on lipid’s phosphates.

gmx editconf -f 2kyv_DMPC_cg_ion.gro -n 2kyv_DMPC_cg_ion.ndx -o DMPC_cg.gro

When prompted, choose Lipid.

Create an index file of the membrane and add a group for phosphates (BFO beads).

echo -e "a BFO \nq" | gmx make_ndx -f DMPC_cg.gro -o DMPC_cg.ndx

Create a position restraint file for phosphate groups in z coordinate.

gmx genrestr -f DMPC_cg.gro -n DMPC_cg.ndx -fc 0 0 100 -o posre_Pz.itp

When prompted, choose BFO.

Edit topol_Lipid_chain_F.itp to include the new position restraints and define the flags POSREZ to switch on these restraints in the input file (.mdp).

Topology before editing

Topology after editing

; Include Position restraint file
#ifdef POSRES
#include "posre_Lipid_chain_F.itp"
#endif




; Include Position restraint file
#ifdef POSRES
#include "posre_Lipid_chain_F.itp"
#endif

#ifdef POSREZ
#include "posre_Pz.itp"
#endif

6.5. Run the simulation

Important

By default in this tutorial we will use input files for GROMACS on GPU (sirah.ff/tutorial/6/GPU). Example input files for using GROMACS on CPU can be found at: sirah.ff/tutorial/6/CPU.

The folder sirah.ff/tutorial/6/GPU/ contains typical input files for energy minimization (em1_CGLIPROT.mdp and em2_CGLIPROT.mdp), equilibration (eq1_CGLIPROT.mdp and eq2_CGLIPROT.mdp) and production (md_CGLIPROT.mdp) runs. Please check carefully the input flags therein.

Make a new folder for the run:

mkdir run; cd run

Energy Minimization of side chains by restraining the backbone:

gmx grompp -f ../sirah.ff/tutorial/6/GPU/em1_CGLIPROT.mdp -p ../topol.top -n ../2kyv_DMPC_cg_ion.ndx -c ../2kyv_DMPC_cg_ion.gro -r ../2kyv_DMPC_cg_ion.gro -o 2kyv_DMPC_cg_em1.tpr
gmx mdrun -deffnm 2kyv_DMPC_cg_em1 &> EM1.log &

Energy Minimization of the whole system:

gmx grompp -f ../sirah.ff/tutorial/6/GPU/em2_CGLIPROT.mdp -p ../topol.top -n ../2kyv_DMPC_cg_ion.ndx -c 2kyv_DMPC_cg_em1.gro -r 2kyv_DMPC_cg_em1.gro -o 2kyv_DMPC_cg_em2.tpr
gmx mdrun -deffnm 2kyv_DMPC_cg_em2 &> EM2.log &

Equilibration 1:

Position restraints are defined in eq1_CGLIPROT.mdp file for protein backbone in xyz and phosphate groups (BFO beads) in z coordinate by setting keywords -DPOSREBB and -DPOSREZ, respectively.

gmx grompp -f ../sirah.ff/tutorial/6/GPU/eq1_CGLIPROT.mdp -p ../topol.top -n ../2kyv_DMPC_cg_ion.ndx -c 2kyv_DMPC_cg_em2.gro -r 2kyv_DMPC_cg_em2.gro -o 2kyv_DMPC_cg_eq1.tpr
gmx mdrun -deffnm 2kyv_DMPC_cg_eq1 &> EQ1.log &

Equilibration 2:

gmx grompp -f ../sirah.ff/tutorial/6/GPU/eq2_CGLIPROT.mdp -p ../topol.top -n ../2kyv_DMPC_cg_ion.ndx -c 2kyv_DMPC_cg_eq1.gro -r 2kyv_DMPC_cg_eq1.gro -o 2kyv_DMPC_cg_eq2.tpr
gmx mdrun -deffnm 2kyv_DMPC_cg_eq2 &> EQ2.log &

Production (1000ns):

gmx grompp -f ../sirah.ff/tutorial/6/GPU/md_CGLIPROT.mdp -p ../topol.top -n ../2kyv_DMPC_cg_ion.ndx -c 2kyv_DMPC_cg_eq2.gro -r 2kyv_DMPC_cg_eq2.gro -o 2kyv_DMPC_cg_md.tpr
gmx mdrun -deffnm 2kyv_DMPC_cg_md &> MD.log &

6.6. Visualizing the simulation

That’s it! Now you can analyze the trajectory.

Process the output trajectory at folder run/ to account for the Periodic Boundary Conditions (PBC):

gmx trjconv -s 2kyv_DMPC_cg_em1.tpr -f 2kyv_DMPC_cg_md.xtc -o 2kyv_DMPC_cg_md_pbc.xtc -n ../2kyv_DMPC_cg_ion.ndx -pbc mol

When prompted, choose System for output.

Now you can check the simulation using VMD:

vmd ../2kyv_DMPC_cg_ion.psf ../2kyv_DMPC_cg_ion.gro 2kyv_DMPC_cg_md_pbc.xtc -e ../sirah.ff/tools/sirah_vmdtk.tcl

Note

The file sirah_vmdtk.tcl is a Tcl script that is part of SIRAH Tools and contains the macros to properly visualize the coarse-grained structures in VMD. Use the command sirah-help in the Tcl/Tk console of VMD to access the manual pages. To learn about SIRAH Tools’ capabilities, you can also go to the SIRAH Tools tutorial.

7. SMD with multiscale solvent

Note

Please report bugs, errors or enhancement requests through Issue Tracker or if you have a question about SIRAH open a New Discussion.

This tutorial shows how to apply the multiscale solvation approach of SIRAH force field in Steered Molecular Dynamics (SMD) simulations. In this tutorial, an immunoglobulin domain known as I10 derived from Titin is enclosed by CG waters referred to as WatFour (WT4). These waters are further embedded within supra-coarse-grained (SCG) molecules called WatElse (WLS), which serve as a representation of bulk water. The main references for this tutorial are: Darré et al. and Machado et al..

Note

We strongly advise you to read and complete Tutorial 2 and Tutorial 3 before starting. We also recommend you to perform the Umbrella Sampling tutorial of GROMACS to get familiar with pulling simulations parameters.

Important

Check the Setting up SIRAH section for download and set up details before starting this tutorial. Since this is Tutorial 7, remember to replace X.X, the files corresponding to this tutorial can be found in: sirah.ff/tutorial/7/

7.1. Setting pulling direction

Caution

The basic idea behind SMD is to apply an external force or mechanical perturbation to a specific part of the biomolecular system under study and observe its response. This external force is applied to one atom or groups of atoms as a “steering” force, guiding or pulling the atoms along a predefined path or direction. Thus, prior to proceeding, we need to remember to prepare the PDB file to specify the pulling direction. In this example, the Z-axis was employed to orient the protein terminus and to determine the box size.

In sirah_[version].ff/tutorial/7/ you will find the I10 domain (PDB code: 4QEG) already aligned along the Z-axis, I10_z.pdb file.

_images/TutorialSMD1.png

Figure 1. I10 domain N-terminal (blue) and C-terminal (red) aligned to the Z-axis from I10_z.pdb.

To get the protein aligned along the Z-axis, we used the following commands in VMD’s Tk Console (Extensions > Tk Console):

#Select the protein
set all [atomselect top "protein"]

#Change protein coordinates to align to its center
$all moveby  [vecinvert [measure center $all]]

#Rotate around x-axis 50 degrees
$all move [transaxis x 50]

#Rotate around y-axis -11 degrees
$all move [transaxis y -11]

#Write PDB file with the protein aligned to z-axis
$all writepdb I10_z.pdb

Note

For setting up your own system you can open your own PDB in VMD and probe different alternatives to get the correct direction.

7.2. Build CG representations

Caution

The mapping to CG requires the correct protonation state of each residue at a given pH. We recommend using the CHARMM-GUI server and use the PDB Reader & Manipulator to prepare your system. An account is required to access any of the CHARMM-GUI Input Generator modules, and it can take up to 24 hours to obtain one.

Other option is the PDB2PQR server and choosing the output naming scheme of AMBER for best compatibility. This server was utilized to generate the PQR file featured in this tutorial. Be aware that modified residues lacking parameters such as: MSE (seleno MET), TPO (phosphorylated THY), SEP (phosphorylated SER) or others are deleted from the PQR file by the server. In that case, mutate the residues to their unmodified form before submitting the structure to the server.

See Tutorial 3 for cautions while preparing and mapping atomistic proteins to SIRAH.

Important

Check the Setting up SIRAH section for download and set up details before starting this tutorial. Since this is Tutorial 7, remember to replace X.X, the files corresponding to this tutorial can be found in: sirah.ff/tutorial/7/

Map the atomistic structure of the I10 domain to its CG representation:

./sirah.ff/tools/CGCONV/cgconv.pl -i sirah.ff/tutorial/7/I10_z.pqr -o I10_cg.pdb

The input file -i I10_z.pqr contains the atomistic representation of 4QEG structure at pH 7.0 and aligned to the Z-axis, while the output -o I10_cg.pdb is its SIRAH CG representation.

Tip

This is the basic usage of the script cgconv.pl, you can learn other capabilities from its help by typing:

./sirah.ff/tools/CGCONV/cgconv.pl -h

Please check both PDB and PQR structures using VMD:

vmd -m sirah.ff/tutorial/7/I10_z.pqr I10_cg.pdb

From now on it is just normal GROMACS stuff!

7.3. PDB to GROMACS format

Use pdb2gmx to convert your PDB file into GROMACS format:

gmx pdb2gmx -f I10_cg.pdb -o I10_cg.gro

When prompted, choose SIRAH force field and then SIRAH solvent models. In this specific case, the charge of the system is -5.000 e; this will be addressed later.

Note

By default charged terminal are used but it is possible to set them neutral with option -ter.

Note

Warning messages about long, triangular or square bonds are fine and expected due to the CG topology of some residues.

Caution

However, missing atom messages are errors which probably trace back to the mapping step. In that case, check your atomistic and mapped structures and do not carry on the simulation until the problem is solved.

7.4. Solvate the system

Danger

Since we are doing a SMD simulation, we need to carefully set the box dimensions to provide enough solvent to the “stretched” protein and keep a good separation between the surfaces of periodic images.

In this system, I10 has 88 amino acids and the maximum extension size of each amino acid is 0.34 nm. Thus, “stretched” protein will be 88 * 0.34 = 29.92 nm (~ 30.0 nm) long. In order to accommodate the pulling, GROMACS stipulates a minimum box size double this value, i.e. 60 nm for the Z-axis. However, for optimal results, it is recommended that the dimensions of the box be 2.5 to 3 times greater than the maximum length of the protein when in its extended conformation. Therefore, for this tutorial the box used is 10 10 90 nm.

_images/TutorialSMD2_1.png

Figure 2. Dimensions of the multiscale solvation box used in this tutorial.

In order to have a multiscale solvent approach using WT4 (CG solvent) and WLS (Supra-CG solvent), two steps are needed to solvate the systems.

First, define the simulation region of the system to be enclosed by WT4 (purple in Figure 2)

gmx editconf -f I10_cg.gro -o I10_cg_box.gro -box 10 10 20 -bt triclinic -c

Note

At this step, if you don’t want to use a multiscale solvent method, the whole box dimension (10 10 90 nm) can be used to add only WT4 molecules.

Add WT4 molecules:

gmx solvate -cp I10_cg_box.gro -cs sirah.ff/wt416.gro -o I10_cg_solv1.gro

Note

Before GROMACS version 5.x, the command gmx solvate was called genbox.

Edit the [ molecules ] section in topol.top to include the number of added WT4 molecules:

Topology before editing

Topology after editing

[ molecules ]
; Compound #mols
Protein_chain_A 1

[ molecules ]
; Compound #mols
Protein_chain_A 1
WT4 6281

Hint

If you forget to read the number of added WT4 molecules from the output of solvate, then use the following command line to get it:

grep -c WP1 I10_cg_solv1.gro

Caution

The number of added WT4 molecules, 6281, may change according to the software version.

Remove misplaced WT4 molecules within 0.3 nm of protein:

echo q | gmx make_ndx -f I10_cg_solv1.gro -o I10_cg_solv1.ndx
gmx grompp -f sirah.ff/tutorial/7/em1_CGPROT.mdp -p topol.top -po delete1.mdp -c I10_cg_solv1.gro -o I10_cg_solv1.tpr -maxwarn 2

Caution

New GROMACS versions may complain about not used macros and the non-neutral charge of the system, aborting the generation of the TPR file by command grompp. We will neutralize the system later, so to overcame these issues, just allow warning messages by adding the following keyword to the grompp command line: -maxwarn 2

gmx select -f I10_cg_solv1.gro -s I10_cg_solv1.tpr -n I10_cg_solv1.ndx -on rm_close_wt4.ndx -select 'not (same residue as (resname WT4 and within 0.3 of group Protein))'
gmx editconf -f I10_cg_solv1.gro -o I10_cg_solv2.gro -n rm_close_wt4.ndx

Edit the [ molecules ] section in topol.top to correct the number of WT4 molecules, 6261.

Hint

If you forget to read the number of added WT4 molecules from the output of solvate, then use the following command line to get it:

grep -c WP1 I10_cg_solv2.gro

Now, we include the second solvent layer of solvent with WLS molecules (green in Figure 2):

gmx editconf -f I10_cg_solv2.gro -o I10_cg_box2.gro -box 10 10 90 -bt triclinic -c

Hint

We can check the final box dimensions with VMD:

vmd I10_cg_box2.gro

In the Tk Console (Extensions > Tk Console) use the command:

pbc box

Add WLS molecules:

gmx solvate -cp I10_cg_box2.gro -cs sirah.ff/wlsbox.gro -o I10_cg_solv3.gro

Note

Before GROMACS version 5.x, the command gmx solvate was called genbox.

Edit the [ molecules ] section in topol.top to include the number of added WLS molecules:

Topology before editing

Topology after editing

[ molecules ]
; Compound #mols
Protein_chain_A 1
WT4 6261

[ molecules ]
; Compound #mols
Protein_chain_A 1
WT4 6261
WLS 4697

Hint

If you forget to read the number of added WLS molecules from the output of solvate, then use the following command line to get it:

grep -c LN1 I10_cg_solv3.gro

Caution

The number of added WLS molecules, 4697, may change according to the software version.

Remove misplaced WLS molecules within 7.8 nm of protein:

echo q | gmx make_ndx -f I10_cg_solv3.gro -o I10_cg_solv3.ndx
gmx grompp -f sirah.ff/tutorial/7/em1_CGPROT.mdp -p topol.top -po delete3.mdp -c I10_cg_solv3.gro -o I10_cg_solv3.tpr -maxwarn 2

Caution

New GROMACS versions may complain about not used macros and the non-neutral charge of the system, aborting the generation of the TPR file by command grompp. We will neutralize the system later, so to overcame this issue, just allow warning messages by adding the following keyword to the grompp command line: -maxwarn 2

gmx select -f I10_cg_solv3.gro -s I10_cg_solv3.tpr -n I10_cg_solv3.ndx -on rm_close_wls.ndx -select 'not (same residue as (resname WLS and within 7.8 of group Protein))'
gmx editconf -f I10_cg_solv3.gro -o I10_cg_solv4.gro -n rm_close_wls.ndx

Edit the [ molecules ] section in topol.top to correct the number of WLS molecules, 4580.

Hint

If you forget to read the number of added WLS molecules from the output of solvate, then use the following command line to get it

grep -c LN1 I10_cg_solv4.gro

Note

Consult sirah.ff/0ISSUES and FAQs for information on known solvation issues.

Add CG counterions and 0.15M NaCl:

gmx grompp -f sirah.ff/tutorial/7/em1_CGPROT.mdp -p topol.top -po delete4.mdp -c I10_cg_solv4.gro -o I10_cg_solv4.tpr -maxwarn 2

Caution

New GROMACS versions may complain about not used macros and the non-neutral charge of the system, aborting the generation of the TPR file by command grompp. We are about to neutralize the system, so to overcame this issue, just allow warning messages by adding the following keyword to the grompp command line: -maxwarn 2

gmx genion -s I10_cg_solv4.tpr -o I10_cg_ion.gro -np 189 -pname NaW -nn 184 -nname ClW

When prompted, choose to substitute WT4 molecules by ions.

Note

The available electrolyte species in SIRAH force field are: Na⁺ (NaW), K⁺ (KW) and Cl⁻ (ClW) which represent solvated ions in solution. One ion pair (e.g., NaW-ClW) each 34 WT4 molecules results in a salt concentration of ~0.15M (see Appendix for details). Counterions were added according to Machado et al..

Edit the [ molecules ] section in topol.top to include the CG ions and the correct number of WT4, WLS, and ions.

Topology before editing

Topology after editing

[ molecules ]
; Compound #mols
Protein_chain_A 1
WT4 6261
WLS 4582


[ molecules ]
; Compound #mols
Protein_chain_A 1
WT4 5888
NaW 189
ClW 184
WLS 4580

Caution

Following the above order is important: the number of WT4 comes first, then the number of ions, and finally WLS.

Before running the simulation it may be a good idea to visualize your molecular system. CG molecules are not recognized by molecular visualizers and will not display correctly. To fix this problem you may generate a PSF file of the system using the script g_top2psf.pl:

./sirah.ff/tools/g_top2psf.pl -i topol.top -o I10_cg_ion.psf

Note

This is the basic usage of the script g_top2psf.pl, you can learn other capabilities from its help:

./sirah.ff/tools/g_top2psf.pl -h

Use VMD to check how the CG system looks like:

vmd I10_cg_ion.psf I10_cg_ion.gro -e sirah.ff/tools/sirah_vmdtk.tcl

Tip

VMD assigns default radius to unknown atom types, the script sirah_vmdtk.tcl sets the right ones, according to the CG representation. It also provides a kit of useful selection macros, coloring methods and backmapping utilities. Use the command sirah_help in the Tcl/Tk console of VMD to access the manual pages. To learn about SIRAH Tools’ capabilities, you can also go to the SIRAH Tools tutorial.

To achive a proper interaction between the protein and solvent, we will perform a equilibration step applying restraints over the protein backbone.

Create an index file including a group for the backbone GN and GO beads:

echo -e "a GN GO\n\nq" | gmx make_ndx -f I10_cg_ion.gro -o I10_cg_ion.ndx

Note

WT4 and CG ions (NaW, ClW) are automatically set to the group SIRAH-Solvent.

Generate restraint files for the backbone GN and GO beads:

gmx genrestr -f I10_cg.gro -n I10_cg_ion.ndx -o bkbres.itp

When prompted, choose the group GN_GO.

gmx genrestr -f I10_cg.gro -n I10_cg_ion.ndx -o bkbres_soft.itp -fc 100 100 100

Here, choose the group GN_GO.

Add the restraints to topol.top:

Topology before editing

Topology after editing

; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif
; Include water topology
#include "../sirah.ff/solv.itp"

#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
; i funct fcx fcy fcz
1 1 1000 1000 1000
#endif








; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif
; Backbone restraints
#ifdef GN_GO
#include "bkbres.itp"
#endif
; Backbone soft restrains
#ifdef GN_GO_SOFT
#include "bkbres_soft.itp"
#endif
; Include water topology
#include "../sirah.ff/solv.itp"
#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
; i funct fcx fcy fcz
1 1 1000 1000 1000
#endif

7.5. Run the simulation

Important

By default in this tutorial we will use input files for GROMACS on GPU (sirah.ff/tutorial/7/).

The folder sirah.ff/tutorial/7/ contains typical input files for energy minimization (em1_CGPROT.mdp, em2_CGPROT.mdp, and em3_CGPROT.mdp), equilibration (eq1_CGPROT.mdp and eq2_CGPROT.mdp), production (md_CGPROT.mdp) and SMD (SMD_Force_CGPROT.mdp and SMD_Velocity_CGPROT.mdp) runs. Please check carefully the input flags therein.

Make a new folder for the run:

mkdir run; cd run

Energy Minimization of side chains and Sirah-solvent by restraining the backbone:

gmx grompp -f ../sirah.ff/tutorial/7/em1_CGPROT.mdp -p ../topol.top -po em1.mdp -n ../I10_cg_ion.ndx -c ../I10_cg_ion.gro -r ../I10_cg_ion.gro -o I10_cg_em1.tpr
gmx mdrun -deffnm I10_cg_em1 &> EM1.log &

Energy Minimization of side chains by restraining the backbone:

gmx grompp -f ../sirah.ff/tutorial/7/em2_CGPROT.mdp -p ../topol.top -po em2.mdp -n ../I10_cg_ion.ndx -c I10_cg_em1.gro -r I10_cg_em1.gro -o I10_cg_em2.tpr
gmx mdrun -deffnm I10_cg_em2 &> EM2.log &

Energy Minimization of whole system:

gmx grompp -f ../sirah.ff/tutorial/7/em3_CGPROT.mdp -p ../topol.top -po em3.mdp -n ../I10_cg_ion.ndx -c I10_cg_em2.gro -o I10_cg_em3.tpr
gmx mdrun -deffnm I10_cg_em3 &> EM3.log &

Solvent equilibration:

gmx grompp -f ../sirah.ff/tutorial/7/eq1_CGPROT.mdp -p ../topol.top -po eq1.mdp -n ../I10_cg_ion.ndx -c I10_cg_em2.gro -r I10_cg_em2.gro -o I10_cg_eq1.tpr
gmx mdrun -deffnm I10_cg_eq1 &> EQ1.log &

Soft equilibration to improve side chain solvation:

gmx grompp -f ../sirah.ff/tutorial/7/eq2_CGPROT.mdp -p ../topol.top -po eq2.mdp -n ../I10_cg_ion.ndx -c I10_cg_eq1.gro -r I10_cg_eq1.gro -o I10_cg_eq2.tpr
gmx mdrun -deffnm I10_cg_eq2 &> EQ2.log &

SMD Force or velocity:

Here, we need to modify the index file to add the “Pull” groups, before running SMD force or velocity simulations.

Copy the I10_cg_ion.ndx with a new name:

cp ../I10_cg_ion.ndx ../I10_cg_ion_pull.ndx

Open the I10_cg_ion_pull.ndx in any text editor and manually add to the end of the I10_cg_ion_pull.ndx file these two new groups:

Groups to add

[ pull1 ]
4
[ pull2 ]
437

In this tutorial we are going to run only the SMD Force simulation:

gmx grompp -f ../sirah.ff/tutorial/7/SMD_Force_CGPROT.mdp -p ../topol.top -po md.mdp -n ../I10_cg_ion_pull.ndx -c I10_cg_eq2.gro -o I10_cg_SMD_F.tpr
gmx mdrun -deffnm I10_cg_SMD_F &> SMD_F.log &

However, you can also run a SMD Velocity simulation:

gmx grompp -f ../sirah.ff/tutorial/7/SMD_Velocity_CGPROT.mdp -p ../topol.top -po md.mdp -n ../I10_cg_ion_pull.ndx -c I10_cg_eq2.gro -o I10_cg_SMD_V.tpr
gmx mdrun -deffnm I10_cg_SMD_V &> SMD_V.log &

Note

It is also possible to use the files you made to run a MD simulation of your system:

Production (1000ns):

gmx grompp -f ../sirah.ff/tutorial/7/md_CGPROT.mdp -p ../topol.top -po md.mdp -n ../I10_cg_ion.ndx -c I10_cg_eq2.gro -o I10_cg_md.tpr
gmx mdrun -deffnm I10_cg_md &> MD.log &

7.6. Visualizing the simulation

That’s it! Now you can analyze the trajectory.

GROMACS automatically creates plot results to the SMD simulations. The position versus simulation time is saved as I10_cg_md_pullx.xvg, and the force versus simulation time is saved as I10_cg_md_pullf.xvg.

You can plot the results using Xmgrace.

xmgrace I10_cg_md_pullx.xvg

For I10_cg_md_pullx.xvg, a plot similar to Figure 3 will appear:

_images/Fuerza_100pN.png

Figure 3. Position vs Time plot created by GROMACS from the SMD simulation.

In addition, you can process the output trajectory at folder run/ to account for the Periodic Boundary Conditions (PBC).

For SMD Force:

gmx trjconv -s I10_cg_em1.tpr -f I10_cg_SMD_F.xtc -o I10_cg_SMD_F_pbc.xtc -n ../I10_cg_ion_pull.ndx -ur compact -center -pbc mol

For SMD Velocity:

gmx trjconv -s I10_cg_em1.tpr -f I10_cg_SMD_V.xtc -o I10_cg_SMD_V_pbc.xtc -n ../I10_cg_ion_pull.ndx -ur compact -center -pbc mol

When prompted, choose Protein for centering and System for output.

Note

If you had also run a MD simulation, you could use the following commands to account for PBC:

gmx trjconv -s I10_cg_em1.tpr -f I10_cg_md.xtc -o I10_cg_md_pbc.xtc -n ../I10_cg_ion.ndx -ur compact -center -pbc mol

Now you can check the simulation using VMD:

vmd ../I10_cg_ion.psf ../I10_cg_ion.gro I10_cg_SMD_F_pbc.xtc -e ../sirah.ff/tools/sirah_vmdtk.tcl

Note

The file sirah_vmdtk.tcl is a Tcl script that is part of SIRAH Tools and contains the macros to properly visualize the coarse-grained structures in VMD. Use the command sirah-help in the Tcl/Tk console of VMD to access the manual pages. To learn about SIRAH Tools’ capabilities, you can also go to the SIRAH Tools tutorial.