#!/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 $InputFile  = ($ARGV[0])? $ARGV[0] : "Si-final.cif";
if(!-f $InputFile) {
	my @cif = glob("*.cif");
	$InputFile = $cif[0];
}
my $InputVASPFile  = ($ARGV[1])? $ARGV[1] : "PARCHG.0008.0001.HOMO.vasp";
if(!-f $InputFile) {
	my @vaspfiles = glob("*.vasp");
	$InputVASPFile = $vaspfiles[0];
}

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

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 $ConvertedVASPFile = $InputVASPFile;
$ConvertedVASPFile    =~ s/\.vasp/\-converted.vasp/;

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 "Lattice conversion\n";
print "Input VASP File    : $InputVASPFile\n";
print "Converted VASP File: $ConvertedVASPFile\n";
print "PresetRule: $PresetRule\n";
print "Direction: $DirectionSelect\n";
print "Check consistency:$CheckConsistency\n";
print "\n";

my $CIF = new CIF();
unless($CIF->Read($InputFile)) {
	print "Error: Can not read $InputFile.\n";
	exit -2;
}
unless($CIF->WriteSimpleCIFFile($SimplifiedCIFFile)) {
	print "Error: Can not write to $SimplifiedCIFFile.\n";
	exit -2;
}

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(), " A3\n";
print "Density of original lattice: ", $Crystal->Density(), " g/cm3\n";
print "\n";

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

$vasp->LatticeConversionForDensityFile($InputVASPFile, $ConvertedVASPFile, $NewCrystal, $T, $tRT);

exit;

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

	my $NewCrystal = $Crystal->ConvertLattice($T, 0, $CheckConsistency);
	$NewCrystal->Symmetrize(1.0e-3);
	my ($a2,$b2,$c2,$alpha2,$beta2,$gamma2) = $NewCrystal->LatticeParameters();
	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";
	print "\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 $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->Symmetrize(1.0e-3);

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

	return ($NewCrystal, $SymCrystal);
}

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