#!/usr/bin/perl -w


$PI = 3.1415926540;
$DEG2RAD = $PI / 180;
$RT = 0.6;

#check if the commonds are correct

if (! defined ($ARGV[0]) || $ARGV[0] =~ /-h|-help|-info/) {
    print "\tUsage: T-Analyst.pl <inputfile> [<outputfile>]\n\tUse -sample for example inputfile\n";
    exit(0);
}
if ($ARGV[0] =~ /-v|-version/) {
    print "\tT-Analyst,\n";
    print "\tVersion 1,\n";
    print "\twritten by Rizi Ai, 2009.\n";
    print "\tFor information contact rai001\@ucr.edu\n";
    exit(0);
}
if ($ARGV[0] =~ /-sample/) {
    print "#\tA sample input file to T-Analyst:\n#\n";
    &print_sample_input();
    exit(0);
}

# Open the input and output files

if (! open (INPUT, $ARGV[0])) {
    print "Enter name of input file to process: ";
    $filename = <STDIN>;
    open (INPUT, $filename) || die "Can't open file\n";
}
if ($ARGV[1]) {    
    $testout=$ARGV[1];
    open OUTPUT,'>', $testout;
} else {
    $testout="$ARGV[0].out";
    open OUTPUT,'>', $testout;
}


@date=localtime(time());
$date[5]+=1900;
$date[4]++;
print OUTPUT "T-Analyst\n\n";
printf OUTPUT "%.2i:%.2i, %.2i.%.2i.%4i\n",$date[2],$date[1],$date[3],$date[4],$date[5];
print OUTPUT "-" x 79, "\n";
print OUTPUT "Parameters:\n\n";

# report.out tells general information
open REPORT,'>', "report.out";
print REPORT "T-Analyst\n\n";
printf REPORT "%.2i:%.2i, %.2i.%.2i.%4i\n",$date[2],$date[1],$date[3],$date[4],$date[5];
print REPORT "-" x 79, "\n";
print REPORT "Parameters:\n\n";
print "\n******************************\n       T-Analyst\n******************************\n";
print "Checking input files ...\n\n";

##############################################################################################

# Read the input file and define default parameters

$use_ptraj=0;
$root_defined=0;
$ptraj_exec='ptraj';
$use_pseudobond=0;
$plp_c=" ";
$plp_calculation=0;
$std_calculation=0;
$std_sort=0;
$plp_sort=0;
$omega_on = 0;
$phi_on = 0;
$psi_on = 0;
$sidechain_on =0;
$cov=0;
$level=0;
$pseudo1=-1;
$pseudo2=-1;

while ($line = <INPUT>) {
    chomp ($line);
    if ($line =~/^\#/) {
	# This is just a comment line in the input
    } elsif ($line =~ /^USE_PTRAJ\s+(\d+)/) {
	$use_ptraj=$1;
	if ($use_ptraj==0) {
	    die("Ptraj must be used for now\n");
	} else {
	    print OUTPUT "Ptraj used to calculate BAT values\n";
	    print REPORT "Ptraj used to calculate BAT values\n";
	}
    } elsif ($line =~ /^RESIDUES\s+(.+)/) {
	print OUTPUT "Solute residues: $1\n";
	print REPORT "Solute residues: $1\n";
	@res0 = split (/,/, $1);
	$i=0;
	foreach $j (@res0) {
	    if ($j =~ /(\d+)-(\d+)/) {
		for ($k=$1 ; $k <= $2 ; $k++) {
		    $res[$i++] = $k;
		}
	    } elsif ($j =~ /\d+/) {
		$res[$i++] = $j;
	    } else {
		die("Residue selection $j not understood\n");
	    }
	}
    } elsif ($line =~ /^ROOT_ATOMS\s+(\d+)\s+(\d+)\s+(\d+)/) {
	$root1 = $1;
	$root2 = $2;
	$root3 = $3;
	$root_defined=1;
	print OUTPUT "Root atoms defined as: $root1, $root2, $root3\n";
	print REPORT "Root atoms defined as: $root1, $root2, $root3\n";
    } elsif ($line =~ /^PTRAJ_EXEC\s+(.+)/) {
	$ptraj_exec = $1;
	print OUTPUT "ptraj executable set to: $ptraj_exec\n";
	print REPORT "ptraj executable set to: $ptraj_exec\n";
    } elsif ($line =~ /^PRMTOP\s+(.+)/) {
	print OUTPUT "Using prmtop file $1\n";
	print REPORT "Using prmtop file $1\n";
	$prmtop = $1;
    } elsif ($line =~ /^CRD\s+(.+)/) {
	print OUTPUT "Using trajectory file $1\n";
	print REPORT "Using trajectory file $1\n";
	$trajectory = $1;
    } elsif ($line =~ /^RESIDUE_PRINT\s+(.+)/) {
	print OUTPUT "Residues printed out: $1\n";
	print REPORT "Residues printed out: $1\n";
	@res_print0 = split (/,/, $1);
	$i=0;
	foreach $j (@res_print0) {
	    if ($j =~ /(\d+)-(\d+)/) {
		for ($k=$1 ; $k <= $2 ; $k++) {
		  if ($k != $res[-1]){
		    $res_print[$i++] = $k;
		  }
		}
	    } elsif ($j =~ /\d+/) {
		if ($j != $res[-1]){
		$res_print[$i++] = $j;
		}
	    } else {
		die("Residue print $j not understood\n");
	    }
	}
    } elsif ($line =~ /^TORSION_SELECT\s+(.+)/) {
	print OUTPUT "Torsion angle selected: $1\n";
	print REPORT "Torsion angle selected: $1 ((1) Omega ,(2) Phi,(3) Psi, (4) Side-chain)\n";
	@angle_select0 = split (/,/, $1);
	$i=0;
	foreach $j (@angle_select0) {
	    if ($j =~ /(\d+)-(\d+)/) {
		for ($k=$1 ; $k <= $2 ; $k++) {
		    $angle_select[$i++] = $k;
		}
	    } elsif ($j =~ /\d+/) {
		$angle_select[$i++] = $j;
	    } else {
		die("Angle selected $j not understood\n");
	    }
	}
	foreach $i(@angle_select){
		if($i==1){
			$omega_on = 1;
		} elsif ($i==2){
			$phi_on = 1;
		} elsif ($i==3) {
			$psi_on = 1;
		} elsif ($i==4) {
			$sidechain_on =1;
		}
	}

    } elsif ($line =~ /^ANGLE_REPRESENTATION\s+(.+)/) {

	print OUTPUT "Angle representation: $1\n";
	print REPORT "Angle representation: $1 ((1) -3.14 ~ 3.14, (2) -180 ~ 180) \n";
	$angle_representation = $1;

    } elsif ($line =~ /^STD_CALCULATION\s+(.+)/) {
 	$std_calculation= $1;
	if($std_calculation ==1){
	print OUTPUT "Standard deviation will be calculated\n";
 	print REPORT "Standard deviation will be calculated \n";
 	}
	else{
	print OUTPUT "No standard deviation will be output\n";
	print REPORT "No standard deviation will be output \n";
	}

   } elsif ($line =~ /^ENTROPY_CALCULATION\s+(.+)/) {
 	$plp_calculation= $1;
	if($plp_calculation ==1){
	$plp_c="plp";
 	print OUTPUT "Configurational entropy will be calculated\n";
 	print REPORT "Configurational entropy will be calculated \n";
 	}
	else{
	print OUTPUT "No configurational entropy will be output\n";
	print REPORT "No configurational entropy will be output \n";
	}

   } elsif ($line =~ /^STD_SORT\s+(.+)/) {
	
	$std_sort = $1;
	if($std_sort == 1){
		print OUTPUT "Sorted standard deviation will be output\n";
		print REPORT "Sorted standard deviation will be output\n";
	} elsif ($std_sort == 2) {
		print OUTPUT "Sorted standard deviation will not be output\n";
		print REPORT "Sorted standard deviation will not be output\n";
	} elsif ($std_sort == 3) {
		print OUTPUT "Both sorted and unsorted standard deviation will be output\n";
		print REPORT "Both sorted and unsorted standard deviation will be output\n";
	}

   } elsif ($line =~ /^ENTROPY_SORT\s+(.+)/) {
	
	$plp_sort = $1;
	if($plp_sort == 1){
		print OUTPUT "Sorted configurational entropy will be output\n";
		print REPORT "Sorted configurational entropy will be output\n";
	} elsif ($plp_sort == 2) {
		print OUTPUT "Sorted configurational entropy will not be output\n";
		print REPORT "Sorted configurational entropy will not be output\n";
	}elsif ($plp_sort == 3) {
		print OUTPUT "Both sorted and unsorted configurational entropy will be output\n";
		print REPORT "Both sorted and unsorted configurational entropy will be output\n";
	}

    } elsif ($line =~ /^COV\s+(.+)/) {
	$cov=$1;
	if($cov == 1){
		print OUTPUT "Covariance matrix will be output\n";
		print REPORT "Covariance matrix will be output\n";

	} elsif ($cov == 2) {
		print OUTPUT "Covariance matrix will not be output\n";
		print REPORT "Covariance matrix will not be output\n";
	}
    } elsif ($line =~ /^PLOT\s+(.+)/) {
	$level=$1;
	if($level == 1){
		print OUTPUT "Heatmap for Correlation matrix will be output\n";
		print REPORT "Heatmap for Correlation matrix will be output\n";

	} elsif ($level == 2) {
		print OUTPUT "Heatmap for Correlation matrix will not be output\n";
		print REPORT "Heatmap for Correlation matrix will not be output\n";
	}

    } elsif ($line =~ /^PSEUDO_BOND\s+(.+)\s+(.+)/) {
	$use_pseudobond = 1;
	$pseudo1 = $1;
	$pseudo2 = $2;
	print OUTPUT "Using pseudobond between atom $1 and $2\n";
	print REPORT "Using pseudobond between atom $1 and $2\n";
    }
}

close INPUT;


# open topology and trajectory files

open (PRMTOP, $prmtop) || die ("Can't open prmtop file $prmtop\n");
open (TRAJ, $trajectory) || die "Can't open trajectory file $trajectory\n";

if (@res==0) {die("No residue selected\n");}
if (@res_print==0) {die("No residue print selected\n");}

if ($root_defined==0) {
	print OUTPUT "Root atoms will be selected automatically\n";
	print REPORT "Root atoms will be selected automatically\n";
}
print OUTPUT "-" x 79, "\n";
print REPORT "-" x 79, "\n";
print OUTPUT "Reading amber prmtop file\n\n";
$line = <PRMTOP>;

if ( ! $line =~ /%VERSION  VERSION_STAMP = .*  DATE = .*/) {
    die ("Error reading input prmtop file\n");
}

while ($line = <PRMTOP>) {
    if ($line =~/FLAG\sPOINTERS/) {
	$num_atoms = &read_pointers();
	last;
    }
}


# # ########################################################################################
# # # Check if .crd matches .top
# # open (TRAJ, $trajectory) || die "Can't open trajectory file $trajectory\n";
# # $crd_line = <TRAJ>;
# # $num_crd=0;
# # while ($crd_line = <TRAJ>) {
# # 	@test=();
# #     	@test = split (" ",$crd_line);
# # 	$num_crd += $#test+1;
# # }
# # $num_frames=$num_crd/$num_atoms;
# # # test whether the number of coordinate divided by number of atoms is an integer
# # if ($num_crd % $num_atoms != 0){
# # 	die("Topology file and coordinate file do not match!");
# # }
# # ############################################################################################

# Read in .top file
while ($line = <PRMTOP>) {
    if ($line =~/FLAG\sATOM_NAME/) {
	@atom_name = ();
	@atom_name = &read_atom_name();
    }
    if ($line =~/FLAG\sRESIDUE_LABEL/) { 
	@res_label=();
	@res_label=&read_res_label(); 
    }	
    if ($line =~/FLAG\sRESIDUE_POINTER/) {
  # @res_pointer is the array carrying all the residue pointers of the protop file
	@res_pointer=();
	@res_pointer=&read_res_pointers();
    }

    if ($line =~/FLAG\sBONDS_INC_HYDROGEN/) {
	&read_bonds();
    }
    if ($line =~/FLAG\sBONDS_WITHOUT_HYDROGEN/) {
	&read_bonds();
	last;
    }
}


#############################################################################################

# Define roots
&build_connectivity_hashes();

if ($root_defined==0) {
    &autoselect_root_atoms();
}

if($use_pseudobond == 1){
	if ($pseudo1 != $root3){
		die ("Please choose the 3rd root as the frist pseudo bond atom!\n");
	}
}
# Construct the BAT tree
# This defines the molecule in terms of its internal bond angle 
# and torsional degrees of freedom as outlined in Killian et al. 
# J. Chem. Phys. 127, 024107, 2007.

# Put the root atoms into the tree

$tree{$root1} = [-1,-1,-1,'n'];
$tree{$root2} = [$root1,-1,-1,'n'];
$tree{$root3} = [$root2,$root1,-1,'n'];

# Check if root1 is terminal
foreach $tmp ( @{ $connect{$root1} } ) {
    if ( ! defined($tree{$tmp})) {
	die ("Error building tree: Root atom 1 is connected to atom $tmp\n");
    }
}

# Deal with root 2 impropers if necessary
foreach $tmp ( @{ $connect{$root2} } ) {
    if ( ! defined($tree{$tmp})) {
	$tree{$tmp} = [$root2, $root1, $root3, 'i'];
	foreach $tmp2 ( @{ $connect{$tmp} } ) {
	    if ( ! defined($tree{$tmp2} ) ) {
		die ("Error building tree: Atom $tmp is connected to non-root atom $tmp2\nAll atoms connected to Root atom 2 must be terminal or root atoms\n");
	    }
	}
    }
}

# Check that root3 has connections
$i=0;
foreach $tmp ( @{ $connect{$root3} } ) {
    if ( ! defined($tree{$tmp})) {
	$i++;
    }
}
if ($i==0) {
    die ("Error building tree: Root atom 3 is a terminal atom\n");
}
@nodes = ();
push @nodes, $root3;

while ($use_pseudobond >= 0) {
    while ($branch_head = pop @nodes) {
	@all_twigs = ();
	@all_twigs = @{ $connect{$branch_head} };
	@twigs=();
	foreach $tmp ( @all_twigs ) {
	    if ( ! defined $tree{$tmp} ) {
		push @twigs, $tmp;
	    }
	}
	foreach $tmp ( @twigs ) {
	    $third=$tree{$branch_head}[0];
	    $fourth=$tree{$branch_head}[1];
	    $tree{$tmp} = [$branch_head, $third, $fourth, 'p'];
	}
	push @nodes, @twigs;

    }

    $use_pseudobond--;
    if ($use_pseudobond >=0) {
	push @nodes, $pseudo1+0;
	push @{ $connect{$pseudo1+0} }, $pseudo2+0;
	push @{ $connect{$pseudo2+0} }, $pseudo1+0;
    }
}

&print_BAT_tree();
############################################################################################

# Calculate BAT values
if ($use_ptraj == 1) {

    @torsionfilenames=();
    # Write an input file for ptraj and run it
    $numcommands=0;
    $times=1;
    open PTRAJIN,'>', "ptraj${times}.in";
    print PTRAJIN "trajin $trajectory\n";

	$i=0;	
	$file_count=0;

    foreach $tmp ( sort {$a <=> $b} keys %tree) {
	$numcommands++;

	if ($tree{$tmp}[0]>0 && $tree{$tmp}[1]>0 && $tree{$tmp}[2]>0) {
	    printf PTRAJIN "dihedral torsion%.4i \@%-4.i \@%-4.i \@%-4.i \@%-4.i out torsion_generated%.4i \n", $tmp, $tmp, $tree{$tmp}[0], $tree{$tmp}[1], $tree{$tmp}[2], $tmp;
 	    push @torsionfilenames, sprintf("torsion_generated%.4i", $tmp);
	}

	# Since ptraj can die on big input, run it separately for every 100 bonds
	if ($numcommands >99) {
	    $file_count++;
	    $numcommands = 0;
	    close PTRAJIN;
	    print OUTPUT "-" x 79, "\n";
	    print OUTPUT "Running ptraj on input file $times\n";
	    print "Running ptraj on input file $times\n";
	    print OUTPUT "PTRAJ output: \n";

	    $ptraj_output = `$ptraj_exec $prmtop ptraj${times}.in 2>&1`;

	    if ( ! defined $ptraj_output) { die ("Error: ptraj produced no output, set PTRAJ_EXEC in the input file or run ptraj by hand\n");}
	    print OUTPUT $ptraj_output;
	    print OUTPUT "-" x 79, "\n";

open FRAME, $torsionfilenames[0];
$num_frames=0;
while ($tmp = <FRAME>) {
$num_frames ++;
}
close FRAME;

	if($file_count == 1){
  	 $file_count = 0;
  	 system("rm -rf *in");
  	 if(@torsionfilenames){
  	   foreach $file(@torsionfilenames){
    	    open FILE, $file;
		$i++;
		$j=0;
		while($i==$root1 or $i==$root2 or $i==$root3){
		for ($j=0;$j<=$num_frames;$j++){
			$torsionlist[$i][$j]=0;
		}
		$i++;
		$j=0;
		}
		while ($tmp = <FILE>) {
			$j++;
			$tmp =~ /\d+.\d+\s+(-?\d+.\d+)/;
			if ($angle_representation == 1){
				$torsionlist[$i][$j]=$1*$DEG2RAD;
			} else{
			$torsionlist[$i][$j]=$1;
			}
		}
		close FILE;
	}
	  system("rm -rf torsion_generated*");
	  @torsionfilenames=();
	}
	}
	$times++;
	open PTRAJIN, '>', "ptraj${times}.in";
	print PTRAJIN "trajin $trajectory\n";
	}
    } 
    close PTRAJIN; 
    $ptraj_output = `$ptraj_exec $prmtop ptraj${times}.in 2>&1`; 
    print OUTPUT "-" x 79, "\n";
    print OUTPUT "PTRAJ output: \n";
    if ( ! defined $ptraj_output) { die ("Error: ptraj produced no output, set PTRAJ_EXEC in the input file or run ptraj by hand\n");}
    print OUTPUT $ptraj_output;
    print OUTPUT "-" x 79, "\n";
} 
close PRMTOP;
close TRAJ;
system("rm -rf *in");

open FRAME, $torsionfilenames[0];
$num_frames=0;
while ($tmp = <FRAME>) {
$num_frames ++;
}
close FRAME;

foreach $file(@torsionfilenames){
    open FILE, $file;
    $i++;
    $j=0;
    while($i==$root1 or $i==$root2 or $i==$root3){
	for ($j=0;$j<=$num_frames;$j++){
		$torsionlist[$i][$j]=0;
	}
	$i++;
	$j=0;
	}
    while ($tmp = <FILE>) {
	$j++;
	$tmp =~ /\d+.\d+\s+(-?\d+.\d+)/;
	if ($angle_representation == 1){
	$torsionlist[$i][$j]=$1*$DEG2RAD;
	}
	else{
	$torsionlist[$i][$j]=$1;
	}
    }
    close FILE;
}
system("rm -rf torsion_generated*");

open TORSIONLIST, '>', "torsions_complete.out";
print TORSIONLIST "#\t";
	for($index=1; $index<=$i; $index++){
		print TORSIONLIST "Torsion $index\t";
	}	
	print TORSIONLIST "\n";


for ($k=1 ; $k<=$j ; $k++) {
    printf TORSIONLIST "$k\t";
    for ($l=1 ; $l<=$i ; $l++) {
	printf TORSIONLIST "%2.8f\t", $torsionlist[$l][$k];
    }
    print TORSIONLIST "\n";
}
close TORSIONLIST;
print OUTPUT "Outputting list of torsion values\n";
print OUTPUT "-" x 79, "\n";

#############################################################################

@omega =();
@phi = ();
@side_chain = ();
@aa_name = ();
@aa_index = ();

&residue_dihedral(\@res_pointer, \@res_label, \@atom_name, \@res_print, \@omega, \@phi,\@side_chain,\@aa_name,\@aa_index,\$omega_on,\$phi_on,\$psi_on,\$sidechain_on);


# Output Omega angle
if ($omega_on == 1){
	open CA_ANGLE, '>', "omega_torsion.out";
        open CA_STD, '>', "omega_std_1.out";

	print CA_ANGLE "#";
	foreach $res(@new_res_print){	
		print CA_ANGLE "\tResidue $res";
	}
	print CA_ANGLE "\n";

	for ($k=1 ; $k<=$j ; $k++) {
	   printf CA_ANGLE "$k\t";
           foreach $l(@omega) {
			printf CA_ANGLE "%2.6f\t", $torsionlist[$l][$k];
			printf CA_STD "%2.6f\t", $torsionlist[$l][$k];
	   }
	   print CA_ANGLE "\n";
	   print CA_STD "\n";
	}
	close CA_ANGLE;
	close CA_STD;
	print OUTPUT "Outputting list of omega dihedral angles\n";
	print REPORT "Outputting list of omega dihedral angles in omega_torsion.out\n";
}


# phi angles
if ($phi_on == 1) {
        open PHI_ANGLE, '>', "phi_torsion.out";
	open PHI_STD, '>', "phi_std_1.out";
        print PHI_ANGLE "#";
	foreach $res(@new_res_print){	
		print PHI_ANGLE "\tResidue $res";
	}

	print PHI_ANGLE "\n";

	for ($k=1 ; $k<=$j ; $k++) {
	   printf PHI_ANGLE "$k\t";
	   foreach $l(@phi){
			printf PHI_ANGLE "%2.6f\t", $torsionlist[$l][$k];
			printf PHI_STD "%2.6f\t", $torsionlist[$l][$k];	
	   }
	   print PHI_ANGLE "\n";
	   print PHI_STD "\n";
	}
	close PHI_ANGLE;
	close PHI_STD;
	print OUTPUT "Outputting list of phi dihedral angles\n";
	print REPORT "Outputting list of phi dihedral angles in phi_torsion.out\n";
}

if ($psi_on == 1){
	open PSI_ANGLE, '>', "psi_torsion.out";
	open PSI_STD, '>', "psi_std_1.out";
	print PSI_ANGLE "#";
	foreach $res(@new_res_print){	
		print PSI_ANGLE "\tResidue $res";
	}
	print PSI_ANGLE "\n";

	for ($k=1 ; $k<=$j ; $k++) {
	   printf PSI_ANGLE "$k\t";

	   foreach $l(@psi) {
		printf PSI_ANGLE "%2.6f\t", $torsionlist[$l][$k];
		printf PSI_STD "%2.6f\t", $torsionlist[$l][$k];	

	}
	print PSI_ANGLE "\n";
	print PSI_STD "\n";
	}
	close PSI_ANGLE;
	close PSI_STD;
	print OUTPUT "Outputting list of psi dihedral angles\n";
	print REPORT "Outputting list of psi dihedral angles in psi_torsion.out\n";
} 

if ($sidechain_on == 1 and $#side_chain != -1){
 	open SIDE_CHAIN, '>', "sidechain_torsion.out";
        open SIDE_CHAIN_STD, '>', "sidechain_std_1.out";
	print SIDE_CHAIN "#";
	for($index=1; $index<=$#side_chain+1; $index++){
		printf SIDE_CHAIN "\tsc $index\t\t";
	}
 	print SIDE_CHAIN "\n";

	for ($k=1 ; $k<=$j ; $k++) {
	   printf SIDE_CHAIN "$k\t";
	   foreach $l(@side_chain) {
		printf SIDE_CHAIN "%2.6f\t", $torsionlist[$l][$k];
		printf SIDE_CHAIN_STD "%2.6f\t", $torsionlist[$l][$k];	
	   }
	print SIDE_CHAIN "\n";
	print SIDE_CHAIN_STD "\n";
	}
	close SIDE_CHAIN;
	close SIDE_CHAIN_STD;
	print OUTPUT "Outputting list of side chain dihedral angles\n";
	print REPORT "Outputting list of side chain dihedral angles in sidechain_torsion.out\n";
}


##########################################################################################
# Read in test.out file

open (TESTOUT, $testout) || die ("Can't open file $testout \n");
$line = <TESTOUT>;

if ( ! $line =~ /% T-Analyst*/) {
    die ("Error reading! \n");
}
close TESTOUT;

@omega_line = ();
@phi_line = ();
@psi_line = ();
@sidechain_line =();
$report_index=-1;
$sc_report_index=-1;

foreach $res(@res_print){
	$is_no_res=0;
 	foreach $no_res (@no_res_print){
	if ($res == $no_res){
		$is_no_res=1;
		$sc_report_index++;
		last;	
	} 
	}
if($is_no_res==0) {
  $report_index++;
  $sc_report_index++;

  open (TESTOUT, $testout) || die ("Can't open file $testout \n");
  $line = <TESTOUT>;
  while ($line = <TESTOUT>) {
	if ($line =~/Bond,\sAngle,\sTorsional Tree\sfor\sresidue\sselection:/) {
		last;
	}
  }
  while ($line = <TESTOUT>) {
	if ($line =~ /--------------------------------------/) {
	  last;
	} else {
	@tmp = split( " ", $line);
		if($psi_on ==1 and defined ($psi[$report_index]) and $tmp[0] == $psi[$report_index]){
			push (@psi_line, $line);

		} elsif ($omega_on == 1 and defined ($omega[$report_index])and $tmp[0] == $omega[$report_index]){
			push (@omega_line, $line);

		} elsif ($phi_on == 1 and defined ($phi[$report_index]) and $tmp[0] == $phi[$report_index]){
			push (@phi_line, $line);
		}

	if($sidechain_on==1){
		for ($i=$aa_index[$sc_report_index][0]; $i<=$aa_index[$sc_report_index][1];$i++){ 
			if( $i!=-1 and $tmp[0] == $side_chain[$i-1]){
				push (@sidechain_line, $line);
			}
		}		
	}
	}
}
}

close TESTOUT;
}

$report_index=-1;
$aa_report_index = -1;
$out_index = 0;
@sidechain_col = ();

foreach $res(@res_print){
	$report_index++;
	$out_index++;
	print REPORT "\nIndex $out_index:\tResidue $res\t$aa_name[$report_index]\n";
	
	$is_no_res=0;
	foreach $no_res (@no_res_print){
	if ($res == $no_res){
		$is_no_res=1;
		print REPORT "Not a normal residue!\n";
		last;
	} 
	}
	if($is_no_res ==0){
	$aa_report_index++;
	if($omega_on == 1 and defined ($omega_line[$aa_report_index])) {
	print REPORT "\@omega\t$omega_line[$aa_report_index]";
	} 
	if($phi_on == 1 and defined ($phi_line[$aa_report_index])){
	print REPORT "\@phi\t$phi_line[$aa_report_index]";
	}
	if($psi_on == 1 and defined ($psi_line[$aa_report_index])){
	print REPORT "\@psi\t$psi_line[$aa_report_index]";
	}
	if($sidechain_on==1){
	print REPORT "\@side chain\n";
	
	for ($i=$aa_index[$report_index][0]; $i<=$aa_index[$report_index][1];$i++){
		if($i==-1){
			print REPORT "\tNo side chain dihedral angles ouput.\n";
		}else{
			print REPORT "\t(sc $i) $sidechain_line[$i-1]";
			push (@sidechain_col, $i);
		}
	}
	}
   
}
}
##########################################################################################
# draw histogram

$cmd = "R --slave <temp.R";
  if ($omega_on == 1){
    @omega_dat=();
    if ($angle_representation == 1){
      @omega_dat=&histogram("omega_torsion.out", -$PI, $PI, 721, "omega");
    } else {
      @omega_dat=&histogram("omega_torsion.out", -180, 180, 721, "omega");
    }
   &angle_correction(\@omega_dat,"omega_torsion.out", "omega",$plp_c,$plp_sort, 721);

	if($std_calculation==1){
        open(RCODE,">temp.R");
        print RCODE "x <-read.table(file=\"omega_std_1.out\")\n";
        print RCODE "y <-sd(x)\n";
        print RCODE "z <- format(y, digits=5, scientific=F)\n";
        print RCODE "write.table(z, file=\"torsion_std_pre.out\", quote=F, sep=\"\t\", row.names = F, col.names = F)\n";

        system("$cmd");
        system("rm -rf temp.R");
	close(RCODE);
	
	&std ("omega");
      }
print "\nOmega torsion angles ...\n";
} 
	

if ($phi_on == 1){
  @phi_dat=();
  if ($angle_representation == 1){
    @phi_dat=&histogram("phi_torsion.out", -$PI, $PI, 361, "phi");
  } else{
    @phi_dat=&histogram("phi_torsion.out", -180, 180, 361, "phi");
  }
  &angle_correction(\@phi_dat,"phi_torsion.out", "phi",$plp_c,$plp_sort, 361);


    if($std_calculation==1){
    open(RCODE,">temp.R");
    print RCODE "x <-read.table(file=\"phi_std_1.out\")\n";
      print RCODE "y <-sd(x)\n";
      print RCODE "z <- format(y, digits=5, scientific=F)\n";
      print RCODE "write.table(z, file=\"torsion_std_pre.out\", quote=F, sep=\"\t\", row.names = F, col.names = F)\n";
        system("$cmd");
        system("rm -rf temp.R");
	close(RCODE);

	&std ("phi");
    }
print "Phi torsion angles ...\n";
} 

if ($psi_on == 1){		
  @psi_dat=();
  if ($angle_representation == 1){
    @psi_dat=&histogram("psi_torsion.out", -$PI, $PI, 361, "psi");
  } else{
    @psi_dat=&histogram("psi_torsion.out", -180, 180, 361, "psi");
  }
  &angle_correction(\@psi_dat,"psi_torsion.out", "psi",$plp_c, $plp_sort,361);

    if($std_calculation==1){
    open(RCODE,">temp.R");
    print RCODE "x <-read.table(file=\"psi_std_1.out\")\n";
      print RCODE "y <-sd(x)\n";
      print RCODE "z <- format(y, digits=5, scientific=F)\n";
      print RCODE "write.table(z, file=\"torsion_std_pre.out\", quote=F, sep=\"\t\", row.names = F, col.names = F)\n";

        system("$cmd");
        system("rm -rf temp.R");
	close(RCODE);
	
	&std ("psi");
    }
print "Psi torsion angles ...\n";
} 
	
if ($sidechain_on == 1 and $#side_chain != -1){ 
  @sidechain_dat=();
  if ($angle_representation == 1){
    @sidechain_dat=&histogram("sidechain_torsion.out", -$PI, $PI, 73, "sidechain");
  } else{
    @sidechain_dat=&histogram("sidechain_torsion.out", -180, 180, 73, "sidechain");
  }
  &angle_correction(\@sidechain_dat,"sidechain_torsion.out", "sidechain",$plp_c,$plp_sort, 73);
	
    if($std_calculation==1){
    open(RCODE,">temp.R");
    print RCODE "x <-read.table(file=\"sidechain_std_1.out\")\n";
      print RCODE "y <-sd(x)\n";
      print RCODE "z <- format(y, digits=5, scientific=F)\n";
      print RCODE "write.table(z, file=\"torsion_std_pre.out\", quote=F, sep=\"\t\", row.names = F, col.names = F)\n";

        system("$cmd");
        system("rm -rf temp.R");
	close(RCODE);
	
	&std ("sidechain");
    }
print "Side-chain torsion angles ...\n\n";
}
	
print "Standard deviation calculation completed\n";
print "Entropy calculations completed\n";

####################################################################################
#covariance matrix

if($cov == 1){

if ($omega_on == 1){
	@ca_cov=();
  open STD, "omega_std_1.out";
 $tt=0;
  while ($line = <STD>) {
    @tmp = split /\t+/, $line;
  $ss=0;
foreach $l(@tmp){
	$ca_cov[$tt][$ss]=$l;
	$ss++;
}
$tt++;
    }
close STD;
$remove="omega_std_1.out";
system("rm -rf $remove");
}

if ($phi_on == 1){
	@phi_cov=();
  open STD, "phi_std_1.out";
$tt=0;
  while ($line = <STD>) {
    @tmp = split /\t+/, $line;
	$ss=0;
   foreach $l(@tmp){
	$phi_cov[$tt][$ss]=$l;
	$ss++;
}
$tt++;
    }
close STD;
$remove="phi_std_1.out";
system("rm -rf $remove");
}

if ($psi_on == 1){
	@psi_cov=();
  open STD, "psi_std_1.out";
$tt=0;
  while ($line = <STD>) {
    @tmp = split /\t+/, $line;
	$ss=0;
   foreach $l(@tmp){
	$psi_cov[$tt][$ss]=$l;
	$ss++;
}
$tt++;
    }
close STD;
$remove="psi_std_1.out";
system("rm -rf $remove");
}

if ($sidechain_on == 1 and $#side_chain != -1){
	@sc_cov=();
  open STD, "sidechain_std_1.out";
$tt=0;
  while ($line = <STD>) {
    @tmp = split /\t+/, $line;
	$ss=0;
  foreach $l(@tmp){
	$sc_cov[$tt][$ss]=$l;
	$ss++;
}
$tt++;
    }
close STD;
$remove="sidechain_std_1.out";
system("rm -rf $remove");
}


@cov=();
@cov_index = ();
$cov_col=0;
$cov_res=0;
$cov_res_sc=0;
$cov_res_all=-1;

foreach $s(@res_print) {
	$skip=0;
	$cov_res_all++;
	foreach $no_res (@no_res_print){
		if ($s == $no_res){
			$skip=1;
		} 
	}
	if($skip==0){
	if ($omega_on == 1){
		#First line in @cov_index is Residue Index
		$cov_index[0][$cov_col]=$s;
		#second line in @cov_index is the Residue Type(1:omega,2:psi,3:psi,4sc)
		$cov_index[1][$cov_col]=1;
		 for ($k=0 ; $k<$num_frames ; $k++) {
			$cov[$k][$cov_col] = $ca_cov[$k][$cov_res];		
		}
	$cov_col++;
	}

	if ($phi_on == 1){
		#First line in @cov_index is Residue Index
		$cov_index[0][$cov_col]=$s;
		#second line in @cov_index is the Residue Type(1:omega,2:psi,3:psi,4sc)
		$cov_index[1][$cov_col]=2;
		 for ($k=0 ; $k<$num_frames ; $k++) {
			$cov[$k][$cov_col] = $phi_cov[$k][$cov_res];		
		}
	$cov_col++;
	}

	if ($psi_on == 1){
		#First line in @cov_index is Residue Index
		$cov_index[0][$cov_col]=$s;
		#second line in @cov_index is the Residue Type(1:omega,2:psi,3:psi,4sc)
		$cov_index[1][$cov_col]=3;
		 for ($k=0 ; $k<$num_frames ; $k++) {
			$cov[$k][$cov_col] = $psi_cov[$k][$cov_res];		
		}
	$cov_col++;
	}

	if ($sidechain_on == 1){
	for ($i=$aa_index[$cov_res_all][0]; $i<=$aa_index[$cov_res_all][1];$i++){
		if($i!=-1){
			#First line in @cov_index is Residue Index
			$cov_index[0][$cov_col]=$s;
			#second line in @cov_index is the Residue, 4:sc)
			$cov_index[1][$cov_col]=4;
			$cov_index[2][$cov_col]=$i;

		 for ($k=0; $k<$num_frames; $k++) {
			$cov[$k][$cov_col] = $sc_cov[$k][$cov_res_sc];	
		 }
			$cov_res_sc++;
			$cov_col++;
		}
	}
}

$cov_res++;
}
}

open COV, '>', "cov.out";
for($k=0; $k<$num_frames ; $k++){
	for($t=0;$t<=$#{$cov[0]};$t++){
	printf COV "%2.6f\t", $cov[$k][$t];
}
print COV "\n";
}
close COV;


open(RCODE,">temp.R");	
print RCODE "x <-read.table(file=\"cov.out\")\n";
print RCODE "c <- cov(x, use=\"pairwise.complete.obs\")\n";
print RCODE "d <- format(c, digits=3, scientific=F)\n";
print RCODE "write.table(d, file=\"cov_pre.out\", quote=F, sep=\"\t\", row.names = F, col.names = F)\n"; 
print RCODE "cor_1 <- cov2cor(c)\n";
print RCODE "cor_2 <- format(cor_1, digits=3, scientific=F)\n";
print RCODE "write.table(cor_2, file=\"cor_pre.out\", quote=F, sep=\"\t\", row.names = F, col.names = F)\n"; 



if ($angle_representation == 2){
		print RCODE "xx <- x*0.0174533\n";
		print RCODE "cc <- cov(xx, use=\"pairwise.complete.obs\")\n";
		print RCODE "cdiag <- diag(cc) \n";
		print RCODE "e <- det(cc)\n";
	} else {
		print RCODE "cdiag <- diag(c) \n";
		print RCODE "e <- det(c)\n";
	}
	print RCODE "f <- 0.3 * (ncol(x)*(1+log(6.2831852))+ log(e))\n";
	print RCODE "ediag <- 1\n";
	print RCODE " for (i in 1 : length(cdiag)){ediag <- ediag*cdiag[i]} \n";
	print RCODE "fdiag <- 0.3 * (ncol(x)*(1+log(6.2831852))+ log(ediag))\n";
	print RCODE "g <- format(f, digits=3, scientific=F)\n";
	print RCODE "gdiag <- format(fdiag, digits=3, scientific=F)\n";
	print RCODE "write.table(g, file=\"QH.out\", quote=F, row.names = \"***\n Quasi-Harmonic Approximation (QH)\n***\n\nConfigurational Entropy Computed by Entire Covariance Matrix\n\n \", col.names = F)\n";
	print RCODE "write.table(gdiag, file=\"QH.out\", quote=F, row.names = \"\n\n\n\nConfigurational Entropy Computed by Diagnoal Elements of Covariance Matrix\n\n\", col.names = F,append = T)\n";


        system("$cmd");
        system("rm -rf temp.R");
	close(RCODE);

open COV_OUT, '>', "cov_mat.out";
open COR_OUT, '>', "cor_mat.out";
open COR_OUT_PRE, '>', "correlation_map_reference.out";
print COR_OUT_PRE "***Specific residue reference of correlation map index***\n\n";
$cor_ii=0;

print COV_OUT "\t\t";
print COR_OUT "\t\t";
for($i=0;$i<=$#{$cov_index[0]};$i++){
	print COV_OUT "Res$cov_index[0][$i]/";
	print COR_OUT "Res$cov_index[0][$i]/";
	$cor_ii++;
	print COR_OUT_PRE "correlation map index $cor_ii, Residue $cov_index[0][$i], ";
	if($cov_index[1][$i]==1){
		print COV_OUT "Omega\t";
		print COR_OUT "Omega\t";
		print COR_OUT_PRE "Omega torsion angle\n";
	} elsif($cov_index[1][$i]==2){
		print COV_OUT "Phi\t";
		print COR_OUT "Phi\t";
		print COR_OUT_PRE "Phi torsion angle\n";
	} elsif($cov_index[1][$i]==3){
		print COV_OUT "Psi\t";
		print COR_OUT "Psi\t";
		print COR_OUT_PRE "Psi torsion angle\n";
	} elsif($cov_index[1][$i]==4){
		print COV_OUT "SC($cov_index[2][$i])\t";
		print COR_OUT "SC($cov_index[2][$i])\t";
		print COR_OUT_PRE "Side-chain torsion angle, side-chain index number:$cov_index[2][$i](refer to \"report.out\")\n";
	} 
}
print COV_OUT "\n";
print COR_OUT "\n";
print COR_OUT_PRE "\n";

open(COV_PRE,"cov_pre.out");
open(COR_PRE,"cor_pre.out");
$i=0;
while ($line = <COV_PRE>){
	print COV_OUT "Res$cov_index[0][$i]/";
	if($cov_index[1][$i]==1){
		print COV_OUT "Omega\t";
	} elsif($cov_index[1][$i]==2){
		print COV_OUT "Phi\t";
	} elsif($cov_index[1][$i]==3){
		print COV_OUT "Psi\t";
	} elsif($cov_index[1][$i]==4){
		print COV_OUT "SC($cov_index[2][$i])\t";
	} 
	print COV_OUT "$line";
	$i++;		
}
$i=0;
while ($line = <COR_PRE>){
	print COR_OUT "Res$cov_index[0][$i]/";
	if($cov_index[1][$i]==1){
		print COR_OUT "Omega\t";
	} elsif($cov_index[1][$i]==2){
		print COR_OUT "Phi\t";
	} elsif($cov_index[1][$i]==3){
		print COR_OUT "Psi\t";
	} elsif($cov_index[1][$i]==4){
		print COR_OUT "SC($cov_index[2][$i])\t";
	} 
	print COR_OUT "$line";
	$i++;		
}
system("rm -rf cov_pre.out");

close COV_PRE;
close COR_PRE;
print COV_OUT "\n\n";
print COR_OUT "\n\n";
close COV_OUT;
close COR_OUT;
close COR_OUT_PRE;
}
system("rm -rf cov.out");

#level plot
if($level == 1){
	print "\nPlotting correlation map...\n";
	open(RCODE,">temp.R");	
	print RCODE "x <-read.table(file=\"cor_pre.out\")\n";
	print RCODE "y <- as.matrix(x)\n";
	print RCODE "library(\"lattice\")\n";
	print RCODE "pdf(file=\"correlation_map.pdf\")\n";
	print RCODE "levelplot(y,region = TRUE, labels=F,col=\"gray40\", at=c(-1.0, -0.9,  -0.75, -0.6,  -0.45, -0.3, -0.15,  0.15, 0.3, 0.45, 0.6, 0.75, 0.9, 1.0),xlab=\"correlation map index\", ylab=\"correlation map index\", main=\"Level Plot: torsion angle correlation map\ngenerated from cor_mat.out\nIndex refer to correlation_map_reference.out\")\n";
	print RCODE "contourplot(y,region = TRUE, labels=F,col=\"gray40\", at=c(-1.0, -0.9,  -0.75, -0.6,  -0.45, -0.3, -0.15,  0.15, 0.3, 0.45, 0.6, 0.75, 0.9, 1.0),xlab=\"correlation map index\", ylab=\"correlation map index\", main=\"Contour Plot: torsion angle correlation map\ngenerated from cor_mat.out\n*** Same as level plot but add outlines for each highlighted area\nIndex refer to correlation_map_reference.out\")\n";
	print RCODE "dev.off()\n";

	system("$cmd");
	system("rm -rf temp.R");
	close(RCODE);
}
system("rm -rf cor_pre.out");
#########################################################################################
# Finalize
print "\n*********************************\nExecution of T-Analyst completed!\n*********************************\n\n";
print OUTPUT "Execution of T-Analyst finished\n";
print OUTPUT "-" x 79, "\n";

close OUTPUT;
close REPORT;
exit(0);

#############################################################################################

#############################################################################################
# Subroutines

sub read_pointers(){
	$line = <PRMTOP>;
	while ( $line = <PRMTOP>) {
		if ($line =~ /FLAG/) {
		last;
		} else {
			@tmp = split /\s+/, $line;
			last;
		}
	}
return $tmp[1];
}


sub read_atom_name () {

    $counter=0;
    @atom_name =();

    $line = <PRMTOP>;
    while ( $line = <PRMTOP>) {
	if ($line =~ /FLAG/) {
	    last;
	} else {
		$tmp[$counter]=$line;
		$counter++;
	} 
	}
for($i=0;$i<=$counter-1;$i++){
    @info = split(" ",$tmp[$i]);
#read in atom name, if the length of atom name is greater than 4, just split them every 4 characters
   for($k=0;$k<=$#info;$k++){
        $length = length($info[$k]);
	if($length<=4){
	push (@atom_name , $info[$k]);
        } elsif($length>4){
		$cycle=int($length/4);
		for($j=0;$j<$cycle;$j++){
			push (@atom_name , substr($info[$k],4*$j ,4));
		}

		if($length%4!=0){
			push (@atom_name , substr($info[$k],4*$cycle));		
		}
           }
        }
}	
return @atom_name;	
}


sub read_res_pointers {
    @res_pointer=();
    $line = <PRMTOP>;
    while ($line = <PRMTOP>) {
	if ($line =~ /FLAG/) {
	    last;
	} else {
	    @tmp = split /\s+/, $line;
	    shift @tmp;
	    $num = push (@res_pointer, @tmp);
	}
    }
    print OUTPUT "$num residue pointers read\n";
    print REPORT "$num residue pointers read\n";

    @atom=();
    foreach $resnum (@res) {
	$first=$res_pointer[$resnum-1];
	if($res_pointer[$resnum]){
 		$last=$res_pointer[$resnum];
        	for ($i=$first ; $i < $last ; $i++) {
		    $atom[$i]=1;
		}
	}else {
		$last = $num_atoms+1;
		for ($i=$first ; $i < $last ; $i++) {
		    $atom[$i]=1;
		}
		last;
	}	
    }

$last_atom = $last-1;
    print OUTPUT "Atom map generated\n";
    print OUTPUT "Atom print map generated\n";
    print REPORT "Atom map generated\n";
    print REPORT "Atom print map generated\n";

 return @res_pointer;
}


sub read_res_label(){
   @res_label=();
    $line = <PRMTOP>;
    while ($line = <PRMTOP>) {
	if ($line =~ /FLAG/) {
	    last;
	} else {
	    @tmp = split /\s+/, $line;
	    $num = push (@res_label, @tmp);
	}
    }
    print OUTPUT "$num residue label read\n";
    print REPORT "$num residue label read\n";
return @res_label;
}


sub read_bonds {
    $line = <PRMTOP>;
    while ( $line = <PRMTOP>) {
	if ($line =~ /FLAG/) {
	    last;
	} else {
	    @tmp = split /\s+/, $line;
	    shift @tmp;
	    $num = push (@bonds, @tmp);
	}
    }
}


sub build_connectivity_hashes() {
    while (defined($a1= shift(@bonds))) {
	$a2 = shift @bonds;
	shift @bonds;
	$a1 = ($a1/3)+1;
	$a2 = ($a2/3)+1;
	if (defined ($atom[$a1]) && defined ($atom[$a2])) {
	    push @{ $connect{$a1} }, $a2;
	    push @{ $connect{$a2} }, $a1;
	} elsif ( ! defined ($atom[$a1]) && ! defined ($atom[$a2]) ) {
	    # not used, probably solvent
	} else {
	    die ("Error: residue selection is not a complete molecule, there is a bond between atoms $a1 and $a2\n");
	}
    }
    print OUTPUT "Connectivity table build\n";
    print REPORT "Connectivity table build\n";
}



sub print_BAT_tree() {
    print OUTPUT "Bond, Angle, Torsional Tree for residue selection:\n";
    foreach $tmp ( sort { $a <=> $b } keys %tree ) {
	printf OUTPUT "%5i %5i %5i %5i   $tree{$tmp}[3]\n",$tmp, $tree{$tmp}[0],$tree{$tmp}[1],$tree{$tmp}[2];
	
}
}


sub autoselect_root_atoms() {
    # Try and use the hydrogen-dihedrals from the prmtop as root atoms
    while ($line = <PRMTOP>) {
	if ($line =~/FLAG\sDIHEDRALS_INC_HYDROGEN/) {
	    $line = <PRMTOP>;
	    while ( $line = <PRMTOP>) {
		if ($line =~ /FLAG/) {
		    last;
		} else {
		    @tmp = split /\s+/, $line;
		    shift @tmp;
		    push (@dihedrals, @tmp);
		}
	    }
	}
    }

    $fail=1;

    while (defined ($dihed0=(shift @dihedrals))) {
	$dihed0=($dihed0/3)+1;
	$dihed1=((shift @dihedrals)/3)+1;
	$dihed2=((shift @dihedrals)/3)+1;
	$dihed3=((shift @dihedrals)/3)+1;
	shift @dihedrals;

	if (defined $atom[$dihed0]) {
	    $fail=0;
	    # Check if dihed0 is terminal
	    foreach $tmp ( @{ $connect{$dihed0} } ) {
		if ($tmp != $dihed1) {
		    $fail = 1;
		}
	    }
	    # Check if dihed1 has only root and terminal connections
	    foreach $tmp ( @{ $connect{$dihed1} } ) {
		if ($tmp != $dihed0 && $tmp != $dihed2) {
		    foreach $tmp2 ( @{ $connect{$tmp} } ) {
			if ( $tmp2 != $dihed1 ) {
			    $fail = 1;
			}
		    }
		}
	    }

	    # Check that root3 has connections
	    $i=0;
	    foreach $tmp ( @{ $connect{$dihed3} } ) {
		if ($tmp != $dihed1) {$i++;}
	    }
	    if ($i==0) { $fail = 1; }
	}
	if ($fail == 0) { 
	    ($root1,$root2,$root3) = ($dihed0,$dihed1,$dihed2);
	    last; 
	}
    }
    if ($fail == 1) {
	die ("Error: Could not autoselect root atoms\n");
    } else {
	print OUTPUT "Atoms $root1, $root2, $root3 selected as Root atoms\n";
	print REPORT "Atoms $root1, $root2, $root3 selected as Root atoms\n";
    }
}



sub std (){
	
	$torsion = $_[0];
	@torsion_std = ();
	$i=0;
	if($torsion eq "sidechain"){
		foreach $col(@sidechain_col){	
			$torsion_std[$i][0] = $col;
			$i++;
		}	
	} else {
		foreach $res(@new_res_print){	
				$torsion_std[$i][0] = $res;
				$i++;
		}
	}
	open(TORSION,"torsion_std_pre.out");
	$j=0;
	while($line = <TORSION>){
		@tmp = split /\t+/, $line;
		$torsion_std[$j][1] = $tmp[0];
		$j++;
	}
	close(TORSION);
	system("rm -rf torsion_std_pre.out");
	$file_std_out =  $torsion."_std.out";
	open TORSION_STD, '>', $file_std_out;
	print TORSION_STD "***********\n$torsion: Standard Deviations\n***********\n\n";
	if ($torsion eq "sidechain"){
		for ($l=0;$l<=$#torsion_std;$l++){
			print TORSION_STD "sc $torsion_std[$l][0]\t$torsion_std[$l][1]";
		}
	} else{
		for ($l=0;$l<=$#torsion_std;$l++){
			print TORSION_STD "Residue $torsion_std[$l][0]\t$torsion_std[$l][1]";
		}
	}
	print TORSION_STD "\n\n";
	close(TORSION_STD);
	
        if($std_sort == 1 or $std_sort == 3){
	@torsion_std_sort = sort { $b->[1] <=> $a->[1] }@torsion_std;
	$file_std_sort=$torsion."_std_sorted.out";
	open TORSION_STD_SORT, '>', $file_std_sort;
	print TORSION_STD_SORT "***********\n$torsion: Sorted Standard Deviations\n***********\n\n";
	if ($torsion eq "sidechain"){
		for ($l=0;$l<=$#torsion_std;$l++){
			print TORSION_STD_SORT "sc $torsion_std_sort[$l][0]\t$torsion_std_sort[$l][1]";
		}
	} else{
		for ($l=0;$l<=$#torsion_std;$l++){
			print TORSION_STD_SORT "Residue $torsion_std_sort[$l][0]\t$torsion_std_sort[$l][1]";
		}
	}
	print TORSION_STD_SORT "\n\n";
	close(TORSION_STD_SORT);
       } 

}



sub histogram (){

$input    = $_[0];
$start    = $_[1];
$end      = $_[2];
$nbins    = $_[3];
$string   = $_[4];

open(TEST,$input);
@info = <TEST>;
@test = split(" ",$info[1]);
$columns = $#test-1;
close(TEST);

open(BATCH,">batch.com");
     $read_block  = "read block \"".$input."\""; 
     print(BATCH "$read_block\n");
     for($i=0;$i<=$columns;$i++){
         $get_histogram  = "histogram(s".$i.",mesh(".$start.",".$end.",".$nbins."), off,off)";
         print(BATCH "$get_histogram\n");
     }
     for($i=0;$i<=$columns;$i++){
         $displace = $columns+1 +$i;
         $move_set = "move g0.s".$displace." to g0.s".$i;
         print(BATCH "$move_set\n");
    }
     $save_all = "saveall \"".$string.".agr\"";
     $get_out  = "exit";
     print(BATCH "$save_all\n");
     print(BATCH "$get_out\n");
     close(BATCH);

# xmgrace command:
$xmgrace = "xmgrace -nxy ".$input." -nosafe -batch batch.com ";
system("$xmgrace");
system("rm -rf batch.com");

$agr_filename=$string . ".agr";
open (STRING, $agr_filename) || die ("Can't open file $agr_filename\n");

@dat = ();
$j=-2;
while ($line = <STRING>) {
	$i=0;
    if ($line =~/\@type\sxy/) {
	  $j+=2;
   	 while ( $line = <STRING>) {
		if ($line =~ /&/) {
		    last;
	} else {
 	  @tmp = split /\s+/, $line;
	  $dat[$i][$j] = $tmp[0];
	  $dat[$i][$j+1] = $tmp[1];
	  $i++;
	} 
	}
    }


}
close(STRING);

 open DAT,'>', "$string.dat";
	for ($s=0;$s<=$i-1;$s++){
		for ($l=0; $l<=$j+1;$l++){
			printf DAT "%15.5g\t", $dat[$s][$l];
		}
	if($s != $i-1){
		print DAT "\n";}
	}
  close(DAT);
  return @dat;
}


sub angle_correction(){
#input_2 file is the _tosion.out file
$input_2 = $_[1];
#the output file name, eg: omega, phi, psi
$output_name = $_[2];

@angle_arr=();
if($output_name eq "omega"){
  @angle_arr=@omega_dat;
  $h=0.5*$DEG2RAD;
} elsif ($output_name eq "phi"){
  @angle_arr=@phi_dat;
  $h=$DEG2RAD;
} elsif ($output_name eq "psi"){
  @angle_arr=@psi_dat;
  $h=$DEG2RAD;
} elsif ($output_name eq "sidechain"){
  @angle_arr=@sidechain_dat;
  $h=5*$DEG2RAD;
} 

$input_plp=$_[3];
$plp_sort = $_[4];
$num_bins = $_[5];

open(INPUT,$input_2);
@angle_modify = <INPUT>;
splice(@angle_modify,0,1);
close(INPUT);


#  @angle_corr is a two dimentional array that carries the INPUT (tosion.out) file
@angle_corr=();
for($i=0;$i<=$#angle_modify;$i++){
    @spl = split (" ",$angle_modify[$i]);
	if(@spl){
	    for($j=0;$j<=$#spl;$j++){
		$angle_corr[$i][$j] = $spl[$j];
    		}
	}
}


#  to shift angle if the probability of the first bin(actually the 2nd bin) is greater than 0.05
# hist_range==0 means the histogram is from -180 to 180
# hist_range==1 means the histogram is from 90 to 450
$hist_range=0;

# @plp_arr is the array that holds the PlnP of each residue
@plp_arr = ();
$n_angle_arr=$#angle_arr;
$n_angle_modify=$#angle_modify;

# $#{$angle_arr[0]} is the number of columns in the first row;
$n_column=$#{$angle_arr[0]};
$col_count=0;

for($j=1; $j<=$n_column+1; $j+=2) {
	$n_below_0=0;
	$n_above_0=0;
   for($i=0;$i<int($n_angle_arr/2);$i++){
	$n_below_0 += $angle_arr[$i][$j];
   }

   for($i=int($n_angle_arr/2);$i<=$n_angle_arr;$i++){
	$n_above_0 += $angle_arr[$i][$j];
   }

# n_total is the total number of frames
	$n_total=$n_below_0 + $n_above_0;

#######################################################################
# Entropy calculation
	$col_count++;

# $prob is the probability of each bin
	$prob=0;
# $plp is the PlnP of each residue
	$plp=0;
# $n_angle_arr = $#angle_arr;
	for($i=0; $i<=$n_angle_arr; $i++){
		$prob=$angle_arr[$i][$j]/($n_total*$h);
# log(0) makes no sense
		if($prob != 0){
			$plp += $prob * log($prob);
		}
	}
	$ST=-$plp*$h*$RT;
	$plp_arr[$col_count-1][0]=$col_count;
	$plp_arr[$col_count-1][1]=$ST;
#####################################################################################

	# Angle correction 
	$bin_123 = $angle_arr[0][$j] +$angle_arr[1][$j]+$angle_arr[2][$j];
	$bin_5 = $angle_arr[0][$j]+$angle_arr[1][$j]+$angle_arr[2][$j]+$angle_arr[3][$j]+$angle_arr[4][$j];

	$l1=$num_bins-1;
	$l2=$num_bins-2;
	$l3=$num_bins-3;
	$l4=$num_bins-4;
	$l5=$num_bins-5;

$bin_last_123 = $angle_arr[$l1][$j]+$angle_arr[$l2][$j]+$angle_arr[$l3][$j];
$bin_last_5 = $angle_arr[$l1][$j]+$angle_arr[$l2][$j]+$angle_arr[$l3][$j]+$angle_arr[$l4][$j]+$angle_arr[$l5][$j];

if($output_name eq "sidechain"){
  if($bin_123>0 and $bin_last_123>0){
	$hist_range=1;
	if($n_below_0 < $n_above_0){
		for ($bin=3;$bin<=$num_bins-3;$bin++){
				if(($angle_arr[$bin][$j]==0) and ($angle_arr[$bin+1][$j]==0) and ($angle_arr[$bin+2][$j]==0)) {
					$bin++;
					last;
				}
			}
		if($bin==$num_bins-2){
			for ($bin=1;$bin<=$num_bins-10;$bin++){
				if((($angle_arr[$bin][$j]+$angle_arr[$bin+1][$j]+$angle_arr[$bin+2][$j]) <($angle_arr[$bin+3][$j]+$angle_arr[$bin+4][$j]+$angle_arr[$bin+5][$j])) and  (($angle_arr[$bin+3][$j]+$angle_arr[$bin+4][$j]+$angle_arr[$bin+5][$j]) <($angle_arr[$bin+6][$j]+$angle_arr[$bin+7][$j]+$angle_arr[$bin+8][$j]))) {$bin+=4;
					last;
				}
			}
		}
			$shift_degree_1=-$PI+$bin*2*$PI/($num_bins-1);
			$shift_degree_2=-180+$bin*360/($num_bins-1);

		for($k=0;$k<=$n_angle_modify;$k++){
			if($angle_representation ==1 and $angle_corr[$k][$col_count]<=$shift_degree_1){
					$angle_corr[$k][$col_count]+=2*$PI;				
			}
			if($angle_representation ==2 and $angle_corr[$k][$col_count]<=$shift_degree_2){
				$angle_corr[$k][$col_count]+=360;
			}
		}
	}else{
		for ($bin_right=$l3;$bin_right>=3;$bin_right--){
				if(($angle_arr[$bin_right][$j]==0) and ($angle_arr[$bin_right-1][$j]==0) and ($angle_arr[$bin_right-2][$j]==0)){
					$bin_right--;
					last;
				}
			}
		if($bin_right == 2){
			for ($bin_right=$l1;$bin_right>=10;$bin_right--){
				if((($angle_arr[$bin_right][$j]+$angle_arr[$bin_right-1][$j]+$angle_arr[$bin_right-2][$j]) <($angle_arr[$bin_right-3][$j]+$angle_arr[$bin_right-4][$j]+$angle_arr[$bin_right-5][$j])) and  (($angle_arr[$bin_right-3][$j]+$angle_arr[$bin_right-4][$j]+$angle_arr[$bin_right-5][$j]) <($angle_arr[$bin_right-6][$j]+$angle_arr[$bin_right-7][$j]+$angle_arr[$bin_right-8][$j]))) {$bin_right-=4;
					last;
				}
			}
		}
		$shift_degree_1=($bin_right-int($num_bins/2))*2*$PI/($num_bins-1);
		$shift_degree_2=($bin_right-int($num_bins/2))*360/($num_bins-1);

		for($k=0;$k<=$n_angle_modify;$k++){		
			if($angle_representation ==1 and $angle_corr[$k][$col_count]>=$shift_degree_1){
  				$angle_corr[$k][$col_count]-=2*$PI;				
			}
			if($angle_representation ==2 and $angle_corr[$k][$col_count]>=$shift_degree_2){	
					$angle_corr[$k][$col_count]-=360;
			}
		}	
	}
  }
  }else{
    if($bin_5>0 and $bin_last_5>0){
	$hist_range=1;
	if($n_below_0 < $n_above_0){
		for ($bin=5;$bin<=$num_bins-10;$bin++){
				if(($angle_arr[$bin][$j]==0) and ($angle_arr[$bin+1][$j]==0) and ($angle_arr[$bin+2][$j]==0) and ($angle_arr[$bin+3][$j]==0) and ($angle_arr[$bin+4][$j]==0) and ($angle_arr[$bin+5][$j]==0) and ($angle_arr[$bin+6][$j]==0) and ($angle_arr[$bin+7][$j]==0) and ($angle_arr[$bin+8][$j]==0)and ($angle_arr[$bin+9][$j]==0)) {
					$bin++;
					last;
				}
			}
		if($bin==$num_bins-9){
			for ($bin=1;$bin<=$num_bins-10;$bin++){
				if((($angle_arr[$bin][$j]+$angle_arr[$bin+1][$j]+$angle_arr[$bin+2][$j]) <($angle_arr[$bin+3][$j]+$angle_arr[$bin+4][$j]+$angle_arr[$bin+5][$j])) and  (($angle_arr[$bin+3][$j]+$angle_arr[$bin+4][$j]+$angle_arr[$bin+5][$j]) <($angle_arr[$bin+6][$j]+$angle_arr[$bin+7][$j]+$angle_arr[$bin+8][$j]))) {$bin+=4;
					last;
				}
			}
		}
			$shift_degree_1=-$PI+$bin*2*$PI/($num_bins-1);
			$shift_degree_2=-180+$bin*360/($num_bins-1);
	
		for($k=0;$k<=$n_angle_modify;$k++){		
			if($angle_representation ==1 and $angle_corr[$k][$col_count]<=$shift_degree_1){
					$angle_corr[$k][$col_count]+=2*$PI;				
			}
			if($angle_representation ==2 and $angle_corr[$k][$col_count]<=$shift_degree_2){
				$angle_corr[$k][$col_count]+=360;
			}		
		}

	}else{
		for ($bin_right=$l5;$bin_right>=10;$bin_right--){
				if(($angle_arr[$bin_right][$j]==0) and ($angle_arr[$bin_right-1][$j]==0) and ($angle_arr[$bin_right-2][$j]==0) and ($angle_arr[$bin_right-3][$j]==0) and ($angle_arr[$bin_right-4][$j]==0) and ($angle_arr[$bin_right-5][$j]==0) and ($angle_arr[$bin_right-6][$j]==0) and ($angle_arr[$bin_right-7][$j]==0) and ($angle_arr[$bin_right-8][$j]==0)and ($angle_arr[$bin_right-9][$j]==0)) {
					$bin_right--;
					last;
				}
			}
		if($bin_right == 9){
			for ($bin_right=$l1;$bin_right>=10;$bin_right--){
				if((($angle_arr[$bin_right][$j]+$angle_arr[$bin_right-1][$j]+$angle_arr[$bin_right-2][$j]) <($angle_arr[$bin_right-3][$j]+$angle_arr[$bin_right-4][$j]+$angle_arr[$bin_right-5][$j])) and  (($angle_arr[$bin_right-3][$j]+$angle_arr[$bin_right-4][$j]+$angle_arr[$bin_right-5][$j]) <($angle_arr[$bin_right-6][$j]+$angle_arr[$bin_right-7][$j]+$angle_arr[$bin_right-8][$j]))) {$bin_right-=4;
					last;
				}
			}
		}
			$shift_degree_1=($bin_right-int($num_bins/2))*2*$PI/($num_bins-1);
			$shift_degree_2=($bin_right-int($num_bins/2))*360/($num_bins-1);

			for($k=0;$k<=$n_angle_modify;$k++){	
				if($angle_representation ==1 and $angle_corr[$k][$col_count]>=$shift_degree_1){
  					$angle_corr[$k][$col_count]-=2*$PI;				
				}
				if($angle_representation ==2 and $angle_corr[$k][$col_count]>=$shift_degree_2){	
						$angle_corr[$k][$col_count]-=360;
				}
			}
	}
   }
}
}

if($input_plp eq "plp"){
	if($plp_sort == 2 or $plp_sort == 3){
		$plp_output_name= $output_name."_entropy.out";
		open PLP,'>', $plp_output_name ;
		print PLP "****************\n$output_name: Configurational Entropy\n****************\n\n";
	if($output_name eq "sidechain"){
		for ($i=0;$i<=$#plp_arr;$i++){
			printf PLP "sc $plp_arr[$i][0]\t";
	 		printf PLP "%2.3f\t", $plp_arr[$i][1];
	 		print PLP "\n";
		}
	} else{
		for ($i=0;$i<=$#plp_arr;$i++){
			printf PLP "Residue $plp_arr[$i][0]\t";
	 		printf PLP "%2.3f\t", $plp_arr[$i][1];
	 		print PLP "\n";
		}
	}
		close PLP;
	}

	if($plp_sort == 1 or $plp_sort == 3){
		@plp_sort = sort { $b->[1] <=> $a->[1] }@plp_arr;
		$plp_output_sorted = $output_name."_entropy_sorted.out";
		open PLPSORT,'>',$plp_output_sorted ;
		print PLPSORT "****************\n$output_name: Sorted Configurational Entropy\n****************\n\n";
		if($output_name eq "sidechain"){
			for ($i=0;$i<=$#plp_arr;$i++){
				printf PLPSORT "sc $plp_sort[$i][0]\t";
				printf PLPSORT "%2.3f\t", $plp_sort[$i][1];
				print PLPSORT "\n";
			}
		} else {
			for ($i=0;$i<=$#plp_arr;$i++){
				printf PLPSORT "Residue $plp_sort[$i][0]\t";
				printf PLPSORT "%2.3f\t", $plp_sort[$i][1];
				print PLPSORT "\n";
			}
		}
	close PLPSORT;
	}
}


if ($hist_range==1){
	$remove=$output_name.".dat";
	system("rm -rf $remove");
	$remove=$output_name.".agr";
	system("rm -rf $remove");
	$remove=$output_name."_torsion.out";
	system("rm -rf $remove");
	$remove=$output_name."_std_1.out";
	system("rm -rf $remove");
	$new_angle_name = $output_name."_torsion.out";
	$new_std_name = $output_name."_std_1.out";

	open ANGLECORR, '>', $new_angle_name;
	open STD, '>', $new_std_name;
	
	print ANGLECORR "#\t";
	if($output_name eq "sidechain"){
 		for ($j=1 ; $j<=$#{$angle_corr[0]} ; $j++) {
			print ANGLECORR "sc $sidechain_col[$j-1]\t\t";
		}
	} else {

		for ($j=1 ; $j<=$#{$angle_corr[0]} ; $j++) {
			print ANGLECORR "Residue $new_res_print[$j-1]\t";
		}
	}
	print ANGLECORR "\n";	


	for ($i=0 ; $i<=$#angle_corr ; $i++) {
		$k=$i+1;
		printf ANGLECORR "$k\t";
	     for ($j=1 ; $j<=$#{$angle_corr[0]} ; $j++) {
		printf ANGLECORR "%2.6f\t", $angle_corr[$i][$j];
		printf STD "%2.6f\t", $angle_corr[$i][$j];
 	   }
  	  print ANGLECORR "\n";
	  print STD "\n";
	}
	close ANGLECORR;
	close STD;

	if($output_name eq "sidechain"){
		if($angle_representation ==1){
 			&histogram($new_angle_name, -3*$PI, 3*$PI, 217, $output_name);
		} else{
 			&histogram($new_angle_name, -540, 540, 217, $output_name);
		}

	}elsif($output_name eq "omega"){
		if($angle_representation ==1){
 			&histogram($new_angle_name, -3*$PI, 3*$PI, 2161, $output_name);
		} else{
 			&histogram($new_angle_name, -540, 540, 2161, $output_name);
		}

	}else{
		if($angle_representation ==1){
 			&histogram($new_angle_name, -3*$PI, 3*$PI, 1081, $output_name);
		} else{
 			&histogram($new_angle_name, -540, 540, 1081, $output_name);
		}
	}
}
}


sub residue_dihedral() {

$num=-1;
$index =1;
print REPORT "Side chain index:\n\n";
@new_res_print = @res_print;
@no_res_print = ();

foreach $res_sel(@res_print){
	if($res_sel == $res[-1]){
		next;
	}
	$num++;
	$res_lab = $res_label[$res_sel-1];
	$res_pointer_begin = $res_pointer[$res_sel-1];
	$res_pointer_end = $res_pointer[$res_sel]-1;

	if($res_lab eq "VAL"){
		push(@aa_name, "VAL");
		for($i=$res_pointer_begin-1; $i<$res_pointer_end-1;$i++){
			if($phi_on == 1 and $atom_name[$i] eq "C"){
				push (@phi, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "CG1"){
				push (@side_chain, $i+1);	
			}
		}

		if($res_sel<($res[-1]-1)){
		  for($i=$res_pointer[$res_sel]-1; $i<$res_pointer[$res_sel+1]-2;$i++){
			if($omega_on == 1 and $atom_name[$i] eq "CA"){
				push (@omega, $i+1);	
			}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
				push (@psi, $i+1);
			}
		  }
 		} elsif ($res_sel==($res[-1]-1)){
			for($i=$res_pointer[$res_sel]-1; $i<$last_atom;$i++){
				if($omega_on == 1 and $atom_name[$i] eq "CA"){
					push (@omega, $i+1);	
				}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
					push (@psi, $i+1);
				}
			}
		}


	if($sidechain_on == 1) {		
		$aa_index[$num][0] = $index;
		$aa_index[$num][1] = $index;
		$index+=1;
	}

	} elsif($res_lab eq "ALA"){
		push(@aa_name, "ALA");
		for($i=$res_pointer_begin-1; $i<=$res_pointer_end-1;$i++){
			if($phi_on == 1 and $atom_name[$i] eq "C"){
				push (@phi, $i+1);	
			} 
		}

		if($res_sel<($res[-1]-1)){
		  for($i=$res_pointer[$res_sel]-1; $i<$res_pointer[$res_sel+1]-2;$i++){
			if($omega_on == 1 and $atom_name[$i] eq "CA"){
				push (@omega, $i+1);	
			}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
				push (@psi, $i+1);
			}
		  }
 		} elsif ($res_sel==($res[-1]-1)){
			for($i=$res_pointer[$res_sel]-1; $i<$last_atom;$i++){
				if($omega_on == 1 and $atom_name[$i] eq "CA"){
					push (@omega, $i+1);	
				}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
					push (@psi, $i+1);
				}
			}
		}


	if($sidechain_on == 1) {	
		$aa_index[$num][0] = -1;
		$aa_index[$num][1] = -1;
	}
	
	} elsif($res_lab eq "PRO"){
		push(@aa_name, "PRO");
		for($i=$res_pointer_begin-1; $i<=$res_pointer_end-1;$i++){
			if($phi_on == 1 and $atom_name[$i] eq "C"){
				push (@phi, $i+1);	
			} 
		}

		if($res_sel<($res[-1]-1)){
		  for($i=$res_pointer[$res_sel]-1; $i<$res_pointer[$res_sel+1]-2;$i++){
			if($omega_on == 1 and $atom_name[$i] eq "CA"){
				push (@omega, $i+1);	
			}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
				push (@psi, $i+1);
			}
		  }
 		} elsif ($res_sel==($res[-1]-1)){
			for($i=$res_pointer[$res_sel]-1; $i<$last_atom;$i++){
				if($omega_on == 1 and $atom_name[$i] eq "CA"){
					push (@omega, $i+1);	
				}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
					push (@psi, $i+1);
				}
			}
		}

	if($sidechain_on == 1) {	
		$aa_index[$num][0] = -1;
		$aa_index[$num][1] = -1;
	}
	
	} elsif($res_lab eq "GLY"){
		push(@aa_name, "GLY");
		for($i=$res_pointer_begin-1; $i<=$res_pointer_end-1;$i++){
			if($phi_on == 1 and $atom_name[$i] eq "C"){
				push (@phi, $i+1);	
			} 
		}


		if($res_sel<($res[-1]-1)){
		  for($i=$res_pointer[$res_sel]-1; $i<$res_pointer[$res_sel+1]-2;$i++){
			if($omega_on == 1 and $atom_name[$i] eq "CA"){
				push (@omega, $i+1);	
			}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
				push (@psi, $i+1);
			}
		  }
 		} elsif ($res_sel==($res[-1]-1)){
			for($i=$res_pointer[$res_sel]-1; $i<$last_atom;$i++){
				if($omega_on == 1 and $atom_name[$i] eq "CA"){
					push (@omega, $i+1);	
				}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
					push (@psi, $i+1);
				}
			}
		}

	if($sidechain_on == 1) {	
		$aa_index[$num][0] = -1;
		$aa_index[$num][1] = -1;
	}
	
	} elsif($res_lab eq "LEU"){
		push(@aa_name, "LEU");
		for($i=$res_pointer_begin-1; $i<=$res_pointer_end-1;$i++){
			if($phi_on == 1 and $atom_name[$i] eq "C"){
				push (@phi, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "CG"){
				push (@side_chain, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "CD1"){
				push (@side_chain, $i+1);	
			}
		}

		if($res_sel<($res[-1]-1)){
		  for($i=$res_pointer[$res_sel]-1; $i<$res_pointer[$res_sel+1]-2;$i++){
			if($omega_on == 1 and $atom_name[$i] eq "CA"){
				push (@omega, $i+1);	
			}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
				push (@psi, $i+1);
			}
		  }
 		} elsif ($res_sel==($res[-1]-1)){
			for($i=$res_pointer[$res_sel]-1; $i<$last_atom;$i++){
				if($omega_on == 1 and $atom_name[$i] eq "CA"){
					push (@omega, $i+1);	
				}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
					push (@psi, $i+1);
				}
			}
		}

	if($sidechain_on == 1) {			
		$aa_index[$num][0] = $index;
		$aa_index[$num][1] = $index+1;

		$index+=2;
	}
	} elsif($res_lab eq "MET"){
		push(@aa_name, "MET");
		for($i=$res_pointer_begin-1; $i<=$res_pointer_end-1;$i++){
			if($phi_on == 1 and $atom_name[$i] eq "C"){
				push (@phi, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "CG"){
				push (@side_chain, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "SD"){
				push (@side_chain, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "CE"){
				push (@side_chain, $i+1);	
			}
		}

		if($res_sel<($res[-1]-1)){
		  for($i=$res_pointer[$res_sel]-1; $i<$res_pointer[$res_sel+1]-2;$i++){
			if($omega_on == 1 and $atom_name[$i] eq "CA"){
				push (@omega, $i+1);	
			}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
				push (@psi, $i+1);
			}
		  }
 		} elsif ($res_sel==($res[-1]-1)){
			for($i=$res_pointer[$res_sel]-1; $i<$last_atom;$i++){
				if($omega_on == 1 and $atom_name[$i] eq "CA"){
					push (@omega, $i+1);	
				}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
					push (@psi, $i+1);
				}
			}
		}

	if($sidechain_on == 1) {	
		$aa_index[$num][0] = $index;
		$aa_index[$num][1] = $index+2;

		$index+=3;
	}
	} elsif($res_lab eq "ILE"){
		push(@aa_name, "ILE");
		for($i=$res_pointer_begin-1; $i<=$res_pointer_end-1;$i++){
			if($phi_on == 1 and $atom_name[$i] eq "C"){
				push (@phi, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "CG1"){
				push (@side_chain, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "CD1"){
				push (@side_chain, $i+1);	
			}
		}

		if($res_sel<($res[-1]-1)){
		  for($i=$res_pointer[$res_sel]-1; $i<$res_pointer[$res_sel+1]-2;$i++){
			if($omega_on == 1 and $atom_name[$i] eq "CA"){
				push (@omega, $i+1);	
			}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
				push (@psi, $i+1);
			}
		  }
 		} elsif ($res_sel==($res[-1]-1)){
			for($i=$res_pointer[$res_sel]-1; $i<$last_atom;$i++){
				if($omega_on == 1 and $atom_name[$i] eq "CA"){
					push (@omega, $i+1);	
				}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
					push (@psi, $i+1);
				}
			}
		}

	if($sidechain_on == 1) {	
		$aa_index[$num][0] = $index;
		$aa_index[$num][1] = $index+1;

		$index+=2;
	}
	} elsif($res_lab eq "SER"){
		push(@aa_name, "SER");
		for($i=$res_pointer_begin-1; $i<=$res_pointer_end-1;$i++){
			if($phi_on == 1 and $atom_name[$i] eq "C"){
				push (@phi, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "OG"){
				push (@side_chain, $i+1);	
			} 
		}

		if($res_sel<($res[-1]-1)){
		  for($i=$res_pointer[$res_sel]-1; $i<$res_pointer[$res_sel+1]-2;$i++){
			if($omega_on == 1 and $atom_name[$i] eq "CA"){
				push (@omega, $i+1);	
			}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
				push (@psi, $i+1);
			}
		  }
 		} elsif ($res_sel==($res[-1]-1)){
			for($i=$res_pointer[$res_sel]-1; $i<$last_atom;$i++){
				if($omega_on == 1 and $atom_name[$i] eq "CA"){
					push (@omega, $i+1);	
				}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
					push (@psi, $i+1);
				}
			}
		}

if($sidechain_on == 1) {	
		$aa_index[$num][0] = $index;
		$aa_index[$num][1] = $index;

		$index+=1;
}
	} elsif($res_lab eq "THR"){
		push(@aa_name, "THR");
		for($i=$res_pointer_begin-1; $i<=$res_pointer_end-1;$i++){
			if($phi_on == 1 and $atom_name[$i] eq "C"){
				push (@phi, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "OG1"){
				push (@side_chain, $i+1);	
			}
		}

		if($res_sel<($res[-1]-1)){
		  for($i=$res_pointer[$res_sel]-1; $i<$res_pointer[$res_sel+1]-2;$i++){
			if($omega_on == 1 and $atom_name[$i] eq "CA"){
				push (@omega, $i+1);	
			}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
				push (@psi, $i+1);
			}
		  }
 		} elsif ($res_sel==($res[-1]-1)){
			for($i=$res_pointer[$res_sel]-1; $i<$last_atom;$i++){
				if($omega_on == 1 and $atom_name[$i] eq "CA"){
					push (@omega, $i+1);	
				}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
					push (@psi, $i+1);
				}
			}
		}


if($sidechain_on == 1) {	

		$aa_index[$num][0] = $index;
		$aa_index[$num][1] = $index;

		$index+=1;
}

	} elsif($res_lab eq "CYS" or $res_lab eq "CYX"){
		push(@aa_name, "CYS");
		for($i=$res_pointer_begin-1; $i<=$res_pointer_end-1;$i++){
			if($phi_on == 1 and $atom_name[$i] eq "C"){
				push (@phi, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "SG"){
				push (@side_chain, $i+1);	
			}
		}

		if($res_sel<($res[-1]-1)){
		  for($i=$res_pointer[$res_sel]-1; $i<$res_pointer[$res_sel+1]-2;$i++){
			if($omega_on == 1 and $atom_name[$i] eq "CA"){
				push (@omega, $i+1);	
			}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
				push (@psi, $i+1);
			}
		  }
 		} elsif ($res_sel==($res[-1]-1)){
			for($i=$res_pointer[$res_sel]-1; $i<$last_atom;$i++){
				if($omega_on == 1 and $atom_name[$i] eq "CA"){
					push (@omega, $i+1);	
				}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
					push (@psi, $i+1);
				}
			}
		}


if($sidechain_on == 1) {	
		$aa_index[$num][0] = $index;
		$aa_index[$num][1] = $index;

		$index+=1;
}

	} elsif($res_lab eq "ASN"){
		push(@aa_name, "ASN");
		for($i=$res_pointer_begin-1; $i<=$res_pointer_end-1;$i++){
			if($phi_on == 1 and $atom_name[$i] eq "C"){
				push (@phi, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "CG"){
				push (@side_chain, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "OD1"){
				push (@side_chain, $i+1);	
			}
		}

		if($res_sel<($res[-1]-1)){
		  for($i=$res_pointer[$res_sel]-1; $i<$res_pointer[$res_sel+1]-2;$i++){
			if($omega_on == 1 and $atom_name[$i] eq "CA"){
				push (@omega, $i+1);	
			}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
				push (@psi, $i+1);
			}
		  }
 		} elsif ($res_sel==($res[-1]-1)){
			for($i=$res_pointer[$res_sel]-1; $i<$last_atom;$i++){
				if($omega_on == 1 and $atom_name[$i] eq "CA"){
					push (@omega, $i+1);	
				}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
					push (@psi, $i+1);
				}
			}
		}


if($sidechain_on == 1) {	
		$aa_index[$num][0] = $index;
		$aa_index[$num][1] = $index+1;

		$index+=2;
}

	}  elsif($res_lab eq "GLN"){
		push(@aa_name, "GLN");
		for($i=$res_pointer_begin-1; $i<=$res_pointer_end-1;$i++){
			if($phi_on == 1 and $atom_name[$i] eq "C"){
				push (@phi, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "CG"){
				push (@side_chain, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "CD"){
				push (@side_chain, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "OE1"){
				push (@side_chain, $i+1);	
			}
		}

		if($res_sel<($res[-1]-1)){
		  for($i=$res_pointer[$res_sel]-1; $i<$res_pointer[$res_sel+1]-2;$i++){
			if($omega_on == 1 and $atom_name[$i] eq "CA"){
				push (@omega, $i+1);	
			}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
				push (@psi, $i+1);
			}
		  }
 		} elsif ($res_sel==($res[-1]-1)){
			for($i=$res_pointer[$res_sel]-1; $i<$last_atom;$i++){
				if($omega_on == 1 and $atom_name[$i] eq "CA"){
					push (@omega, $i+1);	
				}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
					push (@psi, $i+1);
				}
			}
		}


if($sidechain_on == 1) {	
		$aa_index[$num][0] = $index;
		$aa_index[$num][1] = $index+2;

		$index+=3;
}

	} elsif($res_lab eq "PHE"){
		push(@aa_name, "PHE");
		for($i=$res_pointer_begin-1; $i<=$res_pointer_end-1;$i++){
			if($phi_on == 1 and $atom_name[$i] eq "C"){
				push (@phi, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "CG"){
				push (@side_chain, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "CD1"){
				push (@side_chain, $i+1);	
			}
		}

		if($res_sel<($res[-1]-1)){
		  for($i=$res_pointer[$res_sel]-1; $i<$res_pointer[$res_sel+1]-2;$i++){
			if($omega_on == 1 and $atom_name[$i] eq "CA"){
				push (@omega, $i+1);	
			}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
				push (@psi, $i+1);
			}
		  }
 		} elsif ($res_sel==($res[-1]-1)){
			for($i=$res_pointer[$res_sel]-1; $i<$last_atom;$i++){
				if($omega_on == 1 and $atom_name[$i] eq "CA"){
					push (@omega, $i+1);	
				}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
					push (@psi, $i+1);
				}
			}
		}

if($sidechain_on == 1) {	
		$aa_index[$num][0] = $index;
		$aa_index[$num][1] = $index+1;

		$index+=2;
}

	}   elsif($res_lab eq "TYR"){
		push(@aa_name, "TYR");
		for($i=$res_pointer_begin-1; $i<=$res_pointer_end-1;$i++){
			if($phi_on == 1 and $atom_name[$i] eq "C"){
				push (@phi, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "CG"){
				push (@side_chain, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "CD1"){
				push (@side_chain, $i+1);	
			}
		}

		if($res_sel<($res[-1]-1)){
		  for($i=$res_pointer[$res_sel]-1; $i<$res_pointer[$res_sel+1]-2;$i++){
			if($omega_on == 1 and $atom_name[$i] eq "CA"){
				push (@omega, $i+1);	
			}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
				push (@psi, $i+1);
			}
		  }
 		} elsif ($res_sel==($res[-1]-1)){
			for($i=$res_pointer[$res_sel]-1; $i<$last_atom;$i++){
				if($omega_on == 1 and $atom_name[$i] eq "CA"){
					push (@omega, $i+1);	
				}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
					push (@psi, $i+1);
				}
			}
		}

	if($sidechain_on == 1) {	
		$aa_index[$num][0] = $index;
		$aa_index[$num][1] = $index+1;

		$index+=2;
	}

	}  elsif($res_lab eq "LYS"){
		push(@aa_name, "LYS");
		for($i=$res_pointer_begin-1; $i<=$res_pointer_end-1;$i++){
			if($phi_on == 1 and $atom_name[$i] eq "C"){
				push (@phi, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "CG"){
				push (@side_chain, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "CD"){
				push (@side_chain, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "CE"){
				push (@side_chain, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "NZ"){
				push (@side_chain, $i+1);	
			}
		}

		if($res_sel<($res[-1]-1)){
		  for($i=$res_pointer[$res_sel]-1; $i<$res_pointer[$res_sel+1]-2;$i++){
			if($omega_on == 1 and $atom_name[$i] eq "CA"){
				push (@omega, $i+1);	
			}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
				push (@psi, $i+1);
			}
		  }
 		} elsif ($res_sel==($res[-1]-1)){
			for($i=$res_pointer[$res_sel]-1; $i<$last_atom;$i++){
				if($omega_on == 1 and $atom_name[$i] eq "CA"){
					push (@omega, $i+1);	
				}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
					push (@psi, $i+1);
				}
			}
		}

	if($sidechain_on == 1) {	
		$aa_index[$num][0] = $index;
		$aa_index[$num][1] = $index+3;

		$index+=4;
	}
	} elsif($res_lab eq "ARG"){
		push(@aa_name, "ARG");
		for($i=$res_pointer_begin-1; $i<=$res_pointer_end-1;$i++){
			if($phi_on == 1 and $atom_name[$i] eq "C"){
				push (@phi, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "CG"){
				push (@side_chain, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "CD"){
				push (@side_chain, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "NE"){
				push (@side_chain, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "CZ"){
				push (@side_chain, $i+1);	
			}
		}

		if($res_sel<($res[-1]-1)){
		  for($i=$res_pointer[$res_sel]-1; $i<$res_pointer[$res_sel+1]-2;$i++){
			if($omega_on == 1 and $atom_name[$i] eq "CA"){
				push (@omega, $i+1);	
			}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
				push (@psi, $i+1);
			}
		  }
 		} elsif ($res_sel==($res[-1]-1)){
			for($i=$res_pointer[$res_sel]-1; $i<$last_atom;$i++){
				if($omega_on == 1 and $atom_name[$i] eq "CA"){
					push (@omega, $i+1);	
				}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
					push (@psi, $i+1);
				}
			}
		}

if($sidechain_on == 1) {	
		$aa_index[$num][0] = $index;
		$aa_index[$num][1] = $index+3;

		$index+=4;
}
	} elsif($res_lab eq "ASP"){
		push(@aa_name, "ASP");
		for($i=$res_pointer_begin-1; $i<=$res_pointer_end-1;$i++){
			if($phi_on == 1 and $atom_name[$i] eq "C"){
				push (@phi, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "CG"){
				push (@side_chain, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "OD1"){
				push (@side_chain, $i+1);	
			}
		}

		if($res_sel<($res[-1]-1)){
		  for($i=$res_pointer[$res_sel]-1; $i<$res_pointer[$res_sel+1]-2;$i++){
			if($omega_on == 1 and $atom_name[$i] eq "CA"){
				push (@omega, $i+1);	
			}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
				push (@psi, $i+1);
			}
		  }
 		} elsif ($res_sel==($res[-1]-1)){
			for($i=$res_pointer[$res_sel]-1; $i<$last_atom;$i++){
				if($omega_on == 1 and $atom_name[$i] eq "CA"){
					push (@omega, $i+1);	
				}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
					push (@psi, $i+1);
				}
			}
		}

if($sidechain_on == 1) {	
		$aa_index[$num][0] = $index;
		$aa_index[$num][1] = $index+1;

		$index+=2;
}
	} elsif($res_lab eq "GLU"){ 
		push(@aa_name, "GLU");
		for($i=$res_pointer_begin-1; $i<=$res_pointer_end-1;$i++){
			if($phi_on == 1 and $atom_name[$i] eq "C"){
				push (@phi, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "CG"){
				push (@side_chain, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "CD"){
				push (@side_chain, $i+1);
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "OE1"){
				push (@side_chain, $i+1);
			}
		}

		if($res_sel<($res[-1]-1)){
		  for($i=$res_pointer[$res_sel]-1; $i<$res_pointer[$res_sel+1]-2;$i++){
			if($omega_on == 1 and $atom_name[$i] eq "CA"){
				push (@omega, $i+1);	
			}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
				push (@psi, $i+1);
			}
		  }
 		} elsif ($res_sel==($res[-1]-1)){
			for($i=$res_pointer[$res_sel]-1; $i<$last_atom;$i++){
				if($omega_on == 1 and $atom_name[$i] eq "CA"){
					push (@omega, $i+1);	
				}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
					push (@psi, $i+1);
				}
			}
		}

if($sidechain_on == 1) {	
		$aa_index[$num][0] = $index;
		$aa_index[$num][1] = $index+2;

		$index+=3;
}

	} elsif($res_lab eq "TRP"){
		push(@aa_name, "TRP");
		for($i=$res_pointer_begin-1; $i<=$res_pointer_end-1;$i++){
			if($phi_on == 1 and $atom_name[$i] eq "C"){
				push (@phi, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "CG"){
				push (@side_chain, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "CD1"){
				push (@side_chain, $i+1);	
			}
		}

		if($res_sel<($res[-1]-1)){
		  for($i=$res_pointer[$res_sel]-1; $i<$res_pointer[$res_sel+1]-2;$i++){
			if($omega_on == 1 and $atom_name[$i] eq "CA"){
				push (@omega, $i+1);	
			}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
				push (@psi, $i+1);
			}
		  }
 		} elsif ($res_sel==($res[-1]-1)){
			for($i=$res_pointer[$res_sel]-1; $i<$last_atom;$i++){
				if($omega_on == 1 and $atom_name[$i] eq "CA"){
					push (@omega, $i+1);	
				}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
					push (@psi, $i+1);
				}
			}
		}

	if($sidechain_on == 1) {	
		$aa_index[$num][0] = $index;
		$aa_index[$num][1] = $index+1;

		$index+=2;
	}

	}elsif($res_lab eq "HIS" or $res_lab eq "HID" or $res_lab eq "HIE" or $res_lab eq "HIP"){
		push(@aa_name, $res_lab);
		for($i=$res_pointer_begin-1; $i<=$res_pointer_end-1;$i++){
			if($phi_on == 1 and $atom_name[$i] eq "C"){
				push (@phi, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "CG"){
				push (@side_chain, $i+1);	
			} elsif($sidechain_on == 1 and $atom_name[$i] eq "CD2"){
				push (@side_chain, $i+1);	
			}
		}

		if($res_sel<($res[-1]-1)){
		  for($i=$res_pointer[$res_sel]-1; $i<$res_pointer[$res_sel+1]-2;$i++){
			if($omega_on == 1 and $atom_name[$i] eq "CA"){
				push (@omega, $i+1);	
			}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
				push (@psi, $i+1);
			}
		  }
 		} elsif ($res_sel==($res[-1]-1)){
			for($i=$res_pointer[$res_sel]-1; $i<$last_atom;$i++){
				if($omega_on == 1 and $atom_name[$i] eq "CA"){
					push (@omega, $i+1);	
				}elsif($psi_on == 1 and $atom_name[$i] eq "N"){
					push (@psi, $i+1);
				}
			}
		}

if($sidechain_on == 1) {	
		$aa_index[$num][0] = $index;
		$aa_index[$num][1] = $index+1;

		$index+=2;
}
	} else {
		push(@aa_name, $res_lab);
		@new_res_print = grep { $_ != $res_sel } @new_res_print;
		push (@no_res_print, $res_sel);
	if($sidechain_on == 1) {	
		$aa_index[$num][0] = -1;
		$aa_index[$num][1] = -1;
	}
	}       
   }
      print REPORT "\n\n";

}



sub print_sample_input() {
    print << 'EOF'



###################################################################
# Sample T-Analyst input file
###################################################################
# The amber prmtop and trajectory file
#
PRMTOP  HIV.prmtop
#
CRD     HIV.crd
#
###################################################################
# This determines if ptraj will be used to calculate BAT values
# PTRAJ_EXEC should be the path/name of your ptraj executable
# (default 'ptraj')
#
USE_PTRAJ       1

#
#PTRAJ_EXEC     /opt/amber-10/exe/ptraj 
#
###################################################################
# Select which residues make up the solute
# acceptable entries have the form of comma separated
# numbers or number ranges (e.g. 5,7-12,18)
#
RESIDUES        1-198
#
##################################################################
# Select which residue you want to print out
# acceptable entries have the form of comma separated
# numbers or number ranges (e.g. 5,7-12,18)
#
RESIDUE_PRINT     180-198
#
###################################################################
# Select which type(s) of torsion angle(s) you want to print out
# (1) Omega ,(2) Phi,(3) Psi, (4) Side_chain (e.g. 1, 2, 3 or 1-3)
#
TORSION_SELECT    2,3
#
###################################################################
# Select angle representation
# (1) -3.14 ~ 3.14 ,(2) -180 ~ 180 (e.g. 1)
#
ANGLE_REPRESENTATION	2
#
###################################################################
# Standard deviaton calculation
# (1) Yes, (2) No (e.g. 1)
#
STD_CALCULATION     1
#
###################################################################
# Configurational entropy calculation
# (1) Yes, (2) No (e.g. 1)
#
ENTROPY_CALCULATION     1
#
###################################################################
# Standard deviation sorting 
# (1) Yes, (2) No, (3) both files
# (e.g. 1)
#
STD_SORT     3
#
###################################################################
# Entropy sorting 
# (1) Yes, (2) No, (3) both files
# (e.g. 1)
#
ENTROPY_SORT    3
#
###################################################################
# Covariance matrix
# # (1) Yes, (2) No (e.g. 1)
COV	1
###################################################################
# Contour plot for correlation matrix
# please select (1) if you want to view correlation matrix as a contour plot.
#
# # (1) Yes, (2) No (e.g. 1)
PLOT		1
#
###################################################################
# The root atoms to construct the connectivity tree, put AUTO
# instead of three numbers to try automatic root atom selection
# If hand-selected, the root atoms must be connected, root1 must be 
# terminal and root2 must be connected only to root and terminal 
# atoms
#
ROOT_ATOMS       AUTO
#ROOT_ATOMS     43	40 	38
#
###################################################################
# If your solute contains two separate molecules, a pseudo-bond must
# be defined between them. The first atom given must be from the
# THIRD (3rd) root atom. The second pseudo-bond atom should be the 
# first atom of the second part of your protein
#
PSEUDO_BOND	38	1577
#
###################################################################



EOF
}
