#!/usr/bin/perl

BEGIN {
use lib 'd:/Programs/Perl/lib';
use lib '/home/tkamiya/bin/lib';
my $BaseDir = $ENV{'TkPerlDir'};
#print "\nBaseDir: $BaseDir\n";
@INC = ("$BaseDir/lib", "$BaseDir/VNL", "d:/Programs/Perl/lib", @INC);
}

use strict;
#use warnings;

use Sci qw($NA);
use Sci::ChemicalReaction;
use Crystal::CIF;
use Crystal::Crystal;

my $Composition = "Ga2Sn6O15";
#my $Composition = "Ga2Sn6O15He2";
$Composition = $ARGV[0] if($ARGV[0] ne '');
my $Density = 6.0;
$Density = $ARGV[1] if($ARGV[1] ne '');

my $OutFile  = "$Composition.cif";
my $xyzsigma = 0.2; # in A
print "Composition  : $Composition\n";
print "Density      : $Density g/cm3\n";
print "Output CIF   : $OutFile\n";
print "(x,y,z) sigma: $xyzsigma A\n";

#=============================================
# Create ChemicalReaction object
#=============================================
my $R = new ChemicalReaction(); #\@DBFiles, $BaseDBDir); #"Compound" , "Composition");
if(!$R) {
	print "Error: $!: Can not create ChemicalReaction Object.\n";
	exit;
}

#=============================================
# Analyze reaction
#=============================================
$R->AnalyzeAReaction('Reagents', $Composition);
my $pReagentHash = $R->ReagentsHash();
my $pElements    = $R->ElementsArray("Reagents");
#my $pReagentCompounds = $R->pCompounds('Reagents');
#my $pReagentElements  = $R->pElements("Reagents");

#=============================================
# Batch calculation
#=============================================
my %ToArabic = (
	I    => 1,
	II   => 2,
	III  => 3,
	IV   => 4,
	V    => 5,
	VI   => 6,
	VII  => 7,
	VIII => 8,
	IX   => 9,
	X    => 10,
	);

my $nIon = 0;
for(my $i = 0 ; $i < @$pElements ; $i++) {
		my $e = $pElements->[$i];
		my $n = $R->nElement("Reagents", $e);
		$nIon += $n;
}
print "Number of atoms: $nIon\n";

print "\n";
my $wMol        = 0.0;
my $TotalCharge = 0.0;
my (@CName, @AName, @NName);
my ($cC, $cA, $cN) = (0, 0, 0);
my ($nCation, $nNeutral, $nAnion) = (0, 0, 0);
for(my $i = 0 ; $i < @$pElements ; $i++) {
		my $e = $pElements->[$i];
		my $n = $R->nElement("Reagents", $e);
print "$e: $n\n";
		my $pAtomDB = AtomType::GetAtomDBHash($e);
		my $pIonRadius = $pAtomDB->{pIonRadius};
#		my $charge  = $pAtomDB->{charge};
		my $charge        = $pIonRadius->[0];
		$charge = ($charge =~ /-/)? -$charge + 0.0 : $charge + 0.0;
		my $nCoordination = $pIonRadius->[1];
		$nCoordination = $ToArabic{$nCoordination};
		my $ri            = $pIonRadius->[2] + 0.0;

		if($charge !~ /^[\+\-0-9]+$/) {
			print "Error 1: Invalid charge [$charge]\n";
			exit;
		}
		if($charge == 0) {
			$nNeutral += $n;
		}
		elsif($charge > 0) {
			$nCation += $n;
		}
		elsif($charge < 0) {
			$nAnion += $n;
		}
		else {
			print "Error2: Invalid charge [$charge]\n";
			exit;
		}
		$TotalCharge += $charge * $n;
print "tc += $charge * $n\n";

		$pAtomDB    = AtomType::GetAtomDBHash($e, $charge);
print "  Z=$pAtomDB->{Z}: M=$pAtomDB->{mass}\n";
print "    Charge=$charge: nCood=$nCoordination: ri=$ri\n";
		$wMol += $n * $pAtomDB->{mass};
		
		for(my $j = 0 ; $j < $n ; $j++) {
			if($charge == 0) {
				$NName[$cN++] = $e;
			}
			elsif($charge > 0) {
				$CName[$cC++] = $e;
			}
			elsif($charge < 0) {
				$AName[$cA++] = $e;
			}
		}
	}
print "Molecular weight: $wMol\n";
print "Total charge: $TotalCharge\n";
if(abs($TotalCharge) > 1.0e-3) {
	print "Error: Non-zero total charge [$TotalCharge]\n";
	exit;
}
print "Number of ions: cations=$nCation, anions=$nAnion, netral=$nNeutral\n";

print "\n";
my $V = ($wMol / $NA) / $Density * 1.0e24; # in A
my $latt = $V**(1.0/3.0);
print "Volume: $V A^3 for density = $Density g/cm3\n";
print "  Lattice parameter: $latt A\n";

my $v1 = $V / $nIon;
my $a1 = $v1**(1.0/3.0);
print "Volume per atom: $v1 A^3 (a = $a1 A)\n";

print "\n";
my $nCation2 = ($nCation > $nAnion)? $nCation : $nAnion;
my $nAnion2  = ($nCation > $nAnion)? $nCation : $nAnion;
my $nTotal = $nCation2 + $nAnion2 + $nNeutral;
print "Number of ions adjusted: caions=$nCation2, anions=$nAnion2, netral=$nNeutral, total=$nTotal\n";
my $n1 = int($nTotal**(1.0/3.0) + 0.99);
$n1++ if($n1 % 2 != 0);
my $nSite = $n1*$n1*$n1;
print "Number of ions in 1D: $n1  nSite: $nSite\n";
my $R1 = $latt / $n1;
print "Interatomic distance: $R1 A\n";

print "\n";
my (@CV, @AV);
my $nCV = $nSite / 2;
my $nAV = $nSite / 2;
print "$nCation cation sites from $nCV sites:\n";
for(my $i = 0 ; $i < $nCation ; $i++) {
	my $idx = int(rand($nCV));
	my $idx2 = &FindVoid(\@CV, $nSite/2, $idx, $CName[$i]);
print "$i: $CName[$i]: $idx2: seach by $idx in $nCV\n";
	$nCV--;
}

print "\n";
print "$nAnion anion sites from $nAV sites:\n";
for(my $i = 0 ; $i < $nAnion ; $i++) {
	my $idx = int(rand($nAV));
	my $idx2 = &FindVoid(\@AV, $nSite/2, $idx, $AName[$i]);
print "$i: $AName[$i]: $idx/$idx2 in $nAV\n";
	$nAV--;
}

print "\n";
my $nNV = $nCV + $nAV;
print "$nNeutral neutral sites from $nNV sites:\n";
for(my $i = 0 ; $i < 2 ; $i++) {
	my $idx = int(rand($nNV));
	my $idx2 = &FindNeutralVoid(\@CV, \@AV, $nSite/2, $nSite/2, $idx, $NName[$i]);
print "$i: $NName[$i]: $idx/$idx2 in $nNV\n";
	$nNV--;
}

print "\n";
print "Site list (n in 1D: $n1)\n";
my ($iC, $iA) = (0, 0);
for(my $iz = 0 ; $iz < $n1 ; $iz++) {
for(my $iy = 0 ; $iy < $n1 ; $iy++) {
for(my $ix = 0 ; $ix < $n1 ; $ix++) {
	if(($ix+$iy+$iz) % 2 == 0) {
		$CV[$iC] = {Name => $CV[$iC], x => $ix / $n1, y => $iy / $n1, z => $iz / $n1};
#print "iC=$iC: ($CV[$iC]->{x}, $CV[$iC]->{y}, $CV[$iC]->{z})\n";
		$iC++;
	}
	else {
		$AV[$iA] = {Name => $AV[$iA], x => $ix / $n1, y => $iy / $n1, z => $iz / $n1};
#print "iA=$iA: ($AV[$iA]->{x}, $AV[$iA]->{y}, $AV[$iA]->{z})\n";
		$iA++;
	}
}
}
}

my $crystal = new Crystal();
$crystal->SetLatticeParameters($latt, $latt,$latt, 90.0, 90.0, 90.0);

for(my $i = 0 ; $i < @$pElements ; $i++) {
	my $e = $pElements->[$i];
	$crystal->AddAtomType($e, 0, 1.0);
}

for(my $i = 0 ; $i < $nSite/2 ; $i++) {
	next if($CV[$i]->{Name} eq '');

	my $Name = $CV[$i]->{Name};
	my $sx = (1.0 - rand(2.0)) * $xyzsigma/$latt;
	my $sy = (1.0 - rand(2.0)) * $xyzsigma/$latt;
	my $sz = (1.0 - rand(2.0)) * $xyzsigma/$latt;
	my $x    = $CV[$i]->{x} + $sx;
	my $y    = $CV[$i]->{y} + $sy;
	my $z    = $CV[$i]->{z} + $sz;
	print "C$i: $Name ($x, $y, $z)\n";
	my $i1 = $i + 1;
	$crystal->AddAtomSite("$Name$i1", $Name, $x, $y, $z, 1.0);
}

for(my $i = 0 ; $i < $nSite/2 ; $i++) {
	next if($AV[$i]->{Name} eq '');

	my $Name = $AV[$i]->{Name};
	my $sx = (1.0 - rand(2.0)) * $xyzsigma/$latt;
	my $sy = (1.0 - rand(2.0)) * $xyzsigma/$latt;
	my $sz = (1.0 - rand(2.0)) * $xyzsigma/$latt;
	my $x    = $AV[$i]->{x} + $sx;
	my $y    = $AV[$i]->{y} + $sy;
	my $z    = $AV[$i]->{z} + $sz;
	print "A$i: $Name ($x, $y, $z)\n";
	my $i1 = $i + 1;
	$crystal->AddAtomSite("$Name$i1", $Name, $x, $y, $z, 1.0);
}
$crystal->ExpandCoordinates();

print "\n";
my $CIF = new CIF;
if($CIF->CreateCIFFileFromCCrystal($crystal, $OutFile, 0, 'unix')) {
	print "CIF [$OutFile] has been created.\n";
}
else {
	print "Cannot write to [$OutFile].\n";
}

exit;

sub FindVoid
{
	my ($p, $n, $idx, $name) = @_;

	my $c = 0;
	for(my $i = 0 ; $i < $n ; $i++) {
#print "$i($p->[$i])/$n: idx=$idx: c=$c\n";
		if($p->[$i] ne '') {
			next;
		}
		else {
			if($idx == $c) {
				$p->[$i] = $name;
				return $i;
			}
			$c++;
		}
	}
	
	print "Error: Can not find void site\n\n";
exit;
	return undef;
}

sub FindNeutralVoid
{
	my ($pC, $pA, $nC, $nA, $idx, $name) = @_;

	my $c = 0;
	for(my $i = 0 ; $i < $nC ; $i++) {
#print "$i($pC->[$i])/$nC: idx=$idx: c=$c\n";
		if($pC->[$i] ne '') {
			next;
		}
		else {
			if($idx == $c) {
				$pC->[$i] = $name;
				return $i;
			}
			$c++;
		}
	}
	for(my $i = 0 ; $i < $nA ; $i++) {
#print "$i($p->[$i])/$n: idx=$idx: c=$c\n";
		if($pA->[$i] ne '') {
			next;
		}
		else {
			if($idx == $c) {
				$pA->[$i] = $name;
				return $i;
			}
			$c++;
		}
	}

	print "Error: Can not find void site\n\n";
exit;
	return undef;
}