voronota-voromqa

#!/bin/bash

function print_help_and_exit
{
cat >&2 << EOF

'voronota-voromqa' script is an implementation of VoroMQA method using Voronota.

Basic options:
    --input | -i                   string      input structure file in PDB or mmCIF format
    --input-filter-query           string      input atoms filtering query parameters
    --output-atom-scores           string      output text file with atom scores
    --output-atom-scores-pdb       string      output PDB file with atom scores as B-factors
    --output-residue-scores        string      output text file with residue scores
    --output-residue-scores-pdb    string      output PDB file with residue scores as B-factors
    --output-residue-scores-plot   string      output PNG image file with residue scores plot, requires R
    --help | -h                                flag to display help message and exit

Advanced options:
    --cache-dir                    string      path to cache directory
    --smoothing-window             number      residue scores smoothing window size, default is 5
    --atoms-query                  string      atoms query parameters to define selection
    --contacts-query               string      contacts query parameters to define selection
    --output-selected-scores       string      output text file with selected atom scores
    --reference-sequence           string      sequence file or string for residue renumbering
    --output-sequence-alignment    string      output text file with sequence alignment
    --print-header                             flag to print output header
    --multiple-models                          flag to handle multiple models in PDB file
    --score-inter-chain                        flag to output inter-chain interface scores
    --list-all-options                         flag to display list of all command line options and exit

Standard output (one line):
    {input file path} {global score} {number of residues} {number of atoms} [ {selection score} {number of selected atoms} ]

EOF
exit 1
}

function print_all_options_and_exit
{
cat << EOF

All options of 'voronota-voromqa' script, in alphabetical order:

    --add-hydrogens                            string      command to add hydrogens
    --atoms-query                              string      atoms query parameters to define selection
    --cache-dir                                string      path to cache directory
    --contacts-query                           string      contacts query parameters to define selection
    --energy-input-potential                   string      path to custom potential file
    --energy-input-statistics                  string      path to custom energy statistics file
    --energy-mode                              string      energy mode name
    --help                                                 flag to display help message and exit
    --highlight-selection-in-plot                          flag to highlight selection in plot
    --input-filter-query                       string      input atoms filtering query parameters
    --input-is-list                                        flag to treat input as a list of files
    --input-is-structure                                   flag to treat input as a PDB file
    --input                                    string      input structure file in PDB format
    --more-logging                                         flag to enable more logging
    --multiple-energy-modes                    string      list of energy modes
    --multiple-models                                      flag to handle multiple models in PDB file
    --mutation-goal                            string      three-letter residue code
    --mutation-query                           string      atoms query parameters to define residues to mutate
    --neglect-SAS                                          flag to discard solvent-accessible surface
    --no-multiple-energy-modes                             flag to not use multiple energy modes
    --no-split-into-models                                 flag to not split input into separate models
    --no-use-slurm                                         flag to not use Slurm
    --options-for-calculating-contacts         string      custom options for calculating contacts
    --output-atom-depth-values                 string      output file with atom depth values
    --output-atom-energies                     string      output file with atom energy values
    --output-atom-scores-pdb                   string      output PDB file with atom scores as B-factors
    --output-atom-scores                       string      output text file with atom scores
    --output-balls                             string      output file with balls
    --output-contacts-map-svg                  string      output file with contact map image
    --output-contacts                          string      output file with contacts
    --output-cutting-suggestions               string      output file with chain cutting suggestions
    --output-directory-for-global-scores       string      output directory for global scores
    --output-residue-scores-for-CAMEO          string      output file with output for CAMEO
    --output-residue-scores-like-for-CASP      string      output file with output like for CASP
    --output-residue-scores-pdb                string      output PDB file with residue scores as B-factors
    --output-residue-scores-plot               string      output PNG image file with residue scores plot, requires R
    --output-residue-scores                    string      output text file with residue scores
    --output-scores-for-CASP-log               string      output file with logs generated in CASP mode
    --output-scores-for-CASP-pdb               string      output file with structure generated in CASP mode
    --output-scores-for-CASP                   string      output file with scores generated in CASP mode
    --output-screenshot                        string      output file with screenshot
    --output-selected-scores                   string      output text file with selected atom scores
    --output-selected-scores-by-residue        string      output text file with selected scores summed by residue
    --output-sequence-alignment                string      output text file with sequence alignment
    --output-summary-for-potential             string      output file with statistics of contact
    --print-clash-score-of-contacts-selection              flag to print clash score of contacts selection
    --print-energy-of-contacts-selection                   flag to print energy of contacts selection
    --print-header                                         flag to print output header
    --print-horizontally                                   flag to print output horizontally
    --print-vertically                                     flag to print output vertically
    --probe-radius                             number      probe radius for defining contacts
    --processors                               number      maximum number of processors to use
    --rebuild-sidechains                       string      command to rebuild side-chains
    --reference-sequence                       string      sequence file or string for residue renumbering
    --reinterpret-SAS                                      flag to interpret solvent as hydrophobic atoms
    --sbatch-parameters                        string      parameters for sbatch when using slurm
    --score-inter-chain                                    flag to output inter-chain interface scores
    --screenshot-background                    string      screenshot background color
    --screenshot-height                        number      screenshot image height
    --screenshot-width                         number      screenshot image width
    --small-plot                                           flag to generate smaller residue scores plot
    --smoothing-window                         number      residue scores smoothing window size, default is 5
    --split-into-models                        string      path to directory for storing extracted models
    --strip-rotamers                                       flag to remove sidechains
    --surface-craving-atoms-energies-output    string      output text file with surface frustration results
    --surface-craving-atoms-energies-pdb       string      output PDB file with surface frustration results
    --surface-craving-atoms-output-pdb         string      output PDB file with atom scores as B-factors
    --surface-craving-atoms-output             string      output text file with atom scores
    --surface-craving-depth-value              number      depth for surface frustration calculation, default is 3
    --surface-craving-residues-output-pdb      string      output PDB file with residue scores as B-factors
    --surface-craving-residues-output          string      output text file with residue scores
    --w-surface-craving-residues-output-pdb    string      output PDB file with residue scores as B-factors
    --w-surface-craving-residues-output        string      output text file with residue scores
    --tag-peripheral-contacts                              flag to tag peripheral contacts
    --unteaching-dir                           string      path to directory with potential-defining files
    --use-slurm                                string      path to directory to use for Slurm logs

EOF
exit 1
}

function substitute_ids_in_filename
{
	SUBSTITUTE_BASENAME="$(basename "$1")"
	SUBSTITUTE_ENERGY_MODE="$(basename "$2")"
	SUBSTITUTE_TEMPLATE="$3"
	
	echo "$SUBSTITUTE_TEMPLATE" \
	| sed "s|-BASENAME-|$SUBSTITUTE_BASENAME|" \
	| sed "s|-ENERGYMODE-|$SUBSTITUTE_ENERGY_MODE|"
}

readonly ZEROARG=$0
ALLARGS=("$@")

if [[ $ZEROARG == *"/"* ]]
then
	cd "$(dirname ${ZEROARG})"
	export PATH="$(pwd):${PATH}"
	cd - &> /dev/null
fi

export LC_ALL=C

command -v voronota &> /dev/null || { echo >&2 "Error: 'voronota' executable not in binaries path"; exit 1; }
command -v voronota-resources &> /dev/null || { echo >&2 "Error: 'voronota-resources' executable not in binaries path"; exit 1; }

INFILE=""
REFERENCE_SEQUENCE=""
INPUT_FILTER_QUERY_PARAMETERS=""
OUTFILE_SEQUENCE_ALIGNMENT=""
OUTFILE_BALLS=""
OUTFILE_CONTACTS=""
OUTFILE_ATOM_ENERGIES=""
OUTFILE_ATOM_SCORES=""
OUTFILE_ATOM_SCORES_PDB=""
OUTFILE_RESIDUE_SCORES=""
OUTFILE_RESIDUE_SCORES_PDB=""
OUTFILE_ATOM_DEPTHS=""
OUTFILE_CAMEO_RESIDUE_SCORES=""
OUTFILE_CASPLIKE_RESIDUE_SCORES=""
OUTFILE_PLOT=""
OUTFILE_MAP=""
OUTFILE_SCREENSHOT=""
OUTFILE_SCORES_FOR_CASP=""
OUTFILE_SCORES_FOR_CASP_PDB=""
OUTFILE_SCORES_FOR_CASP_LOG=""
OUTFILE_CUTTING_SUGGESTIONS=""
OUTFILE_SUMMARY_FOR_POTENTIAL=""
OUTDIR_GLOBAL_SCORES=""
SMOOTHING_WINDOW="5"
ATOMS_QUERY_PARAMETERS=""
CONTACTS_QUERY_PARAMETERS=""
OUTFILE_SELECTED_ATOM_SCORES=""
OUTFILE_SELECTED_RESIDUE_SCORES=""
CACHE_DIRECTORY=""
MULTIPLE_MODELS_CHAINS_OPTION=""
UNTEACHING_DIRECTORY=""
REBUILD_SIDECHAINS_COMMAND=""
ADD_HYDROGENS_COMMAND=""
SCREENSHOT_WIDTH="500"
SCREENSHOT_HEIGHT="500"
SCREENSHOT_BACKGROUND="black"
PLOT_MODE="normal"
SURFACE_CRAVING_ATOMS_OUTPUT=""
SURFACE_CRAVING_ATOMS_OUTPUT_PDB=""
SURFACE_CRAVING_RESIDUES_OUTPUT=""
SURFACE_CRAVING_RESIDUES_OUTPUT_PDB=""
SURFACE_CRAVING_ATOMS_ENERGIES_OUTPUT=""
SURFACE_CRAVING_ATOMS_ENERGIES_PDB=""
SURFACE_CRAVING_DEPTH_VALUE="3"
WEIGHTED_SURFACE_CRAVING_RESIDUES_OUTPUT=""
WEIGHTED_SURFACE_CRAVING_RESIDUES_OUTPUT_PDB=""
ENERGY_MODE="voromqa_v1"
MULTIPLE_ENERGY_MODES=""
ENERGY_INPUT_POTENTIAL=""
ENERGY_INPUT_STATISTICS=""
PROBE_RADIUS="1.4"
MUTATION_QUERY=""
MUTATION_GOAL=""
OPTIONS_FOR_CALCULATING_CONTACTS=""
MAX_PROCESSORS="1"
USE_SLURM_DIR=""
SBATCH_PARAMETERS=""
SPLIT_INTO_MODELS_DIR=""
TAG_PERIPHERAL_CONTACTS=""
STRIP_ROTAMERS=false
NEGLECT_SAS=false
REINTERPRET_SAS=false
MORE_LOGGING=false
HIGHLIGHT_SELECTION_IN_PLOT=false
PRINT_ENERGY_OF_CONTACTS_SELECTION=false
PRINT_CLASH_SCORE_OF_CONTACTS_SELECTION=false
PRINT_GLOBAL_SCORES_HEADER=false
PRINT_GLOBAL_SCORES_VERTICALLY=false
TREAT_INPUT_AS_LIST_OF_FILES=false
SCORE_INTER_CHAIN=false
LIST_ALL_OPTIONS=false
HELP_MODE=false

while [[ $# > 0 ]]
do
	OPTION="$1"
	OPTARG="$2"
	shift
	case $OPTION in
	-i|--input)
		INFILE="$OPTARG"
		shift
		;;
	--reference-sequence)
		REFERENCE_SEQUENCE="$OPTARG"
		shift
		;;
	--input-filter-query)
		INPUT_FILTER_QUERY_PARAMETERS="$OPTARG"
		shift
		;;
	--output-balls)
		OUTFILE_BALLS="$OPTARG"
		shift
		;;
	--output-contacts)
		OUTFILE_CONTACTS="$OPTARG"
		shift
		;;
	--output-sequence-alignment)
		OUTFILE_SEQUENCE_ALIGNMENT="$OPTARG"
		shift
		;;
	--output-atom-energies)
		OUTFILE_ATOM_ENERGIES="$OPTARG"
		shift
		;;
	--output-atom-scores)
		OUTFILE_ATOM_SCORES="$OPTARG"
		shift
		;;
	--output-atom-scores-pdb)
		OUTFILE_ATOM_SCORES_PDB="$OPTARG"
		shift
		;;
	--output-residue-scores)
		OUTFILE_RESIDUE_SCORES="$OPTARG"
		shift
		;;
	--output-residue-scores-pdb)
		OUTFILE_RESIDUE_SCORES_PDB="$OPTARG"
		shift
		;;
	--output-atom-depth-values)
		OUTFILE_ATOM_DEPTHS="$OPTARG"
		shift
		;;
	--output-residue-scores-for-CAMEO)
		OUTFILE_CAMEO_RESIDUE_SCORES="$OPTARG"
		shift
		;;
	--output-residue-scores-like-for-CASP)
		OUTFILE_CASPLIKE_RESIDUE_SCORES="$OPTARG"
		shift
		;;
	--output-residue-scores-plot)
		OUTFILE_PLOT="$OPTARG"
		shift
		;;
	--output-contacts-map-svg)
		OUTFILE_MAP="$OPTARG"
		shift
		;;
	--output-screenshot)
		OUTFILE_SCREENSHOT="$OPTARG"
		shift
		;;
	--output-scores-for-CASP)
		OUTFILE_SCORES_FOR_CASP="$OPTARG"
		shift
		;;
	--output-scores-for-CASP-pdb)
		OUTFILE_SCORES_FOR_CASP_PDB="$OPTARG"
		shift
		;;
	--output-scores-for-CASP-log)
		OUTFILE_SCORES_FOR_CASP_LOG="$OPTARG"
		shift
		;;
	--output-cutting-suggestions)
		OUTFILE_CUTTING_SUGGESTIONS="$OPTARG"
		shift
		;;
	--output-summary-for-potential)
		OUTFILE_SUMMARY_FOR_POTENTIAL="$OPTARG"
		shift
		;;
	--output-directory-for-global-scores)
		OUTDIR_GLOBAL_SCORES="$OPTARG"
		shift
		;;
	--smoothing-window)
		SMOOTHING_WINDOW="$OPTARG"
		shift
		;;
	--atoms-query)
		ATOMS_QUERY_PARAMETERS="$OPTARG"
		shift
		;;
	--contacts-query)
		CONTACTS_QUERY_PARAMETERS="$OPTARG"
		shift
		;;
	--output-selected-scores)
		OUTFILE_SELECTED_ATOM_SCORES="$OPTARG"
		shift
		;;
	--output-selected-scores-by-residue)
		OUTFILE_SELECTED_RESIDUE_SCORES="$OPTARG"
		shift
		;;
	--cache-dir)
		CACHE_DIRECTORY="$OPTARG"
		shift
		;;
	--unteaching-dir)
		UNTEACHING_DIRECTORY="$OPTARG"
		shift
		;;
	--rebuild-sidechains)
		REBUILD_SIDECHAINS_COMMAND="$OPTARG"
		shift
		;;
	--add-hydrogens)
		ADD_HYDROGENS_COMMAND="$OPTARG"
		shift
		;;
	--screenshot-width)
		SCREENSHOT_WIDTH="$OPTARG"
		shift
		;;
	--screenshot-height)
		SCREENSHOT_HEIGHT="$OPTARG"
		shift
		;;
	--screenshot-background)
		SCREENSHOT_BACKGROUND="$OPTARG"
		shift
		;;
	--surface-craving-atoms-output)
		SURFACE_CRAVING_ATOMS_OUTPUT="$OPTARG"
		shift
		;;
	--surface-craving-atoms-output-pdb)
		SURFACE_CRAVING_ATOMS_OUTPUT_PDB="$OPTARG"
		shift
		;;
	--surface-craving-residues-output)
		SURFACE_CRAVING_RESIDUES_OUTPUT="$OPTARG"
		shift
		;;
	--surface-craving-residues-output-pdb)
		SURFACE_CRAVING_RESIDUES_OUTPUT_PDB="$OPTARG"
		shift
		;;
	--surface-craving-atoms-energies-output)
		SURFACE_CRAVING_ATOMS_ENERGIES_OUTPUT="$OPTARG"
		shift
		;;
	--surface-craving-atoms-energies-pdb)
		SURFACE_CRAVING_ATOMS_ENERGIES_PDB="$OPTARG"
		shift
		;;
	--surface-craving-depth-value)
		SURFACE_CRAVING_DEPTH_VALUE="$OPTARG"
		shift
		;;
	--w-surface-craving-residues-output)
		WEIGHTED_SURFACE_CRAVING_RESIDUES_OUTPUT="$OPTARG"
		shift
		;;
	--w-surface-craving-residues-output-pdb)
		WEIGHTED_SURFACE_CRAVING_RESIDUES_OUTPUT_PDB="$OPTARG"
		shift
		;;
	--energy-mode)
		ENERGY_MODE="$OPTARG"
		shift
		;;
	--multiple-energy-modes)
		MULTIPLE_ENERGY_MODES="$OPTARG"
		shift
		;;
	--no-multiple-energy-modes)
		MULTIPLE_ENERGY_MODES=""
		;;
	--energy-input-potential)
		ENERGY_INPUT_POTENTIAL="$OPTARG"
		shift
		;;
	--energy-input-statistics)
		ENERGY_INPUT_STATISTICS="$OPTARG"
		shift
		;;
	--probe-radius)
		PROBE_RADIUS="$OPTARG"
		shift
		;;
	--mutation-query)
		MUTATION_QUERY="$OPTARG"
		shift
		;;
	--mutation-goal)
		MUTATION_GOAL="$OPTARG"
		shift
		;;
	--options-for-calculating-contacts)
		OPTIONS_FOR_CALCULATING_CONTACTS="$OPTARG"
		shift
		;;
	--processors)
		MAX_PROCESSORS="$OPTARG"
		shift
		;;
	--use-slurm)
		USE_SLURM_DIR="$OPTARG"
		shift
		;;
	--no-use-slurm)
		USE_SLURM_DIR=""
		;;
	--sbatch-parameters)
		SBATCH_PARAMETERS="$OPTARG"
		shift
		;;
	--split-into-models)
		SPLIT_INTO_MODELS_DIR="$OPTARG"
		shift
		;;
	--no-split-into-models)
		SPLIT_INTO_MODELS_DIR=""
		;;
	--tag-peripheral-contacts)
		TAG_PERIPHERAL_CONTACTS="--tag-peripherial"
		;;
	--multiple-models)
		MULTIPLE_MODELS_CHAINS_OPTION="--multimodel-chains"
		;;
	--small-plot)
		PLOT_MODE="small"
		;;
	--strip-rotamers)
		STRIP_ROTAMERS=true
		;;
	--neglect-SAS)
		NEGLECT_SAS=true
		;;
	--reinterpret-SAS)
		REINTERPRET_SAS=true
		;;
	--more-logging)
		MORE_LOGGING=true
		;;
	--highlight-selection-in-plot)
		HIGHLIGHT_SELECTION_IN_PLOT=true
		;;
	--print-energy-of-contacts-selection)
		PRINT_ENERGY_OF_CONTACTS_SELECTION=true
		;;
	--print-clash-score-of-contacts-selection)
		PRINT_CLASH_SCORE_OF_CONTACTS_SELECTION=true
		;;
	--print-header)
		PRINT_GLOBAL_SCORES_HEADER=true
		;;
	--print-vertically)
		PRINT_GLOBAL_SCORES_VERTICALLY=true
		;;
	--print-horizontally)
		PRINT_GLOBAL_SCORES_VERTICALLY=false
		;;
	--input-is-list)
		TREAT_INPUT_AS_LIST_OF_FILES=true
		;;
	--input-is-structure)
		TREAT_INPUT_AS_LIST_OF_FILES=false
		;;
	--score-inter-chain)
		SCORE_INTER_CHAIN=true
		;;
	--list-all-options)
		LIST_ALL_OPTIONS=true
		;;
	-h|--help)
		HELP_MODE=true
		;;
	*)
		echo >&2 "Error: invalid command line option '$OPTION'"
		exit 1
		;;
	esac
done

if $HELP_MODE
then
	print_help_and_exit
fi

if $LIST_ALL_OPTIONS
then
	print_all_options_and_exit
fi

if [ -z "$INFILE" ] && [ ! -t 0 ]
then
	INFILE="-"
fi

if [ -z "$INFILE" ]
then
	echo >&2 "Error: no input"
	exit 1
fi

if [ "$ENERGY_MODE" != "voromqa_v1" ] && [ "$ENERGY_MODE" != "voromqa_v2b" ]
then
	echo >&2 "Error: unsupported energy mode '$ENERGY_MODE'"
	exit 1
fi

if [ -z "$ADD_HYDROGENS_COMMAND" ]
then
	if [ "$ENERGY_MODE" == "voromqa_v2b" ] || [ "$ENERGY_MODE" == "voromqa_v3b" ]
	then
		ADD_HYDROGENS_COMMAND="reduce"
	fi
fi

if [ -z "$TAG_PERIPHERAL_CONTACTS" ]
then
	if [ "$ENERGY_MODE" == "voromqa_v3a" ] || [ "$ENERGY_MODE" == "voromqa_v3b" ] || [ "$ENERGY_MODE" == "voromqa_v5" ]
	then
		TAG_PERIPHERAL_CONTACTS="--tag-peripherial"
	fi
fi

if [ -z "$PROBE_RADIUS" ]
then
	if [ "$ENERGY_MODE" == "voromqa_v4" ] || [ "$ENERGY_MODE" == "voromqa_v5" ]
	then
		PROBE_RADIUS="3.5"
	fi
fi

if [ "$STRIP_ROTAMERS" != "true" ]
then
	if [ "$ENERGY_MODE" == "voromqa_v4" ] || [ "$ENERGY_MODE" == "voromqa_v5" ]
	then
		STRIP_ROTAMERS="true"
	fi
fi

if [ -n "$ENERGY_INPUT_POTENTIAL" ] && [ ! -s "$ENERGY_INPUT_POTENTIAL" ]
then
	echo >&2 "Error: no file '$ENERGY_INPUT_POTENTIAL'"
	exit 1
fi

if [ -n "$ENERGY_INPUT_STATISTICS" ] && [ ! -s "$ENERGY_INPUT_STATISTICS" ]
then
	echo >&2 "Error: no file '$ENERGY_INPUT_STATISTICS'"
	exit 1
fi

if [ -n "$OUTFILE_CAMEO_RESIDUE_SCORES" ] && [ -n "$MULTIPLE_MODELS_CHAINS_OPTION" ]
then
	echo >&2 "Error: '--output-residue-scores-for-CAMEO' and '--multiple-models' cannot be used together"
	exit 1
fi

if [ -n "$REBUILD_SIDECHAINS_COMMAND" ] && [ -n "$MULTIPLE_MODELS_CHAINS_OPTION" ]
then
	echo >&2 "Error: '--rebuild-side-chains' and '--multiple-models' cannot be used together"
	exit 1
fi

if [ -n "$REBUILD_SIDECHAINS_COMMAND" ] && [ -n "$OUTFILE_CAMEO_RESIDUE_SCORES" ]
then
	echo >&2 "Error: '--rebuild-sidechains' and '--output-residue-scores-for-CAMEO' cannot be used together"
	exit 1
fi

if [ -z "$REFERENCE_SEQUENCE" ] && [ -n "$OUTFILE_SEQUENCE_ALIGNMENT" ]
then
	echo >&2 "Error: '--output-sequence-alignment' cannot be used without '--reference-sequence'"
	exit 1
fi

if [ -z "$REFERENCE_SEQUENCE" ] && [ -n "$OUTFILE_SCORES_FOR_CASP" ]
then
	echo >&2 "Error: '--output-scores-for-CASP' cannot be used without '--reference-sequence'"
	exit 1
fi

if [ -z "$OUTFILE_SCORES_FOR_CASP" ] && [ -n "$OUTFILE_SCORES_FOR_CASP_PDB" ]
then
	echo >&2 "Error: '--output-scores-for-CASP-pdb' cannot be used without '--output-scores-for-CASP'"
	exit 1
fi

if [ -z "$OUTFILE_SCORES_FOR_CASP" ] && [ -n "$OUTFILE_SCORES_FOR_CASP_LOG" ]
then
	echo >&2 "Error: '--output-scores-for-CASP-log' cannot be used without '--output-scores-for-CASP'"
	exit 1
fi

if [ -n "$MUTATION_QUERY" ] && [ -z "$REBUILD_SIDECHAINS_COMMAND" ] && [ "$STRIP_ROTAMERS" != "true" ]
then
	echo >&2 "Error: '--mutation-query' cannot be used without '--rebuild-sidechains'"
	exit 1
fi

if [ -n "$MUTATION_QUERY" ] && [ -z "$MUTATION_GOAL" ]
then
	echo >&2 "Error: '--mutation-query' cannot be used without '--mutation-goal'"
	exit 1
fi

if [ -n "$MUTATION_GOAL" ] && [ -z "$MUTATION_QUERY" ]
then
	echo >&2 "Error: '--mutation-goal' cannot be used without '--mutation-query'"
	exit 1
fi

if $SCORE_INTER_CHAIN
then
	CONTACTS_QUERY_PARAMETERS="$CONTACTS_QUERY_PARAMETERS --no-same-chain --no-solvent"
	PRINT_ENERGY_OF_CONTACTS_SELECTION=true
	PRINT_CLASH_SCORE_OF_CONTACTS_SELECTION=true
	PRINT_GLOBAL_SCORES_HEADER=true
fi

MD5SUM_COMMAND="md5sum"
if command -v md5sum &> /dev/null
then
	MD5SUM_COMMAND="md5sum"
else
	MD5SUM_COMMAND="md5"
fi
command -v $MD5SUM_COMMAND &> /dev/null || { echo >&2 "Error: 'md5sum' or 'md5' executable not in binaries path"; exit 1; }


if [ -n "$OUTFILE_PLOT" ]
then
	command -v R &> /dev/null || { echo >&2 "Error: 'R' command, needed for plotting, is not available"; exit 1; }
fi

if [ -n "$OUTFILE_SCREENSHOT" ]
then
	command -v pymol &> /dev/null || { echo >&2 "Error: 'pymol' command, needed for making screenshots, is not available"; exit 1; }
fi

if [ -n "$REBUILD_SIDECHAINS_COMMAND" ]
then
	command -v $REBUILD_SIDECHAINS_COMMAND &> /dev/null || { echo >&2 "Error: '$REBUILD_SIDECHAINS_COMMAND' command, requested for rebuilding side-chains, is not available"; exit 1; }
fi

if [ -n "$ADD_HYDROGENS_COMMAND" ]
then
	command -v $ADD_HYDROGENS_COMMAND &> /dev/null || { echo >&2 "Error: '$ADD_HYDROGENS_COMMAND' command, requested for adding hydrogens, is not available"; exit 1; }
fi

if [ -n "$USE_SLURM_DIR" ]
then
	command -v sbatch &> /dev/null || { echo >&2 "Error: 'sbatch' command, needed for using Slurm, is not available"; exit 1; }
fi

readonly TMPLDIR=$(mktemp -d)
trap "rm -r $TMPLDIR" EXIT

voronota-resources radii > "$TMPLDIR/radii"

if [ ! -s "$TMPLDIR/radii" ]
then
	echo >&2 "Error: failed to get the predefined atomic radii"
	exit 1
fi

if [ -n "$REFERENCE_SEQUENCE" ] && [ ! -s "$REFERENCE_SEQUENCE" ]
then
	if [ "$REFERENCE_SEQUENCE" == "$(echo $REFERENCE_SEQUENCE | tr -dc '[:alpha:]')" ]
	then
		echo "$REFERENCE_SEQUENCE" > "$TMPLDIR/reference_sequence_file"
		REFERENCE_SEQUENCE="$TMPLDIR/reference_sequence_file"
	else
		echo >&2 "Error: invalid reference sequence"
		exit 1
	fi
fi

INFILE_CONTENTS=$INFILE

if [[ "$INFILE" == *"://"* ]]
then
	INFILE_CONTENTS=$TMPLDIR/downloaded_input_file
	curl -s "$INFILE" > $INFILE_CONTENTS
fi

if [ "$INFILE" == "-" ]
then
	INFILE_CONTENTS=$TMPLDIR/downloaded_input_file
	cat > $INFILE_CONTENTS
fi

if [ ! -s "$INFILE_CONTENTS" ]
then
	echo >&2 "Error: input file does not exist or is empty"
	exit 1
fi

INFILE_CONTENTS_EXTRACTED=$INFILE_CONTENTS

if [[ "$INFILE" == *".gz" ]]
then
	INFILE_CONTENTS_EXTRACTED=$TMPLDIR/extracted_input_file
	zcat $INFILE_CONTENTS > $INFILE_CONTENTS_EXTRACTED
	if [ ! -s "$INFILE_CONTENTS_EXTRACTED" ]
	then
		echo >&2 "Error: could not extract gzipped file"
		exit 1
	fi
fi

if $TREAT_INPUT_AS_LIST_OF_FILES
then
	mkdir -p "$TMPLDIR/children_global_scores"
	
	if [ -n "$USE_SLURM_DIR" ]
	then
		mkdir -p "$USE_SLURM_DIR"
		
		cat "$INFILE_CONTENTS_EXTRACTED" \
		| sort \
		| uniq \
		| xargs -L 1 sbatch -o "$USE_SLURM_DIR/slurmjob-%j.out" -e "$USE_SLURM_DIR/slurmjob-%j.err" $SBATCH_PARAMETERS \
		"$ZEROARG" "${ALLARGS[@]}" \
		  --print-header \
		  --print-horizontally \
		  --input-is-structure \
		  --input \
		| egrep '^Submitted batch job ' \
		| awk '{print $4}' \
		> "$TMPLDIR/slurm_job_ids"
		
		if [ ! -s "$TMPLDIR/slurm_job_ids" ]
		then
			echo >&2 "Error: no Slurm jobs submitted"
			exit 1
		fi
		
		REMAINING_SLURM_JOBS="$(squeue | grep -f "$TMPLDIR/slurm_job_ids" | wc -l)"
		while [ "$REMAINING_SLURM_JOBS" -gt "0" ]
		do
			sleep 5
			REMAINING_SLURM_JOBS="$(squeue | grep -f "$TMPLDIR/slurm_job_ids" | wc -l)"
		done
		
		cat "$TMPLDIR/slurm_job_ids" | while read SLURM_JOB_ID
		do
			cat "$USE_SLURM_DIR/slurmjob-${SLURM_JOB_ID}.out" > "$TMPLDIR/children_global_scores/result_${SLURM_JOB_ID}"
		done
	else
		cat "$INFILE_CONTENTS_EXTRACTED" \
		| sort \
		| uniq \
		| xargs -L 1 -P "$MAX_PROCESSORS" "$ZEROARG" "${ALLARGS[@]}" \
		  --print-header \
		  --print-horizontally \
		  --output-directory-for-global-scores "$TMPLDIR/children_global_scores" \
		  --input-is-structure \
		  --input \
		> /dev/null
	fi
	
	find "$TMPLDIR/children_global_scores" -type f -not -empty > "$TMPLDIR/global_score_files"
	
	MAX_NUMBER_OF_FIELDS="$(cat "$TMPLDIR/global_score_files" | xargs -L 1 awk '{print NF}' | awk '{if(MAX=="" || ($1+0)>MAX){MAX=$1}} END {print MAX}')"
	
	{
		cat "$TMPLDIR/global_score_files" | while read GSFILE
		do
			if [ "$(cat "$GSFILE" | tail -1 | awk '{print NF}')" -eq "$MAX_NUMBER_OF_FIELDS" ]
			then
				echo "$GSFILE"
			fi
		done
	} \
	| xargs -L 1 cat \
	| awk '{if(NR==1 || $1!="input") print $0}' \
	> "$TMPLDIR/unsorted_table"
	
	cat "$TMPLDIR/unsorted_table" | head -1
	cat "$TMPLDIR/unsorted_table" | tail -n +2 | sort
	
	exit 0
fi

if [ -n "$SPLIT_INTO_MODELS_DIR" ]
then
	mkdir -p "$SPLIT_INTO_MODELS_DIR"
	
	cat "$INFILE_CONTENTS_EXTRACTED" \
	| voronota x-split-atoms-file --prefix "${SPLIT_INTO_MODELS_DIR}/$(basename $INFILE)__" --postfix '.pdb' \
	| "$ZEROARG" "${ALLARGS[@]}" --input - --input-is-list --no-split-into-models \
	| \
	{
		if [ -n "$OUTDIR_GLOBAL_SCORES" ]
		then
			mkdir -p "$OUTDIR_GLOBAL_SCORES"
			cat | tee "$(mktemp "$OUTDIR_GLOBAL_SCORES/vgs.XXXXXXXXXX")"
		else
			cat
		fi
	}
	
	exit 0
fi

if [ -n "$MULTIPLE_ENERGY_MODES" ]
then
	echo "$MULTIPLE_ENERGY_MODES" \
	| sed 's/\s\+/\n/g' \
	| sed 's/,\+/\n/g' \
	| egrep '\S' \
	| sort \
	| uniq \
	| xargs -L 1 "$ZEROARG" "${ALLARGS[@]}" \
	  --print-header \
	  --print-vertically \
	  --no-multiple-energy-modes \
	  --output-directory-for-global-scores "" \
	  --input "$INFILE_CONTENTS_EXTRACTED" \
	  --energy-mode \
	| egrep -v '^input' \
	| \
	{
		echo "input $INFILE"
		cat
	} \
	| \
	{
		if $PRINT_GLOBAL_SCORES_VERTICALLY
		then
			cat
		else
			cat > "$TMPLDIR/multiple_global_scores"
			cat "$TMPLDIR/multiple_global_scores" | awk '{print $1}' | tr '\n' ' ' | sed 's/\s$/\n/'
			cat "$TMPLDIR/multiple_global_scores" | awk '{print $2}' | tr '\n' ' ' | sed 's/\s$/\n/'
		fi
	} \
	| \
	{
		if [ -n "$OUTDIR_GLOBAL_SCORES" ]
		then
			mkdir -p "$OUTDIR_GLOBAL_SCORES"
			cat | tee "$(mktemp "$OUTDIR_GLOBAL_SCORES/vgs.XXXXXXXXXX")"
		else
			cat
		fi
	}
	
	exit 0
fi

REFERENCE_SEQUENCE="$(substitute_ids_in_filename "$INFILE" "$ENERGY_MODE" "$REFERENCE_SEQUENCE")"
OUTFILE_SEQUENCE_ALIGNMENT="$(substitute_ids_in_filename "$INFILE" "$ENERGY_MODE" "$OUTFILE_SEQUENCE_ALIGNMENT")"
OUTFILE_BALLS="$(substitute_ids_in_filename "$INFILE" "$ENERGY_MODE" "$OUTFILE_BALLS")"
OUTFILE_CONTACTS="$(substitute_ids_in_filename "$INFILE" "$ENERGY_MODE" "$OUTFILE_CONTACTS")"
OUTFILE_ATOM_ENERGIES="$(substitute_ids_in_filename "$INFILE" "$ENERGY_MODE" "$OUTFILE_ATOM_ENERGIES")"
OUTFILE_ATOM_SCORES="$(substitute_ids_in_filename "$INFILE" "$ENERGY_MODE" "$OUTFILE_ATOM_SCORES")"
OUTFILE_ATOM_SCORES_PDB="$(substitute_ids_in_filename "$INFILE" "$ENERGY_MODE" "$OUTFILE_ATOM_SCORES_PDB")"
OUTFILE_RESIDUE_SCORES="$(substitute_ids_in_filename "$INFILE" "$ENERGY_MODE" "$OUTFILE_RESIDUE_SCORES")"
OUTFILE_RESIDUE_SCORES_PDB="$(substitute_ids_in_filename "$INFILE" "$ENERGY_MODE" "$OUTFILE_RESIDUE_SCORES_PDB")"
OUTFILE_ATOM_DEPTHS="$(substitute_ids_in_filename "$INFILE" "$ENERGY_MODE" "$OUTFILE_ATOM_DEPTHS")"
OUTFILE_CAMEO_RESIDUE_SCORES="$(substitute_ids_in_filename "$INFILE" "$ENERGY_MODE" "$OUTFILE_CAMEO_RESIDUE_SCORES")"
OUTFILE_CASPLIKE_RESIDUE_SCORES="$(substitute_ids_in_filename "$INFILE" "$ENERGY_MODE" "$OUTFILE_CASPLIKE_RESIDUE_SCORES")"
OUTFILE_PLOT="$(substitute_ids_in_filename "$INFILE" "$ENERGY_MODE" "$OUTFILE_PLOT")"
OUTFILE_MAP="$(substitute_ids_in_filename "$INFILE" "$ENERGY_MODE" "$OUTFILE_MAP")"
OUTFILE_SCREENSHOT="$(substitute_ids_in_filename "$INFILE" "$ENERGY_MODE" "$OUTFILE_SCREENSHOT")"
OUTFILE_SCORES_FOR_CASP="$(substitute_ids_in_filename "$INFILE" "$ENERGY_MODE" "$OUTFILE_SCORES_FOR_CASP")"
OUTFILE_SCORES_FOR_CASP_PDB="$(substitute_ids_in_filename "$INFILE" "$ENERGY_MODE" "$OUTFILE_SCORES_FOR_CASP_PDB")"
OUTFILE_SCORES_FOR_CASP_LOG="$(substitute_ids_in_filename "$INFILE" "$ENERGY_MODE" "$OUTFILE_SCORES_FOR_CASP_LOG")"
OUTFILE_CUTTING_SUGGESTIONS="$(substitute_ids_in_filename "$INFILE" "$ENERGY_MODE" "$OUTFILE_CUTTING_SUGGESTIONS")"
OUTFILE_SUMMARY_FOR_POTENTIAL="$(substitute_ids_in_filename "$INFILE" "$ENERGY_MODE" "$OUTFILE_SUMMARY_FOR_POTENTIAL")"
OUTFILE_SELECTED_ATOM_SCORES="$(substitute_ids_in_filename "$INFILE" "$ENERGY_MODE" "$OUTFILE_SELECTED_ATOM_SCORES")"
OUTFILE_SELECTED_RESIDUE_SCORES="$(substitute_ids_in_filename "$INFILE" "$ENERGY_MODE" "$OUTFILE_SELECTED_RESIDUE_SCORES")"
SURFACE_CRAVING_ATOMS_OUTPUT="$(substitute_ids_in_filename "$INFILE" "$ENERGY_MODE" "$SURFACE_CRAVING_ATOMS_OUTPUT")"
SURFACE_CRAVING_ATOMS_OUTPUT_PDB="$(substitute_ids_in_filename "$INFILE" "$ENERGY_MODE" "$SURFACE_CRAVING_ATOMS_OUTPUT_PDB")"
SURFACE_CRAVING_RESIDUES_OUTPUT="$(substitute_ids_in_filename "$INFILE" "$ENERGY_MODE" "$SURFACE_CRAVING_RESIDUES_OUTPUT")"
SURFACE_CRAVING_RESIDUES_OUTPUT_PDB="$(substitute_ids_in_filename "$INFILE" "$ENERGY_MODE" "$SURFACE_CRAVING_RESIDUES_OUTPUT_PDB")"
SURFACE_CRAVING_ATOMS_ENERGIES_OUTPUT="$(substitute_ids_in_filename "$INFILE" "$ENERGY_MODE" "$SURFACE_CRAVING_ATOMS_ENERGIES_OUTPUT")"
SURFACE_CRAVING_ATOMS_ENERGIES_PDB="$(substitute_ids_in_filename "$INFILE" "$ENERGY_MODE" "$SURFACE_CRAVING_ATOMS_ENERGIES_PDB")"
WEIGHTED_SURFACE_CRAVING_RESIDUES_OUTPUT="$(substitute_ids_in_filename "$INFILE" "$ENERGY_MODE" "$WEIGHTED_SURFACE_CRAVING_RESIDUES_OUTPUT")"
WEIGHTED_SURFACE_CRAVING_RESIDUES_OUTPUT_PDB="$(substitute_ids_in_filename "$INFILE" "$ENERGY_MODE" "$WEIGHTED_SURFACE_CRAVING_RESIDUES_OUTPUT_PDB")"

cat $INFILE_CONTENTS_EXTRACTED \
| voronota get-balls-from-atoms-file \
  --input-format detect \
  --annotated $MULTIPLE_MODELS_CHAINS_OPTION \
  --radii-file $TMPLDIR/radii \
  --include-heteroatoms \
| voronota query-balls \
  --match 'R<LEU,ALA,GLY,VAL,GLU,SER,LYS,ILE,ASP,THR,ARG,PRO,ASN,PHE,GLN,TYR,HIS,MET,TRP,CYS,MSE>' \
| voronota query-balls \
  --drop-adjuncts \
  --drop-altloc-indicators \
| \
{
	if [ -n "$INPUT_FILTER_QUERY_PARAMETERS" ]
	then
		voronota query-balls $INPUT_FILTER_QUERY_PARAMETERS
	else
		cat
	fi
} \
| \
{
	if [ -s "$REFERENCE_SEQUENCE" ]
	then
		voronota query-balls \
		  --set-ref-seq-num-adjunct "$REFERENCE_SEQUENCE" \
		  --ref-seq-alignment $TMPLDIR/sequence_alignment \
		| voronota query-balls \
		  --renumber-from-adjunct refseq \
		| voronota query-balls \
		  --drop-adjuncts
	else
		cat
	fi
} \
| \
{
	if [ -n "$MUTATION_QUERY" ]
	then
		voronota query-balls $MUTATION_QUERY \
		  --set-tags mutating \
		| sed 's/\(.*R<\)\(...\)\(.*mutating.*\)/\1XXX\3/' \
		| sed "s/R<XXX>/R<$MUTATION_GOAL>/" \
		| sed 's/;mutating;/;/' \
		| sed 's/;mutating//' \
		| sed 's/mutating;//' \
		| sed 's/mutating/./'
	else
		cat
	fi
} \
| \
{
	if [ -n "$REBUILD_SIDECHAINS_COMMAND" ]
	then
		voronota query-balls \
		  --match 'A<CA,C,N,O,OXT>' \
		| voronota write-balls-to-atoms-file \
		  --pdb-output $TMPLDIR/atoms_stripped.pdb \
		> /dev/null

		$REBUILD_SIDECHAINS_COMMAND -h \
		  -i $TMPLDIR/atoms_stripped.pdb \
		  -o $TMPLDIR/atoms_rebuilt.pdb \
		> /dev/null

		cat $TMPLDIR/atoms_rebuilt.pdb \
		| voronota get-balls-from-atoms-file \
		  --annotated \
		  --radii-file $TMPLDIR/radii \
		| voronota query-balls \
		  --drop-adjuncts
	else
		cat
	fi
} \
| \
{
	if [ -n "$ADD_HYDROGENS_COMMAND" ]
	then
		voronota query-balls \
		| voronota write-balls-to-atoms-file \
		  --pdb-output $TMPLDIR/atoms_heavy.pdb \
		> /dev/null

		$ADD_HYDROGENS_COMMAND \
		  -FLIP $TMPLDIR/atoms_heavy.pdb \
		> $TMPLDIR/atoms_heavy_and_hydrogens.pdb \
		2> /dev/null

		cat $TMPLDIR/atoms_heavy_and_hydrogens.pdb \
		| voronota get-balls-from-atoms-file \
		  --annotated $MULTIPLE_MODELS_CHAINS_OPTION \
		  --radii-file $TMPLDIR/radii \
		  --include-heteroatoms \
		  --include-hydrogens \
		| voronota query-balls \
		  --drop-adjuncts
	else
		cat
	fi
} \
| \
{
	if $STRIP_ROTAMERS
	then
		voronota query-balls --match 'A<C,N,O,CA,CB,OXT>'
	else
		cat
	fi
} \
> $TMPLDIR/balls

if $MORE_LOGGING && [ -f "$TMPLDIR/sequence_alignment" ]
then
	cat >&2 "$TMPLDIR/sequence_alignment"
fi

if [ -n "$OUTFILE_SEQUENCE_ALIGNMENT" ] && [ -f "$TMPLDIR/sequence_alignment" ]
then
	cp "$TMPLDIR/sequence_alignment" "$OUTFILE_SEQUENCE_ALIGNMENT"
fi

if [ ! -s "$TMPLDIR/balls" ]
then
	echo >&2 "Error: no atoms in input file"
	exit 1
fi

if [ -n "$OUTFILE_BALLS" ]
then
	mkdir -p $(dirname $OUTFILE_BALLS)
	cp "$TMPLDIR/balls" "$OUTFILE_BALLS"
fi

CACHE_MD5=""
if [ -n "$CACHE_DIRECTORY" ]
then
	CACHE_MD5=$(cat $TMPLDIR/balls <(echo "$PROBE_RADIUS" "$TAG_PERIPHERAL_CONTACTS $OPTIONS_FOR_CALCULATING_CONTACTS") | $MD5SUM_COMMAND | awk '{print $1}')
	if [ -n "$CACHE_MD5" ]
	then
		CACHE_MD5="${CACHE_MD5}.voronota.voromqa"
		if [ -d "$CACHE_DIRECTORY/$CACHE_MD5" ] && [ -s "$CACHE_DIRECTORY/$CACHE_MD5/all_contacts" ] && [ -s "$CACHE_DIRECTORY/$CACHE_MD5/atom_energies" ] && [ -s "$CACHE_DIRECTORY/$CACHE_MD5/depth_values" ]
		then
			cp -r $CACHE_DIRECTORY/$CACHE_MD5 $TMPLDIR/foundation
		fi
	fi
fi

mkdir -p $TMPLDIR/foundation
	
if [ ! -s "$TMPLDIR/foundation/all_contacts" ]
then
	cat $TMPLDIR/balls \
	| voronota calculate-contacts \
	  --annotated \
	  --probe "$PROBE_RADIUS" \
	  --tag-centrality $TAG_PERIPHERAL_CONTACTS $OPTIONS_FOR_CALCULATING_CONTACTS \
	| tee $TMPLDIR/foundation/raw_contacts \
	| voronota query-contacts \
	  --match-min-seq-sep 1 \
	| voronota query-contacts \
	  --match-first 'A<C>' \
	  --match-second 'A<N>' \
	  --match-max-seq-sep 1 \
	  --match-max-dist 1.6 \
	  --invert \
	| voronota query-contacts \
	  --match-min-seq-sep 1 \
	  --match-max-seq-sep 1 \
	  --set-tags 'sep1' \
	| voronota query-contacts \
	  --match-min-seq-sep 2 \
	  --no-solvent \
	  --set-tags 'sep2' \
	> $TMPLDIR/foundation/all_contacts
fi

cat $TMPLDIR/foundation/all_contacts \
| awk '{print $1 " " $2 " " $5 " " $3}' \
| tr ';' '_' \
> $TMPLDIR/contacts

if [ -n "$OUTFILE_CONTACTS" ]
then
	mkdir -p $(dirname $OUTFILE_CONTACTS)
	cp "$TMPLDIR/foundation/all_contacts" "$OUTFILE_CONTACTS"
fi

if [ -n "$OUTFILE_SUMMARY_FOR_POTENTIAL" ]
then
	mkdir -p $(dirname "$OUTFILE_SUMMARY_FOR_POTENTIAL")

	cat $TMPLDIR/balls \
	| voronota query-balls \
	  --chains-summary-output $TMPLDIR/chains_counts \
	> /dev/null

	cat $TMPLDIR/foundation/all_contacts \
	| voronota query-contacts \
	  --match-min-seq-sep 2 \
	  --match-max-seq-sep 6 \
	  --set-tags 'sep1' \
	| sed 's/sep1;sep2/sep1/' \
	| awk '{print $1 " " $2 " " $5 " " $3}' \
	| tr ';' '_' \
	| voronota score-contacts-potential \
	  --multiply-areas $(cat $TMPLDIR/chains_counts | head -1 | awk '{print (24.0/$1*$2)}') \
	> "$OUTFILE_SUMMARY_FOR_POTENTIAL"

	exit 0
fi

mkdir -p $TMPLDIR/resources

{
if [ -n "$ENERGY_INPUT_POTENTIAL" ]
then
	cat "$ENERGY_INPUT_POTENTIAL"
else
	voronota-resources ${ENERGY_MODE}_energy_potential
fi
} > $TMPLDIR/resources/energy_potential

{
if [ -n "$ENERGY_INPUT_STATISTICS" ]
then
	cat "$ENERGY_INPUT_STATISTICS"
else
	voronota-resources ${ENERGY_MODE}_energy_means_and_sds
fi
} > $TMPLDIR/resources/energy_means_and_sds

if [ ! -s "$TMPLDIR/foundation/atom_energies" ]
then
	cat $TMPLDIR/contacts \
	| voronota score-contacts-energy \
	  --potential-file $TMPLDIR/resources/energy_potential \
	  --atom-scores-file $TMPLDIR/foundation/atom_energies \
	> /dev/null
fi

if [ ! -s "$TMPLDIR/foundation/depth_values" ]
then
	cat $TMPLDIR/contacts \
	| voronota x-query-contacts-depth-values \
	> $TMPLDIR/foundation/depth_values
fi

if [ -n "$CACHE_DIRECTORY" ] && [ -n "$CACHE_MD5" ] && [ ! -d "$CACHE_DIRECTORY/$CACHE_MD5" ]
then
	mkdir -p $CACHE_DIRECTORY
	cp -r $TMPLDIR/foundation $CACHE_DIRECTORY/$CACHE_MD5
fi

if $NEGLECT_SAS
then
	cat $TMPLDIR/foundation/all_contacts \
	| awk '{print $1 " " $2 " " $5 " " $3}' \
	| tr ';' '_' \
	| sed 's|c<solvent>\s\+.\s\+\S\+$|c<solvent> . 0.001|' \
	| voronota score-contacts-energy \
	  --potential-file $TMPLDIR/resources/energy_potential \
	  --atom-scores-file $TMPLDIR/foundation/atom_energies \
	> /dev/null
fi

if $REINTERPRET_SAS
then
	cat $TMPLDIR/foundation/all_contacts \
	| awk '{print $1 " " $2 " " $5 " " $3}' \
	| tr ';' '_' \
	| voronota score-contacts-energy \
	  --potential-file <(cat $TMPLDIR/resources/energy_potential | voronota x-simulate-potential-for-membrane-proteins) \
	  --atom-scores-file $TMPLDIR/foundation/atom_energies \
	> /dev/null
fi

if [ -n "$OUTFILE_CUTTING_SUGGESTIONS" ]
then
	mkdir -p $(dirname "$OUTFILE_CUTTING_SUGGESTIONS")
	cat $TMPLDIR/foundation/all_contacts \
	| awk '{print $1 " " $2 " " $5 " " $3}' \
	| tr ';' '_' \
	| voronota x-score-contacts-global-energy-by-cuts \
	  --potential-file $TMPLDIR/resources/energy_potential \
	  --limit-results 10 \
	| column -t \
	> "$OUTFILE_CUTTING_SUGGESTIONS"
fi

if [ -n "$UNTEACHING_DIRECTORY" ] && [ -s "$UNTEACHING_DIRECTORY/contact_areas" ] && [ -s "$UNTEACHING_DIRECTORY/contact_categories" ]
then
	cat $TMPLDIR/balls \
	| voronota query-balls \
	  --chains-summary-output $TMPLDIR/chains_counts \
	> /dev/null

	{
		cat $UNTEACHING_DIRECTORY/contact_areas
		
		cat $TMPLDIR/foundation/all_contacts \
		| voronota query-contacts \
		  --match-min-seq-sep 2 \
		  --match-max-seq-sep 6 \
		  --set-tags 'sep1' \
		| sed 's/sep1;sep2/sep1/' \
		| awk '{print $1 " " $2 " " $5 " " $3}' \
		| tr ';' '_' \
		| voronota score-contacts-potential \
		  --multiply-areas $(cat $TMPLDIR/chains_counts | head -1 | awk '{print (24.0/$1*$2)}') \
		| awk '{print $1 " " $2 " " $3 " -" $4}'
	} \
	| voronota score-contacts-potential \
	  --input-contributions $UNTEACHING_DIRECTORY/contact_categories \
	  --potential-file $TMPLDIR/untaught_potential \
	> /dev/null

	cat $TMPLDIR/foundation/all_contacts \
	| awk '{print $1 " " $2 " " $5 " " $3}' \
	| tr ';' '_' \
	| voronota score-contacts-energy \
	  --potential-file $TMPLDIR/untaught_potential \
	  --atom-scores-file $TMPLDIR/foundation/atom_energies \
	> /dev/null
fi

if [ -n "$OUTFILE_ATOM_ENERGIES" ]
then
	mkdir -p $(dirname "$OUTFILE_ATOM_ENERGIES")
	cp $TMPLDIR/foundation/atom_energies "$OUTFILE_ATOM_ENERGIES"
fi

cat $TMPLDIR/foundation/atom_energies \
| voronota score-contacts-quality \
  --default-mean -0.34 \
  --default-sd 0.19 \
  --means-and-sds-file $TMPLDIR/resources/energy_means_and_sds \
  --external-weights-file $TMPLDIR/foundation/depth_values \
  --smoothing-window $SMOOTHING_WINDOW \
  --atom-scores-file $TMPLDIR/atom_quality_scores \
  --residue-scores-file $TMPLDIR/residue_quality_scores \
> $TMPLDIR/global_quality_score

if [ ! -s "$TMPLDIR/global_quality_score" ] || [ ! -s "$TMPLDIR/atom_quality_scores" ] || [ ! -s "$TMPLDIR/residue_quality_scores" ]
then
	echo >&2 "Error: failed to calculate scores"
	exit 1
fi

if [ -n "$ATOMS_QUERY_PARAMETERS" ] || [ -n "$CONTACTS_QUERY_PARAMETERS" ]
then
	cat $TMPLDIR/balls \
	| voronota query-balls \
	  --set-external-adjuncts $TMPLDIR/foundation/depth_values \
	  --set-external-adjuncts-name depth \
	| voronota query-balls $ATOMS_QUERY_PARAMETERS \
	| awk '{print $1}' \
	> $TMPLDIR/selection_by_atoms_query
	
	cat $TMPLDIR/foundation/all_contacts \
	| voronota query-contacts $CONTACTS_QUERY_PARAMETERS \
	| awk '{print $1 " " $2}' | tr ' ' '\n' | sort | uniq | grep -v solvent \
	> $TMPLDIR/selection_by_contacts_query
	
	join -e EMPTY <(sort -k 1b,1 $TMPLDIR/selection_by_atoms_query) <(sort -k 1b,1 $TMPLDIR/selection_by_contacts_query) \
	| grep -v EMPTY \
	> $TMPLDIR/selection_by_queries
	
	join -e EMPTY <(sort -k 1b,1 $TMPLDIR/foundation/depth_values) <(sort -k 1b,1 $TMPLDIR/selection_by_queries) \
	| grep -v EMPTY \
	| awk '{print $1 " " $2}' \
	> $TMPLDIR/scoring_domain
	
	join -e EMPTY <(sort -k 1b,1 $TMPLDIR/scoring_domain) <(sort -k 1b,1 $TMPLDIR/atom_quality_scores) \
	| grep -v EMPTY \
	> $TMPLDIR/scoring_table
	
	{
		if [ -s "$TMPLDIR/scoring_table" ]
		then
			cat $TMPLDIR/scoring_table | awk '{sum+=$2;wsum+=($2*$3)} END {print (wsum/sum)}'
			cat $TMPLDIR/scoring_table | wc -l
		else
			echo "0"
			echo "0"
		fi
	} > $TMPLDIR/global_quality_score_of_selection
fi

if [ -n "$OUTFILE_ATOM_SCORES" ]
then
	mkdir -p $(dirname "$OUTFILE_ATOM_SCORES")
	cp $TMPLDIR/atom_quality_scores "$OUTFILE_ATOM_SCORES"
fi

if [ -n "$OUTFILE_RESIDUE_SCORES" ]
then
	mkdir -p $(dirname "$OUTFILE_RESIDUE_SCORES")
	cp $TMPLDIR/residue_quality_scores "$OUTFILE_RESIDUE_SCORES"
fi

if [ -n "$OUTFILE_SELECTED_ATOM_SCORES" ] && [ -f "$TMPLDIR/scoring_table" ]
then
	mkdir -p $(dirname "$OUTFILE_SELECTED_ATOM_SCORES")
	cat $TMPLDIR/scoring_table \
	| awk '{print $1 " " $3}' \
	| sort -V \
	> "$OUTFILE_SELECTED_ATOM_SCORES"
fi

if [ -n "$OUTFILE_SELECTED_RESIDUE_SCORES" ] && [ -f "$TMPLDIR/scoring_table" ]
then
	mkdir -p $(dirname "$OUTFILE_SELECTED_RESIDUE_SCORES")
	cat $TMPLDIR/scoring_table \
	| awk '{print $1 " " $3}' \
	|voronota expand-descriptors \
	| awk '{print $1 "__and__" $2 "__and__" $3 "__and__" $6 " " $8}' \
	| awk '{sums[$1]+=$2; counts[$1]++}; END {for(id in sums){print id " " counts[id] " " sums[id]}}' \
	| sed 's/__and__/ /g' \
	| awk '{print "c<" $1 ">r<" $2 ">i<" $3 ">R<" $4 "> " $5 " " $6}' \
	| sed 's/.<\.>//g' \
	| sort -V \
	> "$OUTFILE_SELECTED_RESIDUE_SCORES"
fi

if [ -n "$OUTFILE_ATOM_SCORES_PDB" ]
then
	mkdir -p $(dirname "$OUTFILE_ATOM_SCORES_PDB")
	cat $TMPLDIR/balls \
	| voronota query-balls \
	  --drop-altloc-indicators \
	  --set-adjuncts 'oc=1' \
	  --set-external-adjuncts $TMPLDIR/atom_quality_scores \
	  --set-external-adjuncts-name score \
	| \
	{
		if [ -s "$TMPLDIR/selection_by_queries" ]
		then
			voronota query-balls \
			  --match-external-annotations $TMPLDIR/selection_by_queries \
			  --set-adjuncts 'oc=2'
		else
			cat
		fi
	} \
	| voronota write-balls-to-atoms-file \
	  --pdb-output "$OUTFILE_ATOM_SCORES_PDB" \
	  --pdb-output-b-factor score \
	> /dev/null
fi

if [ -n "$OUTFILE_RESIDUE_SCORES_PDB" ]
then
	mkdir -p $(dirname "$OUTFILE_RESIDUE_SCORES_PDB")
	cat $TMPLDIR/balls \
	| voronota query-balls \
	  --drop-altloc-indicators \
	  --set-adjuncts 'oc=1' \
	  --set-external-adjuncts $TMPLDIR/residue_quality_scores \
	  --set-external-adjuncts-name score \
	| voronota write-balls-to-atoms-file \
	  --pdb-output "$OUTFILE_RESIDUE_SCORES_PDB" \
	  --pdb-output-b-factor score \
	> /dev/null
fi

if [ -n "$OUTFILE_ATOM_DEPTHS" ]
then
	mkdir -p $(dirname $OUTFILE_ATOM_DEPTHS)
	cp $TMPLDIR/foundation/depth_values "$OUTFILE_ATOM_DEPTHS"
fi

if [ -n "$OUTFILE_CAMEO_RESIDUE_SCORES" ]
then
	mkdir -p $(dirname $OUTFILE_CAMEO_RESIDUE_SCORES)
	cat $TMPLDIR/balls \
	| voronota query-balls \
	  --set-external-adjuncts <(awk '{print $1 " " (10*(1-$2))}' $TMPLDIR/residue_quality_scores) \
	  --set-external-adjuncts-name rqscore \
	| voronota write-balls-to-atoms-file \
	  --pdb-output $OUTFILE_CAMEO_RESIDUE_SCORES \
	  --pdb-output-b-factor rqscore \
	  --pdb-output-template $INFILE_CONTENTS_EXTRACTED \
	> /dev/null
fi

if [ -n "$OUTFILE_CASPLIKE_RESIDUE_SCORES" ]
then
	mkdir -p $(dirname "$OUTFILE_CASPLIKE_RESIDUE_SCORES")
	cat $TMPLDIR/balls \
	| voronota query-balls \
	  --match-tags-not 'el=H' \
	| voronota query-balls \
	  --drop-altloc-indicators \
	  --set-adjuncts 'oc=1' \
	  --set-external-adjuncts <(awk 'BEGIN {a_mean=0.310; a_sd=0.097; b_mean=0.487; b_sd=0.203; d0=3.0;} {s=((($2-a_mean)/a_sd)*b_sd+b_mean); if(s<0.01){s=0.01;} if(s>1.0){s=1.0;} print $1 " " (d0*sqrt((1-s)/s)); }' $TMPLDIR/residue_quality_scores) \
	  --set-external-adjuncts-name score \
	| voronota write-balls-to-atoms-file \
	  --pdb-output "$OUTFILE_CASPLIKE_RESIDUE_SCORES" \
	  --pdb-output-b-factor score \
	  --add-chain-terminators \
	> /dev/null
fi

if [ -n "$OUTFILE_PLOT" ]
then

mkdir -p $(dirname $OUTFILE_PLOT)

cat $TMPLDIR/foundation/atom_energies \
| voronota score-contacts-quality \
  --default-mean -0.34 \
  --default-sd 0.19 \
  --means-and-sds-file $TMPLDIR/resources/energy_means_and_sds \
  --external-weights-file $TMPLDIR/foundation/depth_values \
  --smoothing-window "$(echo "0,$SMOOTHING_WINDOW" | sed 's|^0,\(.\+,.\+\).*|\1|')" \
  --residue-scores-file $TMPLDIR/residue_quality_scores_for_plot \
> /dev/null

cat $TMPLDIR/residue_quality_scores_for_plot \
| voronota expand-descriptors \
> $TMPLDIR/table_for_plot

BACK_FILE="none"
if $HIGHLIGHT_SELECTION_IN_PLOT && [ -s "$TMPLDIR/selection_by_queries" ]
then
	cat $TMPLDIR/selection_by_queries \
	| voronota expand-descriptors \
	| awk '{print $1 " " $2}' | sort | uniq \
	> $TMPLDIR/back_table_for_plot
	if [ -s "$TMPLDIR/back_table_for_plot" ]
	then
		BACK_FILE="$TMPLDIR/back_table_for_plot"
	fi
fi

R --slave --vanilla --args "$TMPLDIR/table_for_plot" "$BACK_FILE" "$OUTFILE_PLOT" "$PLOT_MODE" << 'EOF' &> /dev/null
	args=commandArgs(TRUE);
	infile=args[1];
	back_infile=args[2];
	outfile=args[3];
	plot_mode=args[4];
	
	t=read.table(infile, header=FALSE, stringsAsFactors=FALSE);
	chains=sort(union(t$V1, t$V1));
	
	back_t=c();
	if(back_infile!="none")
	{
		back_t=read.table(back_infile, header=FALSE, stringsAsFactors=FALSE);
	}
	
	if(plot_mode=="small")
	{
		png(filename=outfile, width=8*length(chains), height=2.5, units="in", res=80);
		par(mfcol=c(1, length(chains)));
		par(mar=c(2.7,3,1,1));
	} else
	{
		png(filename=outfile, width=10, height=5*length(chains), units="in", res=200);
		par(mfrow=c(length(chains), 1));
	}

	for(chain in chains)
	{
		title=ifelse(chain==".", "VoroMQA local scores for the unnamed chain", paste("VoroMQA local scores for chain ", chain, sep=""));
		
		st=t[which(t$V1==chain),];
		resnums=st$V2;
		v1=st$V8;
		v2=st$V9;
		N=length(v1);
		mask=rep(0, N);
		duplicate=identical(v1, v2);
		
		thick_color=rgb(0.0, 0.0, 0.0);
		thick_color_back=rgb(1.0, 0.0, 0.0);
		thin_color=ifelse(duplicate, thick_color, rgb(0.7, 0.7, 0.7));
		thin_color_back=ifelse(duplicate, thick_color_back, rgb(1.0, 0.7, 0.7));
		
		if(plot_mode=="small")
		{
			plot(c(min(resnums),max(resnums)), c(0, 1), type="n", main="", xlab="", ylab="");
		} else
		{
			plot(c(min(resnums),max(resnums)), c(0, 1), type="n", main=title, xlab="Residue number", ylab="Score");
		}
		
		points(c(min(resnums)-50,max(resnums)+50), c(0.5, 0.5), type="l");
		
		back_resnums=c();
		if(length(back_t)>0)
		{
			back_resnums=back_t$V2[which(back_t$V1==chain)];
		}
		
		resnums_col=rep(thin_color, N);
		resnums_col[which(is.element(resnums, back_resnums))]=thin_color_back;
		i=1;
		while(i<N)
		{
			if((resnums[i+1]-resnums[i])<2)
			{
				points(c(resnums[i], resnums[i+1]), c(v1[i], v1[i+1]), type="l", col=c(resnums_col[i], resnums_col[i+1]), lwd=1);
				mask[i]=1; mask[i+1]=1;
			}
			i=i+1;
		}
		sel=which(mask==0);
		if(length(sel)>0)
		{
			points(resnums[sel], v1[sel], col=resnums_col[sel], pch=16, cex=0.5);
		}
		
		if(duplicate==FALSE)
		{
			resnums_col=rep(thick_color, N);
			resnums_col[which(is.element(resnums, back_resnums))]=thick_color_back;
			i=1;
			while(i<N)
			{
				if((resnums[i+1]-resnums[i])<2)
				{
					points(c(resnums[i], resnums[i+1]), c(v2[i], v2[i+1]), type="l", col=c(resnums_col[i], resnums_col[i+1]), lwd=3);
				}
				i=i+1;
			}
			sel=which(mask==0);
			if(length(sel)>0)
			{
				points(resnums[sel], v2[sel], col=resnums_col[sel], pch=16, cex=0.7);
			}
		}
	}
	
	dev.off();
EOF
fi

if [ -n "$OUTFILE_MAP" ]
then
	mkdir -p $(dirname $OUTFILE_MAP)
	cat $TMPLDIR/foundation/all_contacts \
	| voronota query-contacts \
	  --drop-adjuncts \
	  --set-external-means $TMPLDIR/atom_quality_scores \
	| voronota query-contacts \
	  --no-solvent \
	  --match-adjuncts em=0:1 \
	| sed 's|em=||' \
	| awk '{print $1 " " $2 " " $3 " " $4 " " $5 " em=" ($3*$6)}' \
	| voronota query-contacts \
	  --inter-residue \
	| sed 's|em=||' \
	| awk '{print $1 " " $2 " " $3 " " $4 " " $5 " em=" ($6/$3)}' \
	| voronota x-plot-contacts \
	  --adjunct-gradient em \
	  --adjunct-gradient-blue 1 \
	  --adjunct-gradient-red 0 \
	  --svg-output $TMPLDIR/map.svg \
	> /dev/null
	mv $TMPLDIR/map.svg $OUTFILE_MAP
fi

if [ -n "$OUTFILE_SCORES_FOR_CASP" ] && [ -s "$TMPLDIR/sequence_alignment" ]
then
	mkdir -p $(dirname "$OUTFILE_SCORES_FOR_CASP")
	voronota x-write-qa-scores-in-casp-format \
	  --name "$(basename $INFILE)" \
	  --global-score "$(cat $TMPLDIR/global_quality_score)" \
	  --sequence-length "$(cat $TMPLDIR/sequence_alignment | head -1 | tr -d '-' | tr -d ' ' | tr -d '\n' | tr -d '\t' | wc -c)" \
	  --local-scores "$TMPLDIR/residue_quality_scores" \
	  --wrap_size 20 \
	  --completeness-threshold 0.85 \
	  --converted-local-scores "$TMPLDIR/residue_quality_scores_for_casp" \
	> "$OUTFILE_SCORES_FOR_CASP"
	
	if [ -n "$OUTFILE_SCORES_FOR_CASP_LOG" ]
	then
		mkdir -p $(dirname "$OUTFILE_SCORES_FOR_CASP_LOG")
		paste \
		  "$TMPLDIR/residue_quality_scores" \
		  "$TMPLDIR/residue_quality_scores_for_casp" \
		| column -t \
		> "$OUTFILE_SCORES_FOR_CASP_LOG"
	fi
	
	if [ -n "$OUTFILE_SCORES_FOR_CASP_PDB" ]
	then
		mkdir -p $(dirname "$OUTFILE_SCORES_FOR_CASP_PDB")
		cat $TMPLDIR/balls \
		| voronota query-balls \
		  --drop-altloc-indicators \
		  --set-adjuncts 'oc=1' \
		  --set-external-adjuncts "$TMPLDIR/residue_quality_scores_for_casp" \
		  --set-external-adjuncts-name score \
		| voronota write-balls-to-atoms-file \
		  --pdb-output "$OUTFILE_SCORES_FOR_CASP_PDB" \
		  --pdb-output-b-factor score \
		> /dev/null
	fi
fi

if [ -n "$OUTFILE_SCREENSHOT" ]
then
	cat $TMPLDIR/balls \
	| voronota query-balls \
	  --drop-altloc-indicators \
	  --set-adjuncts 'oc=1' \
	  --set-external-adjuncts "$TMPLDIR/residue_quality_scores" \
	  --set-external-adjuncts-name score \
	| voronota write-balls-to-atoms-file \
	  --pdb-output "$TMPLDIR/structure_for_screenshot.pdb" \
	  --pdb-output-b-factor score \
	> /dev/null
	
	mkdir -p $(dirname "$OUTFILE_SCREENSHOT")

pymol -p -c "$TMPLDIR/structure_for_screenshot.pdb" &> /dev/null << EOF
	set all_states, on
	hide all
	bg_color $SCREENSHOT_BACKGROUND
	show cartoon
	spectrum b, red_white_blue, minimum=0, maximum=1
	orient
	center
	zoom complete=1
	png $OUTFILE_SCREENSHOT, width=$SCREENSHOT_WIDTH, height=$SCREENSHOT_HEIGHT
	quit
EOF

fi

if $PRINT_ENERGY_OF_CONTACTS_SELECTION
then
	{
		if [ -s "$TMPLDIR/selection_by_contacts_query" ]
		then
			cat $TMPLDIR/foundation/all_contacts \
			| voronota query-contacts $CONTACTS_QUERY_PARAMETERS \
			| awk '{print $1 " " $2 " " $5 " " $3}' \
			| tr ';' '_' \
			| voronota score-contacts-energy \
			--potential-file $TMPLDIR/resources/energy_potential \
			| awk '{print $5 " " $2 " " $3 " " $4 " " ($4/($2+$3))}'
		else
			echo "0 0 0 0 0"
		fi
	} > $TMPLDIR/energy_descriptor_of_contacts_selection
fi

if $PRINT_CLASH_SCORE_OF_CONTACTS_SELECTION
then
	CLASH_DIST_HEAVY_HEAVY="2.5"
	CLASH_DIST_HEAVY_HYDROGEN="2.0"
	CLASH_DIST_HYDROGEN_HYDROGEN="1.5"
	HYDROGEN_NAMES="A<H,HA,HA2,HB,HB1,HB2,HD1,HD11,HD2,HD21,HE,HE1,HE2,HE21,HE3,HG,HG1,HG11,HG12,HG2,HG21,HH,HH11,HH2,HZ,HZ1,HZ2,HZ3>"
	
	{
		if [ -s "$TMPLDIR/selection_by_contacts_query" ]
		then
			cat $TMPLDIR/foundation/all_contacts \
			| voronota query-contacts $CONTACTS_QUERY_PARAMETERS \
			| voronota query-contacts --match-min-seq-sep 2 \
			> $TMPLDIR/clash_interest_contacts

			echo \
			"$(cat $TMPLDIR/clash_interest_contacts \
			| awk '{print $1 " " $2}' | tr ' ' '\n' | sort | uniq | wc -l)" \
			"$(cat $TMPLDIR/clash_interest_contacts \
			| voronota query-contacts \
				--match-first-not "$HYDROGEN_NAMES" \
				--match-second-not "$HYDROGEN_NAMES" \
				--match-max-dist "$CLASH_DIST_HEAVY_HEAVY" \
			| awk '{print $1 " " $2}' | tr ' ' '\n' | sort | uniq | wc -l)" \
			"$(cat $TMPLDIR/clash_interest_contacts \
			| voronota query-contacts \
				--match-first-not "$HYDROGEN_NAMES" \
				--match-second "$HYDROGEN_NAMES" \
				--match-max-dist "$CLASH_DIST_HEAVY_HYDROGEN" \
			| awk '{print $1 " " $2}' | tr ' ' '\n' | sort | uniq | wc -l)" \
			"$(cat $TMPLDIR/clash_interest_contacts \
			| voronota query-contacts \
				--match-first "$HYDROGEN_NAMES" \
				--match-second "$HYDROGEN_NAMES" \
				--match-max-dist "$CLASH_DIST_HYDROGEN_HYDROGEN" \
			| awk '{print $1 " " $2}' | tr ' ' '\n' | sort | uniq | wc -l)" \
			| awk '{ if($1>0) {print (($2+$3+$4)/$1)} else {print "0"} }'
		else
			echo "0"
		fi
	} > $TMPLDIR/clash_score_of_contacts_selection
fi

{
	if $PRINT_GLOBAL_SCORES_HEADER
	then
		{
			echo "input ${ENERGY_MODE}_score ${ENERGY_MODE}_residues ${ENERGY_MODE}_atoms"
			[ -s "$TMPLDIR/global_quality_score_of_selection" ] && echo "sel_${ENERGY_MODE}_score sel_${ENERGY_MODE}_atoms"
			[ -s "$TMPLDIR/energy_descriptor_of_contacts_selection" ] && echo "sel_${ENERGY_MODE}_contacts sel_${ENERGY_MODE}_area sel_${ENERGY_MODE}_area_unk sel_${ENERGY_MODE}_energy sel_${ENERGY_MODE}_energy_norm"
			[ -s "$TMPLDIR/clash_score_of_contacts_selection" ] && echo "sel_${ENERGY_MODE}_clash_score"
		} \
		| tr '\n' ' ' \
		| sed 's/\s$/\n/'
	fi
	
	{
		echo $INFILE
		cat $TMPLDIR/global_quality_score
		cat $TMPLDIR/balls | grep 'A<CA>' | wc -l
		cat $TMPLDIR/balls | wc -l
		[ -s "$TMPLDIR/global_quality_score_of_selection" ] && cat $TMPLDIR/global_quality_score_of_selection
		[ -s "$TMPLDIR/energy_descriptor_of_contacts_selection" ] && cat $TMPLDIR/energy_descriptor_of_contacts_selection
		[ -s "$TMPLDIR/clash_score_of_contacts_selection" ] && cat $TMPLDIR/clash_score_of_contacts_selection
	} \
	| tr '\n' ' ' \
	| sed 's/\s$/\n/'
} \
| \
{
	if $PRINT_GLOBAL_SCORES_VERTICALLY
	then
		cat > "$TMPLDIR/htable"
		if [ "$(cat "$TMPLDIR/htable" | wc -l)" == "2" ]
		then
			paste -d ' ' <(cat "$TMPLDIR/htable" | head -1 | tr ' ' '\n') <(cat "$TMPLDIR/htable" | tail -1 | tr ' ' '\n')
		else
			cat "$TMPLDIR/htable" | tr ' ' '\n'
		fi
	else
		cat
	fi
} \
| \
{
	if [ -n "$OUTDIR_GLOBAL_SCORES" ]
	then
		mkdir -p "$OUTDIR_GLOBAL_SCORES"
		cat | tee "$(mktemp "$OUTDIR_GLOBAL_SCORES/vgs.XXXXXXXXXX")"
	else
		cat
	fi
}


if [ -n "$SURFACE_CRAVING_ATOMS_OUTPUT" ] || [ -n "$SURFACE_CRAVING_ATOMS_OUTPUT_PDB" ] || [ -n "$SURFACE_CRAVING_RESIDUES_OUTPUT" ] || [ -n "$SURFACE_CRAVING_RESIDUES_OUTPUT_PDB" ] || [ -n "$SURFACE_CRAVING_ATOMS_ENERGIES_OUTPUT" ] || [ -n "$SURFACE_CRAVING_ATOMS_ENERGIES_PDB" ] || [ -n "$WEIGHTED_SURFACE_CRAVING_RESIDUES_OUTPUT" ] || [ -n "$WEIGHTED_SURFACE_CRAVING_RESIDUES_OUTPUT_PDB" ]
then
	{
		echo "c<solvent>"
		cat $TMPLDIR/foundation/depth_values \
		| egrep ' 1$' \
		| awk '{print $1}'
	} > $TMPLDIR/exposure_ids

	cat $TMPLDIR/foundation/raw_contacts \
	| voronota query-contacts \
	  --match-external-first $TMPLDIR/exposure_ids \
	  --match-external-second $TMPLDIR/exposure_ids \
	| awk '{if($2=="c<solvent>") {print $1 " " $2 " " $5 " " $3} else {print $1 " " $2 " " $5 " 0"}}' \
	| tr ';' '_' \
	| voronota score-contacts-energy \
	  --potential-file $TMPLDIR/resources/energy_potential \
	  --depth "$SURFACE_CRAVING_DEPTH_VALUE" \
	  --atom-scores-file $TMPLDIR/surface_atom_energies \
	> /dev/null
	
	cat $TMPLDIR/surface_atom_energies \
	| voronota score-contacts-quality \
	  --default-mean 0.13 \
	  --default-sd 0.13 \
	  --smoothing-window 0 \
	  --atom-scores-file $TMPLDIR/surface_atom_quality_scores \
	  --residue-scores-file $TMPLDIR/surface_residue_quality_scores \
	> /dev/null
	
	cat $TMPLDIR/surface_atom_quality_scores \
	| awk '{print $1 " " (1-$2)}' \
	> $TMPLDIR/surface_atom_craving_scores
	
	cat $TMPLDIR/surface_residue_quality_scores \
	| awk '{print $1 " " (1-$2)}' \
	> $TMPLDIR/surface_residue_craving_scores
	
	if [ -n "$SURFACE_CRAVING_ATOMS_OUTPUT" ]
	then
		mkdir -p $(dirname "$SURFACE_CRAVING_ATOMS_OUTPUT")
		cp $TMPLDIR/surface_atom_craving_scores "$SURFACE_CRAVING_ATOMS_OUTPUT"
	fi
	
	if [ -n "$SURFACE_CRAVING_RESIDUES_OUTPUT" ]
	then
		mkdir -p $(dirname "$SURFACE_CRAVING_RESIDUES_OUTPUT")
		cp $TMPLDIR/surface_residue_craving_scores "$SURFACE_CRAVING_RESIDUES_OUTPUT"
	fi
	
	if [ -n "$SURFACE_CRAVING_ATOMS_ENERGIES_OUTPUT" ]
	then
		mkdir -p $(dirname "$SURFACE_CRAVING_ATOMS_ENERGIES_OUTPUT")
		cp $TMPLDIR/surface_atom_energies "$SURFACE_CRAVING_ATOMS_ENERGIES_OUTPUT"
	fi
	
	if [ -n "$SURFACE_CRAVING_ATOMS_ENERGIES_PDB" ]
	then
		mkdir -p $(dirname "$SURFACE_CRAVING_ATOMS_ENERGIES_PDB")
		
		cat $TMPLDIR/balls \
		| voronota query-balls \
		  --drop-altloc-indicators \
		  --set-adjuncts 'oc=1' \
		  --set-external-adjuncts <(cat $TMPLDIR/surface_atom_energies | awk '{print $1 " " ($4/100.0)}') \
		  --set-external-adjuncts-name score \
		| voronota write-balls-to-atoms-file \
		  --pdb-output "$SURFACE_CRAVING_ATOMS_ENERGIES_PDB" \
		  --pdb-output-b-factor score \
		> /dev/null
	fi
	
	if [ -n "$SURFACE_CRAVING_ATOMS_OUTPUT_PDB" ]
	then
		mkdir -p $(dirname "$SURFACE_CRAVING_ATOMS_OUTPUT_PDB")
		
		cat $TMPLDIR/balls \
		| voronota query-balls \
		  --drop-altloc-indicators \
		  --set-adjuncts 'oc=1' \
		  --set-external-adjuncts $TMPLDIR/surface_atom_craving_scores \
		  --set-external-adjuncts-name score \
		| voronota write-balls-to-atoms-file \
		  --pdb-output "$SURFACE_CRAVING_ATOMS_OUTPUT_PDB" \
		  --pdb-output-b-factor score \
		> /dev/null
	fi
	
	if [ -n "$SURFACE_CRAVING_RESIDUES_OUTPUT_PDB" ]
	then
		mkdir -p $(dirname "$SURFACE_CRAVING_RESIDUES_OUTPUT_PDB")
		
		cat $TMPLDIR/balls \
		| voronota query-balls \
		  --drop-altloc-indicators \
		  --set-adjuncts 'oc=1' \
		  --set-external-adjuncts $TMPLDIR/surface_residue_craving_scores \
		  --set-external-adjuncts-name score \
		| voronota write-balls-to-atoms-file \
		  --pdb-output "$SURFACE_CRAVING_RESIDUES_OUTPUT_PDB" \
		  --pdb-output-b-factor score \
		> /dev/null
	fi
	
	if [ -n "$WEIGHTED_SURFACE_CRAVING_RESIDUES_OUTPUT" ] || [ -n "$WEIGHTED_SURFACE_CRAVING_RESIDUES_OUTPUT_PDB" ]
	then
		
R --slave --vanilla --args \
  "$TMPLDIR/residue_quality_scores" \
  "$TMPLDIR/surface_residue_craving_scores" \
  "$TMPLDIR/weighted_surface_residue_craving_scores" \
&> /dev/null \
<< 'EOF'
args=commandArgs(TRUE);
infile1=args[1];
infile2=args[2];
outfile=args[3];
t1=read.table(infile1, header=FALSE, stringsAsFactors=FALSE);
t2=read.table(infile2, header=FALSE, stringsAsFactors=FALSE);
t=merge(t1, t2, by="V1");
result=data.frame(residue_id=t$V1, score=(t$V2.x*t$V2.y));
write.table(result, outfile, quote=FALSE, row.names=FALSE, col.names=FALSE);
EOF
		
		if [ -n "$WEIGHTED_SURFACE_CRAVING_RESIDUES_OUTPUT" ]
		then
			mkdir -p $(dirname "$WEIGHTED_SURFACE_CRAVING_RESIDUES_OUTPUT")
			cp $TMPLDIR/weighted_surface_residue_craving_scores "$WEIGHTED_SURFACE_CRAVING_RESIDUES_OUTPUT"
		fi
		
		if [ -n "$WEIGHTED_SURFACE_CRAVING_RESIDUES_OUTPUT_PDB" ]
		then
			mkdir -p $(dirname "$WEIGHTED_SURFACE_CRAVING_RESIDUES_OUTPUT_PDB")
			
			cat $TMPLDIR/balls \
			| voronota query-balls \
			  --drop-altloc-indicators \
			  --set-adjuncts 'oc=1' \
			  --set-external-adjuncts $TMPLDIR/weighted_surface_residue_craving_scores \
			  --set-external-adjuncts-name score \
			| voronota write-balls-to-atoms-file \
			  --pdb-output "$WEIGHTED_SURFACE_CRAVING_RESIDUES_OUTPUT_PDB" \
			  --pdb-output-b-factor score \
			> /dev/null
		fi
	fi
fi