MD
Tutorial
It is possible to use a pdb file as an input file for MD simulations with Gromacs, but it is better to first convert it to the gromos file type (*.gro). Original data in the pdb file is often incomplete, carbon bound hydrogens are generally omitted. The conversion program pdb2gmx will check every residue in the structure file against a database and add all hydrogens. In the conversion process it also creates a topology file, with all the connections between the atoms listed.
pdb2gmx can be used by simply typing it at the prompt. To get a list of available options type:
pdb2gmx -h
This shows a description of the program and the available options. Most of these options are not important here. For the conversion of the structure file 1AKI.pdb only the main options, -f, -p and -o, will be used for structure file input, topology file output and structure file output, respectively:
pdb2gmx -f 1AKI.pdb -p aki.top -o aki.gro
This will first show the description of the program and generate a warning. The program will ask to select a force field:
Select the Force Field: 0: GROMOS96 43a1 force field 1: GROMOS96 43b1 vacuum force field 2: GROMOS96 43a2 force field (improved alkane dihedrals) 3: GROMOS96 45a3 force field (Schuler JCC 2001 22 1205) 4: GROMOS96 53a5 force field (JCC 2004 vol 25 pag 1656) 5: GROMOS96 53a6 force field (JCC 2004 vol 25 pag 1656) 6: OPLS-AA/L all-atom force field (2001 aminoacid dihedrals) 7: [DEPRECATED] Gromacs force field (see manual) 8: [DEPRECATED] Gromacs force field with hydrogens for NMR 9: Encad all-atom force field, using scaled-down vacuum charges 10: Encad all-atom force field, using full solvent charges
Type 2 (for the official distribution forcefield) and press Enter. It will produce a lot of output to the screen. Just have a brief look at it and try to understand what is done.
When the program has finished, list the files in your directory ('ls'). There are four files now, the original pdb file, the generated topology (.top) and gromos (.gro) file and a file called posre.itp, which is created automatically together with the topology. For a description of these files follow this link.
Now have a look at the .gro file ('more aki.gro'). It looks a lot like the original pdb file, containing the same information regarding the positions of the atoms, but the layout is different, hydrogens have been added and units have been converted to nm.
Look at the topology file ('more aki.top'). This file contains the information on the atom names, types, masses and charges, as well as a description of bonds, angles, dihedrals, etc.
The structure is now complete (hydrogens have been added) and a topology file has been created. However, there may be local strain in the protein, due to the generation of the hydrogens, and bad Van der Waals contacts may exist, caused by particles that are too close. The strain has to be removed by energy minimization of the structure. This can be done with the program 'mdrun', which is the MD program. Mdrun uses a single .tpr file as input, which is generated by combining the topology (aki.top), structure (aki.gro) and parameter files (minim.mdp). To generate the .tpr file the program grompp has to be used.
A description of grompp can be obtained by giving the command:
grompp -h
The parameter file minim.mdp has to be copied to the personal subdirectory using the command 'cp' (usage cp [source/]file loc):
cp ~mdcourse/minim.mdp ./
The file can also be found here for people who don't have access to the network.
Check the contents of this file (more minim.mdp) and try to understand what it says.
Use grompp to generate the mdrun input file (.tpr) from the three files mentioned:
grompp -f minim.mdp -c aki.gro -p aki.top -o input.tpr
The file input.tpr has now been created (to verify this, type 'ls'). This is a so called binary file. That means that it is written in computer language and is unreadable to humans. This file contains all information that is needed to perform the run. Have a look at the options for the mdrun program:
mdrun -h
The -s option is obligatory and followed by the name of the run input file. Output is generated by specifying the -o, -c and -e options. This will create a trajectory file, a structure file of the minimized structure and an energy file, with all the energies for every step. Perform the run using the following command (the -v option is used to output all messages to the screen and the option -nice sets the priority level):
mdrun -nice 0 -v -s input.tpr -o minim_traj.trr -c minimized.gro -e minim_ener.edr
The energy minimization may take some time, depending on the CPU in and the load of the computer. The trajectory file is not very important in energy minimizations, but the generated structure file (minimized.gro) will serve as input for the simulation.
During the minimization the potential energy decreases. A plot from the energy over time can be made from the minim_ener.edr file using g_energy. This program also has several options:
g_energy -h
Simply make a plot from the .edr file by executing:
g_energy -f minim_ener.edr -o minim_ener.xvg
This will display something like the following:
Opened minim_ener.edr as single precision energy file
Select the terms you want from the following list
G96Angle Proper-Dih. Improper-Dih. LJ-14 Coulomb-14 LJ-(SR) LJ-(LR) Coulomb-(SR) Coulomb-(LR) RF-excl. Potential Kinetic-En. Total-Energy Temperature Pressure-(bar) Box-X Box-Y Box-Z Volume Density-(SI) pV Vir-XX Vir-XY Vir-XZ Vir-YX Vir-YY Vir-YZ Vir-ZX Vir-ZY Vir-ZZ Pres-XX-(bar) Pres-XY-(bar) Pres-XZ-(bar) Pres-YX-(bar) Pres-YY-(bar) Pres-YZ-(bar) Pres-ZX-(bar) Pres-ZY-(bar) Pres-ZZ-(bar) #Surf*SurfTen Pcoupl-Mu-XX Pcoupl-Mu-YY Pcoupl-Mu-ZZ Mu-X Mu-Y Mu-Z T-System Lamb-System
Select the property you want by typing the name, e.g. Potential, which codes for potential energy and then press return and another return to quit.
The program g_energy produces a .xvg graph, which can be viewed and edited with xmgrace:
xmgrace -nxy minim_ener.xvg
Help on this program can be found at http://plasma-gate.weizmann.ac.il/Grace/doc/UsersGuide.html. The current graph has an error in it, the label for the horizontal axis should be 'step' instead of 'time (ps)'. Correct this by going to 'plot', select Tick labels/tick marks and in Axis label delete Time (ps) and write Step. Then click Accept and Close. To exit xmgrace go to File and select Exit.
Performing a full MD run is essentially the same as an energy
minimization. The structure file generated in the previous run is
needed (minimized.gro), as well as the original topology file
(aki.top).
First, the boundary condition of the system should be defined (box and
solvent).
editconf can be use to create the box :
editconf -h
Use a standard cubic box and be sure that the protein does not see
its images.
This means that the distance from the the box wall should be greater
than half of the cut-off (1.4 nm)
editconf -f minimized.gro -o minimized_box.gro -d 0.75
-bt cubic
Then solvate the protein. The program genbox will do it for you.
genbox -h
Solvent molecule are removed from the box where the distance
between any atom of the protein and the solvent molecules is less than
the sum of Van der Waals radii. As solvent we will use water. The
water model commonly used is SPC (Simple Point Charge) water.
cp ~mdcourse/spc216.gro ./
cp ~mdcourse/spc.itp ./
genbox -cp minimized_box.gro -cs spc216.gro -o minimized_water.gro -p aki.top
The topology file (aki.top) is updated
automatically by genbox: the topology file of the water (spc.itp)
is include and the number of water molecules is added.
Take a look at the new version of aki.top (less aki.top)
Now the boundary conditions are defined. The simulation can be
started. We run a short simulation of 5 ps.
The parameter file (fullmd_sol.mdp)
which can be downloaded here, or copied
using:
cp ~mdcourse/fullmd_sol.mdp ./
This file has information on settings to use for the simulations: the duration, writing of output files, treatment of long-range interactions and thermodynamic quantities, initial velocities and constraints on the protein. Take a look at the file (less fullmd_sol.mdp).
The preprocessor grompp is used again to combine the input files to one run input file (.tpr):
grompp -f fullmd_sol.mdp -c minimized_water.gro -p aki.top -o fullmd.tpr
After grompp has finished the simulation can begin:
mdrun -nice 0 -v -s fullmd.tpr -o md_traj.trr -x md_traj.xtc -c md_final.gro -e md_ener.edr
Most of the options have also been used for the energy minimization, except for -x. This will generate a trajectory file in the .xtc format, which contains information on the trajectories of the atoms in a compact manner (only containing cartesian coordinates).
While the simulation is running, use VMD to look the protein
solvated in water (minimized_water.gro).
To view the system,
select the option 'File' in the main menu. Then select 'New
Molecule'. Type minimized_water.gro in filename section. Finally click
on 'Load'. Try to visualize the solvent accessible surface.
Notice the box wall dimension and the number of water molecules around
the protein.
Now we try to visualize the trajectory of the system. Select again the option 'File' in the main menu and go to 'New Molecule'. Set the 'Load files for' to minimized_water.gro, then type 'md_traj.xtc' in the 'Filename' section and click on 'Load'.
VMD is now reading your trajectory. Use the arrows at the bottom of
the main menu to start, to stop and to change the speed of your animation.
To visualize better your trajectory, to go 'Graphics' select
'Representations'. First type in the 'Selected Atoms': 'protein' and
then enter. Only the protein will be visualized. Then,
click on 'Create Rep' and type in the 'Selected Atoms' : 'water and
within 3 of protein' and then enter. This command will display only
the water molecules around the protein within a distance of 3
Å. Choice the 'orthographic' option in the 'Display' menu
and then zoom using 'Scale mode' in the 'Mouse' menu to see better the
details.
Notice how the protein and the water around moves. Visualize the
hydrogen bond using the 'Draw style' menu
Try to explicitly animate the residues involved in the active site
(Glu 35, Asp 52, Gln 57 and Trp 108). Also try to show the hydrogen
bonds during the animation. VMD has many more features and it is
encouraged to explore the program and the structure by trying
some. Create a nice movie of your simulation, using Draw style
menu
('Selected atoms: Backbone' and 'Draw style Cartoon'). Remember that
you can save your work.
Now the simulation is finished. Take a look what happens to the
total energy.
Use g_energy to check the total energy of your
systems.
g_energy -f md_ener.edr -o md_ener_5ps.xvg
Select 12 and type 0 to quit. It is also possible to select more
than one options, e.g. 10 11 12 0 <enter>
will generate a graph of the potential, kinetic and total energy.
xmgrace -nxy md_ener_5ps.xvg
Notice what happens to the total energy in the first 5 ps.
Go to the next section to know more about the simulation analysis.
©2005 T.A.Wassenaar