#!/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 $ExistedLines = 0;

my $OutFile = "TotalEnergy.csv";
$OutFile    = $ARGV[0] if(defined $ARGV[0]);
my $Dir     = ($ARGV[1] eq '')? '.' : $ARGV[1];
#my $Target  = ($ARGV[1] eq '')? '' : $ARGV[1];

my @OutputTerms = (
	index                  => 'index',
	FileBody               => 'Source',
	TEBEG                  => 'TEBEG',
	TEEND                  => 'TEEND',
	Pexternal              => 'P(ext) [kB]',
	PPullay                => 'P(Pullay) [kB]',
	Ekin                   => 'E(kin) [kB]',
	TMD                    => 'T(MD) [K]',
	Compound               => 'Compound',
	SortedComposition      => 'SortedComposition/Z',
	SumChemicalComposition => 'Composition',
	Type                   => 'Type',
#	FormulaUnit            => 'Z',
#	'Position(16,x)'       => 'x(Se)',
#	'Position(16,y)'       => 'y(Se)',
#	'Position(16,z)'       => 'z(Se)',
	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',
	P11                    => 'P11',
	P22                    => 'P22',
	P33                    => 'P33',
	P23                    => 'P23',
	P31                    => 'P31',
	P12                    => 'P12',
	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',
	);


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

$App->SetOutputMode("console");

my $R = new ChemicalReaction;

my $IsExist = -f $OutFile;
if($IsExist) {
	$ExistedLines = 0;
	my $in = new JFile;
	if(!$in->Open($OutFile, "r")) {
		$App->print("Error: Can not read [$OutFile].\n");
		exit;
	}
	$in->ReadLine();
	while(1) {
		my $line = $in->ReadLine();
		last if(!defined $line);

		$line = Utils::DelSpace($line);
		next if($line eq '');
		$ExistedLines++;
	}
	$in->Close();
}

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

if(!$IsExist) {
	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("\nAdd MD History to [$OutFile]\n");

my $SCFDir = Deps::MakePath($Dir, "SCF", 0);
$App->print("Dir: $Dir\n");
my $ret = &AddVASPSummary($App, $out, $Dir);

$out->Close();

exit;

#===============================================
# Subroutines
#===============================================
sub AddVASPSummary
{
	my ($App, $out, $RootDir) = @_;

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

	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 $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 $SCFOUTCARHash      = $vasp->ReadOUTCARtoHash (Utils::MakePath($SCFDir, 'OUTCAR', '/', 0));
	my $VCRelaxKPOINTSHash = $vasp->ReadKPOINTStoHash(Utils::MakePath($VCRelaxDir, 'KPOINTS', '/', 0));
	my $VCROUTCARHash      = $vasp->ReadOUTCARtoHash (Utils::MakePath($VCRelaxDir, '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));
		$SCFOUTCARHash      = $vasp->ReadOUTCARtoHash (Utils::MakePath($MDDir, 'OUTCAR', '/', 0));
	}
	if(!-d $VCRelaxDir) {
		$VCRelaxKPOINTSHash = $vasp->ReadKPOINTStoHash(Utils::MakePath($MDDir, 'KPOINTS', '/', 0));
		$VCROUTCARHash      = $vasp->ReadOUTCARtoHash (Utils::MakePath($MDDir, 'OUTCAR', '/', 0));
	}

	my $MDOUTCar = Utils::MakePath($MDDir, 'OUTCAR', '/', 0);
	$App->print("  Read MD structures from [$MDOUTCar].\n");
	my @Crystals = $vasp->ReadStructuresFromOUTCAR($MDOUTCar, 0);
	if(@Crystals == 0) {
		$App->print("Error: Structure relax data is not found in [$MDOUTCar].\n");
		return 0;
	}

	my $nCrystal = @Crystals;
$App->print("nCrystal: $nCrystal\n");
	my @AtomSiteList = $Crystals[0]->GetCExpandedAtomSiteList();
	my $nAtomSite    = @AtomSiteList;

	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}//;
	$POSCARHash->{FormulaUnit} = 1 if(defined $POSCARHash->{FormulaUnit});
	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 '');

	for(my $i = 0 ; $i < $nCrystal ; $i++) {
		my $Crystal   = $Crystals[$i];
		my @AtomSite  = $Crystal->GetCAsymmetricAtomSiteList();
		my $nAtomSite = @AtomSite;
		my ($a, $b, $c, $alpha, $beta, $gamma) = $Crystal->LatticeParameters();
		my $Volume  = $Crystal->Volume();
		my $Density = $Crystal->Density();

		my %v = (
		index            => $ExistedLines + $i+1,
		FileBody         => $filebody,
		Compound         => $Compound,
		SortedComposition      => $SortedComposition,
		SumChemicalComposition => $POSCARHash->{SumChemicalComposition},
		FormulaUnit      => $POSCARHash->{FormulaUnit},
		Pexternal        => $Crystal->{Pexternal},
		PPullay          => $Crystal->{PPullay},
		Pkin             => $Crystal->{Pkin},
		P11              => $Crystal->{P11},
		P12              => $Crystal->{P12},
		P13              => $Crystal->{P13},
		P21              => $Crystal->{P21},
		P22              => $Crystal->{P22},
		P23              => $Crystal->{P23},
		P31              => $Crystal->{P31},
		P32              => $Crystal->{P32},
		P33              => $Crystal->{P33},
		Ekin             => $Crystal->{Ekin},
		TMD              => $Crystal->Temperature(),
		TOTEN            => $Crystal->{TOTEN},
		TOTENwithEntropy => $Crystal->{TOTENwithEntropy},
		'Etot/mol(SCF)'  => $Eavrg,
		'EF(SCF)'        => $SCFOUTCARHash->{EF},
		'EF(DOS)'        => $EF,
		VBM              => $VBM,
		CBM              => $CBM,
		Eg               => $Eg,
		Density          => $Density,
		a                => $a,
		b                => $b,
		c                => $c,
		alpha            => $alpha,
		beta             => $beta,
		gamma            => $gamma,
		Volume           => $Volume,
		TEBEG            => $INCARHash->{TEBEG},
		TEEND            => $INCARHash->{TEEND},
		PseudoPotential  => $POTCARHash->{PseudoPotential},
		PREC             => $INCARHash->{PREC},
		'aKmin(SCF)'     => $aKmin,
		'aKmax(SCF)'     => $aKmax,
		'aKmin(VCRelax)' => $aKminVCRelax,
		EDIFF            => $INCARHash->{EDIFF},
		EDIFFG           => $INCARHash->{EDIFFG},
		NBANDS           => $INCARHash->{NBANDS},
		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,
		);

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

		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;
}