#!C:\Perl\bin

$pdbfile=@ARGV[0];
open (PDB,$pdbfile);
@pdblines=<PDB>;
close (PDB);

($fixed, $moving, $interdomain, *atomcount, *atomtype, *resnum, *Axcoord, *Aycoord, *Azcoord, *Bxcoord, *Bycoord, *Bzcoord, $totalnumber) = &informationretrieval(\@pdblines);

(*beginningfixed, *endingfixed, $nsegF) = &domainprocessing($fixed);

(*beginningmoving, *endingmoving, $nsegM) = &domainprocessing($moving);

(*beginninginterdomain, *endinginterdomain) = &interdomainsort($interdomain); 

(*res_to_atom, *natres) = &restoatom(\@resnum, $totalnumber);

(*AMxcord, *AMycord, *AMzcord, *BMxcord, *BMycord, *BMzcord) = &coordinatesorter(\@Axcoord, \@Aycoord, \@Azcoord, \@Bxcoord, \@Bycoord, \@Bzcoord, \@res_to_atom, \@natres, $nsegM, \@beginningmoving, \@endingmoving); 

#(*AFxcord, *AFycord, *AFzcord, *BFxcord, *BFxcord, *BFxcord) = &coordinatesorter(\@Axcoord, \@Aycoord, \@Azcoord, \@Bxcoord, \@Bycoord, \@Bzcoord, \@res_to_atom, \@natres, $nsegF, \@beginningfixed, \@endingfixed); 

#($clashcounterA) = &calculation(\@AMxcord, \@AMycord, \@AMzcord, \@AFxcord, \@AFycord, \@AFzcord);

#($clashcounterB) = &calculation(\@BMxcord, \@BMxcord, \@BMxcord, \@BFxcord, \@BFxcord, \@BFxcord);

#print STDOUT $clashcounterA;

#print STDOUT $clashcounterB;

#####
#SUBROUTINES
#####

sub informationretrieval {

local ($pdblines) = @_;
my @atomcount;
my @atomtype; 
my @resnum;
my @chainID;
my @Axcoord;
my @Aycoord;
my @Azcoord;
my @Bxcoord;
my @Bycoord;
my @Bzcoord;

	#loops through pdb file
	foreach $pdblines(@pdblines) {
		
		#ignores co-ordinate information related to hindge
		if ($pdblines =~ /TER arrow molecule that represents hinge axis for movement from conformer 1 to 2/) {

		last;
		}

	#fixed domain information retrieveal   
	$fixed = $pdblines[7];

	#moving domain information retrieveal
	$moving = $pdblines[10];

	#moving domain information retrieveal
	$interdomainA = $pdblines[19];
	$interdomainB = $pdblines[21];

	# processing of interdomain information
	$interdomainA .= ',';
	$interdomainA =~ s/select//;
	$interdomainA =~ s/\s//; 
	$interdomainA =~ s/-/,/g;
	$interdomainA =~ s/A//g;

	# processing of interdomain information
	$interdomainB =~ s/select//;
	$interdomainB =~ s/\s//; 
	$interdomainB =~ s/-/,/g;
	$interdomainB =~ s/A//g;

	$interdomain = $interdomainA . $interdomainB;

		# identifies atom information and saves atom number, residue number, chainID, X,Y,Z coordinates and atom type in separate arrays
		if ($pdblines =~ /^ATOM/) {
	
		$iat++;
		$atomcount[$iat]=substr($pdblines,7,4);
		$atomtype[$iat]=substr($pdblines,13,3); 
		$resnum[$iat]=substr($pdblines,22,4);
		$chainID[$iat]=substr($pdblines,21,1);
	
			if ($chainID[$iat] eq 'A') {
			
			$Axcoord[$iat] = substr($pdblines,30,8);
			$Aycoord[$iat] = substr($pdblines,38,8);
			$Azcoord[$iat] = substr($pdblines,46,8);

			}

			elsif ($chainID[$iat] eq 'B') {

			$Bxcoord[$iat] = substr($pdblines,30,8);
			$Bycoord[$iat] = substr($pdblines,38,8);
			$Bzcoord[$iat] = substr($pdblines,46,8);

			}
		}
	}

return ($fixed, $moving, $interdomain, \@atomcount, \@atomtype, \@resnum, \@Axcoord, \@Aycoord, \@Azcoord, \@Bxcoord, \@Bycoord, \@Bzcoord, $iat);

}

sub domainprocessing {

local ($domain) = @_;
my @beginningdomain;
my @endingdomain;
@domainnumbers;
my $nseg;

	# processing of domain information
	$domain =~ s/select//;
	$domain =~ s/\s//; 
	$domain =~ s/-/,/g;
	$domain =~ s/A//g;

	# split domain information into an array
	@domainnumbers = split(/,/, $domain);

	$nseg = @domainnumbers/2;

		# divide domain information into beginning and ending regions of the domains 		
		foreach $domainnumbers (@domainnumbers) {
   			
			# if array entry is divisable by 2 (even numbered) then it is put in the beginning array
			if ($i % 2 == 0) {
			$ibeg++;
			$beginningdomain[$ibeg] = $domainnumbers[$i];

			}  
		
			# if array entry is not divisable by 2 (odd numbered) then it is put in the ending array
	   		elsif ($i % 2 != 0) {
			$iend++;
			$endingdomain[$iend] = $domainnumbers[$i];
   			}
		
			# print error if something goes wrong 
   			else {
    			print STDOUT 'Error'; 
   			}

			$i++;
		}

$ibeg = 0;
$iend = 0;
$i = 0;

return (\@beginningdomain, \@endingdomain, $nseg);
}

	
sub interdomainsort {
	
local ($interdomain) = @_;
my @interdomainnumbers;
my @beginninginterdomainunsorted;
my @endinginterdomainunsorted;

# split interdomain information into an array
@interdomainnumbers = split(/,/, $interdomain);

	# divide domain information into beginning and ending regions of the domains 	
	foreach $interdomainnumbers (@interdomainnumbers) {
   	
		# if array entry is divisable by 2 (even numbered) then it is put in the beginning array
		if ($e % 2 == 0) {
		$kbeg++;
		$beginninginterdomainunsorted[$kbeg]= $interdomainnumbers[$e];
		}
  
		# if array entry is not divisable by 2 (odd numbered) then it is put in the ending array
  	 	elsif ($e % 2 != 0) {
		$kend++;
		$endinginterdomainunsorted[$kend]= $interdomainnumbers[$e];
		}
   	
		# print error if something goes wrong
		else {
 		print STDOUT 'Error'; 
 		}

		$e++;
	}

$kbeg = 0;
$kend = 0;
$e = 0;

@beginninginterdomain = sort @beginninginterdomainunsorted; 
@endinginterdomain = sort @endinginterdomainunsorted;

return (\@beginninginterdomain, \@endinginterdomain);

}

sub restoatom {

local ($resnum, $totalnumber)= @_;
my $iatres;
my @natres;
my @res_to_atom;

	#sorting residue to atom number
	$iatres=0;

	$resnum[0]=$resnum[1];

	for ($iat=1; $iat<=$totalnumber; $iat++) {

		#if atom is in the same amino acid
		if ($resnum[$iat] == $resnum[$iat-1]) {
		
		$iatres++;
		$res_to_atom[$resnum[$iat]][$iatres] = $iat;
		}

		#else if atom is in different amino acid
		else {

		$natres[$resnum[$iat-1]] = $iatres;
		$iatres = 1;
		$res_to_atom[$resnum[$iat]][$iatres] = $iat;
		}
	}

return (\@res_to_atom, \@natres);
}

#Sorts through moving coordinate information for calculation
sub coordinatesorter {

local (*Axcoord, *Aycoord, *Azcoord, *Bxcoord, *Bycoord, *Bzcoord, *res_to_atom, *natres, $nseg, *beginning, *ending) = @_;

my @Axcord;
my @Aycord;
my @Azcord;
my @Bxcord;
my @Bycord;
my @Bzcord;
#my $interdomaincounter=1;
my $counter = 0;

	for ($iseg = 1; $iseg <= $nseg; ++$iseg){

		for ($ires = $beginning[$iseg]; $ires <= $ending[$iseg]; ++$ires) {

	        	for ($iat = 1; $iat <= $natres[$ires]; ++$iat) {

				$Axcord->[$counter] = $Axcoord[$res_to_atom[$ires][$iat]];
				$Aycord->[$counter] = $Aycoord[$res_to_atom[$ires][$iat]];
				$Azcord->[$counter] = $Azcoord[$res_to_atom[$ires][$iat]];
				$Bxcord->[$counter] = $Bxcoord[$res_to_atom[$ires][$iat]];
				$Bycord->[$counter] = $Bycoord[$res_to_atom[$ires][$iat]];
				$Bzcord->[$counter] = $Bzcoord[$res_to_atom[$ires][$iat]];

print STDOUT $Axcord[$counter];

			$counter++;

			}
		}
	}

print STDOUT "should print out @Axcord array but does'nt?!", @Axcord;
return (\@Axcord, \@Aycord, \@Azcord, \@Bxcord, \@Bycord, \@Bzcord);

}

#Distance Calculation
sub calculation {			

local ($Mxcoord, $Mycoord, $Mzcoord, $Fxcoord, $Fycoord, $Fzcoord) = @_;

my $fixedcoordsize = @Fxcoord; 
my $movingcoordsize = @Mxcoord;
my $clashcounter;

	for($M=1; $M <= $fixedcoordsize; $M++){
   
		for($F=1; $F <= $movingcoordsize; $F++){

		$distance[$F] = sqrt(($Mxcoord[$F] - $Fxcoord[$M])**2 + ($Mycoord[$F] - $Fycoord[$M])**2 + ($Mzcoord[$F] - $Fzcoord[$M])**2);      
   
			if ($distance[$F] <= 4) {

			$clashcounter++;
			#print STDOUT "$atomtype[$F]";
			#print STDOUT "$resnum[$F]";
			#print STDOUT "$atomcount[$F]\n"; 
			}
    		}
	}

return ($clashcounter);
}