use strict;
BEGIN {
use lib 'c:/Programs/Perl/lib';
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 @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 $CreateAtomCIFs = 0;

my $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);
}

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);
	}
}
exit;

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);
}