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

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

#===============================================
# スクリプト大域変数
#===============================================
my $VASPPath     = "vasp";
my $SGROUPPath   = "sgroup";
my $XCrySDenPath = "xcrysden";
# for klist
my $K   = 1;
# for klist & bxsf
my ($nx, $ny, $nz)  = (30, 30, 30);
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|MakeKPOINTS|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 $VASP = new VASP;
my $INCARFile             = $VASP->GetINCARFileName($DataDir);
my $Title                 = &GetTitle($DataDir);
my $StructPath            = Deps::MakePath($DataDir, "$Title.struct");
my $SymmetrizedStructPath = Deps::MakePath($DataDir, "$Title.sym.struct");
my $KListPath             = $VASP->GetKPOINTSFileName($DataDir);
my $EIGENVALPath          = $VASP->GetEIGENVALFileName($DataDir);
my $BXSFPath              = Deps::MakePath($DataDir, "$Title.bxsf");

my $Crystal = $VASP->ReadStructureFromCARFiles($INCARFile, 1);
unless($Crystal) {
	$App->print("Error: Can not read [$INCARFile].\n");
	return 0;
}

$App->print("INCAR  : [$INCARFile]\n");
$App->print("KPOINTS: [$KListPath]\n");
$App->print("Title: [$Title]\n");
$App->print("Struct File converted from [$StructPath] to [$SymmetrizedStructPath]\n");

my $WIEN2k = new WIEN2k;
if(!$WIEN2k->SaveStructFile($Crystal, $StructPath)) {
	$App->print("Error: Can not write to [$StructPath].\n");
	exit 0;
}
unlink($SymmetrizedStructPath);
my $cmd = "$SGROUPPath -wi $StructPath -wo $SymmetrizedStructPath -set-TOL=0.001";
if(&Execute($cmd)) {
	$App->print("Error: Execute [$cmd] failed.\n");
	exit 0;
}
if(!-f $SymmetrizedStructPath) {
	$App->print("Error: Can not create [$SymmetrizedStructPath].\n");
	exit 0;
}

$Crystal = $WIEN2k->ReadStructFile($SymmetrizedStructPath, 0);
my $SPG = $Crystal->GetCSpaceGroup();
$SPG->SetLatticeSystem(undef, 1, 1.0e-3);
$App->print("Lattice System: ", $Crystal->LatticeSystem(), "\n");
$App->print("Space Group   : ", $Crystal->SPGName(), "\n");

my $ret = 0;
if($Action eq 'MakeKPOINTS') {
	&MakeKPOINTS($Crystal, $KListPath);
	$ret = 1;
}
elsif($Action eq 'CalBands') {
	&CalBands();
}
elsif($Action eq 'MakeBXSF') {
	&MakeBXSF($Crystal, $EIGENVALPath, $BXSFPath);
	$ret = 1;
}
elsif($Action eq 'ViewFS') {
	&ViewFS($BXSFPath);
}
elsif($Action eq 'all') {
	&MakeKPOINTS($Crystal, $KListPath);
	&CalBands();
	&MakeBXSF($Crystal, $EIGENVALPath, $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 MakeKPOINTS
{
	my ($Crystal, $KListPath) = @_;

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

if($Action eq 'all') {
	open(CHECK, ">check.txt") or die "$!: Can not write to [check.txt].\n";
}
	my $count = 0;
	my $BandName = "";
	my @lines;
$App->print("kx range: (0 - $nx)\n");
	for(my $ix = 0 ; $ix <= $nx ; $ix++) {
		for(my $iy = 0 ; $iy <= $ny ; $iy++) {
			my $BandNamePrinted = 0;
			for(my $iz = 0 ; $iz <= $nz ; $iz++) {
				my ($IsConverted, $cix, $ciy, $ciz) = $Crystal->ConvertLatticeIndexByLatticeSystem($ix, $iy, $iz);
				if($UseSymmetry and $IsConverted) {
					next;
				}

				my $x = 1.0 * $K * $ix / $nx;
				my $y = 1.0 * $K * $iy / $ny;
				my $z = 1.0 * $K * $iz / $nz;
				my $weight = 1.0;
				if(!$BandNamePrinted) {
					$BandNamePrinted = 1;
#					$out->printf("%8.4f %8.4f %8.4f %8.4f ! B%02d%02d\n", $x, $y, $z, $weight, $ix+1, $iy+1);
					$lines[$count] = sprintf("%8.4f %8.4f %8.4f %8.4f ! B%02d%02d\n", $x, $y, $z, $weight, $ix+1, $iy+1);
				}
				else {
#					$out->printf("%8.4f %8.4f %8.4f %8.4f\n", $x, $y, $z, $weight);
					$lines[$count] = sprintf("%8.4f %8.4f %8.4f %8.4f\n", $x, $y, $z, $weight);
				}
if($Action eq 'all') {
	$KIndex[$count] = sprintf("%03d%03d%03d", $ix, $iy, $iz);
	printf CHECK "%05d: %2d %2d %2d [%s]\n", $count, $ix, $iy, $iz, $KIndex[$count];
}
				$count++;
			}
		}
	}
if($Action eq 'all') {
	close CHECK;
}

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

	$out->print("$Title\n");
	$out->print("$count\n"); #($nx+1) * ($ny+1) * ($nz+1), "\n");
#	$out->print("Cartesian\n");
	$out->print("Reciprocal\n");

	for(my $i = 0 ; $i < $count ; $i++) {
		$out->print($lines[$i]);
	}
	$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");
	exit 0 if(Execute("$VASPPath"));
}

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

unlink($OutFile);
	$App->print("\nMake Band XSF(BXSF) file\n");
	my $OUTCAR = $VASP->GetOUTCARFileName($EIGENVALPath);
	$App->print("EIGENVAL Path: $EIGENVALPath\n");
	$App->print("BXSF Path    : $BXSFPath\n");

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

	my ($drive, $directory, $filename, $ext, $lastdir, $filebody) = Deps::SplitFilePath($OutFile);

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

	my $EF = $VASP->ReadFermiEnergy($OUTCAR);
#$EFmax += $EF;
#$EFmin += $EF;
#$EFavrg += $EF;
#$EF = 0.0;
print "EF in OUTCAR: $EF eV\n";

	my (@BandEnergy, @BandIndex);
	my $iBand = 0;
	my $nBand = 0;
$App->print("  Read [$EIGENVALPath].\n");
	my $in = new JFile($EIGENVALPath, "rb");
	return undef unless($in);

	my $line;
	for(my $i = 0 ; $i < 5 ; $i++) {
		$line = $in->ReadLine();
	}
	my $SampleName = $line;

	$line = $in->ReadLine();
	my ($nn, $nK, $nLevels) = Utils::Split("\\s+", $line);

	my @EnergyBands;
	for(my $iK = 0 ; $iK < $nK ; $iK++) {
#Skip blank line
		$line = $in->ReadLine();
		$line =~ s/^\s+//;
		$line =~ s/\s+$//;
		if($line eq '') {
			$line = $in->ReadLine();
		}
#print "line2: $line\n";
		my ($kx, $ky, $kz, $w) = Utils::Split("\\s+", $line);
#print "k[$iK]: $kx $ky $kz\n";
		for(my $iL = 0 ; $iL < $nLevels ; $iL++) {
			$line = $in->ReadLine();
#print "l: $line";
			my ($ib, $eup, $edn) = Utils::Split("\\s+", $line);
			$EnergyBands[$iL][$iK][0] = $eup - $EF;
			if(defined $edn) {
#print "Spin Polarized\n";
				$VASP->SpinPolarized(1);
				$EnergyBands[$iL][$iK][1] = $edn - $EF;
			}
			else {
				$VASP->SpinPolarized(0);
			}
		}
	}
	$in->Close();

	for(my $iL = 0 ; $iL < $nLevels ; $iL++) {
print "iL=$iL/$nLevels\n";
		$iBand++;

		my ($emin, $emax) = (1.0e10, -1.0e10);
		my $count = 0;
		my $iK = 0;
		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;
					}

					my $e = $EnergyBands[$iL][$iK][0];
#print "e[$ib][$ix][$iy][$iz][$count]=$e\n";
					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]) {
		Utils::DelSpace($line);
		print "Mismatch Line: [$count] [$line]\n";
		print "   Write: [$KIndex[$count]]\n";
		print "   Read : [$index]\n";
exit;
	}
}

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

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

$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 VASPFS.pl
   #
   # Launch as: xcrysden --bxsf $OutFile
   #
   Fermi Energy: $EFavrg
 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) {
						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);
}