Embedded Files
Step 1 (Preparation of the dimer)
Step 1 (Preparation of the dimer)
We need to download the protein structure that we will work with. We are going to use 7NT3. The protein structure can be downloaded at https://www.rcsb.org/.
The next step is we need to remove any water molecules from the protein. To remove water molecules (“HOH”), we can either remove them using a text editor such as Notepad++ or using the command:
>!grep -v HOH 7NT3.pdb > 7NT3_clean.pdb
PDB file is ready to be input into the first GROMACS module, pdb2gmx. The purpose of pdb2gmx is to generate topology for the molecules, a position restraint file and a post-processed structure file with the command:
> gmx pdb2gmx -f 7nt3_clean.pdb -o 7nt3_processed.gro -water spce
We will be prompted to choose a force field:
Choose 8 at the command prompt, followed by ENTER.
STEP 2 (EXAMINE THE TOPOLOGY)
STEP 2 (EXAMINE THE TOPOLOGY)
Next, we need to look at the topol.top file produced. We can use text editor to examine the file such as Notepad ++. In the file, we will find:
; Include forcefield parameters
#include "charmm27.ff/forcefield.itp"
There are also ions parameters, as follows:
; Include topology for ions
#include "charmm27.ff/ions.itp".
Finally, the [system] directive which gives the name of the system that will be written to output files during simulation. The [molecules] list all the molecules in the system.
[ system ]
; Name
Protein
[ molecules ]
; Compound #mols
Protein_chain_A 1
Protein_chain_B 1
step 3 ( DEFINING THE BOX AND ADDING SOLVENT)
step 3 ( DEFINING THE BOX AND ADDING SOLVENT)
Now, we are going to continue building our system in a simple aqueous system. The first step is to define the box dimensions using the editconf module. We will use a cubic box as a unit cell with the command:
> gmx editconf -f 7nt3_processed.gro -o 7nt3_newbox.gro -c -d 1.0 -bt cubic
The command center the protein in the box , (-c) and places it at least 1nm from the box edges (-d 1.0). The box type is cube.
After that, we will need to fill it with solvent (water) with the command:
> gmx solvate -cp 7nt3_newbox.gro -cs spc216.gro -o 7nt3_solv.gro -p topol.top
Note that, there will be changes in the [molecules]:
[ molecules ]
; Compound #mols
Protein_chain_A 1
Protein_chain_B 1
SOL 34773
STEP 4 (ADDING IONS)
STEP 4 (ADDING IONS)
We already have solvated system that contain a charged protein. The tools for adding ions is called genion. It will replace the water molecules with the ions that user specifies by reading the topology. The input is called a run input file, which has an extension of .tpr which also used to run our first simulation. The tpr file contains all the parameter of all the atoms in the system.
To produce .tpr file, we need an another file which is .mdp file. We will use a .mdp file which is ions.mdp that can be downloaded here : ions.mdp
The file can be edited using Notepad ++
Assemble the .tpr file with command :
> gmx grompp -f ions.mdp -c 7nt3_solv.gro -p topol.top -o ions.tpr
Now, we have ions.tpr file. We will pass this file to genion:
>gmx genion -s ions.tpr -o 7nt3_solv_ions.gro -p topol.top -pname NA -nname CL -neutral
Choose group "13" as we want to replace SOL to ions:
Our [molecules] directive now will change to :
[ molecules ]
; Compound #mols
Protein_chain_A 1
Protein_chain_B 1
SOL 34765
NA 8
STEP 5 (ENERGY MINIMIZATION)
STEP 5 (ENERGY MINIMIZATION)
Next, we must ensure that the system do not have steric clashes or inappropriate geometry. The structure is relaxed through a process called energy minimization (EM). The process of EM is much like addition of ions.
Assembly the binary input using grompp with the file:
Then, run the command:
> gmx grompp -f minim.mdp -c 1AKI_solv_ions.gro -p topol.top -o em.tpr
We are now ready to invoke mdrun to carry out the EM:
> gmx mdrun -v -deffnm em
From that, we will get the following files:
em.log
em.edr
em.trr
em.gro
We will do an analysis. The em.edr file contain all of the energy terms that GROMACS collects during EM. We can analyze that using GROMACS energy module:
> gmx energy -f em.edr -o potential.xvg
At the prompt, we will need to select "10 0" which is Potential (10) and zero (0) terminates input. To plot the data, we will need Xmgrace plotting tool.
Now that our system is at an energy minimum, we can begin real dynamics.
STEP 6 (EQUILIBRATION)
STEP 6 (EQUILIBRATION)
EM ensure that we have a solid starting structure. Nonetheless, to begin the dynamics, we must equilibrate the solvents and ions around protein. If we attempt to unrestrained dynamics here, the system may collapse.
Equilibration is often conducted in two phases. The first phase is conducted under NVT ensemble (constant number of particles, volume, and temperature). Usually, 50-100 ps is sufficient.
The nvt file:
We will call grompp and mdrun just as EM:
> gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr
> gmx mdrun -deffnm nvt
Now, we will analyze the temperature progression, again using energy:
> gmx energy -f nvt.edr -o temperature.xvg
Type "17 0" to choose temperature and press enter.
The previous step, NVT equilibration, stabilized the temperature of the system. Now, we will need to stabilize the pressure and density of the system. It is conducted by NPT ensemble, that the number of particles, pressure and temperature are constant.
The .mdp file can be found here:
We will call grompp and run as we did with NVT. In order to conserve the velocities produced during NVT, we must include this file.
The command:
> gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr
> gmx mdrun -deffnm npt
Then, we will analyze the pressure progression:
> gmx energy -f npt.edr -o pressure.xvg
Type 19 0 to choose pressure and press enter.
Then, we will look at the density as well:
> gmx energy -f npt.edr -o density.xvg
Type "24 0" to choose density.
step 7 (md production)
step 7 (md production)
After completed two equilibration, we will now ready to release the position restraints and run production MD for data collection. We will run the simulation with the script:
This is 1-ns script but in this tutorial, we will run 20-ns of MD simulation.
The command:
> gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md_0_1.tpr
Now, execute MD run
> gmx mdrun -deffnm md_0_1
STEP 8 (ANALYSIS)
STEP 8 (ANALYSIS)
We already run simulation. Now, we will need to do some analysis. The first is trjconv which is used as a post-processing tool to strip out the coordinates. The protein will diffuse through the unit cell, and may appear broken or jump across the other side of the box. To fix such behaviors, issue the following:
> gmx trjconv -s md_0_1.tpr -f md_0_1.xtc -o md_0_1_noPBC.xtc -pbc mol -center
Select 1 ("Protein") as the group to be centered and 0 ("System") for output.
To use rms, use this command:
> gmx rms -s md_0_1.tpr -f md_0_1_noPBC.xtc -o rmsd.xvg -tu ns
Choose 4 which is for both the least-squares fit and the group for RMSD calculation.
If we wish to calculate RMSD relative to the crystal structure, we could issue the following:
> gmx rms -s em.tpr -f md_0_1_noPBC.xtc -o rmsd_xtal.xvg -tu ns
Choose 4 which is for both the least-squares fit and the group for RMSD calculation.
Let's analyze the radius of gyration for lysozyme in our simulation:
> gmx gyrate -s md_0_1.tpr -f md_0_1_noPBC.xtc -o gyrate.xv
Choose group 1 (Protein) for analysis.
Google Sites
Report abuse