use strict;
use lib 'd:/Programs/Perl/lib';

use Utils;
use JFile;
use Crystal::CIF;
use Crystal::Crystal;

my $CIFPath = "aIGZO-Dens608-PBE/VCRelax/a-IGZOm184-Dens608-MDRelaxed-final.cif";
my $OutFile = "CN.csv";
my $MaxR  = 4; # A
my $MaxR2 = 2.5; # A

my $out = new JFile($OutFile, "w");
$out->print("label1,nM,nIn,nGa,Zn,nO\n");

my $cif = new CIF;
$cif->Read($CIFPath);
my $Crystal = $cif->GetCCrystal();
my @latt = $Crystal->LatticeParameters();
print("Latt: ", join(', ', @latt), "\n");

$Crystal->ExpandCoordinates();
my @AtomTypeList = $Crystal->GetCAtomTypeList();
my @AtomSiteList = $Crystal->GetCAsymmetricAtomSiteList();
my @ExpandedAtomSiteList = $Crystal->GetCExpandedAtomSiteList();

my $site1 = $ExpandedAtomSiteList[0];
my $label1       = $site1->Label();
my $type1        = $site1->AtomName();
my ($x1,$y1,$z1) = $site1->Position(1);
my $site2 = $ExpandedAtomSiteList[74];
my $label2       = $site2->Label();
my $type2        = $site2->AtomName();
my ($x2,$y2,$z2) = $site2->Position(1);
#my $d = $Crystal->GetNearestInterAtomicDistance(0.313, 0.618, 0.877, 0.391, 0.461, 0.971);
my $d = $Crystal->GetNearestInterAtomicDistance($x1, $y1, $z1, $x2, $y2, $z2);
print "d=$d: ($x1,$y1,$z1)-($x2,$y2,$z2): $label1-$label2\n";

my $c = 1;
for(my $i = 0 ; $i < @ExpandedAtomSiteList ; $i++) {
	my @List;

	my $site1 = $ExpandedAtomSiteList[$i];
	my $label1       = $site1->Label();
	my $type1        = $site1->AtomName();
	my ($x1,$y1,$z1) = $site1->Position(1);
	for(my $j = 0 ; $j < @ExpandedAtomSiteList ; $j++) {
		next if($i == $j);

		my $site2 = $ExpandedAtomSiteList[$j];
		my $label2       = $site2->Label();
		my $type2        = $site2->AtomName();
		my ($x2,$y2,$z2) = $site2->Position(1);

		my $d = $Crystal->GetNearestInterAtomicDistance($x1, $y1, $z1, $x2, $y2, $z2);
#		my $dx = $latt[0] * ($x1 - $x2);
#		my $dy = $latt[1] * ($y1 - $y2);
#		my $dz = $latt[2] * ($z1 - $z2);
#		my $d2 = $dx*$dx + $dy*$dy + $dz*$dz;
#		$d2 = sqrt($d2);
#print "$d - $d2\n";
#		if(1) {
		if($d < $MaxR) {
#print "deb: $c: $label1 - $label2: $d\n";
#print "deb: $c: $label1 - $label2: ($x1,$y1,$z1)-($x2,$y2,$z2): $d\n";
			my %hash = (
				i      => $i,
				j      => $j,
				c      => $c,
				label1 => $label1,
				type1  => $type1,
				x1     => $x1,
				y1     => $y1,
				z1     => $z1,
				label2 => $label2,
				type2  => $type2,
				x2     => $x2,
				y2     => $y2,
				z2     => $z2,
				d      => $d,
			);
			@List = (@List, \%hash);
#			@List = push(@List, \%hash);
			$c++;
		}
	}
	
	@List = sort { $a->{d} <=> $b->{d}; } @List;
	my ($nIn, $nGa, $nZn, $nO) = (0, 0, 0, 0);
	foreach my $l (@List) {
		if($l->{d} < $MaxR2) {
			$nIn++ if($l->{label2} =~ /In/);
			$nGa++ if($l->{label2} =~ /Ga/);
			$nZn++ if($l->{label2} =~ /Zn/);
			$nO++  if($l->{label2} =~ /O/);
		}
	}
	my $nM = $nIn + $nGa + $nZn;

	print "Atom($i): $label1\n";
	$out->print("$label1,$nM,$nIn,$nGa,$nZn,$nO,");
	foreach my $l (@List) {
#print "l=$l\n";
		print "  $l->{c}: $l->{label2}: $l->{d}\n";
		$out->print("$l->{label2},$l->{d},");
	}
	$out->print("\n");
}


$out->Close();
exit;