#!/bin/csh

# this is the demo molecule
setenv MOL cpeptide

####################
### INTRODUCTION ###
####################
clear 
cat << _EOF_
-----------------------------------------------------------------
-----------------------------------------------------------------
Welcome to the GROMACS demo.

This is a script that takes about 10 min to run.

In this demo we will demonstrate how to simulate Molecular
Dynamics (MD) using the GROMACS software package.

The demo will perform a complete molecular dynamics (MD) simulation
of a small peptide in water. The only input file we need to do this
is a pdb file of a small peptide.

If you have any problems or remarks with respect to this demonstration,
please mail to: gromacs@gromacs.org , or check the resources on
our website http://www.gromacs.org .
-----------------------------------------------------------------
-----------------------------------------------------------------
_EOF_
echo -n "Press <enter>"
set  ans = $<

#########################
### CHECK ENVIRONMENT ###
#########################
clear 
cat << _EOF_
-----------------------------------------------------------------
-----------------------------------------------------------------
Before we you can actually start the GROMACS demo, the programs
must be present in your PATH. This might already be the case if
they are linked to /usr/local/bin. If not, follow the instructions
in the getting started section. If GROMACS is not installed 
properly on your computer, contact your system manager.
-----------------------------------------------------------------
-----------------------------------------------------------------
_EOF_
echo -n "Press <enter>"
set  ans = $<


###############
### PDB2GMX ###
###############
clear
cat << _EOF_
-----------------------------------------------------------------
-----------------------------------------------------------------
Before we can start any simulation we need a molecular toplogy
file. This topology file ( .top extension ) is generated by the
program pdb2gmx. The only input file of the pdb2gmx program is the pdb
file of our peptide ( .pdb extension ). 

Because most pdb files do not contain all hydrogen atoms, the pdb2gmx
program will also add them to our peptide. The output file which
contains the structure of the peptide when hydrogen atoms are added is a
gromos structure file ( .gro extension )  

-----------------------------------------------------------------
-----------------------------------------------------------------
_EOF_

if ( $?DISPLAY ) then
	echo "You seem to have the DISPLAY variable is set, so we will"
        echo "pop up a window with the output of the pdb2gmx program"  
endif
echo -n "Press <enter>"
set  ans = $<


echo "Starting pdb2gmx"
if ( $?DISPLAY ) then 
	xterm -title pdb2gmx -sb -e tail +0f output.pdb2gmx &
endif
echo 0 | pdb2gmx -f ${MOL}.pdb -o ${MOL}.gro -p ${MOL}.top >& ! output.pdb2gmx
echo "pdb2gmx finished"
echo -n "Press <enter>"
set  ans = $<

##############
### GENBOX ###
##############
clear
cat << _EOF_
-----------------------------------------------------------------
-----------------------------------------------------------------
Because a simulation of a peptide in vacua is a bit unrealistic, we
have to solvate our peptide in a box of water. genbox is the program
we use to do this.

The genbox program reads the peptide structure file and an input file
containing the sizes of the desired water box. The output of genbox is
a gromos structure file of a peptide solvated in a box of water. The
genbox program also changes the topology file ( .top extension ) to
include water. First we will use the program editconf to define the
right boxsize for our system.

-----------------------------------------------------------------
-----------------------------------------------------------------
_EOF_

if ( $?DISPLAY ) then
	echo "The output of the genbox program should appear"
        echo "in a separate xterm window"  
endif

echo -n "Press <enter>"
set  ans = $<

echo "Starting editconf and genbox..."
if ( $?DISPLAY ) then 
	xterm -title genbox -sb -e tail +0f output.genbox &
endif
editconf -f ${MOL}.gro -o ${MOL}.gro -d 0.5 >& ! output.genbox 

genbox -cp ${MOL}.gro -cs -o ${MOL}_b4em.gro -p ${MOL}.top >>& ! output.genbox 

echo "editconf and genbox finished"
echo -n "Press <enter>"
set  ans = $<

##############
### GROMPP ###
##############
clear
cat << _EOF_
-----------------------------------------------------------------
-----------------------------------------------------------------
In principle we can start a Molecular Dynamics simulation now. However
it is not very wise to do so, because our system is full of close
contacts. These close contacts are mainly a result of the genbox
program. The added solvent might have some close contacts with the
peptide resulting in very high repulsive energies. If we would start a
Molecular Dynamics (MD) simulation without energy minimisation the
system would not be stable because of these high energies.

The standard procedure to remove these close contacts is
Energy Minimisation (EM). Energy minimisation slightly changes the
coordinates of our system to remove high energies from our system.  

Before we can start the Energy Minimisation we have to preprocess all
the input files with the GROMACS preprocessor named grompp. grompp
preprocesses the topology file (.top), the structure file (.gro) and a
parameter file (.mdp) resulting in a binary topology file (.tpr
extension). This binary topology file contains all information for a 
simulation (in this case an energy minimisation).
-----------------------------------------------------------------
-----------------------------------------------------------------
_EOF_

if ( $?DISPLAY ) then
	echo "The output of the grompp program should appear"  
        echo "in a separate xterm window"
endif

echo -n "Press <enter>"
set  ans = $<

echo generating energy minimisation parameter file...
cat > em.mdp << _EOF_
title               =  ${MOL}
cpp                 =  /lib/cpp
define              =  -DFLEX_SPC
constraints         =  none
integrator          =  steep
dt                  =  0.002    ; ps !
nsteps              =  100
nstlist             =  10
ns_type             =  grid
rlist               =  1.0
rcoulomb            =  1.0
rvdw                =  1.0
;
;       Energy minimizing stuff
;
emtol               =  1000.0
emstep              =  0.01
_EOF_

echo "Starting grompp..."
if ( $?DISPLAY ) then 
	xterm -title grompp -sb -e tail +0f output.grompp_em &
endif
grompp -f em -c ${MOL}_b4em -p ${MOL} -o ${MOL}_em >& ! output.grompp_em

echo "grompp finished"
echo -n "Press <enter>"
set  ans = $<

################
### MDRUN EM ###
################
clear
cat << _EOF_
-----------------------------------------------------------------
-----------------------------------------------------------------
Now the binary topology file is generated, we can start the energy
minimisation (EM). The program which performs the EM is called
mdrun. In fact all simulations are performed by the same program:
mdrun. 

As the Energy Minimisation is running, watch the output of the
program. The first number ( from left to right ) is the number of the
iteration step. The second number is the step size, which gives an
indication of the change in the system. The third number is the
potential energy of the system. This number starts at a high value and
rapidly drops down, and converges, to a large negative value. 
-----------------------------------------------------------------
-----------------------------------------------------------------
_EOF_

if ( $?DISPLAY ) then
	echo "The output of the mdrun program should appear"
        echo "in a separate xterm window"
endif

echo -n "Press <enter>"
set  ans = $<

echo "starting energy minimisation mdrun..."

if ( $?DISPLAY ) then 
	xterm -title mdrun -sb -e tail +0f output.mdrun_em &
endif
mdrun -nice 4 -s ${MOL}_em -o ${MOL}_em -c ${MOL}_b4pr -v >& ! output.mdrun_em 

echo "mdrun finished"
echo -n "Press <enter>"
set  ans = $<

#################
### GROMPP PR ###
#################
clear
cat << _EOF_
-----------------------------------------------------------------
-----------------------------------------------------------------
Once all close contacts are removed from the system, we want to do
molecular dynamics of the water molecules, and keep the peptide
fixed. This is called position restrained (PR) MD.

Position Restrained MD keeps the peptide fixed and lets all water
molecules equilibrate around the peptide in order to fill holes
etc. which were not filled by the genbox program.

We are first going to preprocess the input files to generate the
binary topology. The input files are the topology file, the structure
file ( output of the EM ) a paremeter file, and an index file.

By default our system is split into several groups. In this case we
use two of those groups: Protein and SOL(vent). We use these groups to
put position restraints on all the atoms of the peptide.

The parameter file ( .mdp extension ) contains all information about
the PR-MD like: step size, number of steps, temperature, etc. This
Paramter file also tells GROMACS what kind of simulation should be
performed ( like EM, PR-MD and MD etc. )
-----------------------------------------------------------------
-----------------------------------------------------------------
_EOF_

if ( $?DISPLAY ) then
	echo "The output of the grompp program should appear"  
        echo "in a separate xterm window"
endif

echo -n "Press <enter>"
set  ans = $<

echo "generating parameter file..."
cat > pr.mdp << _EOF_
title               =  ${MOL} position restraining
warnings            =  10
cpp                 =  /lib/cpp
define              =  -DPOSRES
constraints         =  all-bonds
integrator          =  md
dt                  =  0.002	; ps !
nsteps              =  500	; total 1.0 ps.
nstcomm             =  1
nstxout             =  10
nstvout             =  1000
nstfout             =  0
nstlog              =  10
nstenergy           =  10
nstlist             =  10
ns_type             =  grid
rlist               =  1.0
rcoulomb            =  1.0
rvdw                =  1.0
; Berendsen temperature coupling is on in two groups
Tcoupl              =  berendsen
tau_t               =  0.1      	0.1
tc-grps		    =  protein  	sol
ref_t               =  300      	300
; Pressure coupling is not on
Pcoupl              =  no
tau_p               =  0.5
compressibility     =  4.5e-5
ref_p               =  1.0
; Generate velocites is on at 300 K.
gen_vel             =  yes
gen_temp            =  300.0
gen_seed            =  173529
_EOF_


echo "Starting grompp..."
if ( $?DISPLAY ) then 
	xterm -title grompp -sb -e tail +0f output.grompp_pr &
endif
grompp -f pr -c ${MOL}_b4pr -r ${MOL}_b4pr -p ${MOL} -o ${MOL}_pr >& ! output.grompp_pr
echo "grompp finished"

echo -n "Press <enter>"
set  ans = $<

################
### MDRUN PR ###
################
clear
cat << _EOF_
-----------------------------------------------------------------
-----------------------------------------------------------------
Now we start the Position restrained Molecular Dynamics simulation. It
is important to note that in this example the simulated time is too
short (1 ps) to equilibrate our system completely, but that would simple take
too much time. ( about one day ). 

-----------------------------------------------------------------
-----------------------------------------------------------------
_EOF_

if ( $?DISPLAY ) then
	echo "Because your DISPLAY variable is set, I will pop up a" 
	echo "window with the output of the mdrun program"  
endif

echo -n "Press <enter>"
set  ans = $<

echo "starting mdrun..."
if ( $?DISPLAY ) then 
	xterm -title mdrun -sb -e tail +0f output.mdrun_pr &
endif
mdrun -nice 4 -s ${MOL}_pr -o ${MOL}_pr -c ${MOL}_b4md -v >& ! output.mdrun_pr

echo "mdrun finished"
echo -n "Press <enter>"
set  ans = $<

#################
### GROMPP MD ###
#################
clear
cat << _EOF_
-----------------------------------------------------------------
-----------------------------------------------------------------
Now our complete system is finally ready for the actual Molecular
Dynamics simulation. We start again by preprocessing the input files
by the grompp program to generate the binary topology file (.tpb/.tpr
extension).

-----------------------------------------------------------------
-----------------------------------------------------------------
_EOF_

if ( $?DISPLAY ) then
	echo "The output of the grompp program should appear"  
        echo "in a separate xterm window"
endif

echo -n "Press <enter>"
set  ans = $<

echo "generating parameter file..."
cat > md.mdp << _EOF_
title               =  ${MOL} MD
cpp                 =  /lib/cpp
constraints         =  all-bonds
integrator          =  md
dt                  =  0.002	; ps !
nsteps              =  5000	; total 5 ps.
nstcomm             =  1
nstxout             =  50
nstvout             =  0
nstfout             =  0
nstlist             =  10
ns_type             =  grid
rlist               =  1.0
rcoulomb            =  1.0
rvdw                =  1.0
; Berendsen temperature coupling is on in two groups
Tcoupl              =  berendsen
tau_t               =  0.1	     0.1
tc-grps		    =  protein	     sol
ref_t               =  300	     300
; Pressure coupling is not on
Pcoupl              =  no
tau_p               =  0.5
compressibility     =  4.5e-5
ref_p               =  1.0
; Generate velocites is on at 300 K.
gen_vel             =  yes
gen_temp            =  300.0
gen_seed            =  173529
_EOF_

echo "Starting grompp..."
if ( $?DISPLAY ) then 
	xterm -title grompp -sb -e tail +0f output.grompp_md &
endif
grompp -f md -c ${MOL}_b4md  -p ${MOL} -o ${MOL}_md >& ! output.grompp_md

echo "grompp finished"
echo -n "Press <enter>"
set  ans = $<

################
### MDRUN MD ###
################
clear
cat << _EOF_
-----------------------------------------------------------------
-----------------------------------------------------------------
Now we can start the MD simualtion. Watch the number of steps
increasing ( the total number of steps is 2500, for 5 ps ).

-----------------------------------------------------------------
-----------------------------------------------------------------
_EOF_

if ( $?DISPLAY ) then
	echo "The output of the mdrun program should appear"  
        echo "in a separate xterm window"
endif

echo -n "Press <enter>"
set  ans = $<

echo "starting mdrun..."
if ( $?DISPLAY ) then 
	xterm -title mdrun -sb -e tail +0f output.mdrun_md &
endif
mdrun -nice 4 -s ${MOL}_md -o ${MOL}_md -c ${MOL}_after_md -v >& ! output.mdrun_md

echo "mdrun finished"
echo -n "Press <enter>"
set  ans = $<

############
### NGMX ###
############
clear
cat << _EOF_
-----------------------------------------------------------------
-----------------------------------------------------------------
We are finished simulating, and we are going to view the calculated
trajectory. The trajectory file ( .trj extension ) contains all
coordinates, velocities and forces of all the atoms in our system. 

The next program we are going run is ngmx. ngmx is a very simple
trajectory viewer. 

Once the ngmx program has been started you need to click on a few
buttons to view your trajectory.

1. Once the program has been started a dialog box shows up. Click on
the box on the left of the word Protein. ( This means that we want to
view the peptide ). Then Click on the OK Button

2. Now we see the edges of the box with a lines drawing of the peptide
we just simulated. 

3. Select Animation in the Display menu. If you did this correctly. A
dialog box at the bottom of the screen appears. This dialog box is
used to move through your trajectory. 

4. Click on the FastForward button (two triangles pointing to the
right) and watch the peptide moving.
-----------------------------------------------------------------
-----------------------------------------------------------------
_EOF_

if ( $?DISPLAY ) then
	echo Starting Trajectory viewer...
	ngmx -f ${MOL}_md -s ${MOL}_md  &
endif
#last line 




















