#!/bin/bash
#sgi AZ
#/usr/freeware/bin/bash 
################################
# xleap.sh    alias: leapmin
# bash script for running tleap with antechamber, parmchk, running sander
# Stefan Henrich
version="EML/AZ version 13.09.2005"
# new: with reading arguments
# 13.04.2005: cyx-distance instead of cyxbond
# 07.07.2005: settings for EML and AZ
# 11.07.2005: minimization of docking solutions (-gold); default: -noedit
# 13.07.2005: giving an extra water file: -wat
# 23.08.2005: calculating rmsd
################################################################
environment="combine2go_linux"	#environment settings
						#EML_linux
						#AZ_linux
						#combine2go_linux

if [ "$environment" = "EML_linux" ]; then
	echo EML_linux
#	source /sw/mcm/app/prepare/prepare.bash
#	prepare amber8
	sdfdir="/home/henricsn/Combine/tassfun/data4combine/"
    prgdir="/home/henricsn/Combine/tassfun/programs/"  #directory for programs
    paradir="/home/henricsn/Combine/Parameter/" #directory for parameter
    leaprc="leaprc.ff03"  #in directory: /sw/mcm/app/amber/amber8_linux/dat/leap/cmd/
elif [ "$environment" = "AZ_linux" ]; then
	echo AZ_linux
	COMBINEHOME="/home/tassfun/data/combine2go/"
	sdfdir=$COMBINEHOME"sdf/"
	prgdir=$COMBINEHOME"programs/linux/"  #directory for programs
	paradir=$COMBINEHOME"parameter/" #directory for parameter
	leaprc=/home/kemifr/AMBER8/dat/leap/cmd/leaprc.ff02
	topologyh="/usr/software/whatif/whatif99/dbdata/TOPOLOGY.H"
elif [ "$environment" = "combine2go_linux" ]; then
    echo combine2go_linux
	COMBINEHOME="/home/henricsn/combine2go/"
	sdfdir=$COMBINEHOME"sdf/"
    prgdir=$COMBINEHOME"programs/linux/"  #directory for programs
    paradir=$COMBINEHOME"parameter/" #directory for parameter
    leaprc="leaprc.ff03"  #in directory: /sw/mcm/app/amber/amber8_linux/dat/leap/cmd/
	topologyh="/sw/mcm/app/whatif-embl/dbdata/TOPOLOGY.H"
else
	echo give environment settings
fi
################################################################

        pdb="1o5a"          #pdb code 4 letters
        alt=".192B.nc1_v8ff03_restrain_eml"         #additional name for alternative conformations
        alt=""              #no alternative conformations
 ligandname="A04"           #name of ligand in pdbfile 3 letters
     charge="+1"            #total charge of the ligand
 sdf=$sdfdir"A04.sdf"  #sdf file (.sdf) in current directory for setting up bond parameters
  hydroprot="p"             #hydrogens of protein and water are given by:
   hydrolig="p"             #hydrogens of ligand are given by:
                            # a = protonate (amber) (only for protein!)
                            # b = babel
                            # i = insightII
                            # p = pdb-file
                            # w = whatif (only for protein!)
         pH="7.4"           # pH-value for adding hydrogens in insightII 
                            # (just working with hydrolig="i")
##########################
  directory="./"                              #"./"$pdb"/"
#############################################################
# used programs:
# renumb
# cyx-distance
# babel
# insightII
# pdbconv-insight
# /sw/mcm/app/whatif-embl/DO_WHATIF.EML
# pdbconv-whatif2amber
# pdbconv-whatif2amber-new
# antechamber
# parmchk
# tleap
# sander
# ambpdb
# nedit
#############################################################

##################################################################
##################################################################
# DO NOT CHANGE ANYTHING BELOW THIS LINE
# UNLESS YOU KNOW WHAT YOU ARE DOING
##################################################################
##################################################################

#############################################################
# reading arguments from command line
until [ -z "$1" ]; do
	case $1 in
		--help|-h )
			echo bash script for running tleap with antechamber, parmchk, running sander
			echo options:
			echo "-h       help"
			echo "-vers    version ($version)"
			echo "-d       (output) directory ($directory)"
			echo "-pdb     pdb code (base name of output files)"
			echo "-alt     extension of alternative conformation (suffix output	filename)"
			echo "-sdf     sdf file to determ correct bond types in antechamber"
			echo "-n       ligandname (3 characters)"
			echo "-nc      netcharge of ligand ($charge)"
			echo "-prot    name of protein pdb file ($directory\$pdb\$alt\".modif.pdb)\""
			echo "-lig     name of ligand pdb file ($directory\$pdb\".ligand.pdb)\""
			echo "-m       minimize template file (protmin.AZ.template)"
			echo "-min     read external minimize script (protmin.model.in). CAUTION: residue numbers"
			echo "-w       WHATIF for making hydrogens of protein (default: PROTON modul)"
			echo "-whbonds WHATIF for making hydrogens of protein (HBONDS modul) (not working)"
			echo "-i       InsightII for making hydrogens of ligand"
			echo "-notemp  write no temporary files"
			echo "-v       verbose (write all temporary files)"
			echo "         (default: write TMP-directory with important TMP files)"
			echo "-noac    do not run antechamber"
			echo "-nopc    do not run parmchk"
			echo "-noleap  do not run tleap"
			echo "-nos     do not run sander"
			echo "-bres    ambpdb converts residue names of min.pdb in Brookhaven-RESidue-names"
			echo "-noedit  do not start nedit with sander.out (default)"
			echo "-edit    start nedit with sander.out"
			echo "-gold    read single sdf file of GOLD instead of pdb file"
			echo "-wat     extra water pdb-file"
			exit 0
		;;
		--version|-vers )
			echo $version
			exit 0
		;;
		--environment|-e )
			# set to 1 for later testing
			shift
			echo "--environment option doesn't work"
			exit 0
		;;
		--directory|-d )
			shift
			directory=$1
			echo directory $directory
		;;
		--pdb|-pdb )
			shift
			pdb=$1
			echo $pdb
		;;
		--altconf|-alt )
			shift
			alt="."$1
			echo $alt
		;;
		--sdf|-sdf )
			shift
			sdf=$1
			echo $sdf
		;;
		--ligandname|-n )
			shift
			ligandname=$1
			echo $ligandname
		;;
		--netcharge|-nc )
			shift
			charge=$1
			echo $charge
		;;
		--protein|-prot )
			shift
			protein=$1
		;;
		--ligand|-lig )
			shift
			ligand2=$1
		;;
		--minintemplate|-m )
			shift
			minin=$1
		;;
		--minimize|-min )
			shift
			minimize=$1
			externalminimize="yes"
		;;
		--whatif|-w )
			hydroprot="w"
		;;
		--whatifhbonds|-whbonds )
			hydroprot="w"
			hbonds="hbonds"
		;;
		--insight|-i )
			hydrolig="i"
		;;
		--verboseno|-notemp )
			verbose="no"
		;;
		--verbosemax|-v )
			verbose="max"
		;;
		--antechamber|-noac )
			antechamber="no"
		;;
		--parmchk|-nopc )
			parmchk="no"
		;;
		--tleap|-noleap )
			tleap="no"
		;;
		--sander|-nos )
			sander="no"
		;;
		--bres|-bres )
			bres="yes"
		;;
		--noedit|-noedit )
			edit="no"
		;;
		--edit|-edit )
			edit="yes"
		;;
		--gold|-gold )
			gold="yes"
			#ligand2=$sdf
		;;
		--water|-wat )
			shift
			water="yes"
			waterfile=$1
		;;
		-* )
			echo "Unrecognized option: $1"
			exit 1
		;;
		* ) 
			echo no command is given
		;;
	esac 

	shift

	if [ "$#" = "0" ]; then
		break  
	fi
done 
if [ "$environment" = "EML_linux" ]; then
	source /sw/mcm/app/prepare/prepare.bash
	prepare amber8
elif [ "$environment" = "AZ_linux" ]; then
    AMBERHOME=/home/kemifr/AMBER8//
	AMBEREXE=${AMBERHOME}/exe_fix
elif [ "$environment" = "combine2go_linux" ]; then
	source /sw/mcm/app/prepare/prepare.bash
	prepare amber8
else
	echo give environment settings for amber
fi

echo start
# end of reading arguments from command line

#############################################################
##### Input files ##########################################
#protein=${protein:=$directory$pdb".192B.deprH57.modif.pdb"}		
#protein=${protein:=$directory$pdb".192B.modif.pdb"}		
#ligand2=${ligand2:=$directory$pdb".prot.ligand.pdb"} 

														#default (if variables not specified take these default values)
protein=${protein:=$directory$pdb$alt".modif.pdb"}		#default, protein with alternative conformations .modif.pdb
ligand2=${ligand2:=$directory$pdb".ligand.pdb"}      	#default, ligand from pdbconv without sulfate .ligand.pdb
##### Sander Script Template #############################################
minin=${minin:=$paradir"protmin.AZ.template"} 			#template for sander input
#minin=${minin:=$paradir"protmin.restrain.template"} 	#default, template for sander input

##### Output files #########################################
   protein2=$directory$pdb$alt".modifnew.pdb"  #combined protein and ligand file .modifnew.pdb
   xleappdb=$directory$pdb$alt".xleap.pdb"     #pdbfile from xleap   .xleap.pdb
       topo=$directory$pdb$alt".top"           #final topology file
coordinates=$directory$pdb$alt".cor"           #final coordinates
       prep=$directory$pdb$alt".prep"          #atomic charges of ligand
#       prep=$directory$pdb".prep"          #atomic charges of ligand
      param=$directory$pdb$alt".parm"          #force field parameters of ligand
#      param=$directory$pdb".parm"          #force field parameters of ligand
	 minxyz=$directory$pdb$alt".min.xyz"
	 minpdb=$directory$pdb$alt".min.pdb"
	 minbres=$directory$pdb$alt".minbres.pdb"
##### Logfiles ##############################################
   tleaplog=$directory$pdb".tleap.log"    #log from tleap
 convertlog=$directory$pdb".log.tmp"      #log from pdbconv
     minout=$directory$pdb$alt".min.out"  #minimization output
#############################################################
##### Sander Script #############################################
minimize=${minimize:=$directory$pdb$alt".protmin.restrain.in"}       #sander input
externalminimize=${externalminimize:="no"}
##### Parameter #############################################
       gaff=$paradir"gaff.dat"
atomtypegff=$paradir"my_ATOMTYPE_GFF.DEF"
#    myparam=$paradir"my_parm99.dat"
addligparam=$paradir"added_ligand.parm"
 extraparam=$paradir"extra.parm"
  extraprep=$paradir"extra.prep"
   extralib=$paradir"extra.lib"
##### tleap script ##########################################
    script2=$directory$pdb$alt".xleap.script"
##### tempfiles #############################################
    verbose=${verbose:="min"} 	 #write some TMP files into $tmpdir
    tempdir="./TMP."$pdb$alt"/"  #directory for storing tempfiles
      xleap=$directory$pdb$alt".xleap.tmp.pdb"   #tempory pdbfile of first tleap run
     script=$directory$pdb$alt".xleap.tmp.script"
    cyxbond=$directory$pdb$alt".cyxbond.out"
     renumb=$directory$pdb$alt".renumb.pdb"
   protein3=$directory$pdb$alt".modifnew.tmp.pdb"
   protein4=$directory$pdb$alt".whatifout.tmp.pdb"
whatif2amber=$directory$pdb$alt".whatif2amber.pdb"
# xleap/script/cyx-distance/renumb: $directory -> $directory$pdb$alt
#############################################################
#############################################################
echo directory for programs $prgdir
echo directory for parameters $paradir
echo
if [ -e "$convertlog" ]; then
	echo moving old logfile
	echo $convertlog into $directory$pdb".modif.log"
	echo ----------------------------------
	mv -f $convertlog $directory$pdb".modif.log" # rename logfile from convert.sh
fi

echo
echo renumbering of protein residues
echo $prgdir"renumb" $protein $renumb
echo ----------------------------------
if [ "$environment" = "EML_linux" ]; then
	$prgdir"renumb" $protein $renumb
elif [ "$environment" = "AZ_linux" ]; then
	$prgdir"renumb_sgi" $protein $renumb
elif [ "$environment" = "combine2go_linux" ]; then
	$prgdir"renumb" $protein $renumb
else
	echo renumb or renumb_sgi?
fi
	renumb_old=$renumb
### renumb.C (vers. 1.2) writes file CYX-RESIDUES.TMP with residue numbers of CYX residues
### the CYX-RESIDUES.TMP will be later used (after whatif) for renaming CYS bach into CYX

### writing header into $protein2
echo
echo combine protein and ligand pdbfile
echo $renumb and $ligand2 into $protein2
echo ----------------------------------
echo REMARK Input files:  $renumb $ligand2 > $protein2
echo REMARK Output files: $protein2 >> $protein2
echo REMARK date: $(date +%d.%m.%Y) >> $protein2
echo REMARK >> $protein2
echo REMARK Protein with modified ligand.>> $protein2
echo REMARK >> $protein2

#############################################################
#############################################################
# adding hydrogens to ligand
#############################################################

######## babel ########
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# Fehler: Unklar, wie HETATM in ATOM geaendert werden kann, ohne dass Leerstellen verloren gehen.
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if [ "$hydrolig" = "b" ]; then
   prepare babel
   echo
   echo hydrogens of ligand are determined by babel \(all possible hydrogens\!\)
   echo NOT WORKING
	grep "HETATM" $ligand2 > ligand.tmp1 
	while read line  # For as many lines as the input file has...
	do
	  echo "ATOM  "${line:6}
	done < ligand.tmp1 >  ligand.tmp2
#	done < $ligand2 >  ligand.tmp2
   echo babel -ipdb ligand.tmp2 -opdb $directory$pdb".ligand.babel.pdb" -h
   babel -ipdb ligand.tmp2 -opdb $directory$pdb".ligand.babel.pdb" -h
   ligand2=$directory$pdb".ligand.babel.pdb"
   echo $ligand2
   rm -f ligand.tmp1 ligand.tmp2
fi #babel

######## insightII ########
## example for a insightII log-file:
##the log-file can be loaded by: 
##Source_File insight_addH.log
#Builder
#Get Molecule PDB User 1o3p.noh.pdb A_1O3P -Heteroatom -Keep_Alternates -Use_Segids -Keep_All_Frames -Reference_Object
#Color Molecule Atoms A_1O3P Specified By_Atom
#Label Molecule Single_property Atom On A_1O3P Hybrid
##Hydrogens A_1O3P -set_pH -Lone_Pairs Off -Assign_IUPAC_Names
#Hydrogens A_1O3P set_pH 7.4 -Lone_Pairs Off -Assign_IUPAC_Names
###Hybridization Atom Sp3 -Lone_Pairs A_1O3P:246:C15
##Hybridization Atom Sp3 -Lone_Pairs A_1O3P:246:C17
##Delete Atom A_1O3P:246:H4
#Writing molecule A_1O3P in tmp.pdb
if [ "$hydrolig" = "i" ]; then
   echo
   echo hydrogens of ligand are added in insightII
   echo pH value is $pH
   
	#writing log-file for insightII:
	echo > $directory$pdb".insight_addH.log"
#	echo Source_File insight_addH.log >> $directory$pdb".insight.log"
	echo Builder >> $directory$pdb".insight_addH.log"
#	echo Get Molecule PDB User $directory$pdb".ligand.pdb" $ligandname Heteroatom -Keep_Alternates -Use_Segids -Keep_All_Frames -Reference_Object >> $directory$pdb".insight_addH.log"
	echo Get Molecule PDB User $ligand2 $ligandname Heteroatom -Keep_Alternates -Use_Segids -Keep_All_Frames -Reference_Object >> $directory$pdb".insight_addH.log"
	echo Color Molecule Atoms $ligandname Specified By_Atom >> $directory$pdb".insight_addH.log"
	echo Label Molecule Single_property Atom On $ligandname Hybrid >> $directory$pdb".insight_addH.log"
#	echo Hydrogens $ligandname -set_pH -Lone_Pairs Off -Assign_IUPAC_Names >> $directory$pdb".insight_addH.log"
	echo Hydrogens $ligandname set_pH $pH -Lone_Pairs Off -Assign_IUPAC_Names >> $directory$pdb".insight_addH.log"
	echo Put Molecule PDB $ligandname $directory$pdb".insight_addH.pdb" -Transformed -Displayed -Insight_Style -Write_PSF >> $directory$pdb".insight_addH.log"

### Da whatif auch MOL@ files lesen kann, koennte es sinnvoll sein, den Liganden
### als MOL2 file statt als PDB file zu speichern
### lesen von MOL2 in WHATIF:
### Reading TRIPOS MOL2 files (GETML2)
### GETML2 makes WHAT IF prompt you for the name of the Tripos MOL2 coordinate file. The coordinates will be read from this file and added to the SOUP.
	
### start insightII in directory for all complex data, so you need not to start insightII for every new ligand
	echo Please start insightII in directory data\/ \(.\/\) and read the log-file with 
	echo Source_File $directory$pdb".insight_addH.log"
	
	echo Do you want to continue with modified ligand pdb file $pdb".insight_addH.pdb"? \[1\]
	echo 1. $pdb".insight_addH.pdb"
	echo 2. change input file
	echo 3. quit
	read opt

	if [ "$opt" = "3" ]; then
		echo EXIT
		exit 3
	elif [ "$opt" = "2" ]; then
		echo read new file
		read newfile
		echo $newfile
	else
		newfile=$pdb".insight_addH.pdb"
		echo $newfile
	fi
	while [ ! -e $newfile ]; do
		echo 
		echo $newfile does not exist - change file name
		echo 
		read newfile
	done

	newfile2=$pdb".insight_renumbH.pdb"
	$prgdir"pdbconv-insight" $newfile $newfile2 ### rename hydrogens of ligand
	ligand2=$newfile2
	echo
	echo ligand file: $ligand2

fi #insightII

######## hydrogens given in pdb-file ########
if [ "$hydrolig" = "p" ]; then
   echo
   echo hydrogens of ligand are given by pdb-file
fi #pdb
#############################################################


#############################################################
# antechamber
#############################################################
if [ "$antechamber" != "no" ]; then #noantechamber
	echo
	if [ "$gold" = "" ]; then #antechamber - no gold

		if [ "$sdf" = "" ]; then #antechamber - no sdf
			echo antechamber -i $ligand2 -fi pdb -o $prep -fo prepi -c bcc -rn $ligandname -nc $charge -pf y
			if [ "$environment" = "EML_linux" ]; then
				antechamber -i $ligand2 -fi pdb -o $prep -fo prepi -c bcc -rn $ligandname -nc $charge -pf y
			elif [ "$environment" = "AZ_linux" ]; then
				$AMBEREXE/antechamber -i $ligand2 -fi pdb -o $prep -fo prepi -c bcc -rn $ligandname -nc $charge -pf y
			elif [ "$environment" = "combine2go_linux" ]; then
				antechamber -i $ligand2 -fi pdb -o $prep -fo prepi -c bcc -rn $ligandname -nc $charge -pf y
			else
				echo antechamber version?
			fi #antechamber - no sdf
		#antechamber - no sdf

		else #antechamber - sdf
			if [ "$environment" = "EML_linux" ]; then
	#			antechamber -i $ligand2 -fi pdb -a $sdf -fa mdl -ao bond -j 5 -o $prep -fo prepi -c bcc -rn $ligandname -nc $charge -pf y
				#for reading new ATOMTYPE_GFF.DEF
				echo antechamber -i $ligand2 -fi pdb -a $sdf -fa mdl -ao bond -j 5 -o PREP_AC.AC -fo ac -c bcc -rn $ligandname -nc $charge -pf y
				echo atomtype -i PREP_AC.AC -o PREP_TYPE.AC -d $atomtypegff
				echo prepgen -i PREP_TYPE.AC -o $prep -f int -rn $ligandname
				antechamber -i $ligand2 -fi pdb -a $sdf -fa mdl -ao bond -j 5 -o PREP_AC.AC -fo ac -c bcc -rn $ligandname -nc $charge -pf y
				atomtype -i PREP_AC.AC -o PREP_TYPE.AC -d $atomtypegff
				prepgen -i PREP_TYPE.AC -o $prep -f int -rn $ligandname
	#			if [ "$verbose" = "no" ]; then rm -f PREP_AC.AC PREP_TYPE.AC; fi #notemp
				echo
			elif [ "$environment" = "AZ_linux" ]; then
				#for reading new ATOMTYPE_GFF.DEF
				echo antechamber -i $ligand2 -fi pdb -a $sdf -fa mdl -ao bond -j 5 -o PREP_AC.AC -fo ac -c bcc -rn $ligandname -nc $charge -pf y
				echo atomtype -i PREP_AC.AC -o PREP_TYPE.AC -d $atomtypegff
				echo prepgen -i PREP_TYPE.AC -o $prep -f int -rn $ligandname
				$AMBEREXE/antechamber -i $ligand2 -fi pdb -a $sdf -fa mdl -ao bond -j 5 -o PREP_AC.AC -fo ac -c bcc -rn $ligandname -nc $charge -pf y
				$AMBEREXE/atomtype -i PREP_AC.AC -o PREP_TYPE.AC -d $atomtypegff
				$AMBEREXE/prepgen -i PREP_TYPE.AC -o $prep -f int -rn $ligandname
				echo
			elif [ "$environment" = "combine2go_linux" ]; then
				#for reading new ATOMTYPE_GFF.DEF
				echo antechamber -i $ligand2 -fi pdb -a $sdf -fa mdl -ao bond -j 5 -o PREP_AC.AC -fo ac -c bcc -rn $ligandname -nc $charge -pf y
				echo atomtype -i PREP_AC.AC -o PREP_TYPE.AC -d $atomtypegff
				echo prepgen -i PREP_TYPE.AC -o $prep -f int -rn $ligandname
				antechamber -i $ligand2 -fi pdb -a $sdf -fa mdl -ao bond -j 5 -o PREP_AC.AC -fo ac -c bcc -rn $ligandname -nc $charge -pf y
				atomtype -i PREP_AC.AC -o PREP_TYPE.AC -d $atomtypegff
				prepgen -i PREP_TYPE.AC -o $prep -f int -rn $ligandname
				echo
			else
				echo antechamber version?
			fi
		fi #antechamber - sdf
	#antechamber - no gold

	else #antechamber - gold
   		echo SDF file of GOLD is used as input file
		if [ "$environment" = "EML_linux" ]; then
			#for reading new ATOMTYPE_GFF.DEF
			echo antechamber -i $ligand2 -fi pdb -a $sdf -fa mdl -ao bond -j 5 -o PREP_AC.AC -fo ac -c bcc -rn $ligandname -nc $charge -pf y
			echo atomtype -i PREP_AC.AC -o PREP_TYPE.AC -d $atomtypegff
			echo prepgen -i PREP_TYPE.AC -o $prep -f int -rn $ligandname
			antechamber -i $ligand2 -fi pdb -a $sdf -fa mdl -ao bond -j 5 -o PREP_AC.AC -fo ac -c bcc -rn $ligandname -nc $charge -pf y
			atomtype -i PREP_AC.AC -o PREP_TYPE.AC -d $atomtypegff
			prepgen -i PREP_TYPE.AC -o $prep -f int -rn $ligandname
			echo
		elif [ "$environment" = "AZ_linux" ]; then
			#for reading new ATOMTYPE_GFF.DEF
			echo $AMBEREXE/antechamber -i $sdf -fi mdl -o PREP_AC.AC -fo ac -c bcc -rn $ligandname -nc $charge -pf y
			echo atomtype -i PREP_AC.AC -o PREP_TYPE.AC -d $atomtypegff
			echo prepgen -i PREP_TYPE.AC -o $prep -f int -rn $ligandname
			$AMBEREXE/antechamber -i $sdf -fi mdl -o PREP_AC.AC -fo ac -c bcc -rn $ligandname -nc $charge -pf y
			$AMBEREXE/atomtype -i PREP_AC.AC -o PREP_TYPE.AC -d $atomtypegff
			$AMBEREXE/prepgen -i PREP_TYPE.AC -o $prep -f int -rn $ligandname
			echo $AMBEREXE/antechamber -i  PREP_TYPE.AC -fi ac -o ligand.pdb -fo pdb -rn $ligandname -nc $charge -pf y
			$AMBEREXE/antechamber -i  PREP_TYPE.AC -fi ac -o ligand.pdb -fo pdb -rn $ligandname -nc $charge -pf y
			echo
		elif [ "$environment" = "combine2go_linux" ]; then
			#for reading new ATOMTYPE_GFF.DEF
			echo antechamber -i $sdf -fi mdl -o PREP_AC.AC -fo ac -c bcc -rn $ligandname -nc $charge -pf y
			echo atomtype -i PREP_AC.AC -o PREP_TYPE.AC -d $atomtypegff
			echo prepgen -i PREP_TYPE.AC -o $prep -f int -rn $ligandname
			echo antechamber -i  PREP_TYPE.AC -fi ac -o ligand.pdb -fo pdb -rn $ligandname -nc $charge -pf y
			antechamber -i $sdf -fi mdl -o PREP_AC.AC -fo ac -c bcc -rn $ligandname -nc $charge -pf y
			atomtype -i PREP_AC.AC -o PREP_TYPE.AC -d $atomtypegff
			prepgen -i PREP_TYPE.AC -o $prep -f int -rn $ligandname

			#converting PREP_TYPE.AC into pdb file for reading into tleap
			antechamber -i  PREP_TYPE.AC -fi ac -o ligand.pdb -fo pdb -rn $ligandname -nc $charge -pf y
			echo
		else
			echo antechamber version?
		fi
		ligand2=ligand.pdb
		echo new ligand file: $ligand2
   fi #antechamber - gold

fi #noantechamber
#############################################################


#############################################################
# parmchk
#############################################################
if [ "$parmchk" != "no" ]; then #noparmchk
	echo ----------------------------------
	echo
	echo parmchk -i $prep -f prepi -o $param -p $gaff
	echo ----------------------------------
	if [ "$environment" = "EML_linux" ]; then
		parmchk -i $prep -f prepi -o $param -p $gaff
		echo
	elif [ "$environment" = "AZ_linux" ]; then
		$AMBEREXE/parmchk -i $prep -f prepi -o $param -p $gaff
		echo
	elif [ "$environment" = "combine2go_linux" ]; then
		parmchk -i $prep -f prepi -o $param -p $gaff
		echo
	else
		echo parmchk version?
	fi
fi #noparmchk
#############################################################


#############################################################
# adding hydrogens to protein and water
#############################################################
if [ "$hydroprot" = "a" ]; then
   echo
   echo hydrogens of protein are determined by protonate \(amber\)
   echo NOT WORKING
fi
if [ "$hydroprot" = "b" ]; then
   echo
   echo hydrogens of protein are determined by babel \(all possible hydrogens\!\)
   prepare babel
   echo NOT WORKING
fi
if [ "$hydroprot" = "i" ]; then
   echo
   echo hydrogens of protein are given by pdb-file and were determined by insightII
   echo NOT WORKING
fi
if [ "$hydroprot" = "p" ]; then
   echo
   echo hydrogens of protein are given by pdb-file
fi

################################################################
######## whatif ########
if [ "$hydroprot" = "w" ]; then
	echo
	echo hydrogens of protein are determined by whatif
   
	### $protein2 contains header for pdb file
	### $renumb contains renumbered protein
	### --> new protein file: tmpcat2.pdb (tempory file)
	while read line ### change CYX back to CYS
	do
	  neu="$line"
	  neu="${neu/CYX/CYS}"
	  neu="${neu/HID/HIS}"
	  neu="${neu/HIE/HIS}"
	  neu="${neu/HIP/HIS}"
	  neu="${neu/ASH/ASP}"
	  echo "$neu"
	done < $renumb >  tmpcat2.pdb
	echo "rename HIS/ASP/CYX in tmpcat.pdb -> tmpcat2.pdb"
	### --> new protein file: tmpcat2.pdb (tempory file)
	
	echo "cat tmpcat2.pdb $ligand2 > $protein3"
	cat tmpcat2.pdb $ligand2 > $protein3
	if [ "$water" = "yes" ]; then
		echo cat $waterfile >> $protein3
		cat $waterfile >> $protein3
	fi

	echo
	echo whatif input:  $protein3
	echo whatif output: $protein4
	rm -f tmpcat2.pdb
	### --> new protein file: $protein3

if [ "$hbonds" = "hbonds" ]; then
################################################################
### write script for whatif ### old script 
### (TER and chainID is missing in pdb file) ###
	echo getmol $protein3 > $pdb".whatif.script"
	echo >> $pdb".whatif.script"	
	echo HBONDS >> $pdb".whatif.script"
	echo HB2NET >> $pdb".whatif.script"
	echo HB2MAK >> $pdb".whatif.script"
	echo $protein >> $pdb".whatif.script"
	echo $protein4  >> $pdb".whatif.script"
	echo  >> $pdb".whatif.script"
	echo tot  >> $pdb".whatif.script"
	echo 0 >> $pdb".whatif.script"
	echo END >> $pdb".whatif.script"
	echo %HB2LFR >> $pdb".whatif.script"
	echo END >> $pdb".whatif.script"
	echo Y >> $pdb".whatif.script"
################################################################
################################################################
else #proton
### write script for whatif ### new script with keeping TER and chainID
	echo nmr > $pdb".whatif.script"
	echo topolo >> $pdb".whatif.script"
	echo $topologyh >> $pdb".whatif.script"
	echo quit >> $pdb".whatif.script"
	echo getmol $protein3 >> $pdb".whatif.script"
	echo >> $pdb".whatif.script"	
	echo proton >> $pdb".whatif.script"
	echo opthyd >> $pdb".whatif.script"
	echo tot >> $pdb".whatif.script"
	echo 0 >> $pdb".whatif.script"
	echo quit >> $pdb".whatif.script"
	echo soup >> $pdb".whatif.script"
	echo shocys >> $pdb".whatif.script"
	echo makmol >> $pdb".whatif.script"
	echo $protein >> $pdb".whatif.script"
	echo $protein4  >> $pdb".whatif.script"
	echo  >> $pdb".whatif.script"
	echo tot  >> $pdb".whatif.script"
	echo 0 >> $pdb".whatif.script"
	echo END >> $pdb".whatif.script"
	echo %HB2LFR >> $pdb".whatif.script"
	echo END >> $pdb".whatif.script"
	echo Y >> $pdb".whatif.script"
fi #hbonds/proton
################################################################
### start whatif
	rm -f $protein4 #otherwise whatif asks for overwriting the file
if [ "$environment" = "EML_linux" ]; then
	echo
	echo "/sw/mcm/app/whatif-embl/DO_WHATIF.EML < $pdb".whatif.script" > $pdb".whatif.log""
	echo
	prepare whatif-embl
	/sw/mcm/app/whatif-embl/DO_WHATIF.EML < $pdb".whatif.script" > $pdb".whatif.log"
elif [ "$environment" = "AZ_linux" ]; then
	echo
    echo rsh ahps "cd $PWD;whatif < $pdb'.whatif.script'  > $pdb'.whatif.log' "
	echo
    rm -f $pdb'.whatif.log'
	rsh ahps " cd $PWD ; whatif < $pdb'.whatif.script'  >! $pdb'.whatif.log' "
elif [ "$environment" = "combine2go_linux" ]; then
	echo
	echo "/sw/mcm/app/whatif-embl/DO_WHATIF.EML < $pdb".whatif.script" > $pdb".whatif.log""
	echo
	prepare whatif-embl
	/sw/mcm/app/whatif-embl/DO_WHATIF.EML < $pdb".whatif.script" > $pdb".whatif.log"
else
	echo whatif?
fi
	rm -f WHATIF.FIG
	rm -f PDBFILE
	rm -f ALTERR.LOG
	### --> new protein file: $protein4

#	echo Flipped position for residue:
	grep "Flipped position for residue" $pdb".whatif.log"
#	echo There are  2 hydrogen atoms in the following group:
	grep "Group member" $pdb".whatif.log"

# giving S-S-bonds (not working at the moment)	
#	grep "CYS" $pdb".whatif.log" | awk '{if($2=="CYS")}'END'{print }' > tmp
################################################################

	### format whatif output pdb $protein2 for using in amber8
	echo
	echo pdbconv-whatif2amber input:  $protein4
	echo pdbconv-whatif2amber output: $whatif2amber
#	echo pdbconv-whatif2amber $protein4 $whatif2amber $ligandname ### rename hydrogens and residue names for amber
	### rename hydrogens and residue names for amber
	### write only protein into $whatif2amber
	if [ "$hbonds" = "hbonds" ]; then
		pdbconv-whatif2amber $protein4 $whatif2amber $ligandname #for old whatif script
	else
		#for new whatif script
		while read line ### rename oxygens and nitrogens of backbone
		do
		  neu="$line"
		  neu="${neu/"O''"/OXT}"
		  neu="${neu/"O'"/O }"
		  neu="${neu/"HN"/H }"
		  echo "$neu"
		done < $protein4 >  $whatif2amber".tmp"

		if [ "$environment" = "EML_linux" ]; then
			echo $prgdir"pdbconv-whatif2amber-new" $whatif2amber".tmp" $whatif2amber $ligandname #for new whatif script
			$prgdir"pdbconv-whatif2amber-new" $whatif2amber".tmp" $whatif2amber $ligandname #for new whatif script
		elif [ "$environment" = "AZ_linux" ]; then
			echo $prgdir"pdbconv-whatif2amber-new" $whatif2amber".tmp" $whatif2amber $ligandname #for new whatif script
			$prgdir"pdbconv-whatif2amber-new" $whatif2amber".tmp" $whatif2amber $ligandname #for new whatif script
		elif [ "$environment" = "combine2go_linux" ]; then
			echo $prgdir"pdbconv-whatif2amber-new" $whatif2amber".tmp" $whatif2amber $ligandname #for new whatif script
			$prgdir"pdbconv-whatif2amber-new" $whatif2amber".tmp" $whatif2amber $ligandname #for new whatif script
		else
			echo pdbconv-whatif2amber-new?
		fi

		rm -f $whatif2amber".tmp"
	fi #hbonds/proton pdbconv-whatif2amber
	
	renumb=$whatif2amber 
	echo renumb $renumb
	### $renumb is pointing now to $whatif2amber (output of pdbconv-whatif2amber)

	if [ "$verbose" != "max" ]; then #verbose max
		rm -f CYX-RESIDUES.TMP
		rm -f RENUMBTMP.PDB
		rm -f $protein3
		rm -f $protein4
	fi #verbose max

fi #end of whatif

#############################################################
# end of adding hydrogens
#############################################################


#############################################################
### combine renumbered protein file and ligand file
	grep -v "^END" $renumb >> $protein2
	if [ "$gold" = "" ]; then 
		echo "cat $renumb $ligand2 >> $protein2"
		echo "combined protein and ligand file is $protein2"
		cat $ligand2 >> $protein2
	else 
		echo "cat $renumb ligand.pdb \>\> $protein2"
		echo "combined protein and ligand file is $protein2"
   		cat ligand.pdb >> $protein2
	fi
#if [ "$gold" = "yes" ]; then cat ligand.pdb >> $protein2;fi
#############################################################

#exit 1

#############################################################
# tleap
#############################################################
if [ "$tleap" != "no" ]; then #notleap

	echo
	echo writing header into $xleappdb
	echo ----------------------------------
	echo REMARK Input files:  $protein $ligand2 > $xleappdb
	echo REMARK Output files: $protein2 >> $xleappdb
	echo REMARK date: $(date +%d.%m.%Y) >> $xleappdb
	echo REMARK >> $xleappdb
	echo REMARK Protein with modified ligand.>> $xleappdb
	echo REMARK >> $xleappdb

	echo
	echo writing script for tleap
	echo $script2
	#####
	echo logfile $tleaplog  > $script2
	echo logfile leap.log >> $script2
	echo source $leaprc >> $script2
	#echo loadamberparams $myparam >> $script2
	echo loadamberparams $gaff >> $script2
	echo loadamberprep $prep >> $script2
	echo loadamberparams $param >> $script2
	#echo loadamberparams $addligparam >> $script2
	#echo loadamberprep $extraprep >> $script2
	echo loadamberparams $extraparam >> $script2
	echo loadoff $extralib >> $script2
	echo $pdb=loadpdb $protein2 >> $script2
	##### cyx-distance
	echo $prgdir"cyx-distance" $renumb $pdb
		if [ "$environment" = "EML_linux" ]; then
			$prgdir"cyx-distance" $renumb $pdb >> $script2
		elif [ "$environment" = "AZ_linux" ]; then
			$prgdir"cyx-distance" $renumb $pdb >> $script2
		elif [ "$environment" = "combine2go_linux" ]; then
			$prgdir"cyx-distance" $renumb $pdb >> $script2
		else
			echo cyxbond?
		fi #cyxbond environment
	echo savepdb $pdb $xleappdb >> $script2
	echo saveamberparm $pdb $topo $coordinates  >> $script2
	echo quit >> $script2
	echo
	echo ----------------------------------
	cat $script2
	echo ----------------------------------

	##### tleap
	echo
	echo tleap -f $script2
	echo ----------------------------------
		if [ "$environment" = "EML_linux" ]; then
			tleap -f $script2
		elif [ "$environment" = "AZ_linux" ]; then
			$AMBEREXE/tleap -f $script2
		elif [ "$environment" = "combine2go_linux" ]; then
			tleap -f $script2
		else
			echo tleap version?
		fi #tleap environment
	echo "tleap output"
	echo "protein and ligand file is $xleappdb"
	echo "topology $topo"
	echo "coordinates $coordinates"
	echo ----------------------------------
fi #notleap
#############################################################

#############################################################
if [ "$externalminimize" != "yes" ]; then
	### grep last residue of protein and resi numb of ligand
	###### same like grepmin.sh
	firstresi="1"
	last=$(grep "OXT" $xleappdb | awk '{ print $5}')
	lastresi=`expr match "$last" '.*\([0-9]..[0-9]*\)'`
	#lastresi=`expr match "$last" '\([0-9]\)'`  #first OXT
	ligand=$(grep "C1  $ligandname" $xleappdb | awk '{ print $5}')
	while read line  # For as many lines as the input file has...
	do
	  neu="$line"
	  neu="${neu//\$ligand/$ligand}"
	  neu="${neu//\$start/$firstresi}"
	  neu="${neu//\$end/$lastresi}"
	  echo "$neu"
	done < $minin >  $minimize # reading template of minimize script and output into new file
fi #externalminimize
#############################################################

#############################################################
# sander
#############################################################
if [ "$sander" != "no" ]; then #nosander
	echo
	echo running sander
	echo $minimize
	echo in  $topo $coordinates
	echo out $minxyz $minout
	if [ "$environment" = "EML_linux" ]; then
		echo sander -O -i $minimize -o $minout -p $topo -c $coordinates -r $minxyz -ref $coordinates
		sander -O -i $minimize -o $minout -p $topo -c $coordinates -r $minxyz -ref $coordinates 
	elif [ "$environment" = "AZ_linux" ]; then
		echo mpirun -np 1 $AMBEREXE/sander -O -i $minimize -o $minout -p $topo -c $coordinates -r $minxyz -ref $coordinates 
		mpirun -np 1 $AMBEREXE/sander -O -i $minimize -o $minout -p $topo -c $coordinates -r $minxyz -ref $coordinates 
	elif [ "$environment" = "combine2go_linux" ]; then
		if [ "$gold" = "" ]; then #antechamber - no gold
			echo sander -O -i $minimize -o $minout -p $topo -c $coordinates -r $minxyz -ref $coordinates
	 		sander -O -i $minimize -o $minout -p $topo -c $coordinates -r $minxyz -ref $coordinates 
		else
			echo
			echo minimzation in two steps: 
			echo 1. all hydrogens \(minH.restrain.in\) -\> $directory$pdb$alt".minH.pdb"
			echo 2. all atoms \($minimize\) -\> $minpdb
			echo
			echo sander -O -i $paradir"minH.restrain.in" -o $directory$pdb$alt".minH.out" -p $topo -c $coordinates -r $directory$pdb$alt".minH.xyz" -ref $coordinates 
			sander -O -i $paradir"minH.restrain.in" -o $directory$pdb$alt".minH.out" -p $topo -c $coordinates -r $directory$pdb$alt".minH.xyz" -ref $coordinates 
			echo ambpdb -p $topo \< $directory$pdb$alt".minH.xyz" \> $directory$pdb$alt".minH.pdb"
			ambpdb -p $topo < $directory$pdb$alt".minH.xyz" > $directory$pdb$alt".minH.pdb"
			echo sander -O -i $minimize -o $minout -p $topo -c $directory$pdb$alt".minH.xyz" -r $minxyz -ref $directory$pdb$alt".minH.xyz" 
			sander -O -i $minimize -o $minout -p $topo -c $directory$pdb$alt".minH.xyz" -r $minxyz -ref $directory$pdb$alt".minH.xyz" 
			#rm -f $directory$pdb$alt".minH.xyz"
			#rm -f $directory$pdb$alt".minH.pdb"
		fi #gold
	else
		echo sander version?
	fi
	echo ----------------------------------
	echo
	echo running ambpdb
	echo in  $topo $minxyz 
	echo out $minpdb
	if [ "$environment" = "EML_linux" ]; then
		echo ambpdb -p $topo \< $minxyz \> $minpdb
		ambpdb -p $topo < $minxyz > $minpdb
	elif [ "$environment" = "AZ_linux" ]; then
		echo $AMBEREXE/ambpdb -p $topo \< $minxyz \> $minpdb
		$AMBEREXE/ambpdb -p $topo < $minxyz > $minpdb
	elif [ "$environment" = "combine2go_linux" ]; then
		echo ambpdb -p $topo \< $minxyz \> $minpdb
		ambpdb -p $topo < $minxyz > $minpdb
	else
	echo ambpdb version?
	fi
	echo ----------------------------------

	#bres (Brookhaven-residue-names)
	if [ "$bres" = "yes" ]; then
		if [ "$environment" = "EML_linux" ]; then
			echo ambpdb -bres -p $topo \< $minxyz \> $minbres
			ambpdb -bres -p $topo < $minxyz > $minbres
		elif [ "$environment" = "AZ_linux" ]; then
			echo $AMBEREXE/ambpdb -bres -p $topo \< $minxyz \> $minbres
			$AMBEREXE/ambpdb -bres -p $topo < $minxyz > $minbres
		elif [ "$environment" = "combine2go_linux" ]; then
			echo ambpdb -bres -p $topo \< $minxyz \> $minbres
			ambpdb -bres -p $topo < $minxyz > $minbres
		else
		echo ambpdb version?
		fi
	fi #bres
	
	echo loading mimimized molecule in pymol:
	if [ "$bres" = "yes" ]; then
		echo load $minbres
	else
		echo load $minpdb	
	fi

#############################################################
# calculate RMSD
#############################################################
	if [ "$gold" != "yes" ]; then
		echo $ligand2 and $minpdb
		grep $ligandname $minpdb > rmsd.tmp
		echo $prgdir"ligandrmsd" $ligand2 rmsd.tmp
		echo $pdb$alt RMSD $($prgdir"ligandrmsd" $ligand2 rmsd.tmp)
		rm -f rmsd.tmp
	fi #rmsd


fi #nosander

echo
#############################################################
# removing temp files
#############################################################
if [ "$verbose" = "no" ]; then #notemp
	echo no temp files
	echo ----------------------------------
	#antechamber
#	rm -f $topo
#	rm -f $coordinates
#	rm -f $prep
#	rm -f $param
#	if [ -e "$param" ]; then rm -f $param; fi 
#	if [ -e "$prep" ]; then rm -f $prep; fi 


	rm -f $script2 #2nd tleap script
	rm -f $protein2
	rm -f $renumb #renumbered modif.pdb
	if [ -e "$renumb_old" ]; then rm -f $renumb_old; fi
	rm -f $cyxbond
	rm -f $xleap
	rm -f $script
	#tleap
	rm -f $tleaplog
	rm -f leap.log
	rm -f $xleappdb
	rm -f divcon.in
	rm -f divcon.rst
	rm -f divcon.dmx
	rm -f divcon.out
	
	if [ -e "PREP_TYPE.AC" ]; then rm -f PREP_TYPE.AC; fi 
	if [ -e "PREP_AC.AC" ]; then rm -f PREP_AC.AC; fi 
	if [ -e "PREP.INF" ]; then rm -f PREP.INF; fi 
	if [ -e "NEWPDB.PDB" ]; then rm -f NEWPDB.PDB; fi 
	if [ -e "ATOMTYPE.INF" ]; then rm -f ATOMTYPE.INF; fi 

	
	if [ "$sander" != "no" ]; then #nosander
		#sander
		rm -f mdinfo  
		rm -f CYX-RESIDUES.TMP
		if [ "$edit" = "yes" ]; then 
			nedit $minout&
		fi #noedit
		rm -f $minimize
		rm -f $minout
	fi #nosander
	echo
	if [ "$hydroprot" = "w" ]; then
		rm -f $pdb".whatif.script"
		rm -f $pdb".whatif.log"
	fi #hydroprot

	#notemp

else #verbose
	echo
	echo moving tmp files into $tempdir
	echo ----------------------------------
	[ -d "$tempdir" ] || mkdir $tempdir
	mv -f $script2 $tempdir #2nd tleap script
	mv -f $protein2 $tempdir
	mv -f $renumb $tempdir #renumbered modif.pdb
	#tleap
	mv -f $tleaplog $tempdir
	mv -f $xleappdb $tempdir
	#mv -f divcon.in $tempdir
	#mv -f divcon.rst $tempdir
	#mv -f divcon.dmx $tempdir
	#mv -f divcon.out $tempdir

	if [ "$sander" != "no" ]; then #nosander
		#sander
		mv -f mdinfo $tempdir$pdb$alt".mdinfo"  
		mv -f $minout $tempdir
		cp -f $minimize $tempdir
		### check sander output
		nedit $tempdir$minout&
	fi #nosander
	echo
	if [ "$hydroprot" = "w" ]; then
		mv -f $pdb".whatif.script" $tempdir
		mv -f $pdb".whatif.log" $tempdir
	fi #hydroprot
	echo

	#--------------------------------------------
	if [ "$verbose" = "max" ]; then #verbose max
		echo all temp files
		if [ -e "$xleap" ]; then mv -f $xleap $tempdir; fi #1st tleap pdb file
		if [ -e "$cyxbond" ]; then mv -f $cyxbond $tempdir; fi #cyxbond output for 2nd tleap run
		if [ -e "$script" ]; then mv -f $script $tempdir; fi #1st tleap script
		if [ -e "$renumb_old" ]; then mv -f $renumb_old $tempdir; fi
		#mv -f PREP.INF $tempdir
		#mv -f NEWPDB.PDB $tempdir
		#mv -f ATOMTYPE.INF $tempdir
		#mv -f ANTECHAMBER_PREP.AC0 $tempdir
		#mv -f ANTECHAMBER_PREP.AC $tempdir
		#mv -f ANTECHAMBER_AM1BCC_PRE.AC $tempdir
		#mv -f ANTECHAMBER_AM1BCC.AC $tempdir
		#mv -f ANTECHAMBER_BOND_TYPE.AC $tempdir
		#mv -f ANTECHAMBER_AC.AC0 $tempdir
		#mv -f ANTECHAMBER_AC.AC $tempdir

		if [ "$sander" != "no" ]; then #nosander
			mv -f CYX-RESIDUES.TMP $tempdir
		fi #nosander

	#--------------------------------------------
	else #verbose min
		echo some temp files
		if [ -e "$xleap" ]; then rm -f $xleap; fi #1st tleap pdb file
		rm -f $cyxbond #cyxbond output for 2nd tleap run
		rm -f $script #1st tleap script
		if [ -e "$renumb_old" ]; then rm -f $renumb_old; fi

		if [ "$sander" != "no" ]; then #nosander
			rm -f CYX-RESIDUES.TMP
		fi #nosander

	fi #verbose
fi #verbose-notemp
echo $(date)