#!/usr/bin/perl
##!d:/Perl/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 CGI::Carp qw(fatalsToBrowser);
use Cwd;

use Deps;
use Utils;
use MyHTMLApplication;
use Crystal::VASP;
use Sci::ChemicalReaction;

BEGIN {
#	$| = 1;
	open (STDERR, ">&STDOUT");
}

#===============================================
# 大域変数
#===============================================
my $ScriptCharCode     = 'utf8';
my $FileSystemCharCode = 'sjis';

my $nLevel = 3;

my $KeyDir = 'MD'; #'SCF';

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

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

$App->SetOutputMode("console");

$App->AddArgument("--Positions", "--Positions=[indexes;]",       "");
$App->AddArgument("--AtomType",  "--AtomType=[atomtype]",       "");
exit 1 if($App->ReadArgs(1, "sjis", 0) != 1);
my $Args = $App->Args();

my $OutFile = "TotalEnergy.csv";
$OutFile    = $Args->GetGetArg(0) if(defined $Args->GetGetArg(0));

$App->print("Output File: [$OutFile]\n");

my ($pOutputTerms) = &GetConfigArrays($App);
my @OutputTerms  = @$pOutputTerms;

my @Dirs;
#my @Dirs = sort { $a cmp $b; } glob("*");
@Dirs = Utils::SearchFilesRecursive2(".", "*", $nLevel, \@Dirs, 1,
			sub { 
				my ($path) = @_;
				my $KeyDir  = Deps::MakePath($path, $KeyDir, 0);
#				print "dir: $path [$KeyDir]\n"; 

				return 0 if(!-d $path);
				if(-d $KeyDir) {
					print "dir: $path [$KeyDir]\n"; 
					return 1;
				}
#				return 1;# if(-d $SCFDir);
				return 0;
			},
			SearchDir => 1,
#			Debug     => 1,
			);
#exit;
#for(my $i= 0 ; $i < @Dirs ; $i++) {
#	print "$i: $Dirs[$i]\n";
#}
#exit;

my $R = new ChemicalReaction;

my $out = new JFile;
if(!$out->Open($OutFile, "w")) {
	$App->print("Error: Can not write to [$OutFile].\n");
	exit;
}

for(my $i = 0 ; $i < @OutputTerms ; $i += 2) {
	if($i == 0) {
		$out->print($OutputTerms[$i+1]);
	}
	else {
		$out->print(',' . $OutputTerms[$i+1]);
	}
}
$out->print("\n");

print("\nMakeSummaryHTML\n");
#print "nDirs = ", scalar @Dirs, "\n";
for(my $i = 0 ; $i < @Dirs ; $i++) {
#print "dir[$i]: $Dirs[$i]\n";
	next if(!-d $Dirs[$i]);
	next if($Dirs[$i] =~ /junk/i);
	
	my $KeyDir = Deps::MakePath($Dirs[$i], $KeyDir, 0);
	next if(!-d $KeyDir);

	$App->print("\n$i: $Dirs[$i]\n");
	my $ret = &MakeVASPSummary($App, $out, $Dirs[$i]);
}

$out->Close();

exit;

#===============================================
# Subroutines
#===============================================
my $Positions;
my $AtomType;

my $SiteNameThreshold = 0.05;
my @SiteNames = (
	['Cs', 0.0,  0.0],
	['Cs', 0.5,  0.5],
	['Cs', 0.0,  1.0],
	['Cs', 1.0,  0.0],
	['Cs', 1.0,  1.0],
	['F',  0.0,  0.5],
	['F',  0.5,  0.0],
	['F',  1.0,  0.5],
	['F',  0.5,  1.0],
	['P4', 0.25, 0.25],
	['P4', 0.25, 0.75],
	['P4', 0.75, 0.25],
	['P4', 0.75, 0.75],
	['P2', 0.0,  0.25],
	['P2', 0.0,  0.75],
	['P2', 0.25, 0.0],
	['P2', 0.25, 0.5],
	['P2', 0.25, 1.0],
	['P2', 0.75, 0.0],
	['P2', 0.75, 0.5],
	['P2', 0.75, 1.0],
	['P2', 1.0,  0.25],
	['P2', 1.0,  0.75],
	);

sub GetConfigArrays
{
	my ($App) = @_;

	$Positions = $Args->GetGetArg("Positions");
	$AtomType  = $Args->GetGetArg("AtomType");

	$App->print("Positions  : [$Positions]\n");
	$App->print("Atom Type  : [$AtomType]\n");

my @OutputTerms = (
	FileBody               => 'Source',
	Compound               => 'Compound',
	SortedComposition      => 'SortedComposition/Z',
	SumChemicalComposition => 'Composition',
	Type                   => 'Type',
	FormulaUnit            => 'Z',
	NELECT                 => 'NELECT',
#	AdAtom                 => 'AdAtom',
#	SiteName0              => 'SiteName0',
#	SiteName               => 'SiteName',
#	"Position0($Positions,x)" => "x0($AtomType)",
#	"Position0($Positions,y)" => "y0($AtomType)",
#	"Position0($Positions,z)" => "z0($AtomType)",
#	"Position($Positions,x)"  => "x($AtomType)",
#	"Position($Positions,y)"  => "y($AtomType)",
#	"Position($Positions,z)"  => "z($AtomType)",
	TOTEN                  => 'Etot',
	'Etot/mol(SCF)'        => 'Emol',
	'EF(SCF)'              => 'EF(SCF)',
	'EF(DOS)'              => 'EF(DOS)',
	VBM                    => 'EVBM',
	CBM                    => 'ECBM',
	Eg                     => 'Eg',
	Density                => 'Density',
	a                      => 'a',
	b                      => 'b',
	c                      => 'c',
	Volume                 => 'V',
	PseudoPotential        => 'PP',
	PREC                   => 'PREC',
	'aKmin(SCF)'           => 'aK(SCF)',
	'aKmin(VCRelax)'       => 'aK(VCRelax)',
	EDIFF                  => 'EDIFF',
	EDIFFG                 => 'EDIFFG',
	NBANDS                 => 'NBANDS',
	ISPIN                  => 'ISPIN',
	NUPDOWN                => 'NUPDOWN',
	LDAU                   => 'LDAU',
	LDAUU                  => 'U',
	LDAUJ                  => 'J',
	ISMEAR                 => 'ISMEAR',
	SIGMA                  => 'SIGMA',
	Status                 => 'Status',
	Ref                    => 'Ref',
	);

	return (\@OutputTerms);
}

sub MakeVASPSummary
{
	my ($App, $out, $RootDir) = @_;

	my ($drive, $directory, $filename, $ext1, $lastdir, $filebody) = Deps::SplitFilePath($RootDir);

	my $OriginalDir = Deps::MakePath($RootDir, "Original", 0);
	my $MDDir       = Deps::MakePath($RootDir, "MD",       0);
	my $VCRelaxDir  = Deps::MakePath($RootDir, "VCRelax",  0);
	my $SCFDir      = Deps::MakePath($RootDir, "SCF",      0);
	my $DOSDir      = Deps::MakePath($RootDir, "DOS",      0);

	my $vasp = new VASP;
	my $OriginalPOSCARHash = $vasp->ReadPOSCARtoHash( Utils::MakePath($OriginalDir, 'POSCAR',  '/', 0));
	my $VCRelaxKPOINTSHash = $vasp->ReadKPOINTStoHash(Utils::MakePath($VCRelaxDir,  'KPOINTS', '/', 0));
	my $INCARHash          = $vasp->ReadINCARtoHash(  Utils::MakePath($SCFDir,      'INCAR',   '/', 0));
	my $KPOINTSHash        = $vasp->ReadKPOINTStoHash(Utils::MakePath($SCFDir,      'KPOINTS', '/', 0));
	my $POTCARHash         = $vasp->ReadPOTCARtoHash( Utils::MakePath($SCFDir,      'POTCAR',  '/', 0));
	my $POSCARHash         = $vasp->ReadPOSCARtoHash( Utils::MakePath($SCFDir,      'POSCAR',  '/', 0));
	my $MDOUTCARHash       = $vasp->ReadOUTCARtoHash( Utils::MakePath($MDDir,       'OUTCAR',  '/', 0));
	my $VCROUTCARHash      = $vasp->ReadOUTCARtoHash( Utils::MakePath($VCRelaxDir,  'OUTCAR',  '/', 0));
	my $SCFOUTCARHash      = $vasp->ReadOUTCARtoHash( Utils::MakePath($SCFDir,      'OUTCAR',  '/', 0));
	if(!-d $SCFDir) {
		$INCARHash          = $vasp->ReadINCARtoHash(  Utils::MakePath($MDDir,      'INCAR',   '/', 0));
		$KPOINTSHash        = $vasp->ReadKPOINTStoHash(Utils::MakePath($MDDir,      'KPOINTS', '/', 0));
		$POTCARHash         = $vasp->ReadPOTCARtoHash( Utils::MakePath($MDDir,      'POTCAR',  '/', 0));
		$POSCARHash         = $vasp->ReadPOSCARtoHash( Utils::MakePath($MDDir,      'POSCAR',  '/', 0));
	}

	my $Status = '';
	if(-d $VCRelaxDir and $VCROUTCARHash->{RequiredAccuracyMessage} =~ /^\s*$/) {
		$Status = 'Relaxation not converged(1)';
	}
	elsif(-d $VCRelaxDir and $VCROUTCARHash->{AbortMessage} !~ /is\s+reached/) {
		$Status = 'Relaxation not converged(2)';
	}
	elsif(-d $VCRelaxDir and $VCROUTCARHash->{TotCPUTime} eq '') {
		$Status = 'VCR aborted(1)';
	}
	elsif($SCFOUTCARHash->{AbortMessage} =~ /^\s*$/) {
		$Status = 'SCF not converged(1)';
	}
	elsif($SCFOUTCARHash->{TotCPUTime} eq '') {
		$Status = 'SCF aborted(1)';
	}
	my ($EF, $VBM, $CBM, $Eg) = VASP->new()->ReadEgFromDOSOUTCAR(Deps::MakePath($DOSDir, 'OUTCAR', 0));

	$App->print("ISPIN=$INCARHash->{ISPIN}\n");
	$App->print("PREC=$INCARHash->{PREC}\n");
	$App->print("EDIFF=$INCARHash->{EDIFF}\n");
	$App->print("EDIFFG=$INCARHash->{EDIFFG}\n");
	$App->print("LDAU=$INCARHash->{LDAU}\n");
	$App->print("LDAUTYPE=$INCARHash->{LDAUTYPE}\n");
	$App->print("LDAUU=$INCARHash->{LDAUU}\n");
	$App->print("LDAUJ=$INCARHash->{LDAUJ}\n");
	$App->print("ISMEAR=$INCARHash->{ISMEAR}\n");
	$App->print("SIGMA=$INCARHash->{SIGMA}\n");
	$POTCARHash->{PseudoPotential} =~ s/ .*$//s;
	$App->print("PseudoPotential=$POTCARHash->{PseudoPotential}\n");

	$VCRelaxKPOINTSHash->{KPOINTS} =~ s/^[^\n]*\n[^\n]*\n[^\n]*\n//m;
	$VCRelaxKPOINTSHash->{KPOINTS} =~ s/\n.*$//m;
	Utils::DelSpace($VCRelaxKPOINTSHash->{KPOINTS});
	my @nKVCRelax = Utils::Split("\\s+", $VCRelaxKPOINTSHash->{KPOINTS});

	$KPOINTSHash->{KPOINTS} =~ s/^[^\n]*\n[^\n]*\n[^\n]*\n//m;
	$KPOINTSHash->{KPOINTS} =~ s/\n.*$//m;
	Utils::DelSpace($KPOINTSHash->{KPOINTS});
	my @nK = Utils::Split("\\s+", $KPOINTSHash->{KPOINTS});

	$App->print("KPOINTS=$KPOINTSHash->{KPOINTS}\n");
	$App->print("ISPIN=$INCARHash->{ISPIN}\n");
	$App->print("Etot=$SCFOUTCARHash->{TOTEN}\n");
	$App->print("EF(SCF)=$SCFOUTCARHash->{EF}\n");
	$App->print("EF(DOS)=$EF\n");
	$App->print("EVBM(DOS)=$VBM\n");
	$App->print("ECBM(DOS)=$CBM\n");
	$App->print("Eg(DOS)=$Eg\n");
	$App->print("SumChemicalComposition=$POSCARHash->{SumChemicalComposition}\n");
	$App->print("ChemicalComposition=$POSCARHash->{ChemicalComposition}\n");
	$App->print("FormulaUnit=$POSCARHash->{FormulaUnit}\n");

	$POSCARHash->{Latt} =~ s/\n/,/g;
	my @latt = Utils::Split(",+", $POSCARHash->{Latt});

	my @aKVCRelax = ($latt[0]*$nKVCRelax[0], $latt[1]*$nKVCRelax[1], $latt[2]*$nKVCRelax[2]);
	my ($aKminVCRelax, $aKmaxVCRelax) = Utils::CalMinMax(\@aKVCRelax);
	$aKminVCRelax = sprintf("%4.1f", $aKminVCRelax);

	my @aK = ($latt[0]*$nK[0], $latt[1]*$nK[1], $latt[2]*$nK[2]);
	my ($aKmin, $aKmax) = Utils::CalMinMax(\@aK);
	$aKmin = sprintf("%4.1f", $aKmin);

	$App->print("Latt=$POSCARHash->{Latt}\n");
	$App->print("aKMax=$aKmin-$aKmax\n");
	$App->print("Volume=$POSCARHash->{Volume}\n");
	$App->print("Density=$POSCARHash->{Density}\n");

	my $Ref = Deps::MakePath(cwd(), $RootDir, 0);
	$Ref =~ s/$ENV{HOME}//;
	my $Eavrg    = $SCFOUTCARHash->{TOTEN} / $POSCARHash->{FormulaUnit};
	my $Compound = $POSCARHash->{ChemicalComposition};
#	my $Compuond = $POSCARHash->{SumChemicalComposition};
	my $SortedComposition = $R->SortChemicalFormula($POSCARHash->{ChemicalComposition});
	$Compound =~ s/\s//g;
	$Status .= " / SCFWarning:$SCFOUTCARHash->{Warning}" if($SCFOUTCARHash->{Warning} ne '');
	$Status .= " / VCRWarning:$VCROUTCARHash->{Warning}" if($VCROUTCARHash->{Warning} ne '');
	
	my %v = (
		FileBody         => $filebody,
		Compound         => $Compound,
		SortedComposition      => $SortedComposition,
		SumChemicalComposition => $POSCARHash->{SumChemicalComposition},
		FormulaUnit      => $POSCARHash->{FormulaUnit},
		TOTEN            => $SCFOUTCARHash->{TOTEN},
		'Etot/mol(SCF)'  => $Eavrg,
		'EF(SCF)'        => $SCFOUTCARHash->{EF},
		'EF(DOS)'        => $EF,
		VBM              => $VBM,
		CBM              => $CBM,
		Eg               => $Eg,
		Density          => $POSCARHash->{Density},
		a                => $latt[0],
		b                => $latt[1],
		c                => $latt[2],
		alpha            => $latt[3],
		beta             => $latt[4],
		gamma            => $latt[5],
		Volume           => $POSCARHash->{Volume},
		PseudoPotential  => $POTCARHash->{PseudoPotential},
		PREC             => $INCARHash->{PREC},
		'aKmin(SCF)'     => $aKmin,
		'aKmax(SCF)'     => $aKmax,
		'aKmin(VCRelax)' => $aKminVCRelax,
		EDIFF            => $INCARHash->{EDIFF},
		EDIFFG           => $INCARHash->{EDIFFG},
		NBANDS           => $INCARHash->{NBANDS},
		NELECT           => $SCFOUTCARHash->{NELECT},
		ISPIN            => $INCARHash->{ISPIN},
		NUPDOWN          => $INCARHash->{NUPDOWN},
		LDAU             => $INCARHash->{LDAU},
		LDAUU            => $INCARHash->{LDAUU},
		LDAUJ            => $INCARHash->{LDAUJ},
		ISMEAR           => $INCARHash->{ISMEAR},
		SIGMA            => $INCARHash->{SIGMA},
		Status           => $Status,
		Ref              => $Ref,
		);

	my $pAtomSite0 = $OriginalPOSCARHash->{pAtomSite};
	for(my $i = 0 ; $i < @$pAtomSite0 ; $i++) {
		my ($x, $y, $z) = $pAtomSite0->[$i]->Position();
		$v{"Position0($i,x)"} = $x;
		$v{"Position0($i,y)"} = $y;
		$v{"Position0($i,z)"} = $z;
	}

	my $pAtomSite = $POSCARHash->{pAtomSite};
	for(my $i = 0 ; $i < @$pAtomSite ; $i++) {
		my ($x, $y, $z) = $pAtomSite->[$i]->Position();
		$v{"Position($i,x)"} = $x;
		$v{"Position($i,y)"} = $y;
		$v{"Position($i,z)"} = $z;
	}

	my ($AdAtom) = ($filebody =~ /^([A-Z][a-z]?)on/);
	$v{AdAtom} = $AdAtom;

	my ($x0, $y0, $z0) = ($v{"Position0($Positions,x)"}, $v{"Position0($Positions,y)"}, $v{"Position0($Positions,z)"});
	my $minr0 = 1.0e9;
	my $pmin0;
	for(my $i = 0 ; $i < @SiteNames ; $i++) {
		my $p = $SiteNames[$i];
		my ($name, $xc, $yc, $zc) = ($p->[0], $p->[1], $p->[2], $p->[3]);
		
		my $r2 = ($x0 - $xc)*($x0 - $xc) + ($y0 - $yc)*($y0 - $yc); # + ($z0 - $zc)*($z0 - $zc);
		my $r  = sqrt($r2);

		if($r < $minr0) {
			$pmin0 = $p;
			$minr0 = $r;
		}
	}
	if($minr0 > $SiteNameThreshold) {
		$v{SiteName0} = "near $pmin0->[0]";
	}
	else {
		$v{SiteName0} = "$pmin0->[0]";
	}

	my ($x1, $y1, $z1) = ($v{"Position($Positions,x)"},  $v{"Position($Positions,y)"},  $v{"Position($Positions,z)"});
	my $minr1 = 1.0e9;
	my $pmin1;
	for(my $i = 0 ; $i < @SiteNames ; $i++) {
		my $p = $SiteNames[$i];
		my ($name, $xc, $yc, $zc) = ($p->[0], $p->[1], $p->[2], $p->[3]);
		
		my $r2 = ($x1 - $xc)*($x1 - $xc) + ($y1 - $yc)*($y1 - $yc); # + ($z1 - $zc)*($z1 - $zc);
		my $r  = sqrt($r2);

		if($r < $minr1) {
			$pmin1 = $p;
			$minr1 = $r;
		}
	}
	if($minr1 > $SiteNameThreshold) {
		$v{SiteName} = "near $pmin1->[0]";
	}
	else {
		$v{SiteName} = "$pmin1->[0]";
	}

#========================================================
# Add to CSV
#========================================================
	for(my $i = 0 ; $i < @OutputTerms ; $i += 2) {
		my $key = $OutputTerms[$i];
		if($i == 0) {
			$out->print($v{$key});
		}
		else {
			$out->print(',' . $v{$key});
		}
	}
	$out->print("\n");

	return 1;
}