share/tutor/gmxdemo/demo

   1 #!/bin/csh
   2
   3 # this is the demo molecule
   4 setenv MOL cpeptide
   5
   6 ####################
   7 ### INTRODUCTION ###
   8 ####################
   9 clear
  10 cat << _EOF_
  11 -----------------------------------------------------------------
  12 -----------------------------------------------------------------
  13 Welcome to the GROMACS demo.
  14
  15 This is a script that takes about 10 min to run.
  16
  17 In this demo we will demonstrate how to simulate Molecular
  18 Dynamics (MD) using the GROMACS software package.
  19
  20 The demo will perform a complete molecular dynamics (MD) simulation
  21 of a small peptide in water. The only input file we need to do this
  22 is a pdb file of a small peptide.
  23
  24 If you have any problems or remarks with respect to this demonstration,
  25 please mail to: gromacs@gromacs.org , or check the resources on
  26 our website http://www.gromacs.org .
  27 -----------------------------------------------------------------
  28 -----------------------------------------------------------------
  29 _EOF_
  30 echo -n "Press <enter>"
  31 set  ans = $<
  32
  33 #########################
  34 ### CHECK ENVIRONMENT ###
  35 #########################
  36 clear
  37 cat << _EOF_
  38 -----------------------------------------------------------------
  39 -----------------------------------------------------------------
  40 Before we you can actually start the GROMACS demo, the programs
  41 must be present in your PATH. This might already be the case if
  42 they are linked to /usr/local/bin. If not, follow the instructions
  43 in the getting started section. If GROMACS is not installed
  44 properly on your computer, contact your system manager.
  45 -----------------------------------------------------------------
  46 -----------------------------------------------------------------
  47 _EOF_
  48 echo -n "Press <enter>"
  49 set  ans = $<
  50
  51
  52 ###############
  53 ### PDB2GMX ###
  54 ###############
  55 clear
  56 cat << _EOF_
  57 -----------------------------------------------------------------
  58 -----------------------------------------------------------------
  59 Before we can start any simulation we need a molecular toplogy
  60 file. This topology file ( .top extension ) is generated by the
  61 program pdb2gmx. The only input file of the pdb2gmx program is the pdb
  62 file of our peptide ( .pdb extension ).
  63
  64 Because most pdb files do not contain all hydrogen atoms, the pdb2gmx
  65 program will also add them to our peptide. The output file which
  66 contains the structure of the peptide when hydrogen atoms are added is a
  67 gromos structure file ( .gro extension )
  68
  69 -----------------------------------------------------------------
  70 -----------------------------------------------------------------
  71 _EOF_
  72
  73 if ( $?DISPLAY ) then
  74         echo "You seem to have the DISPLAY variable is set, so we will"
  75         echo "pop up a window with the output of the pdb2gmx program"
  76 endif
  77 echo -n "Press <enter>"
  78 set  ans = $<
  79
  80
  81 echo "Starting pdb2gmx"
  82 if ( $?DISPLAY ) then
  83         xterm -title pdb2gmx -sb -e tail +0f output.pdb2gmx &
  84 endif
  85 pdb2gmx -f ${MOL}.pdb -o ${MOL}.gro -p ${MOL}.top >& ! output.pdb2gmx << KOKO
  86 1
  87 1
  88 KOKO
  89
  90 echo "pdb2gmx finished"
  91 echo -n "Press <enter>"
  92 set  ans = $<
  93
  94 ##############
  95 ### GENBOX ###
  96 ##############
  97 clear
  98 cat << _EOF_
  99 -----------------------------------------------------------------
 100 -----------------------------------------------------------------
 101 Because a simulation of a peptide in vacua is a bit unrealistic, we
 102 have to solvate our peptide in a box of water. genbox is the program
 103 we use to do this.
 104
 105 The genbox program reads the peptide structure file and an input file
 106 containing the sizes of the desired water box. The output of genbox is
 107 a gromos structure file of a peptide solvated in a box of water. The
 108 genbox program also changes the topology file ( .top extension ) to
 109 include water. First we will use the program editconf to define the
 110 right boxsize for our system.
 111
 112 -----------------------------------------------------------------
 113 -----------------------------------------------------------------
 114 _EOF_
 115
 116 if ( $?DISPLAY ) then
 117         echo "The output of the genbox program should appear"
 118         echo "in a separate xterm window"
 119 endif
 120
 121 echo -n "Press <enter>"
 122 set  ans = $<
 123
 124 echo "Starting editconf and genbox..."
 125 if ( $?DISPLAY ) then
 126         xterm -title genbox -sb -e tail +0f output.genbox &
 127 endif
 128 editconf -f ${MOL}.gro -o ${MOL}.gro -d 0.5 >& ! output.genbox
 129
 130 genbox -cp ${MOL}.gro -cs -o ${MOL}_b4em.gro -p ${MOL}.top >>& ! output.genbox
 131
 132 echo "editconf and genbox finished"
 133 echo -n "Press <enter>"
 134 set  ans = $<
 135
 136 ##############
 137 ### GROMPP ###
 138 ##############
 139 clear
 140 cat << _EOF_
 141 -----------------------------------------------------------------
 142 -----------------------------------------------------------------
 143 In principle we can start a Molecular Dynamics simulation now. However
 144 it is not very wise to do so, because our system is full of close
 145 contacts. These close contacts are mainly a result of the genbox
 146 program. The added solvent might have some close contacts with the
 147 peptide resulting in very high repulsive energies. If we would start a
 148 Molecular Dynamics (MD) simulation without energy minimisation the
 149 system would not be stable because of these high energies.
 150
 151 The standard procedure to remove these close contacts is
 152 Energy Minimisation (EM). Energy minimisation slightly changes the
 153 coordinates of our system to remove high energies from our system.
 154
 155 Before we can start the Energy Minimisation we have to preprocess all
 156 the input files with the GROMACS preprocessor named grompp. grompp
 157 preprocesses the topology file (.top), the structure file (.gro) and a
 158 parameter file (.mdp) resulting in a binary topology file (.tpr
 159 extension). This binary topology file contains all information for a
 160 simulation (in this case an energy minimisation).
 161 -----------------------------------------------------------------
 162 -----------------------------------------------------------------
 163 _EOF_
 164
 165 if ( $?DISPLAY ) then
 166         echo "The output of the grompp program should appear"
 167         echo "in a separate xterm window"
 168 endif
 169
 170 echo -n "Press <enter>"
 171 set  ans = $<
 172
 173 echo generating energy minimisation parameter file...
 174 cat > em.mdp << _EOF_
 175 title               =  ${MOL}
 176 cpp                 =  /usr/bin/cpp
 177 define              =  -DFLEX_SPC
 178 constraints         =  none
 179 integrator          =  steep
 180 dt                  =  0.002    ; ps !
 181 nsteps              =  100
 182 nstlist             =  10
 183 ns_type             =  grid
 184 rlist               =  1.0
 185 rcoulomb            =  1.0
 186 rvdw                =  1.0
 187 ;
 188 ;       Energy minimizing stuff
 189 ;
 190 emtol               =  1000.0
 191 emstep              =  0.01
 192 _EOF_
 193
 194 echo "Starting grompp..."
 195 if ( $?DISPLAY ) then
 196         xterm -title grompp -sb -e tail +0f output.grompp_em &
 197 endif
 198 grompp -f em -c ${MOL}_b4em -p ${MOL} -o ${MOL}_em >& ! output.grompp_em
 199
 200 echo "grompp finished"
 201 echo -n "Press <enter>"
 202 set  ans = $<
 203
 204 ################
 205 ### MDRUN EM ###
 206 ################
 207 clear
 208 cat << _EOF_
 209 -----------------------------------------------------------------
 210 -----------------------------------------------------------------
 211 Now the binary topology file is generated, we can start the energy
 212 minimisation (EM). The program which performs the EM is called
 213 mdrun. In fact all simulations are performed by the same program:
 214 mdrun.
 215
 216 As the Energy Minimisation is running, watch the output of the
 217 program. The first number ( from left to right ) is the number of the
 218 iteration step. The second number is the step size, which gives an
 219 indication of the change in the system. The third number is the
 220 potential energy of the system. This number starts at a high value and
 221 rapidly drops down, and converges, to a large negative value.
 222 -----------------------------------------------------------------
 223 -----------------------------------------------------------------
 224 _EOF_
 225
 226 if ( $?DISPLAY ) then
 227         echo "The output of the mdrun program should appear"
 228         echo "in a separate xterm window"
 229 endif
 230
 231 echo -n "Press <enter>"
 232 set  ans = $<
 233
 234 echo "starting energy minimisation mdrun..."
 235
 236 if ( $?DISPLAY ) then
 237         xterm -title mdrun -sb -e tail +0f output.mdrun_em &
 238 endif
 239 mdrun -nice 4 -s ${MOL}_em -o ${MOL}_em -c ${MOL}_b4pr -v >& ! output.mdrun_em
 240
 241 echo "mdrun finished"
 242 echo -n "Press <enter>"
 243 set  ans = $<
 244
 245 #################
 246 ### GROMPP PR ###
 247 #################
 248 clear
 249 cat << _EOF_
 250 -----------------------------------------------------------------
 251 -----------------------------------------------------------------
 252 Once all close contacts are removed from the system, we want to do
 253 molecular dynamics of the water molecules, and keep the peptide
 254 fixed. This is called position restrained (PR) MD.
 255
 256 Position Restrained MD keeps the peptide fixed and lets all water
 257 molecules equilibrate around the peptide in order to fill holes
 258 etc. which were not filled by the genbox program.
 259
 260 We are first going to preprocess the input files to generate the
 261 binary topology. The input files are the topology file, the structure
 262 file ( output of the EM ) a paremeter file, and an index file.
 263
 264 By default our system is split into several groups. In this case we
 265 use two of those groups: Protein and SOL(vent). We use these groups to
 266 put position restraints on all the atoms of the peptide.
 267
 268 The parameter file ( .mdp extension ) contains all information about
 269 the PR-MD like: step size, number of steps, temperature, etc. This
 270 Paramter file also tells GROMACS what kind of simulation should be
 271 performed ( like EM, PR-MD and MD etc. )
 272 -----------------------------------------------------------------
 273 -----------------------------------------------------------------
 274 _EOF_
 275
 276 if ( $?DISPLAY ) then
 277         echo "The output of the grompp program should appear"
 278         echo "in a separate xterm window"
 279 endif
 280
 281 echo -n "Press <enter>"
 282 set  ans = $<
 283
 284 echo "generating parameter file..."
 285 cat > pr.mdp << _EOF_
 286 title               =  ${MOL} position restraining
 287 cpp                 =  /usr/bin/cpp
 288 define              =  -DPOSRES
 289 constraints         =  all-bonds
 290 integrator          =  md
 291 dt                  =  0.002    ; ps !
 292 nsteps              =  500      ; total 1.0 ps.
 293 nstcomm             =  1
 294 nstxout             =  10
 295 nstvout             =  1000
 296 nstfout             =  0
 297 nstlog              =  10
 298 nstenergy           =  10
 299 nstlist             =  10
 300 ns_type             =  grid
 301 rlist               =  1.0
 302 rcoulomb            =  1.0
 303 rvdw                =  1.0
 304 ; Berendsen temperature coupling is on in two groups
 305 Tcoupl              =  berendsen
 306 tau_t               =  0.1              0.1
 307 tc-grps             =  protein          sol
 308 ref_t               =  300              300
 309 ; Pressure coupling is not on
 310 Pcoupl              =  no
 311 tau_p               =  0.5
 312 compressibility     =  4.5e-5
 313 ref_p               =  1.0
 314 ; Generate velocites is on at 300 K.
 315 gen_vel             =  yes
 316 gen_temp            =  300.0
 317 gen_seed            =  173529
 318 _EOF_
 319
 320
 321 echo "Starting grompp..."
 322 if ( $?DISPLAY ) then
 323         xterm -title grompp -sb -e tail +0f output.grompp_pr &
 324 endif
 325 grompp -f pr -c ${MOL}_b4pr -r ${MOL}_b4pr -p ${MOL} -o ${MOL}_pr >& ! output.grompp_pr
 326 echo "grompp finished"
 327
 328 echo -n "Press <enter>"
 329 set  ans = $<
 330
 331 ################
 332 ### MDRUN PR ###
 333 ################
 334 clear
 335 cat << _EOF_
 336 -----------------------------------------------------------------
 337 -----------------------------------------------------------------
 338 Now we start the Position restrained Molecular Dynamics simulation. It
 339 is important to note that in this example the simulated time is too
 340 short (1 ps) to equilibrate our system completely, but that would simple take
 341 too much time. ( about one day ).
 342
 343 -----------------------------------------------------------------
 344 -----------------------------------------------------------------
 345 _EOF_
 346
 347 if ( $?DISPLAY ) then
 348         echo "Because your DISPLAY variable is set, I will pop up a"
 349         echo "window with the output of the mdrun program"
 350 endif
 351
 352 echo -n "Press <enter>"
 353 set  ans = $<
 354
 355 echo "starting mdrun..."
 356 if ( $?DISPLAY ) then
 357         xterm -title mdrun -sb -e tail +0f output.mdrun_pr &
 358 endif
 359 mdrun -nice 4 -s ${MOL}_pr -o ${MOL}_pr -c ${MOL}_b4md -v >& ! output.mdrun_pr
 360
 361 echo "mdrun finished"
 362 echo -n "Press <enter>"
 363 set  ans = $<
 364
 365 #################
 366 ### GROMPP MD ###
 367 #################
 368 clear
 369 cat << _EOF_
 370 -----------------------------------------------------------------
 371 -----------------------------------------------------------------
 372 Now our complete system is finally ready for the actual Molecular
 373 Dynamics simulation. We start again by preprocessing the input files
 374 by the grompp program to generate the binary topology file (.tpb/.tpr
 375 extension).
 376
 377 -----------------------------------------------------------------
 378 -----------------------------------------------------------------
 379 _EOF_
 380
 381 if ( $?DISPLAY ) then
 382         echo "The output of the grompp program should appear"
 383         echo "in a separate xterm window"
 384 endif
 385
 386 echo -n "Press <enter>"
 387 set  ans = $<
 388
 389 echo "generating parameter file..."
 390 cat > md.mdp << _EOF_
 391 title               =  ${MOL} MD
 392 cpp                 =  /usr/bin/cpp
 393 constraints         =  all-bonds
 394 integrator          =  md
 395 dt                  =  0.002    ; ps !
 396 nsteps              =  5000     ; total 5 ps.
 397 nstcomm             =  1
 398 nstxout             =  50
 399 nstvout             =  0
 400 nstfout             =  0
 401 nstlist             =  10
 402 ns_type             =  grid
 403 rlist               =  1.0
 404 rcoulomb            =  1.0
 405 rvdw                =  1.0
 406 ; Berendsen temperature coupling is on in two groups
 407 Tcoupl              =  berendsen
 408 tau_t               =  0.1           0.1
 409 tc-grps             =  protein       sol
 410 ref_t               =  300           300
 411 ; Pressure coupling is not on
 412 Pcoupl              =  no
 413 tau_p               =  0.5
 414 compressibility     =  4.5e-5
 415 ref_p               =  1.0
 416 ; Generate velocites is on at 300 K.
 417 gen_vel             =  yes
 418 gen_temp            =  300.0
 419 gen_seed            =  173529
 420 _EOF_
 421
 422 echo "Starting grompp..."
 423 if ( $?DISPLAY ) then
 424         xterm -title grompp -sb -e tail +0f output.grompp_md &
 425 endif
 426 grompp -f md -c ${MOL}_b4md  -p ${MOL} -o ${MOL}_md >& ! output.grompp_md
 427
 428 echo "grompp finished"
 429 echo -n "Press <enter>"
 430 set  ans = $<
 431
 432 ################
 433 ### MDRUN MD ###
 434 ################
 435 clear
 436 cat << _EOF_
 437 -----------------------------------------------------------------
 438 -----------------------------------------------------------------
 439 Now we can start the MD simualtion. Watch the number of steps
 440 increasing ( the total number of steps is 2500, for 5 ps ).
 441
 442 -----------------------------------------------------------------
 443 -----------------------------------------------------------------
 444 _EOF_
 445
 446 if ( $?DISPLAY ) then
 447         echo "The output of the mdrun program should appear"
 448         echo "in a separate xterm window"
 449 endif
 450
 451 echo -n "Press <enter>"
 452 set  ans = $<
 453
 454 echo "starting mdrun..."
 455 if ( $?DISPLAY ) then
 456         xterm -title mdrun -sb -e tail +0f output.mdrun_md &
 457 endif
 458 mdrun -nice 4 -s ${MOL}_md -o ${MOL}_md -c ${MOL}_after_md -v >& ! output.mdrun_md
 459
 460 echo "mdrun finished"
 461 echo -n "Press <enter>"
 462 set  ans = $<
 463
 464 ############
 465 ### NGMX ###
 466 ############
 467 clear
 468 cat << _EOF_
 469 -----------------------------------------------------------------
 470 -----------------------------------------------------------------
 471 We are finished simulating, and we are going to view the calculated
 472 trajectory. The trajectory file ( .trj extension ) contains all
 473 coordinates, velocities and forces of all the atoms in our system.
 474
 475 The next program we are going run is ngmx. ngmx is a very simple
 476 trajectory viewer.
 477
 478 Once the ngmx program has been started you need to click on a few
 479 buttons to view your trajectory.
 480
 481 1. Once the program has been started a dialog box shows up. Click on
 482 the box on the left of the word Protein. ( This means that we want to
 483 view the peptide ). Then Click on the OK Button
 484
 485 2. Now we see the edges of the box with a lines drawing of the peptide
 486 we just simulated.
 487
 488 3. Select Animation in the Display menu. If you did this correctly. A
 489 dialog box at the bottom of the screen appears. This dialog box is
 490 used to move through your trajectory.
 491
 492 4. Click on the FastForward button (two triangles pointing to the
 493 right) and watch the peptide moving.
 494 -----------------------------------------------------------------
 495 -----------------------------------------------------------------
 496 _EOF_
 497
 498 if ( $?DISPLAY ) then
 499         echo Starting Trajectory viewer...
 500         ngmx -f ${MOL}_md -s ${MOL}_md  &
 501 endif
 502 #last line
 503
 504
 505
 506
 507
 508
 509
 510
 511
 512
 513
 514
 515
 516
 517
 518
 519
 520
 521
 522