#!/usr/bin/perl

use Math::Trig;
use Math::Round;
use Math::VectorReal;
use Math::MatrixReal;

#=== This program is described on.. ?

# --- Program defaults are here!
$centralAtomName = "Cu";
$ligandAtomName = "O";
$coordination = 4;
$minNNDistAngstroms = 1.5;
$maxNNDistAngstroms = 2.3;
@cellsPerSupercell = (10, 10, 7);


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


$eyeFile = $ARGV[0];
@supercellEdge = (0, 0, 0);

print "\nThis is local_cmco.pl. Hit enter at a prompt to use default values.\n";

print "\nEnter the number of unit cells in each direction (default: @cellsPerSupercell): ";

$userInput = <STDIN>;

chomp($userInput);
if ($userInput ne "") { @cellsPerSupercell = split(" ",$userInput); }

print "Using supercell size = @cellsPerSupercell\n";


print "\nEnter the symbol of the central atom (default: $centralAtomName): ";

$userInput = <STDIN>;
chomp($userInput);
if ($userInput ne "") { $centralAtomName = $userInput; }

print "Using central atom = $centralAtomName\n";


print "\nEnter the symbol of the ligand atom (default: $ligandAtomName): ";

$userInput = <STDIN>;
chomp($userInput);
if ($userInput ne "") { $ligandAtomName = $userInput; }

print "\nUsing ligand atom = $ligandAtomName\n";

print "\nEnter the coordination number (default: $coordination): ";

$userInput = <STDIN>;
chomp($userInput);
if ($userInput ne "") { $coordination = $userInput; }

print "\nUsing coordination = $coordination\n";

print "\nEnter the minimum and maximum coordination distances in A (default: $minNNDistAngstroms $maxNNDistAngstroms): ";

$userInput = <STDIN>;
chomp($userInput);
if ($userInput ne "") {
	@nnDistInput = split(" ",$userInput);
	$minNNDistAngstroms = $nnDistInput[0];
	$maxNNDistAngstroms = $nnDistInput[1];
}

print "\nUsing min distance = $minNNDistAngstroms A and max distance = $maxNNDistAngstroms A\n";


open (READIT, $eyeFile) || die("** Can't read file! ** \nScript use: perl ab_split.pl [input].cfg\n");
@eyeIn = <READIT>;
close (READIT);


print "Done reading in file $eyeFile \n";


$i=0;
foreach (@eyeIn) {
	chomp();
	if (/Number of particles \=/) {
		$totalAtoms = trim($');
	}

	if (/H0\(1,1\) = /) {
		$supercellEdge[0] = trim($');		
	}
	
	if (/H0\(2,2\) = /) {
		$supercellEdge[1] = trim($');		
	}

	if (/H0\(3,3\) = /) {
		$supercellEdge[2] = trim($');		
	}

	if (!/^\w/) {
		push(@atomIn, $eyeIn[$i]);
	}

	$i++;
	
}

foreach (@supercellEdge) {
	chop(); chop;
}

print "I found $totalAtoms atoms\n";
print "Supercell dimensions are @supercellEdge Angstroms\n\n";

$maxNNDist2Angstroms = $maxNNDistAngstroms*$maxNNDistAngstroms;
$minNNDist2Angstroms = $minNNDistAngstroms*$minNNDistAngstroms;
$maxXNNDist = $maxNNDistAngstroms/$supercellEdge[0];
$maxYNNDist = $maxNNDistAngstroms/$supercellEdge[1];
$maxZNNDist = $maxNNDistAngstroms/$supercellEdge[2];

print "Squared distances are $minNNDist2 and $maxNNDist2 \n";

$lastAtomType = "XX";
$i=0;
$j=0;
$k=0;
foreach (@atomIn) {
	my @atomLine = split(/\s+/, $atomIn[$i]);
	if ($atomLine[2] eq $lastAtomType) {
		$atom[$i][0] = $k;
		$j++;
	} else {
		if ($j > 0 ) {
			push(@totalAtomsOfType, $j);
			print "Total number of $lastAtomType atoms is $j\n";
			$k++;
		}
		print "Found new atom type: $atomLine[2]\n";

		if ($centralAtomName eq $atomLine[2]) {
			$centralAtomType = $k;
			print "** Central atom is now $centralAtomName (type $k)\n";
		}
		if ($ligandAtomName eq $atomLine[2]) {
			$ligandAtomType = $k;
			print "** Ligand atom is now $ligandAtomName (type $k)\n";
		}
		$atom[$i][0] = $k;
		$lastAtomType = $atomLine[2];
		$j=1;
	}
	$atom[$i][1] = $atomLine[3];
	$atom[$i][2] = $atomLine[4];
	$atom[$i][3] = $atomLine[5];

	$i++;
}
print "Total number of $lastAtomType atoms is $j\n\n";
push(@totalAtomsOfType, $j);


print "Atomic positions are loaded into memory!\n";

#for this script we could load the central and ligand arrays directly, but for general purpose let's split up the two processes.
$j=0; $k=0;
for ($i = 0; $i < $#atom; $i++) {
	if ($atom[$i][0] == $centralAtomType) {
		$centralAtom[$j][0] = $atom[$i][1];
		$centralAtom[$j][1] = $atom[$i][2];
		$centralAtom[$j][2] = $atom[$i][3];
		$j++;
	} #this used to be an elsif, but that fails if center/ligand are the same
	if ($atom[$i][0] == $ligandAtomType) {
		$ligandAtom[$k][0] = $atom[$i][1];
		$ligandAtom[$k][1] = $atom[$i][2];
		$ligandAtom[$k][2] = $atom[$i][3];
		$k++;
	}
}

$coordGood = 0;
$coordBad = 0;

print "done with central/ligand arrays. central: $totalAtomsOfType[$centralAtomType] ligand: $totalAtomsOfType[$ligandAtomType]\n";

print "finding nearest neighbors...\n";

# search over all central atoms and give a list of central atom xyz and corresponding ligand xyz
for ($i=0; $i < $totalAtomsOfType[$centralAtomType]; $i++) {

	
	$numberNN = 0;
	$ligandLoopCount = 0;

	# search over all oxygens to find NN pairs, while correcting for PBC
	for ($j=0; $j < $totalAtomsOfType[$ligandAtomType]; $j++) {

		$xDist = $ligandAtom[$j][0]- $centralAtom[$i][0];
		$yDist = $ligandAtom[$j][1]- $centralAtom[$i][1];
		$zDist = $ligandAtom[$j][2]- $centralAtom[$i][2];

		if ($xDist > 0.5) {
			$xDist = $xDist - 1;
		} elsif ($xDist < -0.5) {
			$xDist = $xDist + 1;
		}

		if (abs($xDist) > $maxXNNDist) {
			next;
		}

		if ($yDist > 0.5) {
			$yDist = $yDist - 1;
		} elsif ($yDist < -0.5) {
			$yDist = $yDist + 1;
		}

		if (abs($yDist) > $maxYNNDist) {
			next;
		}

		if ($zDist > 0.5) {
			$zDist = $zDist - 1;
		} elsif ($zDist < -0.5) {
			$zDist = $zDist + 1;
		}

		if (abs($zDist) > $maxZNNDist) {
			next;
		}

		#old version (supercell units) assumed that the 3 directions are on the same scale, i.e. **CUBIC**
		#check bvs and angle scripts!

		#$vectDist2 = $xDist*$xDist + $yDist*$yDist + $zDist*$zDist;

		$xDistAngstroms = $xDist*$supercellEdge[0];
		$yDistAngstroms = $yDist*$supercellEdge[1];
		$zDistAngstroms = $zDist*$supercellEdge[2];

		$vectDist2Angstroms = $xDistAngstroms*$xDistAngstroms + $yDistAngstroms*$yDistAngstroms + $zDistAngstroms*$zDistAngstroms;

		if (($vectDist2Angstroms < $maxNNDist2Angstroms) && ($vectDist2Angstroms > $minNNDist2Angstroms)) {
			$neighborVector[$i][$numberNN] = vector($xDistAngstroms, $yDistAngstroms, $zDistAngstroms);
			$neighbors[$numberNN] = $j;
			$numberNN++;
		}

	}

	if ($numberNN == $coordination) {

		$coordGood++;
	} else {
		print "bad atom $i has $numberNN neighbors\n";
		push(@coordBadCentralAtoms,$i);
		$coordBad++;
		$coordBadSum += $numberNN;
	}

	# calculate bond valence sums for each polyhedron, for each oxidation state [n/a here]
	
	# calculate bond angles for the NN pairs [n/a here]

	

}


print "\nPolyhedra found for $coordGood central atoms, failed for $coordBad central atoms. ";
print $coordGood/($coordGood+$coordBad)*100 . "\%\n";
if ($coordBad > 0) { print "Average coordination for failed central atoms is " . $coordBadSum/$coordBad . "\n"; }


if ($centralAtomName eq "Cu") {
	# calculate bond valence sums for each polyhedron, for each oxidation state
	# bvs values from Brese and O'Keeffe, Acta Cryst. B47, p192 (1991),
	# although Cu+ is given as 1.600 in I. D. Brown, J. Solid State Chem. v90 p 155 (1991)
	$minValence = 1;
	$maxValence = 3;
	$R1Center[1] = 1.593; 
	$R1Center[2] = 1.679;
	$R1Center[3] = 1.739;

} elsif ($centralAtomName eq "Mn") {
	$minValence = 2;
	$maxValence = 4;
	$R1Center[2] = 1.790; 
	$R1Center[3] = 1.760;
	$R1Center[4] = 1.753;

} elsif ($centralAtomName eq "Mg") {
	$minValence = 2;
	$maxValence = 2;
	$R1Center[2] = 1.693; 

} elsif ($centralAtomName eq "Cr") {
	$minValence = 3;
	$maxValence = 4;
	$R1Center[3] = 1.724; 
	$R1Center[4] = 1.794;  

} elsif ($centralAtomName eq "Bi") {
	$minValence = 3;
	$maxValence = 5;
	$R1Center[3] = 2.09; 
	$R1Center[5] = 2.06;

} elsif ($centralAtomName eq "Ti") {
	$minValence = 3;
	$maxValence = 4;
	$R1Center[3] = 1.791;
	$R1Center[4] = 1.815;  

} elsif ($centralAtomName eq "Ru") {
	$minValence = 4;
	$maxValence = 4;
	$R1Center[4] = 1.834;  

} else {
	print "\n** I don't have BVS parameters for central atom type $centralAtomName. Exiting.\n";
	exit;
}

for ($i=0; $i < $totalAtomsOfType[$centralAtomType]; $i++) {

	$skipThisBVS = 0;
	if ($coordBad > 0) {
		for ($j=0; $j < $coordBad; $j++) {
			if ($i == $coordBadCentralAtoms[$j]) {
				$skipThisBVS = 1;
			}
		}
	} else {
		$skipThisBVS = 0;
	}

	if ($skipThisBVS == 1) {
		next;
	} 


	#note that 1/b = 1/0.37 = 2.7027

	for ($j = $minValence; $j <= $maxValence; $j++) {
		$BVS[$j] = 0;
	}

	# loop over allowed valence states
	for ($j = $minValence; $j <= $maxValence; $j++) {

		for ($k = 0; $k < $coordination; $k++) {
			$vector = $neighborVector[$i][$k]->length;
			$BVS[$j] += exp(($R1Center[$j] - $vector)*2.7027);
		}

		$BVSDiff[$j] = abs($j - $BVS[$j]);
	}

	# if there are multiple valence choices find the best match
	# otherwise just output the single BVS
	if ($maxValence > $minValence) {
		for ($j = $minValence; $j < $maxValence; $j++) {
			if ($BVSDiff[$j] < $BVSDiff[$j+1]) {
				$BVSMin = $BVSDiff[$j];
				$BVSPicked = $BVS[$j];
				$valence = $j;
			} else {
				$BVSMin = $BVSDiff[$j+1];
				$BVSPicked = $BVS[$j+1];
				$valence = $j+1;
			}
		}
	} else {
		$BVSPicked = $BVS[$minValence];
		$valence = $minValence;
	}

	push(@BVSum,$BVSPicked);
	push(@BValence,$valence);
	

	# calculate bond angles for the NN pairs
	for ($j = 0; $j < $coordination -1; $j++) {

		for ($k = $j+1; $k < $coordination; $k++) {


			$dotP = $neighborVector[$i][$j] . $neighborVector[$i][$k];

			$lengthsP = $neighborVector[$i][$j]->length * $neighborVector[$i][$k]->length;

			$angle = 180 * acos($dotP / $lengthsP) / 3.14159;

			push(@angles,$angle);

		}

	}

}

print "\n\n";


@angleHist = histogram(\@angles,0,180,0.25);

@BVSumHist = histogram(\@BVSum,0,5,0.05);

@BValenceHist = histogram(\@BValence,0,($maxValence+1),1);

$centralAtomName =~ tr/[A-Z]/[a-z]/;

$newExtension = "_$centralAtomName" . "_angles.dat";

$eyeFile =~ s/\.cfg/$newExtension/;

print "writing output:\nxmgrace $eyeFile\n";
open (WRITEANGLES, ">$eyeFile") || die("** Can't read file! ** \nScript use: perl ab_split.pl [output].cfg\n");

print WRITEANGLES "\#$centralAtomName coordination $coordination angles\n";
print WRITEANGLES join("\n",@angleHist);

close (WRITEANGLES);


$eyeFile =~ s/angles/bvs/;
print "writing output:\nxmgrace $eyeFile\n";
open (WRITEBVS, ">$eyeFile") || die("** Can't read file! ** \nScript use: perl ab_split.pl [output].cfg\n");

print WRITEBVS "\#$centralAtomName coordination $coordination bv state\n";
print WRITEBVS join("\n",@BValenceHist);

print WRITEBVS "\n&newxy\;\n\#$centralAtomName coordination $coordination bv sum\n";
print WRITEBVS join("\n",@BVSumHist);

close (WRITEBVS);


# report
print "Done! \n\n";


sub mod { return $_[0] - int ($_[0] / $_[1]) * $_[1]; };

sub histogram {
	my(@HistOut);
	my($entry);
	# usage: histogram(@list, min, max, binsize)
	# returns an array with "BIN (tab) COUNTS" from min to max
	# uses nearest rounding, not floor or ceiling
	for ($i = 0; $i < (($_[2] - $_[1])/$_[3]); $i++) {
		$HistOut[$i] = 0;
	}
	
	foreach $entry (@{$_[0]}) {

		if ($entry > $_[2]) {
			print "--- a histogram value is outside the bin range! things may break\n";
		}

		$entry = nearest($_[3],$entry);
		$HistOut[nearest(1,(($entry-$_[1])/$_[3]))]++;
	}

	$i = $_[1];
	foreach $entry (@HistOut) {
		$entry = nearest($_[3],$i) . "\t$entry";
		$i += $_[3];
	}
	return @HistOut;
}


sub trim($)
{
	my $foo = shift;
	$foo =~ s/^\s+//;
	$foo =~ s/\s+$//;
	return $foo;
}