use strict;
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 Crystal::Crystal;
use Crystal::CIF;
use Sci::ChemicalReaction;

my $a = 15.0;
my $CreateAtomCIFs = 0;

my @molecules = qw();
#my @molecules = qw(CH4 SiH4 GeH4 SnH4);
#my @molecules = qw(BH3 NH3 PH3 AsH3);
#my @molecules = qw(H2O CO2 NO2 SO2 H2S);
#my @molecules = qw(CO NO SO);
#my @molecules = qw(H2 N2 O2 F2 Cl2 Br2 I2 HF HCl HBr HI);
#BH3 CH4 NH3 H2O SiH4 PH3 H2S GeH4 AsH3 H2Se);
my $MetalCIFTemplate = "{Structure}-Template.cif";
#my @metals = ();
my @metals = qw(Li Be B C Na Mg Al Si K Ca Sc Ti Ni Cu Zn Ga Ge Rb Sr Y Zr Nb Mo Rh Pd Ag Cd In Sn Sb Te Cs Ba La Hf Ta Au Hg Tl Pd Bi Po);
my $ABCIFTemplate = "{Structure}-Template.cif";
my @ABs = ();
#my @ABs = qw(LiF LiCl LiBr LiI NaF NaCl NaBr NaI KF KCl KBr KI RbF RbCl RbBr RbI CsF CsCl CsBr CsI 
#			BeO BeS BeSe BeTe MgO MgS MgSe MgTe CaO CaS CaSe CaTe SrO SrS SrSe SrTe BaO BaS BaSe BaTe 
#			NiAs ZnO ZnS ZnSe ZnTe GaN GaP GaAs GaSb InN InP InAs InSb AlN AlP AlAs AlSb);

my $CIFFile;

if($CreateAtomCIFs) {
	for(my $i = 1 ; ; $i++) {
		my $AtomName = AtomType::GetAtomName($i);
		last if($AtomName eq '');

		print("$i: $AtomName\n");
		$CIFFile = "${AtomName}-atom.cif";
		&CreateMoleculeCIF($AtomName, $a, $CIFFile);
	}
}

for(my $i = 0 ; $i < @molecules ; $i++) {
	print("$i: $molecules[$i]\n");

	my (@a) = Utils::FindAll("[A-Z][a-z]*\\d*", $molecules[$i]);
	if(@a >= 2 or $a[0] =~ /\d/) {
		$CIFFile = "$molecules[$i]-molecule.cif";
	}
	&CreateMoleculeCIF($molecules[$i], $a, $CIFFile);
}

for(my $i = 0 ; $i < @metals ; $i++) {
	print("$i: $metals[$i]\n");

	my $AtomNum1 = AtomType::GetAtomicNumber($metals[$i]);
	my ($filepath, $AtomJName, $AtomEName,
		$AtomicMass, $FoundYear, $AtomicRadius1,
		$MeltingPoint, $BoilingPoint, $Density,
		$HeatCapacity, $IonizationPotential, $ElectronAffinity)
				 = AtomType::GetIonInfFromPeriodicTable2($AtomNum1);
	my $d = $AtomicRadius1 + $AtomicRadius1;

	my $Template = $MetalCIFTemplate;
	$Template =~ s/\{Structure\}/FCC/;
	$CIFFile = "$metals[$i]-FCC-metal.cif";
	my $a = $d * sqrt(2);
	&CreateMetalCIF($metals[$i], $Template, $a, $a, $a, $CIFFile);

	my $Template = $MetalCIFTemplate;
	$Template =~ s/\{Structure\}/BCC/;
	$CIFFile = "$metals[$i]-BCC-metal.cif";
	my $a = $d * 2.0 / sqrt(3.0);
	&CreateMetalCIF($metals[$i], $Template, $a, $a, $a, $CIFFile);

	my $Template = $MetalCIFTemplate;
	$Template =~ s/\{Structure\}/HCP/;
	$CIFFile = "$metals[$i]-HCP-metal.cif";
	my $a = $d;
	my $c = $d * 2.0 / 3.0 * sqrt(6.0);
	&CreateMetalCIF($metals[$i], $Template, $a, $a, $c, $CIFFile);

	my $Template = $MetalCIFTemplate;
	$Template =~ s/\{Structure\}/SC/;
	$CIFFile = "$metals[$i]-SC-metal.cif";
	my $a = $d;
	&CreateMetalCIF($metals[$i], $Template, $a, $a, $a, $CIFFile);

	my $Template = $MetalCIFTemplate;
	$Template =~ s/\{Structure\}/Diamond/;
	$CIFFile = "$metals[$i]-Diamond-metal.cif";
	my $a = $d * 4.0 / sqrt(3.0);
	&CreateMetalCIF($metals[$i], $Template, $a, $a, $a, $CIFFile);
}

for(my $i = 0 ; $i < @ABs ; $i++) {
	print("$i: $ABs[$i]\n");

	my ($atom1, $atom2) = ($ABs[$i] =~ /^([A-Z][a-z]?)([A-Z][a-z]?)$/);
	my $AtomNum1 = AtomType::GetAtomicNumber($atom1);
	my ($filepath, $AtomJName, $AtomEName,
		$AtomicMass, $FoundYear, $AtomicRadius1,
		$MeltingPoint, $BoilingPoint, $Density,
		$HeatCapacity, $IonizationPotential, $ElectronAffinity)
				 = AtomType::GetIonInfFromPeriodicTable2($AtomNum1);
	my $AtomNum2 = AtomType::GetAtomicNumber($atom2);
	my ($filepath, $AtomJName, $AtomEName,
		$AtomicMass, $FoundYear, $AtomicRadius2,
		$MeltingPoint, $BoilingPoint, $Density,
		$HeatCapacity, $IonizationPotential, $ElectronAffinity)
				 = AtomType::GetIonInfFromPeriodicTable2($AtomNum2);
	$AtomicRadius1 = 0.8 if($AtomicRadius1 == 0.0);
	$AtomicRadius2 = 0.8 if($AtomicRadius2 == 0.0);
	my $d = $AtomicRadius1 + $AtomicRadius2;

	my $key = "NaCl";
	my $Template = $MetalCIFTemplate;
	$Template =~ s/\{Structure\}/$key/;
	$CIFFile = "$ABs[$i]-$key-AB.cif";
	my $a = $d * 2.0;
	&CreateABCIF($ABs[$i], $Template, $a, $a, $a, $CIFFile);

	my $key = "CsCl";
	my $Template = $MetalCIFTemplate;
	$Template =~ s/\{Structure\}/$key/;
	$CIFFile = "$ABs[$i]-$key-AB.cif";
	my $a = $d * 2.0 / sqrt(3);
	&CreateABCIF($ABs[$i], $Template, $a, $a, $a, $CIFFile);

	my $key = "NiAs";
	my $Template = $MetalCIFTemplate;
	$Template =~ s/\{Structure\}/$key/;
	$CIFFile = "$ABs[$i]-$key-AB.cif";
	my $a = $d * 1.48;
	my $c = $d * 2.06;
	&CreateABCIF($ABs[$i], $Template, $a, $a, $c, $CIFFile);

	my $key = "ZB";
	my $Template = $MetalCIFTemplate;
	$Template =~ s/\{Structure\}/$key/;
	$CIFFile = "$ABs[$i]-$key-AB.cif";
	my $a = $d * 4.0 / sqrt(3.0);
	&CreateABCIF($ABs[$i], $Template, $a, $a, $a, $CIFFile);

	my $key = "WZ";
	my $Template = $MetalCIFTemplate;
	$Template =~ s/\{Structure\}/$key/;
	$CIFFile = "$ABs[$i]-$key-AB.cif";
	my $a = $d * 1.64;
	my $c = $d * 2.66;
	&CreateABCIF($ABs[$i], $Template, $a, $a, $c, $CIFFile);
}

exit;

sub CreateABCIF
{
	my ($atom, $Template, $alat, $blat, $clat, $path) = @_;
	
	my ($atom1, $atom2) = ($atom =~ /^([A-Z][a-z]?)([A-Z][a-z]?)$/);
	if(!$atom2) {
		print "Error: Compound [$atom] is not AB type.\n";
		return 0;
	}

	my $in = new JFile;
	$in->Open($Template, "r") or die "$!: Can not read [$Template].\n";
	my $out = new JFile;
	$out->Open($path, "w") or die "$!: Can not write to [$path].\n";
	while(1) {
		my $line = $in->ReadLine();
		last if(!defined $line);
		
		$line =~ s/\{atom1\}/$atom1/g;
		$line =~ s/\{atom2\}/$atom2/g;
		$line =~ s/\{a\}/$alat/g;
		$line =~ s/\{b\}/$blat/g;
		$line =~ s/\{c\}/$clat/g;
		
		$out->print($line);
	}
	$out->Close();
	$in->Close();
}

sub CreateMetalCIF
{
	my ($atom, $Template, $alat, $blat, $clat, $path) = @_;
	
	my $in = new JFile;
	$in->Open($Template, "r") or die "$!: Can not read [$Template].\n";
	my $out = new JFile;
	$out->Open($path, "w") or die "$!: Can not write to [$path].\n";
	while(1) {
		my $line = $in->ReadLine();
		last if(!defined $line);
		
		$line =~ s/\{atom\}/$atom/g;
		$line =~ s/\{a\}/$alat/g;
		$line =~ s/\{b\}/$blat/g;
		$line =~ s/\{c\}/$clat/g;
		
		$out->print($line);
	}
	$out->Close();
	$in->Close();
}

sub CreateMoleculeCIF
{
	my ($molecule, $a, $path) = @_;

	my $crystal = new Crystal;
	$crystal->SetLatticeParameters($a, $a, $a, 90.0, 90.0, 90.0);
	$crystal->SetSpaceGroup("P 1");

	my $d = 3.0;
	my $dx = $d / 2.0 / $a;
	my (@a) = Utils::FindAll("[A-Z][a-z]*\\d*", $molecule);
	if($a[0] =~ /^([A-Z][a-z]?)(\d+)$/) {
		my $atom = $1;
		my $AtomNum = AtomType::GetAtomicNumber($atom);
		my ($filepath, $AtomJName, $AtomEName,
			$AtomicMass, $FoundYear, $AtomicRadius,
			$MeltingPoint, $BoilingPoint, $Density,
			$HeatCapacity, $IonizationPotential, $ElectronAffinity)
				 = AtomType::GetIonInfFromPeriodicTable2($AtomNum);
		if($AtomicRadius <= 0) {
			print "Error: Could not get Atomic Radius(1) for [$atom].\n";
		}
		my $dx = $AtomicRadius / $a;

		$crystal->AddAtomSite("${atom}1", $atom, 0.5-$dx, 0.5, 0.5);
		$crystal->AddAtomSite("${atom}2", $atom, 0.5+$dx, 0.5, 0.5);
	}
	elsif($molecule =~ /^([A-Z][a-z]?)2([A-Z][a-z]?)$/) {
		my ($atom0, $atom1) = ($1, $2);
		my $AtomNum0 = AtomType::GetAtomicNumber($atom0);
		my $AtomNum1 = AtomType::GetAtomicNumber($atom1);
		my ($filepath, $AtomJName, $AtomEName,
			$AtomicMass, $FoundYear, $AtomicRadius0,
			$MeltingPoint, $BoilingPoint, $Density,
			$HeatCapacity, $IonizationPotential, $ElectronAffinity)
				 = AtomType::GetIonInfFromPeriodicTable2($AtomNum0);
		my ($filepath, $AtomJName, $AtomEName,
			$AtomicMass, $FoundYear, $AtomicRadius1,
			$MeltingPoint, $BoilingPoint, $Density,
			$HeatCapacity, $IonizationPotential, $ElectronAffinity)
				 = AtomType::GetIonInfFromPeriodicTable2($AtomNum1);
		if($AtomicRadius0 <= 0) {
			print "Error: Could not get Atomic Radius(1) for [$atom0].\n";
		}
		if($AtomicRadius1 <= 0) {
			print "Error: Could not get Atomic Radius(1) for [$atom1].\n";
		}
		my $dx0 = $AtomicRadius0 / $a;
		my $dx1 = $AtomicRadius1 / $a;
		$dx0 += $dx1;

		$crystal->AddAtomSite("${atom1}2", $atom1, 0.5, 0.5, 0.5);
		$crystal->AddAtomSite("${atom0}1", $atom0, 0.5-$dx0*0.9, 0.5, 0.5+$dx0*0.1);
		$crystal->AddAtomSite("${atom0}2", $atom0, 0.5+$dx0*0.9, 0.5, 0.5+$dx0*0.1);
	}
	elsif($molecule =~ /^([A-Z][a-z]?)([A-Z][a-z]?)2$/) {
		my ($atom0, $atom1) = ($1, $2);
		my $AtomNum0 = AtomType::GetAtomicNumber($atom0);
		my $AtomNum1 = AtomType::GetAtomicNumber($atom1);
		my ($filepath, $AtomJName, $AtomEName,
			$AtomicMass, $FoundYear, $AtomicRadius0,
			$MeltingPoint, $BoilingPoint, $Density,
			$HeatCapacity, $IonizationPotential, $ElectronAffinity)
				 = AtomType::GetIonInfFromPeriodicTable2($AtomNum0);
		my ($filepath, $AtomJName, $AtomEName,
			$AtomicMass, $FoundYear, $AtomicRadius1,
			$MeltingPoint, $BoilingPoint, $Density,
			$HeatCapacity, $IonizationPotential, $ElectronAffinity)
				 = AtomType::GetIonInfFromPeriodicTable2($AtomNum1);
		if($AtomicRadius0 <= 0) {
			print "Error: Could not get Atomic Radius(1) for [$atom0].\n";
		}
		if($AtomicRadius1 <= 0) {
			print "Error: Could not get Atomic Radius(1) for [$atom1].\n";
		}
		my $dx0 = $AtomicRadius0 / $a;
		my $dx1 = $AtomicRadius1 / $a;
		$dx0 += $dx1;

		$crystal->AddAtomSite("${atom1}2", $atom0, 0.5, 0.5, 0.5);
		$crystal->AddAtomSite("${atom0}1", $atom1, 0.5-$dx0*0.9, 0.5, 0.5+$dx0*0.1);
		$crystal->AddAtomSite("${atom0}2", $atom1, 0.5+$dx0*0.9, 0.5, 0.5+$dx0*0.1);
	}
	elsif($molecule =~ /^([A-Z][a-z]?)([A-Z][a-z]?)3$/) {
		my ($atom0, $atom1) = ($1, $2);
		my $AtomNum0 = AtomType::GetAtomicNumber($atom0);
		my $AtomNum1 = AtomType::GetAtomicNumber($atom1);
		my ($filepath, $AtomJName, $AtomEName,
			$AtomicMass, $FoundYear, $AtomicRadius0,
			$MeltingPoint, $BoilingPoint, $Density,
			$HeatCapacity, $IonizationPotential, $ElectronAffinity)
				 = AtomType::GetIonInfFromPeriodicTable2($AtomNum0);
		my ($filepath, $AtomJName, $AtomEName,
			$AtomicMass, $FoundYear, $AtomicRadius1,
			$MeltingPoint, $BoilingPoint, $Density,
			$HeatCapacity, $IonizationPotential, $ElectronAffinity)
				 = AtomType::GetIonInfFromPeriodicTable2($AtomNum1);
		if($AtomicRadius0 <= 0) {
			print "Error: Could not get Atomic Radius(1) for [$atom0].\n";
		}
		if($AtomicRadius1 <= 0) {
			print "Error: Could not get Atomic Radius(1) for [$atom1].\n";
		}
		my $dx0 = $AtomicRadius0 / $a;
		my $dx1 = $AtomicRadius1 / $a;
		$dx0 += $dx1;

		$crystal->AddAtomSite("${atom0}1", $atom0, 0.5, 0.5, 0.5+$dx0*0.2);
		$crystal->AddAtomSite("${atom1}1", $atom1, 0.5+$dx0*0.7, 0.5, 0.5);
		$crystal->AddAtomSite("${atom1}2", $atom1, 0.5-$dx0*0.5, 0.5+$dx0*0.6, 0.5);
		$crystal->AddAtomSite("${atom1}3", $atom1, 0.5-$dx0*0.5, 0.5-$dx0*0.6, 0.5);
	}
	elsif($molecule =~ /^([A-Z][a-z]?)([A-Z][a-z]?)4$/) {
		my ($atom0, $atom1) = ($1, $2);
		my $AtomNum0 = AtomType::GetAtomicNumber($atom0);
		my $AtomNum1 = AtomType::GetAtomicNumber($atom1);
		my ($filepath, $AtomJName, $AtomEName,
			$AtomicMass, $FoundYear, $AtomicRadius0,
			$MeltingPoint, $BoilingPoint, $Density,
			$HeatCapacity, $IonizationPotential, $ElectronAffinity)
				 = AtomType::GetIonInfFromPeriodicTable2($AtomNum0);
		my ($filepath, $AtomJName, $AtomEName,
			$AtomicMass, $FoundYear, $AtomicRadius1,
			$MeltingPoint, $BoilingPoint, $Density,
			$HeatCapacity, $IonizationPotential, $ElectronAffinity)
				 = AtomType::GetIonInfFromPeriodicTable2($AtomNum1);
		if($AtomicRadius0 <= 0) {
			print "Error: Could not get Atomic Radius(1) for [$atom0].\n";
		}
		if($AtomicRadius1 <= 0) {
			print "Error: Could not get Atomic Radius(1) for [$atom1].\n";
		}
		my $dx0 = $AtomicRadius0 / $a;
		my $dx1 = $AtomicRadius1 / $a;
		$dx0 += $dx1;

		$crystal->AddAtomSite("${atom0}1", $atom0, 0.5, 0.5, 0.5);
		$crystal->AddAtomSite("${atom1}1", $atom1, 0.5+$dx0*0.8, 0.5,          0.5+$dx0*0.5);
		$crystal->AddAtomSite("${atom1}2", $atom1, 0.5,          0.5+$dx0*0.8, 0.5-$dx0*0.5);
		$crystal->AddAtomSite("${atom1}3", $atom1, 0.5-$dx0*0.8, 0.5,          0.5+$dx0*0.5);
		$crystal->AddAtomSite("${atom1}3", $atom1, 0.5,          0.5-$dx0*0.8, 0.5-$dx0*0.5);
	}
	elsif(@a == 2) {
		my $AtomNum0 = AtomType::GetAtomicNumber($a[0]);
		my $AtomNum1 = AtomType::GetAtomicNumber($a[1]);
		my ($filepath, $AtomJName, $AtomEName,
			$AtomicMass, $FoundYear, $AtomicRadius0,
			$MeltingPoint, $BoilingPoint, $Density,
			$HeatCapacity, $IonizationPotential, $ElectronAffinity)
				 = AtomType::GetIonInfFromPeriodicTable2($AtomNum0);
		my ($filepath, $AtomJName, $AtomEName,
			$AtomicMass, $FoundYear, $AtomicRadius1,
			$MeltingPoint, $BoilingPoint, $Density,
			$HeatCapacity, $IonizationPotential, $ElectronAffinity)
				 = AtomType::GetIonInfFromPeriodicTable2($AtomNum1);
		if($AtomicRadius0 <= 0) {
			print "Error: Could not get Atomic Radius(1) for [$a[0]].\n";
		}
		if($AtomicRadius1 <= 0) {
			print "Error: Could not get Atomic Radius(1) for [$a[1]].\n";
		}
		my $dx0 = $AtomicRadius0 / $a;
		my $dx1 = $AtomicRadius1 / $a;

		$crystal->AddAtomSite("$a[0]1", $a[0], 0.5-$dx0, 0.5, 0.5);
		$crystal->AddAtomSite("$a[1]1", $a[1], 0.5+$dx1, 0.5, 0.5);
	}
	else {
		$crystal->AddAtomSite("$a[0]", $a[0], 0.5, 0.5, 0.5);
	}
	$crystal->ExpandCoordinates();

	return CIF->new()->CreateCIFFileFromCCrystal($crystal, $path);
}