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

#===============================================
# 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 = $App->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 $ret;
	($ret, @KIndex) = $WIEN2k->MakeKListForFS($KListPath, "check.txt", $Crystal, $Title, $UseSymmetry, $nx, $ny, $nz, $div, $emin, $emax, $K, 1);
	return $ret;
}

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

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

	if(!XCrySDen->new()->WriteBXSFFile($OutFile, $PrimCrystal, 
				$Title, $nBand, $nx, $ny, $nz, $EF, \@BandIndex, \@BandEnergy, "WIEN2kFS.pl", $OneByOne)) {
		die "$!: Can not write to [$OutFile]\n";
	}
}