#!/bin/bash
################################################################
# calc_interaction.sh
#
# bash script for calculating interaction energies with AMBER8 
# Stefan Henrich		
# 07.07.2005
# 16.09.2005 changes similar to calc_inter_model.sh
# 20.09.2005 
# EML/AZ version
################################################################
environment="combine2go_linux"	#environment settings
						#EML_linux
						#AZ_linux
						#combine2go_linux
if [ "$environment" = "EML_linux" ]; then
	echo EML_linux
	COMBINEHOME="/home/henricsn/combine2go/programs/"
	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/
	#ligdir=$COMBINEHOME"data/urokinase/pdb/ligand/"          #directory with prep/parm files of ligand
	#directory=$COMBINEHOME"data/urokinase/prepared/" 
	ligdir=$COMBINEHOME"data/urokinase/prepared/"
elif [ "$environment" = "AZ_linux" ]; then
	echo combine2go_AZ
	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"  #in directory: /sw/mcm/app/amber/amber8_linux/dat/leap/cmd/
	#ligdir=$COMBINEHOME"data/urokinase/pdb/ligand/"		  #directory with prep/parm files of ligand
	#directory=$COMBINEHOME"data/urokinase/prepared/" 
	ligdir=$COMBINEHOME"data/urokinase/prepared/"
elif [ "$environment" = "combine2go_linux" ]; then
	echo combine2go_linux
	COMBINEHOME="/home/henricsn/combine2go/"
	sdfdir=$COMBINEHOME"sdf/"
	sdf=$sdfdir"all_300805_3D.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/
	ligdir=$COMBINEHOME"data/urokinase/pdb/ligand/rename/"          #directory with pdb files of ligand
	#directory=$COMBINEHOME"data/urokinase/prepared/" 
#	ligdir="../"
	targetdir="./" 
	directory="../" #source directory with top/cor/parm/prep files
else
	echo give environment settings
fi
################################################################
# COMBINEHOME		project directory
# sdfdir			directory for main library SD file
# sdf				main library SD file in sdfdir
# prgdir			program directory
# paradir			parameter directory (minimization template, anal.in, special ligand parameter)
# leaprc			run command for tleap
# ligdir			directory with pdb file of ligand for calculating prep
# targetdir			directory for storing the output (./)
# directory			directory of prep, parm, min.xyz, cor, and top files
# getcutpdb			pdb files (eg.: *.cr.pdb) in targetdir, 
#					where the interaction energies should be calculated for
# prop_activity		flag of SD file of activity (eg. uPA)
################################################################
prop_activity=uPA ### gave flag of SD file of activity (eg. uPA)

#directory="/home/henricsn/Combine/tassfun/urokinase/pdb/" 
directory=${directory:="./"}  #source directory with top/cor files
#ligdir="/home/henricsn/Combine/tassfun/urokinase/pdb/ligand/"
#ligdir=${ligdir:="../"}  #root directory of subdirectories with prep/parm-files
targetdir=${targetdir:="./"} # output directory for aout

################################################################
##### Input files ##########################################
### read all modified min.pdb-files (eg. .cr.pdb or .cut.pdb) of targetdir 
### if no min.xyz exists or the sequence was cut
#  getminpdb=$(find $targetdir -name "*.min.pdb")
#  getcutpdb=$(find $targetdir -name "*.cut.pdb")
  getcutpdb=$(find $targetdir -maxdepth 1 -name "*.cr.pdb")
#  getcutpdb=$targetdir"1c5y1.A21A.cr.pdb"    #for reading a single file
#  getcutpdb=$targetdir"1f92.0.cr.pdb"    #for reading a single file

################################################################
##### Output files ##########################################
#     pdbout="y"			# y, write output pdb ($xleappdb, .xcut.pdb)
#activitytext="n"		# write new activity.txt (overwrites old one!!!)
#desolvation="y"			# y, calculate desolvation free energy with uhbd
   golpeout="output"
   golpedat=$targetdir$golpeout".dat"
 golpenames=$targetdir$golpeout".dat.VarNames"
golpevartyp=$targetdir$golpeout".dat.VarType"
##################################################################
##################################################################
# DO NOT CHANGE ANYTHING BELOW THIS LINE
# UNLESS YOU KNOW WHAT YOU ARE DOING
##################################################################
##################################################################

################################################################
#source /sw/mcm/app/prepare/prepare.bash
#prepare amber8
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
#############################################################
# log-file
##################################################################
echo $(date +%d.%m.%Y) > interaction.log
#echo $(zdump EST) > interaction.log
echo source directory $directory >> interaction.log
echo target directory $targetdir >> interaction.log
echo  >> interaction.log
#echo uPA without water >> interaction.log
##################################################################
### write golpe file ($golpedat)
molnumber=1
##################################################################
FILES="$getcutpdb"
echo 

for file in $FILES
do
   if [ ! -e "$file" ]
   then
       echo "$file does not exist."; echo
	   continue
   fi

################################################################
#### get pdb and suffix
#### neu von calc_inter_model.sh
	### (./A14.0.cor or ./1gj9.192A.cr.pdb)
#	path=$(dirname $file)/ #folder of file
	base=`basename $file `   #file without path (A14.0.cor or 1gj9.192A.cr.pdb)
	name=`expr "$base" : '\([0-9a-zA-Z][0-9a-zA-Z]*\)' ` #base name (A14 or 1gj9)
	pdb=$name
	suffix="${base#*.}" #suffix (0.cor or 192A.cr.pdb)
	counter=`expr "$suffix" : '\([0-9a-zA-Z][0-9a-zA-Z]*\)' ` #counter (0 or 192A)
	alt="."$counter
	echo $base $name $suffix $counter
	###grep residue name and sequence number of ligand pdb file
	line=$(grep "C1" $file|uniq -w6)
	residue=${line:17:3}
	lig_number=${line:20:6}
	echo residue $residue sequence number of ligand $lig_number
	count=$file #used for residue number of ligand

##################################################################
##################################################################
$prgdir"pdb2seq" $file > sequence.cut
$prgdir"pdb2seq" $directory$pdb$alt.min.pdb > sequence.ori
oricut=$(cmp -s sequence.ori sequence.cut) 
		#compare sequence of min.pdb and cut.pdb
		#if the sequences are not identical ($? -ne 0) and min.xyz does not exist then
		#run tleap 
#rm -f sequence.cut sequence.ori #remove sequence files
##################################################################
if [[ $? -ne 0 || ! -e "$directory$pdb$alt.min.xyz" || ! -e "$directory$pdb$alt.top" ]]; then #cor top

# Bedingung noch falsch
#if [[ $? -ne 0 || ! -e "$directory$pdb$alt.min.xyz" || ! -e "$directory$pdb$alt.top"  ||  ! -e "$targetdir$pdb$alt.cor" || ! -e "$targetdir$pdb$alt.top"]] ; then #cor top
echo creating amber-script
#exit
##################################################################
########### Input files ##########################################
##################################################################
#	prep=$targetdir$pdb".prep"          #atomic charges of ligand
#	param=$targetdir$pdb".parm"          #force field parameters of ligand
    prep=$directory$pdb".prep"          #atomic charges of ligand
    param=$directory$pdb".parm"          #force field parameters of ligand
##### Output files #########################################
	xleappdb=$targetdir$pdb$alt".xcut.pdb"     #pdbfile from xleap   .xleap.pdb
	topo=$targetdir$pdb$alt".top"           #final topology file
	coordinates=$targetdir$pdb$alt".cor"           #final coordinates
##### Parameter #############################################
#    paradir="/home/henricsn/Combine/Parameter/" #directory for parameter
#     leaprc="leaprc.ff03"  #in directory: /sw/mcm/app/amber/amber8_linux/dat/leap/cmd/
	gaff=$paradir"gaff.dat"
#    myparam=$paradir"my_parm99.dat"
	addligparam=$paradir"added_ligand.parm"
	extraparam=$paradir"extra.parm"
	extraprep=$paradir"extra.prep"
	extralib=$paradir"extra.lib"
	atomtypegff=$paradir"my_ATOMTYPE_GFF.DEF"
##### Programs #############################################
#     prgdir="/home/henricsn/Combine/tassfun/programs/"  #directory for programs
##### tleap script ##########################################
	script=$targetdir$pdb$alt".xleap.script"
##### tempfiles #############################################
	renumb=$targetdir$pdb$alt".renumb.pdb"
##################################################################
##################################################################

if [ ! -e "$prep" ] && [ ! -e "$targetdir$pdb.prep" ]; then #prep
	echo run antechamber
	echo ---------------
	if [ -e "$targetdir$pdb.prep" ]; then prep=$targetdir$pdb".prep";fi          #atomic charges of ligand

#	residue=$(grep $pdb $table | awk '{ print $2 } ')
#	if [ ! -e "$directory$residue.sdf" ]; then 
#		ID=$(grep $pdb $table | awk '{ print $1 } ')
#		#echo $ID > tmp.label
#		#echo $prgdir"select_sdf" -labelfile tmp.label -property_name ID \< $sdf \> $directory$residue.sdf
#		#$prgdir"select_sdf" -labelfile tmp.label -property_name ID < $sdf > $directory$residue.sdf
#		echo $prgdir"grepsdtag.linux" -t ID $ID \< $sdf \> $path$residue.sdf
#		$prgdir"grepsdtag.linux" -t ID $ID < $sdf > $path$residue.sdf
#                echo $directory$residue.sdf
#		echo
#		echo netcharge $netcharge
#		rm -f tmp.label
#	fi #select_sdf
#	netcharge=$(grep $pdb $table | awk '{ print $3 } ')
#	ligandname=$(grep $pdb $table | awk '{ print $2 } ')
localsdf=$targetdir$residue.sdf
echo $prgdir"grepsdtag.linux" -t residue $residue \< $sdf \>  $localsdf
$prgdir"grepsdtag.linux" -t residue $residue < $sdf >  $localsdf
#$directory$residue.sdf
#ligandname
#	echo netcharge $netcharge residue $ligandname
		### calculate netcharge
		charge_numb=0;netcharge=0;charge=0;column=0
		chargeline=$(grep "M  CHG" $localsdf)
		charged_atoms=`echo $chargeline| awk '{ print $3 } '`
		until [ "$charged_atoms" -le "$charge_numb" ]; do
			column=$( echo "5 + $charge_numb * 2" |bc -l)
			#variable column of bash will be used in awk
			charge=`echo $chargeline| awk -v column="$column" '{ total = $column } END {print $column}'`
			netcharge=$( echo "$netcharge + $charge" |bc -l)
			let charge_numb+=1
		done
echo residue $residue localsdf $localsdf netcharge $netcharge

###############
	###grep residue name and sequence number of ligand pdb file
################	

	ligdir=${ligdir:=$directory} #if no ligand dir is given than take $directory
	#alternative conformation of ligand is specified in protein name by _A or _B and in .A.prot.ligand.align.pdb or .B.prot.ligand.align.pdb
#	if [[ "$base" == ????.*_A.*.pdb ]]; then
#		ligand2=$ligdir$pdb".A.prot.ligand.align.pdb"
#	elif [[ "$base" == ????.*_B.*.pdb ]]; then
#		ligand2=$ligdir$pdb".B.prot.ligand.align.pdb"
#	else
#		ligand2=$ligdir$pdb".prot.ligand.align.pdb"
#	fi
	if [ -e "$ligdir$name.0.pdb" ]; then ligand2=$ligdir$name.0.pdb
	else ligand2=$ligdir$name.1.pdb;fi

	echo antechamber -i $ligand2 -fi pdb -a $directory$residue.sdf -fa mdl -ao bond -j 5 -o $prep -fo prepi -c bcc -rn $ligandname -nc $netcharge -pf y
	if [ "$environment" = "EML_linux" ]; then
#		antechamber -i $ligand2 -fi pdb -a $directory$residue.sdf -fa mdl -ao bond -j 5 -o $prep -fo prepi -c bcc -rn $ligandname -nc $netcharge -pf y	
		#for reading new ATOMTYPE_GFF.DEF
		echo antechamber -i $ligand2 -fi pdb -a $directory$residue.sdf -fa mdl -ao bond -j 5 -o PREP_AC.AC -fo ac -c bcc -rn $ligandname -nc $netcharge -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 $netcharge -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
		#antechamber -i $ligand2 -fi pdb -a $directory$residue.sdf -fa mdl -ao bond -j 5 -o $prep -fo prepi -c bcc -rn $ligandname -nc $netcharge -pf y	
	    #for reading new ATOMTYPE_GFF.DEF
		echo $AMBEREXE/antechamber -i $ligand2 -fi pdb -a $directory$residue.sdf -fa mdl -ao bond -j 5 -o PREP_AC.AC -fo ac -c bcc -rn $ligandname -nc $netcharge -pf y
		echo $AMBEREXE/atomtype -i PREP_AC.AC -o PREP_TYPE.AC -d $atomtypegff
		echo $AMBEREXE/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 $netcharge -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
#		antechamber -i $ligand2 -fi pdb -a $directory$residue.sdf -fa mdl -ao bond -j 5 -o $prep -fo prepi -c bcc -rn $ligandname -nc $netcharge -pf y	
		#for reading new ATOMTYPE_GFF.DEF
		echo antechamber -i $ligand2 -fi pdb -a $localsdf -fa mdl -ao bond -j 5 -o PREP_AC.AC -fo ac -c bcc -rn $residue -nc $netcharge -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 $residue
		antechamber -i $ligand2 -fi pdb -a $localsdf -fa mdl -ao bond -j 5 -o PREP_AC.AC -fo ac -c bcc -rn $residue -nc $netcharge -pf y
		atomtype -i PREP_AC.AC -o PREP_TYPE.AC -d $atomtypegff
		prepgen -i PREP_TYPE.AC -o $prep -f int -rn $residue
		echo
	else
		echo antechamber version?
	fi #antechamber
	
fi #prep
#exit 1

#if [ ! -e "$param" ]; then #param
if [ ! -e "$param" ] && [ ! -e "$targetdir$pdb.parm" ]; then #param
	if [ -e "$targetdir$pdb.parm" ]; then param=$targetdir$pdb".parm";fi          #force field parameters of ligand
	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
	elif [ "$environment" = "combine2go_linux" ]; then
		parmchk -i $prep -f prepi -o $param -p $gaff
		echo
	else
		echo parmchk version?
	fi
fi #param
	if [ -e "divcon.rst" ]; then rm -f divcon.rst; fi 
	if [ -e "divcon.dmx" ]; then rm -f divcon.dmx; fi 
	if [ -e "divcon.in" ]; then rm -f divcon.in; fi 
	if [ -e "divcon.out" ]; then rm -f divcon.out; fi 
	if [ -e "ANTECHAMBER_BOND_TYPE.AC0" ]; then rm -f ANTECHAMBER_BOND_TYPE.AC0; fi 
	if [ -e "ANTECHAMBER_BOND_TYPE.AC" ]; then rm -f ANTECHAMBER_BOND_TYPE.AC; fi 
	if [ -e "ANTECHAMBER_AC.AC0" ]; then rm -f ANTECHAMBER_AC.AC0; fi 
	if [ -e "ANTECHAMBER_AC.AC" ]; then rm -f ANTECHAMBER_AC.AC; fi 
	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 
#exit 2
##################################################################
##################################################################

echo
echo $prgdir"renumb" $file $renumb
$prgdir"renumb" $file $renumb

##################################################################
### write tleap script
##################################################################
echo writing script for tleap
echo $script
echo -------------------------------------------------------------
#####
echo source $leaprc > $script
#echo loadamberparams $myparam >> $script
echo loadamberparams $gaff >> $script
echo loadamberprep $extraprep >> $script
echo loadamberparams $extraparam >> $script
echo loadamberprep $prep >> $script
echo loadamberparams $param >> $script
echo loadamberparams $addligparam >> $script
echo loadoff $extralib >> $script
echo $pdb=loadpdb $renumb >> $script

echo
echo $prgdir"cyx-distance" $renumb $pdb
		if [ "$environment" = "EML_linux" ]; then
			$prgdir"cyx-distance" $renumb $pdb >> $script
		elif [ "$environment" = "AZ_linux" ]; then
			$prgdir"cyx-distance" $renumb $pdb >> $script
		elif [ "$environment" = "combine2go_linux" ]; then
			$prgdir"cyx-distance" $renumb $pdb >> $script
		else
			echo cyx-distance?
		fi #cyxbond environment

if [ "$pdbout" = "y" ]; then
echo savepdb $pdb $xleappdb >> $script
fi
echo saveamberparm $pdb $topo $coordinates  >> $script
echo quit >> $script
##################################################################
### end of tleap script
##################################################################

echo
echo running tleap with $script
echo -------------------------------------------------------------
#tleap -f $script
	if [ "$environment" = "EML_linux" ]; then
		tleap -f $script
	elif [ "$environment" = "AZ_linux" ]; then
		$AMBEREXE/tleap -f $script
	elif [ "$environment" = "combine2go_linux" ]; then
		tleap -f $script
	else
		echo tleap version?
	fi #tleap environment

#rm -f $renumb
rm -f $script
rm -f leap.log
##################################################################
### end of tleap
##################################################################
count=$renumb #used for residue number of ligand
fi #cor top
##################################################################


##################################################################
### anal module of amber
##################################################################
#	ligandname=$(grep $pdb $table | awk '{ print $2 } ')
#	ligand=$(grep "C1  $ligandname" $file | awk '{ print $5}')
	while read line  # For as many lines as the input file has...
	do
	  neu="$line"
	  neu="${neu//999/$lig_number}"
	  echo "$neu"
	done < $paradir"anal.in" >  $targetdir"anal.in" # reading template of minimize script and output into new file
	analin=$targetdir"anal.in"
	analin=${analin:=$paradir"anal.in"}
##################################################################
	if [ "$environment" = "EML_linux" ]; then
		echo anal -O -i $analin -o $targetdir$pdb$alt.outa -p $directory$pdb$alt.top -c $directory$pdb$alt.cor
		anal -O -i $analin -o $targetdir$pdb$alt.outa -p $directory$pdb$alt.top -c $directory$pdb$alt.cor
	elif [ "$environment" = "AZ_linux" ]; then
		echo $AMBEREXE/anal -O -i $analin -o $targetdir$pdb$alt.outa -p $directory$pdb$alt.top -c $directory$pdb$alt.cor
		$AMBEREXE/anal -O -i $analin -o $targetdir$pdb$alt.outa -p $directory$pdb$alt.top -c $directory$pdb$alt.cor
	elif [ "$environment" = "combine2go_linux" ]; then
		echo anal -O -i $analin -o $targetdir$pdb$alt.outa -p $topo -c $coordinates
#		anal -O -i $analin -o $targetdir$pdb$alt.outa -p $topo -c $coordinates
		anal -O -i $analin -o $targetdir$pdb$alt.outa -p $topo -c $coordinates
	else
		echo tleap version?
	fi #tleap environment
##################################################################
### end of anal module of amber
##################################################################


##################################################################
### write golpe file ($golpedat)
##################################################################
#	residue=$(grep $pdb $table | awk '{ print $2 } ')
if [ $molnumber == 1 ]; then
	#title
	echo Golpe dat for COMBINE analysis, $(date +"%D %T") > $golpedat
	lig_number=$[$(grep "C1  $residue" $count | awk '{ print $5 } ')]
	prot_components=$[lig_number-1]
	lig_components=1
	#number of variables (x+y)
	#ligand*protein*(vdw+ele)+(sum_vdw+sum_ele)+(999999.999+999999.999)+activity+deltaG
	echo $[($lig_components*$prot_components)*2+2+2+2] >> $golpedat 
	#number of objects
	echo $getcutpdb|wc| awk '{ print $2 } ' >> $golpedat
	##################################################################
	#write .dat.VarNames
	#Reste zaehlen
	echo $[($lig_components*$prot_components)*2+2+2+2] >  $golpenames 
#	echo $(awk '{if($1=="ATOM" && $3=="CA"){num++}}END{print (num*2)}' $file)  >  $golpenames
	#variable names as residue names and numbers
	awk  '($1=="ATOM" && $3=="CA") {print "v"$4_$5}' $file >>  $golpenames
	awk  '($1=="ATOM" && $3=="CA") {print "e"$4_$5}' $file >>  $golpenames
	echo "sum_vdw"  >>  $golpenames
	echo "sum_elec" >>  $golpenames
	echo "dummy"    >>  $golpenames
	echo "dummy"    >>  $golpenames
	echo "activity" >>  $golpenames
	echo "deltaG"   >>  $golpenames
	##################################################################
	#write .dat.VarType
	echo 7 > $golpevartyp #number of blocks

	echo 1 >> $golpevartyp
	echo $prot_components >> $golpevartyp
	echo 2 >> $golpevartyp
	echo " vdw" >> $golpevartyp
	echo 1 >> $golpevartyp #x variable

	echo $[prot_components+1] >> $golpevartyp
	echo $[prot_components*2] >> $golpevartyp
	echo 2 >> $golpevartyp
	echo " elec" >> $golpevartyp
	echo 1 >> $golpevartyp #x variable

	echo $[prot_components*2+1] >> $golpevartyp
	echo $[prot_components*2+1] >> $golpevartyp
	echo 2 >> $golpevartyp
	echo " sum_vdw" >> $golpevartyp
	echo 3 >> $golpevartyp #not used

	echo $[prot_components*2+2] >> $golpevartyp
	echo $[prot_components*2+2] >> $golpevartyp
	echo 2 >> $golpevartyp
	echo " sum_elec" >> $golpevartyp
	echo 3 >> $golpevartyp #not used

	echo $[prot_components*2+3] >> $golpevartyp
	echo $[prot_components*2+4] >> $golpevartyp
	echo 2 >> $golpevartyp
	echo " dummy" >> $golpevartyp
	echo 3 >> $golpevartyp #not used

	echo $[prot_components*2+5] >> $golpevartyp
	echo $[prot_components*2+5] >> $golpevartyp
	echo 2 >> $golpevartyp
	echo " activity" >> $golpevartyp
	echo 3 >> $golpevartyp #y varible not used

	echo $[prot_components*2+6] >> $golpevartyp
	echo $[prot_components*2+6] >> $golpevartyp
	echo 2 >> $golpevartyp
	echo " deltaG" >> $golpevartyp
	echo 2 >> $golpevartyp #y varible
	##################################################################
fi
#sequential number of molecule
echo $molnumber >> $golpedat
#name of molecule
echo $residue_$pdb$alt >> $golpedat
#x-variables
$prgdir"anal2golpe" $targetdir$pdb$alt.outa $prot_components $lig_components >> $golpedat
echo anal2golpe $targetdir$pdb$alt.outa $prot_components $lig_components
#y-variables
#activity=$(grep $pdb $table | awk '{ print $6 } ')
activity=$(grepsdtag.linux -t residue $residue < $sdf | extract_prop_sdf -property_name $prop_activity)
if [ "$activity" = "undef" ]; then activity=888888.888;fi
#dummy lines
echo 999999.999 >> $golpedat
echo 999999.999 >> $golpedat
#activity
echo $activity >> $golpedat
#deltaG
#kJ/mol
if [ "$activity" = "888888.888" ]; then 
	echo "$activity" >> $golpedat
	echo -------------------------------------------------
	echo "No $prop_activity activity is given for $residue"
	echo -------------------------------------------------
	echo "No $prop_activity activity is given for $residue" >> err.log
else
	echo "298*0.008314*l($activity*10^-6)" | bc -l >> $golpedat
fi # no activity
##################################################################
### end of write golpe file ($golpedat)
##################################################################

##################################################################
### calculate desolvation free energy with uhbd
##################################################################
if [ "$desolvation" = "y" ]; then
	$prgdir"amber7toqcd" $directory$pdb$alt.top $directory$pdb$alt.cor $directory$pdb$alt.qcd
	grep $residue $directory$pdb$alt.qcd > $directory$pdb$alt-ligand.qcd
	grep -v $residue $directory$pdb$alt.qcd > $directory$pdb$alt-protein.qcd

#	prepare uhbd
#	uhbd < desolv.inp > desolv.out

fi #desolvation
##################################################################
### end of calculation of desolvation free energy with uhbd
##################################################################

molnumber=$[molnumber+1]
rm -f $renumb
echo +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
echo
done