#!/usr/bin/perl
use strict ;

#
# Ahoy mates
#

my $file_orig = "elec_starting.dat";

open( FILE_ORIG, "< $file_orig" || die "$! $file_orig" ) ;
my @file_orig = <FILE_ORIG> ;
close( FILE_ORIG ) ;

my ( $atom_pairs_ref ) = read_into_hashes_arrays ( @file_orig );

my @global_atom_pairs = @$atom_pairs_ref;         # [ atom_a, atom_b, distance ]
                                                  # atom_a = surrounding, atom_b = our residue
my $global_all_res = 0;
my $global_flipped_res = 0;
my @global_flipped_residues;
my @global_b_res_atom_num;

control_subroutine( @global_atom_pairs );
write_macro( \@global_flipped_residues, \@global_b_res_atom_num );


foreach my $elem ( @global_flipped_residues ){
    print "< $elem > \n";
}
print "total res: $global_all_res  / flipped res: $global_flipped_res \n";



# --------------------------------------------------------------------------------

sub write_macro
{
    my ( @flipped_residues ) = @{$_[0]} ;
    my ( @b_res_atom_num )   = @{$_[1]} ;

    my @elec_flip ;
    my $file_com = "hbond_donor_aceptor_flips.com" ;

    open ( FILE_COM, "> $file_com" || die "$! $file_com" ) ;

	push @elec_flip, "key active_save select active end\n";
	push @elec_flip, "key active select .n all end\n";
    my $nn = 0 ;
    for ( my $i=0; $i < @flipped_residues; $i++ )
    {
	$nn++ ;
	push @elec_flip, "key active select active .o by resid atom numb $b_res_atom_num[$i] end \n";
	#print "$i: $flipped_residues[$i]   \t $b_res_atom_num[$i] \n";
    }
	push @elec_flip, "< ?MAIN_CMDS_NEW:auto_stuff_side.cmds WORK_SEGM FLIP \n";
	push @elec_flip, "key active select active_save end \n";
	push @elec_flip, "key passive select by bond active end\n";
	push @elec_flip, "key active_save drop\n";

	push @elec_flip, "return \n";


    if ( $nn > 0 ) 
    { print FILE_COM @elec_flip ; }
    else
    { print FILE_COM "return\n" ; }

#    foreach ( @elec_flip ){
#	print FILE_COM $_ ;
#    }
    close ( FILE_COM ) ;
}

sub control_subroutine
{

    my @atom_pairs =  @_ ;
    my %b_pairs;
    my %my_resi_list ;
    my ( $resid, $count ) ; 



    for ( $count = 0; $count < @atom_pairs; $count++ )
    {
	my @a =  @{$atom_pairs[$count][0]} ;
	my @b =  @{$atom_pairs[$count][1]} ;

	$resid = join(" ", @b[2..4] ) ;
	print "< $resid > ", @a, " <->", @b, " \n";
	
	if ( !exists( $my_resi_list{$resid} ) ){
	    @{$my_resi_list{$resid}} = ($count) ;
	}
	else {
	    push( @{$my_resi_list{$resid}}, $count ) ;
	}
    }

    foreach $resid ( keys %my_resi_list )
    {
	$global_all_res++;
	
	#test 1B57 obakrat hoce flip narest!!
	#$resid = "HIS 1 129";

	#my @test = split(" ", $resid);
	#if ( $test[0] ne "HIS"){ next; }
	#test end


	#print " <", $resid, "> ", join(" - ", @{$my_resi_list{$resid}}), "\n" ;
	print "< $resid > ", @{$my_resi_list{$resid}}, " \n";

	my @residue = split(" ", $resid);
	my $residue_name = $residue[0];
	my $b_res_atom_num;


	# flip
	my $score_P1_THR = 0;
	my $score_P2_THR = 0;
	my $score_P1_ASN = 0;
	my $score_P2_ASN = 0;
	my $score_P1_GLN = 0;
	my $score_P2_GLN = 0;
	my $score_P1_HIS = 0;
	my $score_P2_HIS = 0;
	my $score_P3_HIS = 0;
	my $score_P4_HIS = 0;
	my $score_P5_HIS = 0;
	my $score_P6_HIS = 0;

	# THR
	my $score_0_OG1  = 0;
	my $score_AD_OG1 = 0;
	my $score_0_CG2  = 0;
	my $score_AD_CG2 = 0;

	# ASN
	my $score_A_OD1 = 0;
	my $score_D_OD1 = 0;
	my $score_A_ND2 = 0;
	my $score_D_ND2 = 0;

	# GLN
	my $score_A_OE1 = 0;
	my $score_D_OE1 = 0;
	my $score_A_NE2 = 0;
	my $score_D_NE2 = 0;

	# HIS
	my $score_0_CE1 = 0;
	my $score_A_CE1 = 0;
	my $score_D_CE1 = 0;
	my $score_0_NE2 = 0;
	my $score_A_NE2 = 0;
	my $score_D_NE2 = 0;
	my $score_0_CD2 = 0;
	my $score_A_CD2 = 0;
	my $score_D_CD2 = 0;
	my $score_0_ND1 = 0;
	my $score_A_ND1 = 0;
	my $score_D_ND1 = 0;

	my %bs = ();
	my %bs_list = ();
	my $nn = 0 ;
#	my $b_atom ; my $a_atom; my $b_res; 
#	my $distance ;
#	foreach $count (  @{$my_resi_list{$resid}} )
#	{
#	    my $a_atom = $atom_pairs[$count][0];
#	    my $b_atom = $atom_pairs[$count][1];
#	    $b_res = $atom_pairs[$count][1][2];
#	    $distance = $atom_pairs[$count][2];
#	    $b_res_atom_num = $atom_pairs[$count][1][0];    #last atom of res
#	    $bs{$b_atom}{$distance} = $count ;
#	    $bs_list{$b_atom} = "done" ;
#	    print \$a_atom, \$b_atom, " done\n" ;
#	}	    
#	my @b_atom_list = () ;
#	my $key ;
#	foreach $key ( keys %bs_list ) { push (@b_atom_list, $key ) ; }
#
#	my %bs1 ;
#
#	foreach $b_atom ( @b_atom_list )
#	{
#	    print $b_atom, "\n" ;
#	    foreach $distance ( keys %{$bs{$b_atom}} ) 
#	    {
#	    $count = $bs{$b_atom}{$distance} ;
#	    $bs1{$distance} = $count;
#	    }
#
#	    foreach $distance ( sort keys %bs1 ) 
#	    {
#		$count = $bs1{$distance} ;
#		print "$b_atom => $distance $count \n" ;
#	    }
#	}
#	exit ;
#
#	my @sorted = reverse sort { $b <=> $a } ( keys @{$bs{$key} ); 

	foreach $count ( @{$my_resi_list{$resid}} )
	{

	    my $a_atom = $atom_pairs[$count][0];
	    my $b_atom = $atom_pairs[$count][1];
	    my $b_res = $atom_pairs[$count][1][2];
	    my $distance = $atom_pairs[$count][2];
	    if ( $distance <= 0.1 ) { next ; }
	    $b_res_atom_num = $atom_pairs[$count][1][0];    #last atom of res
	    get_atom_properties ( $a_atom );   # adds A, D or 0
	    print join(" ", @{$a_atom}), " ++ ", join(" ", @{$b_atom}), "done\n" ;

	    my $a = @{$a_atom}[5];
	    my $b_atom_position = @{$b_atom}[6];  # 0AD

#	    print " property: <$a> <$b_res>\n";
# if no a_atom property found, skip pair
	    if ( $a eq "none"){ next; }


	# for THR
	    if ( $b_res eq "THR" ){
		
		print "    ", join(" ", @{$atom_pairs[$count][0]}, "\t:  $distance  : ", @{$atom_pairs[$count][1]} ), "\n";

		if ( ${$b_atom}[1] eq "OG1" ){

		    $score_0_OG1 = $score_0_OG1 + sum_ab( join( "_", $a, "0"), \$distance );
		    $score_AD_OG1 = $score_AD_OG1 + sum_ab( join( "_", $a, "AD"), \$distance );
		}
		if ( ${$b_atom}[1] eq "CG2" ){
		
		    $score_0_CG2 = $score_0_CG2 + sum_ab( join( "_", $a, "0"), \$distance );
		    $score_AD_CG2 = $score_AD_CG2 + sum_ab( join( "_", $a, "AD"), \$distance );
		}

		print "    ", "0: $score_0_OG1    \t$score_0_CG2 \n";
		print "    ", "AD: $score_AD_OG1    \t$score_AD_CG2\n\n";
	    }

	# for ASN
	    if ( $b_res eq "ASN" ){

		#print "    ", join(" ", @{$atom_pairs[$count][0]}, "\t:  $distance  : ", @{$atom_pairs[$count][1]} ), "\n";

		if ( @{$b_atom}[1] eq "OD1" )
		{
		    $score_A_OD1 = $score_A_OD1 + sum_ab( join( "_", $a, "A"), \$distance );
		    $score_D_OD1 = $score_D_OD1 + sum_ab( join( "_", $a, "D"), \$distance );
		    print " OD1 $score_A_OD1 $score_D_OD1 \n" ;
		}
		if ( @{$b_atom}[1] eq "ND2" ){
		
		    $score_A_ND2 = $score_A_ND2 + sum_ab( join( "_", $a, "A"), \$distance );
		    $score_D_ND2 = $score_D_ND2 + sum_ab( join( "_", $a, "D"), \$distance );
		    print " ND1 $score_A_ND2 $score_D_ND2 \n" ;
		}

		print "    ", "A: $score_A_OD1    \t$score_A_ND2 \n";
		print "    ", "D: $score_D_OD1    \t$score_D_ND2\n\n";
	    }


	#for GLN
	    if ( $b_res eq "GLN" ){

		#print "    ", join(" ", @{$atom_pairs[$count][0]}, "\t:  $distance  : ", @{$atom_pairs[$count][1]} ), "\n";

		if ( @{$b_atom}[1] eq "OE1" ){

		    $score_A_OE1 = $score_A_OE1 + sum_ab( join( "_", $a, "A"), \$distance );
		    $score_D_OE1 = $score_D_OE1 + sum_ab( join( "_", $a, "D"), \$distance );
		}
		if ( @{$b_atom}[1] eq "NE2" ){
		
		    $score_A_NE2 = $score_A_NE2 + sum_ab( join( "_", $a, "A"), \$distance );
		    $score_D_NE2 = $score_D_NE2 + sum_ab( join( "_", $a, "D"), \$distance );
		}

		print "    ", "A: $score_A_OE1    \t$score_A_NE2 \n";
		print "    ", "D: $score_D_OE1    \t$score_D_NE2\n\n";
	    }
	    
	#for HIS
	    if ( $b_res eq "HIS" ){

		print "    ", join(" ", @{$atom_pairs[$count][0]}, "\t:  $distance  : ", @{$atom_pairs[$count][1]} ), "\n";

		if ( @{$b_atom}[1] eq "CE1" ){

		    $score_0_CE1 = $score_0_CE1 + sum_ab( join( "_", $a, "0"), \$distance );
		    $score_A_CE1 = $score_A_CE1 + sum_ab( join( "_", $a, "A"), \$distance );
		    $score_D_CE1 = $score_D_CE1 + sum_ab( join( "_", $a, "D"), \$distance );
		}
		if ( @{$b_atom}[1] eq "NE2" ){
		
		    $score_0_NE2 = $score_0_NE2 + sum_ab( join( "_", $a, "0"), \$distance );
		    $score_A_NE2 = $score_A_NE2 + sum_ab( join( "_", $a, "A"), \$distance );
		    $score_D_NE2 = $score_D_NE2 + sum_ab( join( "_", $a, "D"), \$distance );
		}
		if ( @{$b_atom}[1] eq "CD2" ){
		
		    $score_0_CD2 = $score_0_CD2 + sum_ab( join( "_", $a, "0"), \$distance );
		    $score_A_CD2 = $score_A_CD2 + sum_ab( join( "_", $a, "A"), \$distance );
		    $score_D_CD2 = $score_D_CD2 + sum_ab( join( "_", $a, "D"), \$distance );
		}
		if ( @{$b_atom}[1] eq "ND1" ){
		
		    $score_0_ND1 = $score_0_ND1 + sum_ab( join( "_", $a, "0"), \$distance );
		    $score_A_ND1 = $score_A_ND1 + sum_ab( join( "_", $a, "A"), \$distance );
		    $score_D_ND1 = $score_D_ND1 + sum_ab( join( "_", $a, "D"), \$distance );
		}

#		print "    ", "0: $score_0_CE1    \t$score_0_NE2   \t$score_0_CD2   \t$score_0_ND1 \n";
#		print "    ", "A: $score_A_CE1    \t$score_A_NE2   \t$score_A_CD2   \t$score_A_ND1 \n";
#		print "    ", "D: $score_D_CE1    \t$score_D_NE2   \t$score_D_CD2   \t$score_D_ND1 \n";
	    }
	    
	}
	

	if ( $residue_name eq "THR" ){
	$score_P1_THR = $score_AD_OG1 + $score_0_CG2 ;
	$score_P2_THR = $score_AD_CG2 + $score_0_OG1 ;

	print " THR scores: $score_P1_THR  $score_P2_THR \n" ;
	if ( $score_P1_THR < $score_P2_THR ){ 
	    print "NAREDI FLIP!!!";  

	    $global_flipped_res++;
	    push @global_flipped_residues, $resid;
	    push @global_b_res_atom_num, $b_res_atom_num;
	}
	if ( $score_P1_THR == $score_P2_THR ){
	    #print "ENAKA!!!";
	}
	#print "    P1 = $score_P1_THR  \tP2 = $score_P2_THR\n\n";
	}

	
	if ( $residue_name eq "ASN" )
	{
	    $score_P1_ASN = $score_A_OD1 + $score_D_ND2 ;
	    $score_P2_ASN = $score_A_ND2 + $score_D_OD1 ;
	    print "$score_P1_ASN = $score_A_OD1 + $score_D_ND2 \n";
	    print "$score_P2_ASN = $score_A_ND2 + $score_D_OD1 \n";

	    if ( $score_P1_ASN < $score_P2_ASN )
	    {
		print "NAREDI FLIP!!!"; 
		
		$global_flipped_res++;
		push @global_flipped_residues, $resid;
		push @global_b_res_atom_num, $b_res_atom_num;
	    }
	    if ( $score_P1_ASN == $score_P2_ASN ){
		print "ENAKA!!!";
	    }
	    print "    P1 = $score_P1_ASN  \tP2 = $score_P2_ASN\n\n";
	}


	if ( $residue_name eq "GLN" ){
	$score_P1_GLN = $score_A_OE1 + $score_D_NE2 ;
	$score_P2_GLN = $score_A_NE2 + $score_D_OE1 ;

	if ( $score_P1_GLN < $score_P2_GLN ){
	    #print "NAREDI FLIP!!!";

	    $global_flipped_res++;
	    push @global_flipped_residues, $resid;
	    push @global_b_res_atom_num, $b_res_atom_num;
	}
	if ( $score_P1_GLN == $score_P2_GLN ){
	    #print "ENAKA!!!";
	}
	#print "    P1 = $score_P1_GLN  \tP2 = $score_P2_GLN\n\n";
	}
	

	if ( $residue_name eq "HIS" ){
	    $score_P1_HIS = $score_0_CE1 + $score_D_NE2 + $score_0_CD2 + $score_D_ND1 ;
	    $score_P2_HIS = $score_0_CE1 + $score_D_NE2 + $score_0_CD2 + $score_A_ND1 ;
	    $score_P3_HIS = $score_0_CE1 + $score_A_NE2 + $score_0_CD2 + $score_D_ND1 ;
	    
	    $score_P4_HIS = $score_0_NE2 + $score_D_CE1 + $score_0_ND1 + $score_D_CD2 ;
	    $score_P5_HIS = $score_0_NE2 + $score_D_CE1 + $score_0_ND1 + $score_A_CD2 ;
	    $score_P6_HIS = $score_0_NE2 + $score_A_CE1 + $score_0_ND1 + $score_D_CD2 ;
	    
	    my %not_sorted = ( "P1" => $score_P1_HIS,
			       "P2" => $score_P2_HIS,
			       "P3" => $score_P3_HIS,
			       "P4" => $score_P4_HIS,
			       "P5" => $score_P5_HIS,
			       "P6" => $score_P6_HIS,
		);
	    
	    my @sorted = sort { $not_sorted{$a} <=> $not_sorted{$b} } keys %not_sorted ;
    	    
	    
	    foreach my $elem ( @sorted ){
		print "    $elem:  $not_sorted{$elem}\n";
	    }
#	    print 
#		"P1 => $score_P1_HIS\n",
#		"P2 => $score_P2_HIS\n",
#		"P3 => $score_P3_HIS\n",
#		"P4 => $score_P4_HIS\n",
#		"P5 => $score_P5_HIS\n",
#		"P6 => $score_P6_HIS\n" ;
	    
	    if ( $sorted[5] ne "P1" and $sorted[5] ne "P2" and $sorted[5] ne "P3" ) {
		#print "NAREDI FLIP!!!"; 

		$global_flipped_res++;
		push @global_flipped_residues, $resid;
		push @global_b_res_atom_num, $b_res_atom_num;
	    }
	}
	#exit;
	#print "\n";
    }
}



sub sum_ab
{
    my ( $atom_pair ) = $_[0];
    my ( $distance )  = ${$_[1]};



    my $score = 0;
    my %scoring_matrix = ("0_0"   =>  0,
			  "0_A"   => -0.1,
			  "0_AD"  => -0.1,
			  "0_D"   => -0.1,
	                  "A_A"   => -1,
			  "AD_A" =>  2,
			  "AD_AD" =>  2,
			  "AD_D" =>  2,
			  "A_0"   => -0.1,
			  "A_0D"  => -0.1,
			  "A_AD"  =>  2,
			  "A_D"   =>  2,
			  "A_D0"  => -0.1,
			  "A_DA"  =>  2,
			  "A_DD"  =>  2,
			  "DD_0" =>  -4,
			  "DD_A" =>  4,
			  "DD_D" =>  -4,
			  "D_0"   => -0.1,
			  "D_A"   =>  2,
			  "D_AD"  =>  2,
			  "D_D"   => -1,
			  "D_DA"  =>  2,
       );

    $score = ( $scoring_matrix{$atom_pair} / $distance**2 );
    print "    sum_ab:   ", $atom_pair, "   ", $distance, "\t score:  $score\n";

   return ( $score );
}



sub read_into_hashes_arrays
{

    my @file = @_;
    my $line;
    my $count = 0;
    my @words;

    my @pairs;
    my %pairs_hash;


    shift @file;
    shift @file;
    shift @file;

    foreach $line ( @file )
    {
	print $line;
	chomp ( $line ) ;

	@words = split(" ", $line) ;

	# ommits any C atoms not belonging to THR or HIS
	if ( (substr($words[2], 0, 3) ne "THR" and substr($words[2], 0, 3) ne "HIS")
	     and ($words[1] =~ /[C]/ or $words[7] =~ /[C]/) ){ next; }	


	# substr reads only first 3 characters of residue name
	my @b_atom = ( $words[0], $words[1], substr($words[2], 0, 3), $words[3], $words[4]  );
	my @a_atom = ( $words[6], $words[7], substr($words[8], 0, 3), $words[9], $words[10] );
	my $distance = $words[11];

	#@pairs = ( [ atom_a, atom_b, distance ], ... );
	push @{$pairs[$count]}, [ @a_atom ];
	push @{$pairs[$count]}, [ @b_atom ];
	push @{$pairs[$count]}, $distance;
	$count++;

    }

    return (\@pairs);
}




sub get_atom_properties
{
    my ( $atom ) = @_;

    my $atom_name = @{$atom}[1];
    my $res_name = @{$atom}[2];

    #  main chain
    my @rules_main = (
	[ "H"  , "D" ],
	[ "N"  , "D" ],
	[ "CA" , "0" ],
	[ "C"  , "0" ],
	[ "O"  , "A" ],
	[ "CB" , "0" ],
#
	[ "OXT", "A" ],
        [ "OH2", "AD"],
	[ "CG" , "0" ],
	[ "CD" , "0" ],
	[ "CG1", "0" ],
	[ "CG2", "0" ],
	[ "CE1", "0" ],
	[ "CE2", "0" ],
	[ "CD1", "0" ],
	[ "CD2", "0" ],
	[ "CZ" , "0" ],
	[ "CZ2", "0" ],
	[ "CE" , "0" ],
	);

    #  side chain (res, atom name, atom properties)
    my @rules_side = (
	[ "HIS", "ND1", "AD" ],
	[ "HIS", "CE1", "0"  ],
	[ "HIS", "NE2", "AD" ],
	[ "HIS", "CD2", "0"  ],
	[ "ASN", "OD1", "A"  ],
	[ "ASN", "ND2", "D"  ],
	[ "GLN", "OE1", "A"  ],
	[ "GLN", "NE2", "D"  ],
	[ "THR", "OG1", "AD" ],
	[ "THR", "CG2", "0"  ],
	[ "ARG", "NE" , "D"  ],
	[ "ARG", "NH1", "D"  ],
	[ "ARG", "NH2", "D"  ],
	[ "ASP", "OD1", "A"  ],
	[ "ASP", "OD2", "A"  ],
	[ "GLU", "OE1", "A"  ],
	[ "GLU", "OE2", "A"  ],
	[ "LYS", "NZ" , "D"  ],
	[ "SER", "OG" , "AD" ],
	[ "TRP", "NE1", "D"  ],
	[ "TYR", "OH" , "AD" ],
	[ "ASN", "CG" , "0"  ],
	[ "ASP", "CG" , "0"  ],
	[ "GLN", "CG" , "0"  ],
	[ "LYS", "CE" , "0"  ],
	[ "LEU", "CD1", "0"  ],
	[ "H2O", "OH2", "AD" ],
	[ "HOH", "O", "AD" ],
#
	[ "CYS", "SG" , "AD" ],
	[ "MET", "SD" , "AD" ],
	[ "ILE", "CG1", "0"  ],
	[ "ILE", "CG2", "0"  ],
	[ "ZN", "ZN", "DD"  ],
	[ "NA", "BA", "DD"  ],
	[ "MG", "MG", "DD"  ],
	[ "CAL", "CA", "DD"  ],
	[ "CL", "CL", "AA"  ],
	[ "BR", "BR", "AA"  ],
	[ "SO4", "O1", "AA"  ],
	[ "SO4", "O2", "AA"  ],
	[ "SO4", "O3", "AA"  ],
	[ "SO4", "O4", "AA"  ],
	[ "PO4", "O1", "AA"  ],
	[ "PO4", "O2", "AA"  ],
	[ "PO4", "O3", "AA"  ],
	[ "PO4", "O4", "AA"  ],
	);
    

	my $find_atom_property = "no";

#  side chain
    	for ( my $i=0; $i<@rules_side; $i++){

	    
	    if ( $res_name eq $rules_side[$i][0] and $atom_name eq $rules_side[$i][1] ){

		push @{$atom}, $rules_side[$i][2] ;
		$find_atom_property = "yes";
#		print $i, "  ", $res_name, " ", $atom_name, " -- ", "@{$rules_side[$i]}", "\n";
		last;
	    }
	}


#  main chain atom
	if ( $find_atom_property eq "no" ){    
	    for ( my $i=0; $i<@rules_main; $i++){    
#		print $i, "  ", $res_name, " ", $atom_name, " -- ", "@{$rules_side[$i]}", "\n";
		if ( $atom_name eq $rules_main[$i][0] ){

		    push @{$atom}, $rules_main[$i][1];
		    $find_atom_property = "yes";
		    last;
		}
	    }
	}

	if ( $find_atom_property eq "no" ){
	    print "in get_atom_properties()\tfind_atom_property for @{$atom} not found!\n\n";
	    push @{$atom}, "none";
	}

    return ;
}
