GROMACS run script

#!/bin/bash
#SBATCH
#SBATCH -J 2D6-2F9Q_52
#SBATCH -t 24:00:00
#SBATCH -p gpu
#SBATCH -N 12
#SBATCH -n 12
#SBATCH --mail-user=parker.dewaal09@kzoo.edu
#SBATCH --mail-type=begin  # email me when the job starts
#SBATCH --mail-type=end    # email me when the job finishes

# Script designed to test a single node with 16 CPU gromacs run
# written by Parker de Waal 5/30/2013
# Thanks to Dirk Colbry for generiously sharing his scripts
# Stampede settings
export OMP_NUM_THREADS=16
module load gromacs/4.6.1

# GROMACS settings can be adjusted here 

pdb_name=2F9Q_A_52          # Name of initial pdb file
forcefield=gromos53a6     # 53a6 for heme containing proteins
water=spc                 # simple point change is recommended for 53a6

#make working directory        
working=${WORK}/${SLURM_JOB_NAME}/${pdb_name}_output/${SLURM_JOB_ID}
mkdir -p $working
cd ${working}
#cd ${SLURM_SUBMIT_DIR}
cp $HOME/GROMACS_dev/* ${working}

#mkdir ${pdb_name}
#cp ./* ./${pdb_name}
#echo "it worked"
#cd ${pdb_name}

#convert pdb files to topology (top) and coordinate files (gro)
pdb2gmx -ff ${forcefield} -water ${water} -f ${pdb_name}.pdb -o protein_processed.gro -p protein.top -ignh
if [ ! $? == 0 ]
then
	echo "ERROR - pdb2gmx returnd with an error"
	exit 1
fi

# Ensure that the coordinate file does not include nans Thanks Dirk!
if [ ! `grep -e nan protein_processed.gro | wc -l` = 0 ]
then
	echo "ERROR - nans found in protein_processed.gro file"
	exit 2
fi

# Box definition and solvate
editconf -f protein_processed.gro -o protein_newbox.gro -d 1 -bt cubic
genbox -cp protein_newbox.gro -cs spc216.gro -o protein_solv.gro -p protein.top
if [ ! $? == 0 ]
then
	echo "ERROR - editconf or genbox failed :("
	exit 3
fi

#3 Adding Ions
grompp -f ions.mdp -c protein_solv.gro -p protein.top -o ions.tpr
echo 15 | genion -s ions.tpr -o protein_solv_ions.gro -p protein.top -pname NA -nname CL -conc 0.15 -neutral
echo "Ions added and saved?"

if [ ! $? == 0 ]
then
	echo "ERROR - ions failed"
	exit 4
fi

#minimization
grompp -f minim.mdp -c protein_solv_ions.gro -p protein.top -o em.tpr
mdrun_gpu -v -pd -deffnm em
echo 12 0 | g_energy -f em.edr -o potential.xvg
echo -e "First round Energy Minimization saved as potential.xvg"
if [ ! $? == 0 ]
then
        echo "ERROR - Energy minimization"
        exit 5
fi

#pr
#grompp -v -f pr.mdp -c em.gro -p protein.top -o PR.tpr
#mdrun -v -pd -deffnm PR
#echo 10 0 | g_energy -f PR.edr -o relaxed_potential.xvg
#if [ ! $? == 0 ]
#then
#        echo "ERROR - PR"
#        exit 6
#fi
#exit

#5. Establish EQ within the system
grompp -f nvt.mdp -c em.gro -p protein.top -o nvt.tpr
ibrun tacc_affinity mdrun_mpi_gpu -v -deffnm nvt
if [ ! $? == 0 ]
then
	echo "ERROR - nvt"
	exit 7
fi

echo 14 0 | g_energy -f nvt.edr -o temperature.xvg
echo -e "Temp Equilibrium established saved as temperature.xvg"

grompp -f npt.mdp -c nvt.gro -t nvt.cpt -p protein.top -o npt.tpr
ibrun tacc_affinity mdrun_mpi_gpu -v -deffnm npt

if [ ! $? == 0 ]
then
	echo "ERROR - npt failed"
	exit 8
fi
echo 16 22 0 | g_energy -f npt.edr -o pressure_density.xvg

# Production
grompp -f md.mdp -c npt.gro -t npt.cpt -p protein.top -o ${pdb_name}.tpr
ibrun tacc_affinity mdrun_mpi_gpu -v -deffnm ${pdb_name}
if [ ! $? == 0 ]
then
    echo "ERROR - mdrun returnd with an error"
	exit=9
fi


exit $ex
exit

ions.mdp

title        	= Ion addition	; Title of run

; The following line tell the program the standard locations where to find certain files
cpp		= /lib/cpp	; Preprocessor

; Define can be used to control processes
define          = -DFLEXIBLE

; Parameters describing what to do, when to stop and what to save
integrator	= steep		; Algorithm (steep = steepest descent minimization)
emtol		= 1.0		; Stop minimization when the maximum force < 1.0 kJ/mol
nsteps		= 50		; Maximum number of (minimization) steps to perform default 500,
nstenergy	= 1		; Write energies to disk every nstenergy steps
energygrps	= System	; Which energy group(s) to write to disk

; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
ns_type		= grid		; Method to determine neighbor list (simple, grid)
coulombtype	= PME		; Treatment of long range electrostatic interactions
rcoulomb	= 1.4		; long range electrostatic cut-off
rvdw		= 1.4		; long range Van der Waals cut-off
rlist		= 1.4		; Cut-off for making neighbor list (short range forces)
constraints	= none		; Bond types to replace by constraints
pbc		= xyz		; Periodic Boundary Conditions (yes/no)

md.mdp

; RUN CONTROL PARAMETERS
integrator         		= md
tinit 				= 0 			; Starting time
dt 					= 0.002 		; 2 femtosecond time step for integration
nsteps 				= 25000000		; Make it 50 ns
cutoff-scheme       = Verlet
; OUTPUT CONTROL OPTIONS
nstxout 			= 250000 		; Writing full precision coordinates every 0.5 ns
nstvout 			= 250000 		; Writing velocities every 0.5 ns
nstlog 				= 1000 			; Writing to the log file every ps
nstenergy 			= 1000 			; Writing out energy information every ps
nstxtcout 			= 1000 			; Writing coordinates every ps
energygrps			= Protein Non-Protein
; NEIGHBORSEARCHING PARAMETERS
nstlist 			= 5
ns-type				= Grid
pbc 				= xyz
rlist 				= 1.0
; OPTIONS FOR ELECTROSTATICS AND VDW
coulombtype 		= PME
pme_order 			= 4 			; cubic interpolation
fourierspacing 		= 0.16 			; grid spacing for FFT
rcoulomb 			= 1.0
vdw-type 			= Cut-off
rvdw 				= 1.0
; Dispersion correction
DispCorr 			= EnerPres 		; account for vdw cut-off
; Temperature coupling 
Tcoupl 				= v-rescale
tc-grps 			= Protein Non-Protein
tau_t 				= 0.1 0.1
ref_t 				= 300 300
; Pressure coupling 
Pcoupl 				= Berendsen 
Pcoupltype 			= Isotropic
tau_p 				= 2.0
compressibility 	= 4.5e-5
ref_p 				= 1.0
; GENERATE VELOCITIES FOR STARTUP RUN
gen_vel 			= no
; OPTIONS FOR BONDS 
constraints 		= all-bonds
constraint-algorithm = lincs
continuation 		= yes ; Restarting after NPT without position restraints
lincs-order			= 4
lincs-iter 			= 1
lincs-warnangle 	= 30

minim.mdp

; Define can be used to control processes
define         		= -DFLEXIBLE 		; Use flexible water model

; Parameters describing what to do, when to stop and what to save
integrator 		= steep 			; 
emtol 			= 250.0			; Stop minimization when the maximum force < 1.0 kJ/mol
nsteps 			= 1000 				; Maximum number of (minimization) steps to perform
nstenergy 		= 1 				; Write energies to disk every nstenergy steps
energygrps 		= System 			; Which energy group(s) to write to disk
; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist 		= 1 				; Frequency to update the neighbor list
ns_type 		= grid 				; Method to determine neighbor list (simple, grid)
coulombtype 	        = PME 				; Treatment of long range electrostatic interactions
rlist 			= 1.0 				; Cut-off for making neighbor list (short range forces)
rcoulomb 		= 1.0 				; long range electrostatic cut-off
rvdw 			= 1.0 				; long range Van der Waals cut-off
constraints 	        = none 				; Bond types to replace by constraints
pbc 			= xyz 				; Periodic Boundary Conditions (yes/no)
cutoff-scheme = Verlet

npt.mdp

; VARIOUS PREPROCESSING OPTIONS
define         			= -DPOSRES
; RUN CONTROL PARAMETERS
integrator 			= md
dt 					= 0.002 					; time step (in ps)
nsteps 				= 2500000 					; number of steps 
; OUTPUT CONTROL OPTIONS
nstxout 			= 1000 						; save coordinates every ps
nstvout 			= 1000						; save velocities every ps
nstenergy 			= 1000						; save energies every ps
nstlog 				= 1000						; update log file every ps
energygrps			= Protein Non-Protein
; NEIGHBORSEARCHING PARAMETERS
nstlist 			= 5
ns_type 			= grid
pbc 				= xyz
rlist 				= 1.0
; OPTIONS FOR ELECTROSTATICS AND VDW
coulombtype 		        = PME 						; Particle Mesh Ewald for long-rangeelectrostatics
pme_order 			= 4 						; cubic interpolation
fourierspacing 		        = 0.16 						; grid spacing for FFT
rcoulomb 			= 1.0
vdw-type 			= Cut-off
rvdw 				= 1.0
; Temperature coupling
tcoupl 				= v-rescale 				; Couple temperature to external heat bath according to velocity re-scale method
tc-grps 			= Protein Non-Protein 		; Use separate heat baths for Protein and Non-Protein groups
tau_t 				= 0.1 0.1 					; Coupling time constant, controlling strength of coupling
ref_t 				= 300 300 					; Temperature of heat bath
; Dispersion correction
DispCorr 			= EnerPres 					; account for vdw cut-off
; Pressure coupling is off
Pcoupl 				= Berendsen
Pcoupltype 			= Isotropic
tau_p 				= 2.0
compressibility 	= 4.5e-5
ref_p 				= 1.0
refcoord_scaling 	= com				
; GENERATE VELOCITIES FOR STARTUP RUN
gen_vel 			= no
; OPTIONS FOR BONDS 
constraints 		= all-bonds 				; All bonds will be treated as constraints (fixed length)
continuation 		= yes 						; second dynamics run
constraint_algorithm = lincs 					; holonomic constraints
lincs_iter 			= 1 						; accuracy of LINCS
lincs_order 		= 4 						; also related to accuracy
cutoff-scheme       = Verlet

nvt.mdp

; Run parameters
integrator        = md		; leap-frog integrator
nsteps		= 50000		; 2 * 50000 = 100 ps 50000
dt		= 0.002		; 2 fs
; Output control
nstxout		= 1000		; save coordinates every 2 ps
nstvout		= 1000		; save velocities every 2 ps
nstenergy	= 1000		; save energies every 2 ps
nstlog		= 1000		; update log file every 2 ps
; Bond parameters
continuation	= no		; first dynamics run
constraint_algorithm = lincs	; holonomic constraints 
constraints	= all-bonds	; all bonds (even heavy atom-H bonds) constrained
lincs_iter	= 1		; accuracy of LINCS
lincs_order	= 4		; also related to accuracy
; Neighborsearching
ns_type		= grid		; search neighboring grid cells
nstlist		= 5		; 10 fs
rlist		= 1.4		; short-range neighborlist cutoff (in nm)
rcoulomb	= 1.4		; short-range electrostatic cutoff (in nm)
rvdw		= 1.4		; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype	= PME		; Particle Mesh Ewald for long-range electrostatics
pme_order	= 4		; cubic interpolation
fourierspacing	= 0.16		; grid spacing for FFT
; Temperature coupling is on
tcoupl		= V-rescale	; modified Berendsen thermostat
tc-grps		= Protein Non-Protein	; two coupling groups - more accurate
tau_t		= 0.1	0.1	; time constant, in ps
ref_t		= 300 	300	; reference temperature, one for each group, in K
; Pressure coupling is off
pcoupl		= no 		; no pressure coupling in NVT
; Periodic boundary conditions
pbc		= xyz		; 3-D PBC
; Dispersion correction
DispCorr	= EnerPres	; account for cut-off vdW scheme
; Velocity generation
gen_vel		= yes		; assign velocities from Maxwell distribution
gen_temp	= 300		; temperature for Maxwell distribution
gen_seed	= -1		; generate a random seed
cutoff-scheme        = Verlet