#!/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 Deps;
use Utils;
use JFile;

use Crystal::CIF;
use Crystal::Crystal;
use Crystal::WIEN2k;
use Crystal::VASP;
use Sci::Algorism;

my $SGROUPPath   = ($^O =~ /^MSWin/)? "D:/Programs/Perl/WIEN2k/sgroup.exe" : "sgroup";

my $InputFile         = "Si-final.cif";
my ($drive, $dir, $filename, $ext, $lastdir, $filebody) = Deps::SplitFilePath($InputFile);
my $SampleName    = $filebody;

my $vasp = new VASP;

my $SimplifiedCIFFile = $InputFile;
$SimplifiedCIFFile =~ s/\.cif/-simple\.cif/i;
my $ConvertedFile    = $InputFile;
$ConvertedFile =~ s/\.cif/-converted\.cif/i;
my $StructPath    = $InputFile;
$StructPath       =~ s/\.cif/\.struct/i;
my $SymStructPath = $InputFile;
$SymStructPath    =~ s/\.cif/\-Symmetrized.struct/i;
my $SymCIFPath    = $InputFile;
$SymCIFPath       =~ s/\.cif/\-Symmetrized.cif/i;

my $InputVASPFile     = "PARCHG.0008.0001.HOMO.vasp";
my $ConvertedVASPFile = $InputVASPFile;
$ConvertedVASPFile =~ s/\.vasp/\-converted.vasp/;
print "Input VASP File    : $InputVASPFile\n";
print "Converted VASP File: $ConvertedVASPFile\n";

print "<H1>Lattice conversion</H1>\n";

my $CheckConsistency = 1; #&GetArg('CheckConsistency');
my $ParametersSelect = "ai"; #ai, xyz, a*i, hkl
my $DirectionSelect  = "OriginalToConverted"; #OriginalToConverted, ConvertedlToOriginal
my $PresetRule = "PrimFCC"; #RhombHex, HexRhomb, HexOrtho, FCCPrim, BCCPrim, ACenterPrim, BCenterPrim, CCenterPrim
print "PresetRule: $PresetRule\n";
print "Direction: $DirectionSelect\n";
print "Check consistency:$CheckConsistency\n";

my $CIF = new CIF();
unless($CIF->Read($InputFile)) {
	print "<H2>Error: Can not read $InputFile.</H2>\n";
	return -2;
}
if($CIF->WriteSimpleCIFFile($SimplifiedCIFFile)) {
}

my $Crystal = $CIF->GetCCrystal();
$Crystal->ExpandCoordinates();
my ($a,$b,$c,$alpha,$beta,$gamma) = $Crystal->LatticeParameters();
print "Lattice parameters:\n";
print "a=$a b=$b c=$c  alpha=$alpha beta=$beta gamma=$gamma\n";
print "Volume of original lattice: ", $Crystal->Volume(), " A<sup>3</sup>\n";
print "Density of original lattice: ", $Crystal->Density(), " g/cm<sup>3</sup>\n";
print "\n";

my ($g11, $g12, $g13, $g21, $g22, $g23, $g31, $g32, $g33) = $Crystal->Metrics();
print "Real space metrics conversion:\n";
printf "     |%10.4f %10.4f %10.4f|\n", $g11, $g12, $g13;
printf "     |%10.4f %10.4f %10.4f|\n", $g21, $g22, $g23;
printf "     |%10.4f %10.4f %10.4f|\n", $g31, $g32, $g33;

#基本ベクトルの変換行列
my $sxyz = new Math::MatrixReal(3, 1);
my $shkl = new Math::MatrixReal(3, 1);
	my $a     = '';#0; #&GetArg('a'); $a = eval($a);
	my $b     = 0; #&GetArg('b'); $b = eval($b);
	my $c     = 0; #&GetArg('c'); $c = eval($c);
	my $alpha = 0; #&GetArg('alpha'); $alpha = eval($alpha);
	my $beta  = 0; #&GetArg('beta');  $beta  = eval($beta);
	my $gamma = 0; #&GetArg('gamma'); $gamma = eval($gamma);

	my $v = 0; #&GetArg('sx'); $v = eval($v);
	$sxyz->assign(1, 1, $v);
	$v = 0; #&GetArg('sy'); $v = eval($v);
	$sxyz->assign(2, 1, $v);
	$v = 0; #&GetArg('sz'); $v = eval($v);
	$sxyz->assign(3, 1, $v);

	$v = 0; #&GetArg('sh'); $v = eval($v);
	$shkl->assign(1, 1, $v);
	$v = 0; #&GetArg('sk'); $v = eval($v);
	$shkl->assign(2, 1, $v);
	$v = 0; #&GetArg('sl'); $v = eval($v);
	$shkl->assign(3, 1, $v);

my ($sT, $T, $tRT, $RT, $tR) = Crystal->new()->GetLatticeConversionMatrix($PresetRule, $DirectionSelect, $ParametersSelect);
&PrintMatrix("Input matrix: (tij)\n", $sT);
#実空間べクトル(ai)の変換行列
&PrintMatrix("Real space vector: (a'i) = (Tij)(ai)\n", $T);
#実空間座標(x,y,z)の変換行列
&PrintMatrix("Real space coordinate conversion matrix (x') = ((Tij)t)-1(x)\n", $tRT);
#逆空間ベクトル(a*i)の変換行列: $tRT
&PrintMatrix("Reciprocal space vector conversion matrix (a*') = ((Tij)t)-1(a*)\n", $RT);
#逆空間座標(h k l)の変換行列: $T
&PrintMatrix("Reciprocal space coordinate conversion matrix (h') = (Tij)(h)\n", $tR);

my ($NewCrystal, $SymCrystal) = &LatticeConversion($sT, $T, $tRT, $RT, $tR, $Crystal);
$NewCrystal->Symmetrize(1.0e-3);

my $pDensity = $vasp->ReadDensityFileToHash($InputVASPFile);
print "Title: [$pDensity->{Title}]\n";

my ($a11,$a12,$a13,$a21,$a22,$a23,$a31,$a32,$a33) = $NewCrystal->LatticeVectors();
$pDensity->{Base} = $a11;
($pDensity->{a11}, $pDensity->{a12}, $pDensity->{a13}) = ($a11/$a11, $a12/$a11, $a13/$a11);
($pDensity->{a21}, $pDensity->{a22}, $pDensity->{a23}) = ($a21/$a11, $a22/$a11, $a23/$a11);
($pDensity->{a31}, $pDensity->{a32}, $pDensity->{a33}) = ($a31/$a11, $a32/$a11, $a33/$a11);

my @AtomType = $NewCrystal->GetCAtomTypeList();
my @site = $NewCrystal->GetCAsymmetricAtomSiteList();
my %nAtoms;
for(my $i = 0 ; $i < @site ; $i++) {
	my $type = lc $site[$i]->AtomNameOnly();
	$nAtoms{$type}++;
#print "i=$i: $type\n";
}
my @nAtoms;
my @pos;
my $ic = 0;
for(my $i = 0 ; $i < @AtomType ; $i++) {
	my $name = lc $AtomType[$i]->AtomNameOnly();
	$nAtoms[$i] = $nAtoms{$name};
#print "i=$i: $name\n";
	for(my $is = 0 ; $is < @site ; $is++) {
		my $name1 = lc $site[$is]->AtomNameOnly();
		next if($name ne $name1);

		my ($x, $y, $z) = $site[$is]->Position();
		$pos[$ic] = [$x, $y, $z];
		$ic++;
	}
}
$pDensity->{pnAtoms}  = \@nAtoms;
$pDensity->{pAtomPos} = \@pos;

my $Tinv = $tRT->inverse();
&PrintMatrix("FCC to Primitive coordinate matrix:\n", $Tinv);
&PrintConvertedCoordinate($Tinv,   0,   0,   0);
&PrintConvertedCoordinate($Tinv, 0.9,   0,   0);
&PrintConvertedCoordinate($Tinv, 0.9, 0.9,   0);
&PrintConvertedCoordinate($Tinv, 0.9, 0.9, 0.9);
&PrintConvertedCoordinate($Tinv, 0.5, 0.5, 0.5);

my $pPrimValues = $pDensity->{pValues};
my @ConvValues;
my $v = new Math::MatrixReal(3, 1);
for(my $ix = 0 ; $ix < $pDensity->{nx} ; $ix++) {
	for(my $iy = 0 ; $iy < $pDensity->{ny} ; $iy++) {
		for(my $iz = 0 ; $iz < $pDensity->{nz} ; $iz++) {
			my $x = $ix / $pDensity->{nx};
			my $y = $iy / $pDensity->{ny};
			my $z = $iz / $pDensity->{nz};
			$v->assign(1, 1, $x);
			$v->assign(2, 1, $y);
			$v->assign(3, 1, $z);
			my $v2 = $T * $v;
			my ($x2, $y2, $z2) = ($v2->element(1,1), $v2->element(2,1), $v2->element(3,1));
my $v = Algorism::InterpolateValuesInCube($pDensity->{nx}, $pDensity->{ny},$pDensity->{nz}, $pPrimValues, $x2, $y2, $z2);
			$ConvValues[$ix][$iy][$iz] = $v;
#			printf "$ix,$iy,$iz: (%g,%g,%g) => (%g,%g,%g)\n", $x, $y, $z, $x2, $y2, $z2;
#			print "($ix,$iy,$iz) => ($ix2,$iy2,$iz2)\n";
		}
	}
}
$pDensity->{pValues} = \@ConvValues;

$vasp->SaveDensityFileFromHash($ConvertedVASPFile, $pDensity);

exit;

sub PrintConvertedCoordinate
{
	my ($T, $x, $y, $z) = @_;
	my $v = new Math::MatrixReal(3, 1);
	$v->assign(1, 1, $x);
	$v->assign(2, 1, $y);
	$v->assign(3, 1, $z);
	my $v2 = $T * $v;
	my ($x2, $y2, $z2) = ($v2->element(1,1), $v2->element(2,1), $v2->element(3,1));
	$x2 = Utils::Reduce01($x2);
	$y2 = Utils::Reduce01($y2);
	$z2 = Utils::Reduce01($z2);
	printf "(%g,%g,%g) => (%g,%g,%g)\n", $v->element(1,1),  $v->element(2,1),  $v->element(3,1), 
										$x2, $y2, $z2;
}

sub PrintMatrix
{
	my ($message, $T) = @_;
	print $message if($message ne '');
	printf "     |%10.4f %10.4f %10.4f|\n", $T->element(1,1), $T->element(1,2), $T->element(1,3);
	printf "     |%10.4f %10.4f %10.4f|\n", $T->element(2,1), $T->element(2,2), $T->element(2,3);
	printf "     |%10.4f %10.4f %10.4f|\n", $T->element(3,1), $T->element(3,2), $T->element(3,3);
}

sub LatticeConversion
{
	my ($sT, $T, $tRT, $RT, $tR, $Crystal) = @_;

#座標(sxyz)の変換: tRT
	if($sxyz->element(1,1) ne '') {
		my $a1 = 
		my $txyz = $tRT * $sxyz;
		print "Real space coordinate conversion:\n";
		printf "(%10.4f, %10.4f, %10.4f) => (%10.4f, %10.4f, %10.4f)\n",
			$sxyz->element(1,1), $sxyz->element(2,1), $sxyz->element(3,1),
			$txyz->element(1,1), $txyz->element(2,1), $txyz->element(3,1);
		print "\n";
	}

#逆空間座標(shkl)の変換: tR
	if($shkl->element(1,1) ne '') {
		my $thkl = $tR * $shkl;
		print "Reciprocal space coordinate conversion:\n";
		printf "(%10.4f, %10.4f, %10.4f) => (%10.4f, %10.4f, %10.4f)\n",
			$shkl->element(1,1), $shkl->element(2,1), $shkl->element(3,1),
			$thkl->element(1,1), $thkl->element(2,1), $thkl->element(3,1);
		print "\n";
	}

#格子定数の変換
	if($a ne '') {
		print "Original lattice parameters:\n";
		printf "a=%12.6f b=%12.6f c=%12.6f  alpha=%8.4f beta=%8.4f gamma=%8.4f\n",
			$a, $b, $c, $alpha, $beta, $gamma;
		my $Crystal = new Crystal;
		$Crystal->SetLatticeParameters($a,$b,$c,$alpha,$beta,$gamma);

		my $NewCrystal = $Crystal->ConvertLattice($T, 1);

		my ($a2, $b2, $c2, $alpha2, $beta2, $gamma2) = $NewCrystal->LatticeParameters();
		print "Converted lattice parameters:\n";
		printf "a'=%12.6f b'=%12.6f c'=%12.6f  alpha'=%8.4f beta'=%8.4f gamma'=%8.4f\n",
			$a2, $b2, $c2, $alpha2, $beta2, $gamma2;
		print "\n";

		my ($g11, $g12, $g13, $g21, $g22, $g23, $g31, $g32, $g33) = $Crystal->Metrics();

		print "Real space metrics conversion:\n";
		printf "     |%10.4f %10.4f %10.4f|\n", $g11, $g12, $g13;
		printf "     |%10.4f %10.4f %10.4f|\n", $g21, $g22, $g23;
		printf "     |%10.4f %10.4f %10.4f|\n", $g31, $g32, $g33;

		($g11, $g12, $g13, $g21, $g22, $g23, $g31, $g32, $g33) = $NewCrystal->Metrics();
		printf "     |%10.4f %10.4f %10.4f|\n", $g11, $g12, $g13;
		printf "     |%10.4f %10.4f %10.4f|\n", $g21, $g22, $g23;
		printf "     |%10.4f %10.4f %10.4f|\n", $g31, $g32, $g33;
	}

	my $NewCrystal = $Crystal->ConvertLattice($T, 0, $CheckConsistency);
	my ($a2,$b2,$c2,$alpha2,$beta2,$gamma2) = $NewCrystal->LatticeParameters();
	($g11, $g12, $g13, $g21, $g22, $g23, $g31, $g32, $g33) = $NewCrystal->Metrics();
	printf "     |%10.4f %10.4f %10.4f|\n", $g11, $g12, $g13;
	printf "     |%10.4f %10.4f %10.4f|\n", $g21, $g22, $g23;
	printf "     |%10.4f %10.4f %10.4f|\n", $g31, $g32, $g33;

	print "Converted Lattice parameters:\n";
	print "a'=$a2 b'=$b2 c'=$c2  alpha'=$alpha2 beta'=$beta2 gamma'=$gamma2\n";
	print "Volume of converted lattice: ", $NewCrystal->Volume(), " A3\n";
	print "Density of converted lattice: ", $NewCrystal->Density(), " g/cm3\n";

	my $NewCIF = new CIF();
	$NewCIF->SetFileName($ConvertedFile);
	$NewCIF->SetCrystalName($Crystal->CrystalName());
	$NewCIF->SetLatticeParameters($a2, $b2, $c2, $alpha2, $beta2, $gamma2);
	$NewCIF->SetSpaceGroup("P 1", 1);
	$NewCIF->SetVolume($NewCrystal->Volume());
	$NewCIF->AddSymmetryOperation("x,y,z");
	my @ExpandedAtomSiteList = $NewCrystal->GetCExpandedAtomSiteList();
	my $nTotalExpandedAtomSite = $NewCrystal->nTotalExpandedAtomSite();
	for(my $i = 0 ; $i < $nTotalExpandedAtomSite ; $i++) {
		my $atom      = $ExpandedAtomSiteList[$i];
		my $label     = $atom->Label();
		my $atomname  = $atom->AtomNameOnly();
		my $charge    = $atom->Charge();
		my ($x,$y,$z) = $atom->Position();
		my $occ       = $atom->Occupancy();

		$NewCIF->AddAsymmetricAtomSite($label, $atomname, $x, $y, $z,$occ);
	}
	$NewCIF->FillCIFData();

	unless($NewCIF->WriteSimpleCIFFile($ConvertedFile)) {
		print "Error: Could not write to ($ConvertedFile).\n";
		return -1;
	}

	if(abs($Crystal->Density() - $NewCrystal->Density()) / $Crystal->Density() > 0.05) {
		print "<H1>Causion!!.</H1>\n";
		print "<H2>Densities of the original and converted cells are different.</H2>\n";
		print "<H2>The conversion rule employed would not be compatible with the original cell.</H2>\n";
		die;
	}

	if($PresetRule eq 'HexOrtho' or 
	   $PresetRule eq 'FCCPrim' or
	   $PresetRule eq 'BCCPrim' or
	   $PresetRule eq 'PrimFCC' or
	   $PresetRule eq 'PrimBCC'
	   ) {
		print "<H2>Find symmetry for the converted lattice</H2>";

		my ($a0,$b0,$c0,$alpha0,$beta0,$gamma0) = $NewCrystal->LatticeParameters();
		if($PresetRule eq 'HexOrtho') {
			$NewCrystal->SetLatticeParameters($a0, $a0*1.3, $c0,$alpha0,$beta0,$gamma0);
		}
		elsif($PresetRule eq 'FCCPrim' or $PresetRule eq 'BCCPrim') {
			$NewCrystal->SetLatticeParameters($a0,$b0, $c0, 65,65,65);
		}

		my $WIEN2k2 = new WIEN2k;
		$WIEN2k2->SaveSymmetrizedStructFileFromCrystal($NewCrystal, $StructPath, $SymStructPath, $SGROUPPath);#, $tol, $IsPrint);
		my $SymCrystal = $WIEN2k2->ReadStructFile($SymStructPath, 0);
		unless($SymCrystal) {
			print("Error: Can not read [$SymStructPath].\n");
			return 0;
		}
#		$SymCrystal->SetLatticeSystem('', 1, 1.0e-3);
		$SymCrystal->Symmetrize(1.0e-3);

		my $SPG = $SymCrystal->GetCSpaceGroup();
		my $SPGName       = $SPG->SPGName();
		my $LatticeSystem = $SPG->LatticeSystem();
print "SPGName in Symmetrized Struct file: $SPGName\n";
print "LatticeSystem in Symmetrized Struct file: $LatticeSystem\n";
		if($SPGName =~ /^R/i) {
			$SPG->SetSPGName("$SPGName R");
			my ($a,$b,$c,$alpha,$beta,$gamma) = $SymCrystal->LatticeParameters();
			print("  Hexagonal axis: a=$a   c=$c\n");
			my $ar = sqrt(3.0*$a*$a + $c*$c) / 3.0;
			my $alphah = 1.5 / sqrt(3.0 + ($c/$a)*($c/$a));
			$alphah = 2.0 * Sci::asin($alphah) * Sci::todeg();
			$alphah = $alpha0 if(abs($alphah -65.0) < 1.0e-3);
			print("  Rhombohedral axis: a=$ar   alpha=$alphah\n");
			$SymCrystal->SetLatticeParameters($ar,$ar,$ar,$alphah,$alphah,$alphah);
		}
		elsif($SPGName =~ /^[PFI].*3/i) {
			$SPG->SetSPGName("$SPGName R");
			my ($a,$b,$c,$alpha,$beta,$gamma) = $SymCrystal->LatticeParameters();
			my $ar = ($a+$b+$c) / 3.0;
			$SymCrystal->SetLatticeParameters($ar,$ar,$ar,90,90,90);
		}

		my $CIF2 = new CIF;
		if(!$CIF2->CreateCIFFileFromCCrystal($SymCrystal, $SymCIFPath, 0, 'unix')) {
			print("Error: Can not create [$SymCIFPath].\n");
			return 0;
		}
	}

	return ($NewCrystal, $SymCrystal);
}