#!/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 warnings;
use File::Path;
use File::Basename;
use File::Find;

use Deps;
use Utils;
use JFile;

use MyApplication;

use Crystal::CIF;
use Crystal::Crystal;
use Crystal::SpaceGroup;
use Crystal::AtomType;
use Crystal::VASP;

#===============================================
# デバッグ関係変数
#===============================================
#$PrintLevelが大きいほど、情報が詳しくなる
my $PrintLevel = 0;

#===============================================
# 文字コード関係変数
#===============================================
# sjis, euc, jis, noconv
my $PrintCharCode      = Deps::PrintCharCode();
my $OSCharCode         = Deps::OSCharCode();
my $FileSystemCharCode = Deps::FileSystemCharCode();

my $LF        = Deps::LF();
my $DirSep    = Deps::DirSep();
my $RegDirSep = Deps::RegDirSep();

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

#$App->SetLF("<br>\n");
#$App->SetPrintCharCode("sjis");
#$App->SetDebug($Debug);
$App->SetDeleteHTMLFlag(1);

#===============================================
# スクリプト大域変数
#===============================================
my $RootDir = '.';

my ($nData, $pTE, $pDOS) = &MakeDOSCSVfromPRODOSCARs($RootDir);

exit;

sub MakeDOSCSVfromPRODOSCARs
{
	my ($RootDir) = @_;

	$App->print("\n<b>Make lm-decomposed DOS CSV file from PROCAR/DOSCAR:</b>\n");

	my $VASP = new VASP();
	my $CSVPath     = $VASP->GetFileName($RootDir, 'DOSlm.csv');
	$App->print("CSV           : $CSVPath\n");

#	my $Width = $Args->GetGetArg("Width");
#	$Width = 0 if(!defined $Width);
#	$App->print("Width: $Width eV\n");
#	my $IgnoreZero = $Args->GetGetArg("IgnoreZero");
#	$IgnoreZero = 0 if(!defined $IgnoreZero);
#	$App->print("$IgnoreZero: $IgnoreZero\n");

	&ReadDOSCAR($RootDir);

	$App->print("\n<b>Make lm-decomposed DOS CSV file from PROCAR/DOSCAR: finished</b>\n");
}

sub ReadDOSCAR
{
	my ($RootDir) = @_;

	my $VASP = new VASP();

	my $PROCARPath  = $VASP->GetFileName($RootDir, 'PROCAR');
print "ReadDOSCAR: Read [$PROCARPath].\n";
	my @Orbitals = &GetOrbitalsFromPROCAR($PROCARPath);
	if(@Orbitals == 0) {
		return -1;
	}

	my $POSCAR = $VASP->GetPOSCARFileName($RootDir);
print "ReadDOSCAR: Read [$POSCAR].\n";
	my $Crystal = $VASP->ReadStructureFromCARFiles($POSCAR, 1);

	my @AtomTypeList = $Crystal->GetCAtomTypeList();
	my $nAtomType = @AtomTypeList;
	my %AtomNum;
	for(my $i = 0 ; $i < $nAtomType ; $i++) {
		my $name = $AtomTypeList[$i]->AtomNameOnly();;
		$AtomNum{$name} = 0 if(!defined $AtomNum{$name});
		my $n = ++$AtomNum{$name};
		$AtomTypeList[$i]->SetAtomName("$name$n");
#print "$i: $name$n\n";
	}
	my @ExpandedAtomSiteList = $Crystal->GetCExpandedAtomSiteListByOutputMode();
#	my @M = $Crystal->GetCnMultiplicityExpandedAtomSiteList();
	my $nTotalExpandedAtomSite = @ExpandedAtomSiteList;
	my @AtomNames;
	for(my $i = 0 ; $i < $nTotalExpandedAtomSite ; $i++) {
		my $atomsite = $ExpandedAtomSiteList[$i];
		$AtomNames[$i] = $atomsite->AtomName(); #AtomNameOnly();
#print "Site[$i]: $AtomNames[$i] ($M[$i])\n";
print "Site[$i]: $AtomNames[$i]\n";
	}

	my $OUTCAR = $VASP->GetOUTCARFileName($RootDir);
	my $EF     = $VASP->ReadFermiEnergy($OUTCAR);
print "EF: $EF eV\n";

	my $DOSCARPath  = $VASP->GetFileName($RootDir, 'DOSCAR');
print "ReadDOSCAR: Read [$DOSCARPath].\n";

	my $in = new JFile($DOSCARPath, 'r');
	if(!$in) {
		print "Error in ReadDOSCAR: $!: Can not read [$DOSCARPath].\n";
		return -1;
	}

	my $line;
	for(my $i = 0 ; $i < 5 ; $i++) {
		$line = $in->ReadLine();
	}
# BaZn2As2-I4mmm [PAW_PBE] (dos)          
	my $SampleName = $line;
	$SampleName =~ s/^\s*//;
	$SampleName =~ s/\s*$//;
print "Sample [$SampleName]\n";

#     15.00000000    -35.00000000 5000      3.50770154      1.00000000
	$line = $in->ReadLine();
	my ($E1, $E0, $nData, $EF, $w) = Utils::Split("\\s+", $line);
print "nData=$nData\n";

	my (@TE, %DOS);
	for(my $i = 0 ; $i < $nData ; $i++) {
		$line = $in->ReadLine();
		my ($e, $tot, $ne) = Utils::Split("\\s+", $line);
		$TE[$i]              = $e;
		$DOS{Total}[$i]        += $tot;
		$DOS{Interstitial}[$i] += $tot;
		$DOS{NE}[$i]           += $ne;
#print "  $i: $e\t$up\t$dn\n";
	}

	my $iAtom = 0;
	for(my $i = 0 ; $i < @ExpandedAtomSiteList ; $i++) {
#	while(1) {
		$line = $in->ReadLine();
		my ($E1, $E0, $nData, $EF, $w) = Utils::Split("\\s+", $line);
		my $atomsite = $ExpandedAtomSiteList[$iAtom];
print "iAtom=$iAtom\n";
		my $AtomName = $atomsite->AtomNameOnly();
print "Site[$iAtom]: $AtomName: nData=$nData\n";

		for(my $i = 0 ; $i < $nData ; $i++) {
			$line = $in->ReadLine();
			my ($e, @a) = Utils::Split("\\s+", $line);
			last if(!$line or @a <= 1);
if($i == 0) {
	Utils::DelSpace($line);
	my $nOrb = scalar @a;
	print "  nO=$nOrb [$line]\n";
}
			for(my $j = 0 ; $j < @a ; $j++) {
				my $orb = $Orbitals[$j];
#print("$j: $orb  ");
				$DOS{Interstitial}[$i]     -= $a[$j];
				$DOS{$AtomName}{Total}[$i] += $a[$j];
				$DOS{$AtomName}{$orb}[$i]  += $a[$j];
			}
		}
		last if($in->eof());

		$iAtom++;
	}

	$in->Close();

	my $CSVPath  = $VASP->GetFileName($RootDir, 'DOSlm.csv');
print "ReadDOSCAR: Write to [$CSVPath].\n";
	my $csv = new JFile($CSVPath, 'w');
	if(!$csv) {
		print "Error in MakeDOSCAR:: $!: Can not write to [$CSVPath].\n";
		return -2;
	}

	$csv->print("E(eV),Total,Interstitial,NE");
	my $nAtomType = @AtomTypeList;
	for(my $i = 0 ; $i < $nAtomType ; $i++) {
		my $atom     = $AtomTypeList[$i];
		my $atomname = $atom->AtomNameOnly();
		$csv->print(",$atomname(Total)");
		for(my $j = 0 ; $j < @Orbitals ; $j++) {
			my $orb = $Orbitals[$j];
			$csv->print(",$atomname($orb)");
		}
	}
	$csv->print("\n");
	
	for(my $i = 0 ; $i < $nData ; $i++) {
		$csv->print("$TE[$i],$DOS{Total}[$i],$DOS{Interstitial}[$i],$DOS{TotalNE}[$i]");
		for(my $it = 0 ; $it < $nAtomType ; $it++) {
			my $atom     = $AtomTypeList[$it];
			my $atomname = $atom->AtomNameOnly();
			$csv->print(",$DOS{$atomname}{Total}[$i]");
			for(my $j = 0 ; $j < @Orbitals ; $j++) {
				my $orb = $Orbitals[$j];
				$csv->print(",$DOS{$atomname}{$orb}[$i]");
			}
		}
		$csv->print("\n");
	}

	$csv->Close();
	
	return ($nData, \@TE, \%DOS);
}


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

	my $in = new JFile($path, 'r');
	if(!$in) {
		print "Error in GetOrbitalsFromPROCAR: $!: Can not read [$path].\n";
		return -1;
	}
	if(!$in->SkipTo("band\\s+\\d")) {
		print "Error in GetOrbitalsFromPROCAR: Can not find 'band' label in [$path].\n";
		return -2;
	}
	$in->ReadLine();
	my $line = $in->ReadLine();
	my ($ion, @a) = Utils::Split("\\s+", $line);
	my $na = @a;
	if($a[$na-1] eq 'tot') {
		delete $a[$na-1];
	}
	
	$in->Close();
	
	return @a;
}