#!/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);

#===============================================
# スクリプト大域変数
#===============================================
my $DataDir = ($ARGV[0])? $ARGV[0] : "/home/tkamiya/SrGeO3-Wien2k/";
$DataDir = "." if($DataDir eq '');
$DataDir = Utils::ReduceDirectory($DataDir);
print "DataDir: $DataDir\n";

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

my $Title                 = &GetTitle($DataDir);
my $INCARFile             = "INCAR";
my $KPOINTSFile           = $VASP->GetKPOINTSFileName($INCARFile);
my $StructPath            = Deps::MakePath($DataDir, "$Title.struct");
my $InEnergyPath          = Deps::MakePath($DataDir, "SrGeO3-Wien2k.energy.scissors");

#===============================================
# メインルーチン
#===============================================

my $Crystal = $WIEN2k->ReadStructFile($StructPath, 0);#, tol => 1.0e-2);
my $LatticeSystem = $Crystal->LatticeSystem();
print "Lattice system in WIEN2k: $LatticeSystem\n";

my $EF = $WIEN2k->ReadFermiEnergy($StructPath, 1);
print "EF=$EF eV\n";

print "\nRead [$InEnergyPath].\n";
my ($pBandInf, $pBands) = &ReadWIENEnergy($InEnergyPath, $Crystal);

&SaveVASPKPOINTSFromBandHash($KPOINTSFile, $Crystal, $pBandInf, $pBands);


exit;

sub SaveVASPKPOINTSFromBandHash
{
	my ($path, $Crystal, $pBandInf, $pBands) = @_;
	
	my $out = new JFile;
	$out->Open($path, "w") or die "$!: Can not write to [$path].\n";
	$out->print("Automatically generated mesh\n");
	my $nk = $pBandInf->{nk};
	$out->print("$nk\n");
	$out->print("Reciprocal lattice\n");
	for(my $i = 0 ; $i < $nk ; $i++) {
		my $inf = $pBands->[$i];
		$out->printf(" %19.14f %19.14f %19.14f %23f\n",
			$inf->{kx}, $inf->{ky}, $inf->{kz}, $inf->{Weight});
	}
	
	$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 ReadWIENEnergy
{
	my ($fname, $Crystal) = @_;

	my $in = new JFile($fname, "r") or die "$!: Can not read [$fname]\n";
	my @header;

	my @AtomSiteList = $Crystal->GetCAsymmetricAtomSiteList();
	my $nSite = @AtomSiteList;
	for(my $i = 0 ; $i < $nSite*2 ; $i++) {
		$header[$i] = $in->ReadLine();
	}

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

# 0.714285714286E-01 0.714285714286E-01 0.714285714286E-01         1   411    37  8.0   
		my ($kx, $ky, $kz, $iband) = Utils::Split("\\s+", substr($line, 0, 58));
		my ($iband, $nlevels, $w) = Utils::Split("\\s+", substr($line, 68));
		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,
			Weight  => $w,
			pLevels => \@levels
			);
		$Bands[$nk] = \%KInf;
		$nk++;
print "$nk: k=($kx,$ky,$kz): nL=$nlevels: w=$w\n";
#exit;
	}
print "nk=$nk\n";
	$in->Close();

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

	$BandInf{pBands} = \@Bands;

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