#!/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;
use Sci qw($RyToeV);

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

#===============================================
# スクリプト大域変数
#===============================================
my $VASPPath     = "vasp";
my $SGROUPPath   = ($^O =~ /^MSWin/)? "D:/Programs/Perl/WIEN2k/sgroup.exe" : "sgroup";
my $XCrySDenPath = "xcrysden";
my $tol          = 1.0e-3;

my $Denominator   = 1000;

#==========================================
# コマンドラインオプション読み込み
#==========================================
$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 $DataDir = $App->Args()->GetGetArg(0);
$DataDir = "." if($DataDir eq '');
$DataDir = Utils::ReduceDirectory($DataDir);
my $INCARFile = "INCAR";

my $VASP   = new VASP;
my $WIEN2k = new WIEN2k;

my $Title                 = &GetTitle($DataDir);
my $KListPath             = $VASP->GetKPOINTSFileName($DataDir);
my $EIGENVALPath          = $VASP->GetEIGENVALFileName($DataDir);
my $OUTCAR                = $VASP->GetOUTCARFileName($EIGENVALPath);
my $StructPath            = Deps::MakePath($DataDir, "$Title.nosym.struct");
my $SymmetrizedStructPath = Deps::MakePath($DataDir, "$Title.struct");
my $InEnergyPath          = "LaFeOP.energy";
my $EnergyPath            = Deps::MakePath($DataDir, "$Title.energy");
my $IntransPath           = Deps::MakePath($DataDir, "$Title.intrans");
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;
}

my $LatticeSystem = $Crystal->LatticeSystem();
print "Lattice system in VASP: $LatticeSystem\n";

if(1) {
print "\Make [$StructPath][$SymmetrizedStructPath].\n";
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-2);
$LatticeSystem = $Crystal->LatticeSystem();
#my $SPG = $crystal->GetCSpaceGroup();
#$LatticeSystem = SetLatticeSystem(undef, 1, 1.0e-2);
print "Lattice system in WIEN2k: $LatticeSystem\n";

my ($pBandInf, $pBands);
if(1) {
my $EF = $VASP->ReadFermiEnergy($OUTCAR);
$App->print("EF in OUTCAR: $EF eV\n");
#my ($nn, $nK, $nLevels, $pEnergyBands) = $VASP->ReadEIGENVALtoArray($EIGENVALPath, $EF);

print "\nRead [$EIGENVALPath].\n";
($pBandInf, $pBands) = &ReadVASPEIGENVAL($EIGENVALPath);
my $nK = @$pBands;
print "nK=$nK\n";
my @BandArray;
for(my $i = 0 ; $i < $nK ; $i++) {
	my $pKInf = $pBands->[$i];
print "$i: pKInf: $pKInf\n";
	my $ikx = int($pKInf->{kx}*$Denominator+0.01);
	my $iky = int($pKInf->{ky}*$Denominator+0.01);
	my $ikz = int($pKInf->{kz}*$Denominator+0.01);
	&AddBandToArray(\@BandArray, $ikx, $iky, $ikz, $pKInf);
	if($LatticeSystem eq 'cubic') {
		&AddBandToArray(\@BandArray, $ikx, $ikz, $iky, $pKInf);
		&AddBandToArray(\@BandArray, $iky, $ikx, $ikz, $pKInf);
		&AddBandToArray(\@BandArray, $iky, $ikz, $ikx, $pKInf);
		&AddBandToArray(\@BandArray, $ikz, $ikx, $iky, $pKInf);
		&AddBandToArray(\@BandArray, $ikz, $iky, $ikx, $pKInf);
	}
	elsif($LatticeSystem eq 'tetragonal') {
		&AddBandToArray(\@BandArray, $iky, $ikx, $ikz, $pKInf);
	}
	else {
print "Lattice system [$LatticeSystem] is not supported\n";
	}
}
}

if(0) {
print "\nRead [$InEnergyPath].\n";
($pBandInf, $pBands) = &ReadWIENEnergy($InEnergyPath);
}

&SaveWIENEnergy($EnergyPath, $pBandInf, $pBands, $Crystal);

&CreateIntrans($IntransPath);

exit;

sub CreateIntrans
{
	my ($path) = @_;
	
	open(OUT, ">$path") or die "Error: Can not create [$path].\n";
	print OUT<<EOT;
WIEN                      # Format of DOS
0 1 0 0.0                 # iskip (not presently used) idebug setgap shiftgap
0.54051 0.0005 0.4 240  # Fermilevel (Ry), energygrid, energy span around Fermilevel, number of electrons
CALC                      # CALC (calculate expansion coeff), NOCALC read from file
5                         # lpfac, number of latt-points per k-point
BOLTZ                     # run mode (only BOLTZ is supported)
0.15                     # (efcut) energy range of chemical potential
800 50                # Tmax, temperature grid
-1.                        # energyrange of bands given individual DOS output sig_xxx and dos_xxx (xxx is band number)
TETRA
EOT
	close(OUT);
}

sub SaveWIENEnergy {
	my ($path, $pBandInf, $pBands, $Crystal) = @_;

print "\nCreate [$path].\n";
	
	my $out = new JFile;
	$out->Open($path, "w") or die "Error: Can not write to [$path].\n";
	my @AtomSiteList = $Crystal->GetCAsymmetricAtomSiteList();
	my $nSite = @AtomSiteList;
print "nSite=$nSite\n";
	for(my $i = 0 ; $i < $nSite ; $i++) {
		$out->print(
			 "200.30000200.30000200.33500  0.30000  0.30000  0.30000  0.30000  0.30000  0.30000  0.30000  0.30000  0.30000  0.30000  0.00000\n"
			."  0.30000  0.30000  0.33500999.00000997.00000 -3.10650997.00000999.00000999.00000999.00000999.00000999.00000\n");
	}

	my $nK = @$pBands;
#print "nK=$nK\n";
	my @BandArray;
	for(my $i = 0 ; $i < $nK ; $i++) {
		my $pKInf = $pBands->[$i];
#print "$i: pKInf: $pKInf\n";
		my $mul = 1.0;
		my ($IsConverted, $kx, $ky, $kz) = $Crystal->ConvertLatticeIndexByLatticeSystem(
				$pKInf->{kx}, $pKInf->{ky}, $pKInf->{kz}, 1);
		$out->printf(" %18.11e %18.11e %18.11e%-5s %4d %5d%6d%5.1g\n", 
				$kx, $ky, $kz, '', 1, 713, $pKInf->{nLevels}, $mul);

		my $pLevels = $pKInf->{pLevels};
		for(my $il = 0 ; $il < $pKInf->{nLevels} ; $il++) {
			$out->printf("     %7d %23.15e\n", $il+1, $pLevels->[$il] / $RyToeV);
		}
	}
	$out->Close();
}

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

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

sub AddBandToArray
{
	my ($pBandArray, $ikx, $iky, $ikz, $pKInf) = @_;

	$pBandArray->[$ikx][$iky][$ikz] = $pKInf;
print "ik: $ikx, $iky, $ikz\n";
}


sub ReadWIENEnergy
{
	my ($fname) = @_;

	my $in = new JFile($fname, "r") or die "$!: Can not read [$fname]\n";
	my @header;
	for(my $i = 0 ; $i < 8 ; $i++) {
		$header[$i] = $in->ReadLine();
	}

	my @Bands;
	my $nk = 0;
	while(1) {
		my $line = $in->ReadLine();
#print "line: [$line]\n";
		last if(!defined $line);

		my ($kx, $ky, $kz, $iband, $nlevels, $f) = Utils::Split("\\s+", $line);
		last if(!defined $kz);

		$kx += 0.0;
		$ky += 0.0;
#print "kz=$kz => ";
		$kz =~ s/(E[+-]\d\d).*?$/$1/;
#print "$kz\n";
#exit;
		$kz += 0.0;
		

		my @levels;
		for(my $i = 0 ; $i < $nlevels ; $i++) {
			my $pos = $in->tell();
			my $line = $in->ReadLine();
			if($line !~ /^     /) {
				$in->seek($pos, 0);
				last;
			}
#print "l: $line\n";
			my $idx;
			($idx,  $levels[$i]) = Utils::Split("\\s+", $line);
			$levels[$i] += 0.0;
#print "  $nlevels: $idx  $e Ry\n";
		}
		my %KInf = (
			nLevels => $nlevels,
			kx      => $kx,
			ky      => $ky,
			kz      => $kz,
			pLevels => \@levels
			);
		$Bands[$nk] = \%KInf;
		$nk++;
print "$nk: k=($kx,$ky,$kz): nL=$nlevels\n";
#exit;
	}
print "nk=$nk\n";
	$in->Close();

	my %BandInf = (
		pHeader => \@header,
		nk      => $nk,
		);

	$BandInf{pBands} = \@Bands;

	return (\%BandInf, \@Bands);
}

sub ReadVASPEIGENVAL
{
	my ($fname) = @_;

	my $in = new JFile($fname, "r") or die "$!: Can not read [$fname]\n";
	my @header;
	for(my $i = 0 ; $i < 5 ; $i++) {
		$header[$i] = $in->ReadLine();
	}
	my ($na, $nk, $nlevels) = Utils::Split("\\s+", $in->ReadLine());
print "na, nk, nlevels: $na, $nk, $nlevels\n";

	my %BandInf = (
		pHeader => \@header,
		na      => $na,
		nk      => $nk,
		nLevels => $nlevels,
		);

	my $iSkip = 1;
	my @Bands;
	my $c = 0;
	for(my $ik = 0 ; $ik < $nk ; $ik++) {
		next if($ik % $iSkip != 0);

		my $blank  = $in->ReadLine();
		my $header = $in->ReadLine();
		my ($kx, $ky, $kz, $w) = Utils::Split("\\s+", $header);
		my @levels;
		for(my $il = 0 ; $il < $nlevels ; $il++) {
			my $ilevel;
			($ilevel, $levels[$il]) = Utils::Split("\\s+", $in->ReadLine());
		}
		my %KInf = (
			nLevels => $nlevels,
			Blank   => $blank,
			Header  => $header,
			kx      => $kx,
			ky      => $ky,
			kz      => $kz,
			w       => $w,
			pLevels => \@levels
			);
		$Bands[$c] = \%KInf;
		$c++;
	}
	$in->Close();
	
	$BandInf{pBands} = \@Bands;

	return (\%BandInf, \@Bands);
}