#!/usr/bin/perl

#
# SpinPolarize計算に未対応
# エネルギー準位毎にFSを出力している。
#正しくは、バンド毎に出さないといけない
#

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     = "mpirun -n 16 -machinefile hosts vasp_std";
my $SGROUPPath   = "sgroup";
my $XCrySDenPath = "xcrysden";
my $tol          = 1.0e-3;
# for klist
my $K   = 1;
# for klist & bxsf
#my ($nx, $ny, $nz)  = (20, 20, 20);
#my ($nx, $ny, $nz)  = (30, 30, 30);
#my ($nx, $ny, $nz)  = (40, 40, 40);
#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("--EFmin",       "--EFmin=val (Def: -1.0)", "-1.0");
$App->AddArgument("--EFmax",       "--EFmax=val (Def:  1.0)", "1.0");
$App->AddArgument("--aKProduct",   "--aKProduct=val (Def: 2.0)", "2.0");
$App->AddArgument("--UseSymmetry", "--UseSymmetry=[0|1] (Def: 1)", "1");
$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 $EFmin     = $App->Args()->GetGetArg('EFmin');
   $EFmin     = -1.0 if(!defined $EFmin);
my $EFmax     = $App->Args()->GetGetArg('EFmax');
   $EFmax     =  1.0 if(!defined $EFmax);
my $aKProduct = $App->Args()->GetGetArg('aKProduct');
   $aKProduct = 2.0 if(!defined $aKProduct);
my ($nx, $ny, $nz, $NKREDX, $NKREDY, $NKREDZ) = $VASP->AnalyzeaKProduct($Crystal, $aKProduct);
my $UseSymmetry = $App->Args()->GetGetArg('UseSymmetry');
   $UseSymmetry = 1 if(!defined $UseSymmetry);

$App->print("EF range   : $EFmin - $EFmax\n");
$App->print("UseSymmetry: $UseSymmetry\n");
$App->print("aKProduct  : $aKProduct [nk = ($nx, $ny, $nz)]\n");

my $Function = "FS";
if(!$VASP->ModifyINCARFile($Function, $INCARFile, "${INCARFile}.$Function")) {
	$App->print("Error: Can not write to [${INCARFile}.$Function]\n");
	exit 0;
}
Utils::CopyFile("${INCARFile}.$Function", $INCARFile);

my $WIEN2k = new WIEN2k;
unlink($SymmetrizedStructPath);
if(!$WIEN2k->SaveSymmetrizedStructFileFromCrystal(
		$Crystal, $StructPath, $SymmetrizedStructPath, $SGROUPPath, $tol, 1)) {
	$App->print("Error: Can not create [$SymmetrizedStructPath].\n");
	exit 0;
}

$Crystal = $WIEN2k->ReadStructFile($SymmetrizedStructPath, 0, tol => 1.0e-3);
#$Crystal->SetLatticeSystem(undef, 1, 1.0e-3, 1.0e-2);
$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 = $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 MakeKPOINTS
{
	my ($Crystal, $KListPath) = @_;

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

	my $ret;
	($ret, @KIndex) = $VASP->MakeKPOINTSForFS($KListPath, "check.txt", $Crystal, $Title, $UseSymmetry, $nx, $ny, $nz, $K, 1);
	return $ret;
}

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

	$App->print("\nCalculate bands\n");
	exit 0 if($App->Execute("$VASPPath"));
}

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

unlink($OutFile);
	$App->print("\nMake Band XSF(BXSF) file\n");
	my $OUTCAR = $VASP->GetOUTCARFileName($EIGENVALPath);
	my $pHash  = $VASP->ReadINCARtoHash($EIGENVALPath);
	my $ISPIN  = $pHash->{ISPIN};
	my $EF     = $VASP->ReadFermiEnergy($OUTCAR);
	my $EFavrg = ($EFmin + $EFmax) / 2.0;

	$App->print("EIGENVAL Path: $EIGENVALPath\n");
	$App->print("BXSF Path    : $BXSFPath\n");
	$App->print("ISPIN          : $ISPIN\n");
	$App->print("EF in OUTCAR   : $EF eV\n");
	$App->print("Search EF range: $EFmin - $EFmax eV\n");
	$App->print("Average EF     : $EFavrg\n");

	my $PrimCrystal = $Crystal->GetPrimitiveCrystal(0, 1);

$App->print("  Read [$EIGENVALPath].\n");
	my ($nn, $nK, $nLevels, $pEnergyBands) = $VASP->ReadEIGENVALtoArray($EIGENVALPath, $EF);
	my @EnergyBands = @$pEnergyBands;
	return undef if(@EnergyBands == 0);

	my (@BandEnergy, @BandIndex);
	my $iBand = 0;
	my $nBand = 0;
	my $line;
	for(my $is = 0 ; $is < $ISPIN ; $is++) {
	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++) {
#print "i=($ix,$iy,$iz)\n";
#print "cs: ", $Crystal->LatticeSystem(), "\n";
					my ($IsConverted, $cix, $ciy, $ciz) = $Crystal->ConvertLatticeIndexByLatticeSystem($ix, $iy, $iz);
#print("VASPFS.pl: $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][$is];
#print "e[$iL][$ix][$iy][$iz][$count]=$e\n";
					if(!defined $e) {
						$App->print("Error: Wrong format [iL=$iL  iK=$iK (ix,iy,iz)=($ix,$iy,$iz)/($cix,$ciy,$ciz)]\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++;
				}
			}
		}
#print "Examin... e=($emin - $emax) EF=($EFmin - $EFmax)\n";
		next if($emin > $EFmax or $emax < $EFmin);
#print "Passed: e=($emin - $emax) EF=($EFmin - $EFmax)\n";

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

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