#!/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 Utils;
use MyApplication;
use Crystal::WIEN2k;

#===============================================
# Applicationオブジェクト作成
#===============================================
my $App = new MyApplication;
exit if($App->Initialize() < 0);
$App->SetDeleteHTMLFlag(1);

#===============================================
# スクリプト大域変数
#===============================================
my $XCrySDenPath = "xcrysden";
# for klist
my $K   = 1;
# for klist & bxsf
my $div             = 60;
my ($nx, $ny, $nz)  = (20, 20, 20);
my ($emin, $emax)   = (-1.5, 1.5);
my $UseSymmetry     = 1;

# for bxsf
my ($EFmin, $EFmax) = (-1.0, 1.0);
my $Eoffset         = 0.0;
my $OneByOne        = 1;

my @KIndex;

#==========================================
# コマンドラインオプション読み込み
#==========================================
$App->AddArgument("--Action",    "--Action=[all|MakeKList|MakeINSP|CalBands|MakeBXSF|ViewFS]", '');
$App->AddArgument("--DebugMode", "--DebugMode: Set DebugMode", '');
exit 1 if($App->ReadArgs(0) != 1);
my %Args = $App->Args()->GetArgHash();

#==========================================
# メイン関数スタート
#==========================================
my $Debug = $Args{DebugMode};
$App->SetDebug($Debug);
my $Action  = $Args{Action};
my $DataDir = $App->Args()->GetGetArg(0);
$DataDir = "." if($DataDir eq '');
$DataDir = Utils::ReduceDirectory($DataDir);
my $StructFile = &GetStructPath($DataDir);
if(!-e $StructFile) {
	$App->print("Error: StructFile [$StructFile] could not obtained.\n");
	exit 0;
}
my $Title     = &GetTitle($StructFile);
my $KListPath = Deps::ReplaceExtension($StructFile, ".klist");
my $BXSFPath  = Deps::ReplaceExtension($StructFile, ".bxsf");
my $PrefPath  = Deps::MakePath($DataDir, ".prefs-scf");

my $WIEN2k = new WIEN2k;
my $Crystal = $WIEN2k->ReadStructFile($StructFile, 0);
unless($Crystal) {
	$App->print("Error: Can not read [$StructFile].\n");
	return 0;
}

my $config = new IniFile();
$config->ReadAll($PrefPath);

$App->print("Struct: $StructFile\n");
$App->print("Lattice System: ", $Crystal->LatticeSystem(), "\n");
$App->print("Space Group   : ", $Crystal->SPGName(), "\n");
if($config->{spin} eq '') {
	$App->print("Spin polarized: no\n");
}
else {
	$App->print("Spin polarized: yes\n");
}

my $ret = 0;
if($Action eq 'MakeKList') {
	&MakeKList($Crystal, $KListPath);
	$ret = 1;
}
elsif($Action eq 'MakeINSP') {
	if(!$WIEN2k->SaveINSPFile($StructFile)) {
		my $INSPPath = WIEN2k::GetINSPFile($StructFile);
		$App->print("Error: $!: Can not write to [$INSPPath].\n");
		exit 0;
	}
}
elsif($Action eq 'CalBands') {
	&CalBands($config);
}
elsif($Action eq 'MakeBXSF') {
	&MakeBXSF($Crystal, $config, $BXSFPath);
	$ret = 1;
}
elsif($Action eq 'ViewFS') {
	&ViewFS($BXSFPath);
}
elsif($Action eq 'all') {
	&MakeKList($Crystal, $KListPath);
	if(!$WIEN2k->SaveINSPFile($StructFile)) {
		my $INSPPath = WIEN2k::GetINSPFile($StructFile);
		$App->print("Error: $!: Can not write to [$INSPPath].\n");
		exit 0;
	}
	&CalBands($config);
	&MakeBXSF($Crystal, $config, $BXSFPath);
	&ViewFS($BXSFPath);
}
else {
	$App->print("Error: Invald Action: [$Action]\n");
	$App->Usage();
}

exit $ret;

#===============================================
# スクリプト終了
#===============================================

#==========================================
# &Subroutines
#==========================================
sub ViewFS
{
	my ($BXSFPath) = @_;
	my $cmd = "$XCrySDenPath --bxsf $BXSFPath";
	my $ret = Execute($cmd);
	if($ret) {
		$App->print("Error: Can not execute [$cmd] with ret=$ret.\n");
		exit 0;
	}
}

sub GetTitle
{
	my ($path) = @_;

	my ($drive, $directory, $filename, $ext, $lastdir, $filebody) = Deps::SplitFilePath($path);
	my $Title   = $filebody;
	$Title =~ s/\s/_/g;
	return $Title;
}

sub GetStructPath
{
	my ($DataDir) = @_;

	my ($drive, $directory, $filename, $ext, $lastdir, $filebody) = Deps::SplitFilePath($DataDir);
	my $StructFile = Deps::MakePath($DataDir, "$filename.struct", 0);
	if(!-e $StructFile) {
		my $fmask = Deps::MakePath($DataDir, "*.struct");
		my @files = glob($fmask);
		if($files[0] =~ /setrmt\.struct$/) {
			$StructFile = $files[1];
		}
		else {
			$StructFile = $files[0];
		}
		unless($StructFile) {
			print("Struct file is not found from Data Directory [$DataDir].\n");
			return undef;
		}
	}
	return $StructFile;
}

sub MakeKList
{
	my ($Crystal, $KListPath) = @_;

	$App->print("\nMake .klist file\n");
	$App->print("KList : $KListPath\n");

$App->print("Save to [$KListPath]\n");
	my $out = new JFile;
	if(!$out->Open($KListPath, "w")) {
		$App->print("Error: $!: Can not write to [$KListPath]\n");
		exit 0;
	}

	my $ra = $Crystal->GetMatrixConventionalToPrimitiveReciprocalCell();
if(1) {
	$App->print("Primitive RLattice vector: ($ra->[0][0], $ra->[0][1], $ra->[0][2])\n");
	$App->print("                           ($ra->[1][0], $ra->[1][1], $ra->[1][2])\n");
	$App->print("                           ($ra->[2][0], $ra->[2][1], $ra->[2][2])\n");
}

if($Action eq 'all') {
	open(CHECK, ">check.txt") or die "$!: Can not write to [check.txt].\n";
}
	my $count = 1;
	my $BandName = "";
$App->print("kx range: (0 - $nx)\n");
	for(my $ix = 0 ; $ix <= $nx ; $ix++) {
#		my $iy0 = ($UseSymmetry)? $Crystal->GetYRangeForIteration($ix) : 0;
#$App->print("ky range: $ix-($iy0 - $ny)\n");
		my $iy0 = 0;
		for(my $iy = $iy0 ; $iy <= $ny ; $iy++) {
#			my $iz0 = ($UseSymmetry)? $Crystal->GetZRangeForIteration($ix, $iy) : 0;
#$App->print("kz range: $ix-$iy-($iz0 - $nz)\n");
			my $iz0 = 0;
			for(my $iz = $iz0 ; $iz <= $nz ; $iz++) {
				if($iz == $iz0) {
					$BandName = sprintf("B%02d%02d", $ix+1, $iy+1);
				}
				else {
					$BandName = "";
				}
				my ($IsConverted, $cix, $ciy, $ciz) = $Crystal->ConvertLatticeIndexByLatticeSystem($ix, $iy, $iz);
				if($UseSymmetry and $IsConverted) {
					next;
				}

				my $px = $K * $div * $ix / $nx;
				my $py = $K * $div * $iy / $ny;
				my $pz = $K * $div * $iz / $nz;
				my $x = $ra->[0][0] * $px + $ra->[1][0] * $py + $ra->[2][0] * $pz;
				my $y = $ra->[0][1] * $px + $ra->[1][1] * $py + $ra->[2][1] * $pz;
				my $z = $ra->[0][2] * $px + $ra->[1][2] * $py + $ra->[2][2] * $pz;
#				my $x = ($ra->[0][0] + $ra->[1][0] + $ra->[2][0]) * $K * $div * $ix / $nx;
#				my $y = ($ra->[0][1] + $ra->[1][1] + $ra->[2][1]) * $K * $div * $iy / $ny;
#				my $z = ($ra->[0][2] + $ra->[1][2] + $ra->[2][2]) * $K * $div * $iz / $nz;
#				my $x = $K * $div * $ix / $nx;
#				my $y = $K * $div * $iy / $ny;
#				my $z = $K * $div * $iz / $nz;
				if($count == 1) {
					$out->printf("%5s       %3d  %3d  %3d  %3d%4.1f %4.1f %4.1f\n",
						$BandName, $x, $y, $z, $div, 1.0, $emin, $emax);
				}
				else {
					$out->printf("%5s       %3d  %3d  %3d  %3d%4.1f\n",
						$BandName, $x, $y, $z, $div, 1.0);
				}
if($Action eq 'all') {
	$KIndex[$count-1] = sprintf("%03d%03d%03d", $ix, $iy, $iz);
	printf CHECK "%05d: %2d %2d %2d\n", $count, $ix, $iy, $iz;
}
				$count++;
			}
		}
	}
if($Action eq 'all') {
	close CHECK;
}
	$out->printf("END\n\n");
	$out->Close();
}

sub Execute
{
	my ($cmd) = @_;
	$App->print("  Execute [$cmd]...\n");
	my $ret = system($cmd);
#print "$ret\n";
	if($ret) {
		$App->print("    Error: execute [$cmd] failed with ret=$ret\n");
	}
	return $ret;
}

sub CalBands
{
	my ($config) = @_;

	$App->print("\nCalculate bands\n");
	my $opt = "";
	if($config->{spin} eq '') {
		exit 0 if(Execute("x lapw1"));
		exit 0 if(Execute("x spaghetti"));
	}
	else {
		exit 0 if(Execute("x lapw1 -up"));
		exit 0 if(Execute("x lapw1 -dn"));
		exit 0 if(Execute("x spaghetti -up"));
		exit 0 if(Execute("x spaghetti -dn"));
	}
}

sub MakeBXSF
{
	my ($Crystal, $config, $OutFile) = @_;

	$App->print("\nMake Band XSF(BXSF) file\n");
	$App->print("BXSF Path: $BXSFPath\n");

	my $EF = ($EFmin + $EFmax) / 2.0;
	$App->print("Search EF range: $EFmin - $EFmax eV\n");
	$App->print("Average EF: $EF\n");

	my ($drive, $directory, $filename, $ext, $lastdir, $filebody) = Deps::SplitFilePath($OutFile);
	my @InFile;
	if($config->{spin} eq '') {
		$InFile[0]  = "$filebody.spaghetti_ene";
	}
	else {
		$InFile[0]  = "$filebody.spaghettiup_ene";
		$InFile[1]  = "$filebody.spaghettiup_ene";
	}
	if(!-e $InFile[0]) {
		$App->print("Error: Can not read [$InFile[0]]\n");
		exit;
	}

	my $T = new Math::MatrixReal(3, 3);
	my $matrix = $Crystal->GetMatrixConventionalToPrimitiveCell();
	$T = Math::MatrixReal->new_from_cols($matrix);
	$T->transpose($T);
	my $PrimCrystal = $Crystal->ConvertLattice($T, 0, 1);
#実空間べクトル(ai)の変換行列
if(0) {
	print "Conversion matrix for real space vector: (a'i) = (Tij)(ai)\n";
		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);
}
if(0) {
	my ($a11, $a12, $a13, $a21, $a22, $a23, $a31, $a32, $a33) = $Crystal->LatticeVectors();
	print "Lattice vector: ($a11, $a12, $a13)\n";
	print "                ($a21, $a22, $a23)\n";
	print "                ($a31, $a32, $a33)\n";
	($a11, $a12, $a13, $a21, $a22, $a23, $a31, $a32, $a33) = $Crystal->ReciprocalLatticeVectors();
	print "RLattice vector: ($a11, $a12, $a13)\n";
	print "                 ($a21, $a22, $a23)\n";
	print "                 ($a31, $a32, $a33)\n";

	my ($a11, $a12, $a13, $a21, $a22, $a23, $a31, $a32, $a33) = $PrimCrystal->LatticeVectors();
	print "Lattice vector: ($a11, $a12, $a13)\n";
	print "                ($a21, $a22, $a23)\n";
	print "                ($a31, $a32, $a33)\n";
	($a11, $a12, $a13, $a21, $a22, $a23, $a31, $a32, $a33) = $PrimCrystal->ReciprocalLatticeVectors();
	print "RLattice vector: ($a11, $a12, $a13)\n";
	print "                 ($a21, $a22, $a23)\n";
	print "                 ($a31, $a32, $a33)\n";
}

	my $iBand = 0;
	my $nBand = 0;
	my @BandEnergy;
	my @BandIndex;
	foreach my $InFile (@InFile) {
$App->print("  Read [$InFile].\n");
		open(IN, "$InFile")    or die "$!: Can not read [$InFile]\n";

		while(1) {
			$iBand++;
			my $line = <IN>;
			last if(!defined $line);
			$line = Utils::DelSpace($line);
			last if($line eq '');
$App->print("Index line: $line\n");
			my $count = 1;

			my ($emin, $emax) = (1.0e10, -1.0e10);
			my @Energy = ([], [], []);
			for(my $ix = 0 ; $ix <= $nx ; $ix++) {
				for(my $iy = 0 ; $iy <= $ny ; $iy++) {
					for(my $iz = 0 ; $iz <= $nz ; $iz++) {
						my ($IsConverted, $cix, $ciy, $ciz) = $Crystal->ConvertLatticeIndexByLatticeSystem($ix, $iy, $iz);
#print("$ix-$iy-$iz: $IsConverted, $cix, $ciy, $ciz\n");
						if($UseSymmetry and $IsConverted) {
#print "Conv: $ix,$iy,$iz <= $cix, $ciy, $ciz [$Energy[$cix][$ciy][$ciz]]\n";
							$Energy[$ix][$iy][$iz] = $Energy[$cix][$ciy][$ciz];
							next;
						}

						$line = Utils::DelSpace(<IN>);
						my ($kx, $ky, $kz, $dk, $e) = Utils::Split("\\s+", $line);
						if(!defined $e) {
							$App->print("Error: Wrong format\n");
							$App->print("  line: $line\n");
							exit;
						}
if($Action eq 'all') {
	my $index = sprintf("%03d%03d%03d", $ix, $iy, $iz);
	if($index ne $KIndex[$count-1]) {
		Utils::DelSpace($line);
		print "Mismatch Line: [$count] [$line]\n";
		print "   Write: [$KIndex[$count-1]]\n";
		print "   Read : [$index]\n";
exit;
	}
#print "OK: $count: $ix, $iy, $iz\n";
}

						$Energy[$ix][$iy][$iz] = $e;
						$emin = $e if($e < $emin);
						$emax = $e if($e > $emax);
						$count++;
					}
				}
			}
			next if($emin > $EFmax or $emax < $EFmin);

$App->print("  To be recorded: e=($emin - $emax)\n");
			$BandEnergy[$nBand] = \@Energy;
			$BandIndex[$nBand]  = $iBand;
			$nBand++;
		}
		close(IN);
	}

$App->print("Save to [$OutFile]\n");
	open(OUT, ">$OutFile") or die "$!: Can not write to [$OutFile]\n";

	my ($nx1, $ny1, $nz1) = ($nx+1, $ny+1, $nz+1);
	print OUT <<EOT;
 BEGIN_INFO
   #
   # This is a Band-XCRYSDEN-Structure-File
   # Produced by WIEN2kFS.pl
   #
   # Launch as: xcrysden --bxsf $OutFile
   #
   Fermi Energy: $EF
 END_INFO

 BEGIN_BLOCK_BANDGRID_3D
   $Title
   BEGIN_BANDGRID_3D_fermi
     $nBand
     $nx1 $ny1 $nz1
EOT
	printf(OUT "    %12.4g %12.4g %12.4g\n", 0.0, 0.0, 0.0);
	printf(OUT "    %12.4g %12.4g %12.4g\n", $PrimCrystal->{Ra11}, $PrimCrystal->{Ra12}, $PrimCrystal->{Ra13});
	printf(OUT "    %12.4g %12.4g %12.4g\n", $PrimCrystal->{Ra21}, $PrimCrystal->{Ra22}, $PrimCrystal->{Ra23});
	printf(OUT "    %12.4g %12.4g %12.4g\n", $PrimCrystal->{Ra31}, $PrimCrystal->{Ra32}, $PrimCrystal->{Ra33});

	for(my $ib = 0 ; $ib < $nBand ; $ib++) {
		printf(OUT "   BAND: %3d\n", $BandIndex[$ib]);
		my $ie = 0;
		for(my $ix = 0 ; $ix <= $nx ; $ix++) {
			for(my $iy = 0 ; $iy <= $ny ; $iy++) {
				for(my $iz = 0 ; $iz <= $nz ; $iz++) {
					printf(OUT "    %8.5f", $BandEnergy[$ib][$ix][$iy][$iz]);
					print(OUT "\n") if($OneByOne);
					$ie++;
					if($ie % 10 == 0) {
#					if($ie % 10 == 0 and $iz != $nz) {
						print(OUT "\n") if(!$OneByOne);
					}
				}
				print(OUT "\n");
			}
			if($OneByOne) {
			}
			else {
				print(OUT "\n\n");
			}
			$ie = 0;
		}
	}
	print OUT <<EOT;
   END_BANDGRID_3D
 END_BLOCK_BANDGRID_3D
EOT

	close(OUT);
}