| blob | 427a435f9e4af27276443e5d6ad79c33dbfb59c7 |
1 #!/bin/csh
3 # this is the demo molecule
4 setenv MOL cpeptide
6 ####################
7 ### INTRODUCTION ###
8 ####################
9 clear
11 -----------------------------------------------------------------
12 -----------------------------------------------------------------
13 This is the GROMACS demo. ( This demo takes about 15 to 30 min. )
15 In this demo we will demonstrate you how to simulate Molecular
16 Dynamics (MD) using the GROMACS software package.
18 This demo will perform a complete molecular dynamics (MD) simulation
19 of a small peptide in water. The only input file we need to do this is
20 a pdb file of this peptide.
22 If you have a problem or remark with respect to this demonstration,
23 please mail to: gromacs@chem.rug.nl
24 -----------------------------------------------------------------
25 -----------------------------------------------------------------
26 _EOF_
30 #########################
31 ### CHECK ENVIRONMENT ###
32 #########################
33 clear
35 -----------------------------------------------------------------
36 -----------------------------------------------------------------
37 Before we you can actually start the GROMACS demo, I am going to check
38 if GROMACS is installed properly on your system. If GROMACS is not
39 installed properly on your computer, I strongly recommend you to
40 contact your system manager.
41 -----------------------------------------------------------------
42 -----------------------------------------------------------------
43 _EOF_
48 ###############
49 ### PDB2GMX ###
50 ###############
51 clear
53 -----------------------------------------------------------------
54 -----------------------------------------------------------------
55 Before we can start any simulation we need a molecular toplogy
56 file. This topology file ( .top extension ) is generated by the
57 program pdb2gmx. The only input file of the pdb2gmx program is the pdb
58 file of our peptide ( .pdb extension ).
60 Because most pdb files do not contain all hydrogen atoms, the pdb2gmx
61 program will also add them to our peptide. The output file which
62 contains the structure of the peptide when hydrogen atoms are added is a
63 gromos structure file ( .gro extension )
65 -----------------------------------------------------------------
66 -----------------------------------------------------------------
67 _EOF_
72 endif
80 endif
86 ##############
87 ### GENBOX ###
88 ##############
89 clear
91 -----------------------------------------------------------------
92 -----------------------------------------------------------------
93 Because a simulation of a peptide in vacua is a bit unrealistic, we
94 have to solvate our peptide in a box of water. genbox is the program
95 we use to do this.
97 The genbox program reads the peptide structure file and an input file
98 containing the sizes of the desired water box. The output of genbox is
99 a gromos structure file of a peptide solvated in a box of water. The
100 genbox program also changes the topology file ( .top extension ) to
101 include water. First we will use the program editconf to define the
102 right boxsize for our system.
104 -----------------------------------------------------------------
105 -----------------------------------------------------------------
106 _EOF_
111 endif
119 endif
132 ##############
133 ### GROMPP ###
134 ##############
135 clear
137 -----------------------------------------------------------------
138 -----------------------------------------------------------------
139 In principle we can start a Molecular Dynamics simulation now. However
140 it is not very wise to do so, because our system is full of close
141 contacts. These close contacts are mainly a result of the genbox
142 program. The added solvent might have some close contacts with the
143 peptide resulting in very high repulsive energies. If we would start a
144 Molecular Dynamics (MD) simulation without energy minimisation the
145 system would not be stable because of these high energies.
147 The standard procedure to remove these close contacts is
148 Energy Minimisation (EM). Energy minimisation slightly changes the
149 coordinates of our system to remove high energies from our system.
151 Before we can start the Energy Minimisation we have to preprocess all
152 the input files with the GROMACS preprocessor named grompp. grompp
153 preprocesses the topology file (.top), the structure file (.gro) and a
154 parameter file (.mdp) resulting in a binary topology file (.tpb/.tpr
155 extension). This binary topology file contains all information for a
156 simulation (in this case an energy minimisation).
157 -----------------------------------------------------------------
158 -----------------------------------------------------------------
159 _EOF_
164 endif
172 cpp = /lib/cpp
173 define = -DFLEX_SPC
174 constraints = none
175 integrator = steep
176 tinit = 0.0
177 dt = 0.002 ; ps !
178 nsteps = 100
179 nstlist = 10
180 ns_type = grid
181 deltagrid = 2
182 rlist = 1.0
183 rcoulomb = 1.0
184 rvdw = 1.0
185 ;
186 ; Energy minimizing stuff
187 ;
188 emtol = 1000.0
189 emstep = 0.01
190 _EOF_
195 endif
202 ################
203 ### MDRUN EM ###
204 ################
205 clear
207 -----------------------------------------------------------------
208 -----------------------------------------------------------------
209 Now the binary topology file is generated, we can start the energy
210 minimisation (EM). The program which performs the EM is called
211 mdrun. In fact all simulations are performed by the same program:
212 mdrun.
214 As the Energy Minimisation is running, watch the output of the
215 program. The first number ( from left to right ) is the number of the
216 iteration step. The second number is the step size, which gives an
217 indication of the change in the system. The third number is the
218 potential energy of the system. This number starts at a high value and
219 rapidly drops down, and converges, to a large negative value.
220 -----------------------------------------------------------------
221 -----------------------------------------------------------------
222 _EOF_
227 endif
235 endif
242 #################
243 ### GROMPP PR ###
244 #################
245 clear
247 -----------------------------------------------------------------
248 -----------------------------------------------------------------
249 Once all close contacts are removed from the system, we want to do
250 molecular dynamics of the water molecules, and keep the peptide
251 fixed. This is called position restrained (PR) MD.
253 Position Restrained MD keeps the peptide fixed and lets all water
254 molecules equilibrate around the peptide in order to fill holes
255 etc. which were not filled by the genbox program.
257 We are first going to preprocess the input files to generate the
258 binary topology. The input files are the topology file, the structure
259 file ( output of the EM ) a paremeter file, and an index file.
261 By default our system is split into several groups. In this case we
262 use two of those groups: Protein and SOL(vent). We use these groups to
263 put position restraints on all the atoms of the peptide.
265 The parameter file ( .mdp extension ) contains all information about
266 the PR-MD like: step size, number of steps, temperature, etc. This
267 Paramter file also tells GROMACS what kind of simulation should be
268 performed ( like EM, PR-MD and MD etc. )
269 -----------------------------------------------------------------
270 -----------------------------------------------------------------
271 _EOF_
276 endif
284 warnings = 10
285 cpp = /lib/cpp
286 define = -DPOSRE
287 constraints = all-bonds
288 integrator = md
289 tinit = 0.0
290 dt = 0.002 ; ps !
291 nsteps = 100 ; total 0.2 ps.
292 nstcomm = 1
293 nstxout = 10
294 nstvout = 1000
295 nstfout = 0
296 nstlog = 10
297 nstenergy = 10
298 nstlist = 10
299 ns_type = grid
300 deltagrid = 2
301 rlist = 1.0
302 rcoulomb = 1.0
303 rvdw = 1.0
304 box = rectangular
305 ; Temperature coupling is on in two groups
306 Tcoupl = yes
307 tau_t = 0.1 0.1
308 tc-grps = protein sol
309 ref_t = 300 300
310 ; Pressure coupling is not on
311 Pcoupl = no
312 tau_p = 0.5
313 compressibility = 4.5e-5
314 ref_p = 1.0
315 ; Generate velocites is on at 300 K.
316 gen_vel = yes
317 gen_temp = 300.0
318 gen_seed = 173529
319 _EOF_
325 endif
332 ################
333 ### MDRUN PR ###
334 ################
335 clear
337 -----------------------------------------------------------------
338 -----------------------------------------------------------------
339 Now we start the Position restrained Molecular Dynamics simulation. It
340 is important to note that in this example the simulated time is too
341 short to equilibrate our system completely, but that would simple take
342 too much time. ( about one day ).
344 -----------------------------------------------------------------
345 -----------------------------------------------------------------
346 _EOF_
351 endif
359 endif
366 #################
367 ### GROMPP MD ###
368 #################
369 clear
371 -----------------------------------------------------------------
372 -----------------------------------------------------------------
373 Now our complete system is finally ready for the actual Molecular
374 Dynamics simulation. We start again by preprocessing the input files
375 by the grompp program to generate the binary topology file (.tpb/.tpr
376 extension).
378 -----------------------------------------------------------------
379 -----------------------------------------------------------------
380 _EOF_
385 endif
393 cpp = /lib/cpp
394 constraints = all-bonds
395 integrator = md
396 tinit = 0.0
397 dt = 0.002 ; ps !
398 nsteps = 250 ; total 0.5 ps.
399 nstcomm = 1
400 nstxout = 1
401 nstvout = 0
402 nstfout = 0
403 nstlist = 10
404 ns_type = grid
405 deltagrid = 2
406 rlist = 1.0
407 rcoulomb = 1.0
408 rvdw = 1.0
409 box = rectangular
410 ; Temperature coupling is on in two groups
411 Tcoupl = yes
412 tau_t = 0.1 0.1
413 tc-grps = protein sol
414 ref_t = 300 300
415 ; Pressure coupling is not on
416 Pcoupl = no
417 tau_p = 0.5
418 compressibility = 4.5e-5
419 ref_p = 1.0
420 ; Generate velocites is on at 300 K.
421 gen_vel = yes
422 gen_temp = 300.0
423 gen_seed = 173529
424 _EOF_
429 endif
436 ################
437 ### MDRUN MD ###
438 ################
439 clear
441 -----------------------------------------------------------------
442 -----------------------------------------------------------------
443 Now we can start the MD simualtion. Watch the number of steps
444 increasing ( the total number of steps is 250 ).
446 -----------------------------------------------------------------
447 -----------------------------------------------------------------
448 _EOF_
453 endif
461 endif
468 ############
469 ### NGMX ###
470 ############
471 clear
473 -----------------------------------------------------------------
474 -----------------------------------------------------------------
475 We are finished simulating, and we are going to view the calculated
476 trajectory. The trajectory file ( .trj extension ) contains all
477 coordinates, velocities and forces of all the atoms in our system.
479 The next program we are going run is ngmx. ngmx is a very simple
480 trajectory viewer.
482 Once the ngmx program has been started you need to click on a few
483 buttons to view your trajectory.
485 1. Once the program has been started a dialog box shows up. Click on
486 the box on the left of the word Protein. ( This means that we want to
487 view the peptide ). Then Click on the OK Button
489 2. Now we see the edges of the box with a lines drawing of the peptide
490 we just simulated.
492 3. Select Animation in the Display menu. If you did this correctly. A
493 dialog box at the bottom of the screen appears. This dialog box is
494 used to move through your trajectory.
496 4. Click on the FastForward button (two triangles pointing to the
497 right) and watch the peptide moving.
498 -----------------------------------------------------------------
499 -----------------------------------------------------------------
500 _EOF_
503 echo Starting Trajectory viewer...
505 endif
506 #last line