#!/usr/bin/perl

use strict;

use lib 'c:/Programs/Perl/lib';
use lib 'd:/Programs/Perl/lib';
use lib "$ENV{TkPerlDir}/lib";

use Deps;
use Utils;
use JFile;
use CSV;
use Sci qw($e $me $pi $kB $h $hbar $KToeV $JToeV erfc);
use Sci::Algorism;
use Crystal::VASP;

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

#===============================================
# 大域変数
#===============================================

# Which Ne used to determin 0 K EF
#   NELECT     : Use NELECT in OUTCAR
#   GivenEF    : Use EF in OUTCAR
#   SpecifiedEF: Use EF specified by $SpecifiedEF;
#my $NeSelect = 'NELECT';
my $NeSelect = 'GivenEF';
#my $NeSelect = 'SpecifiedEF';
my $SpecifiedEF = 0.1; # eV


my ($Emin, $Emax, $nE) = (-10.0, +10.0, 1001);
my $Estep = ($Emax - $Emin) / ($nE - 1);

my ($Tmin, $Tmax, $Tstep) = (0.0, 1200.0, 100.0);
my $gEPS = 1.0e-4;

my ($dEFmin, $dEFmax, $EFEPS, $NeEPS) = (-1.0, 1.0, 1.0e-6, 1.0e-4);

# EPS of Ne error to determine EF
my $NeEPS = 1.0e-4; # eV

# Tempearture to output f(E,T)
my $Temp = 1000.0; # K

my $CARDir = ".";
$CARDir = $ARGV[0] if(defined $ARGV[0]);
$CARDir = 'ZnO' if($CARDir eq '.');

my $INCAR     = Deps::MakePath($CARDir, 'INCAR',  0);
my $POSCAR    = Deps::MakePath($CARDir, 'POSCAR', 0);
my $POTCAR    = Deps::MakePath($CARDir, 'POTCAR', 0);
my $OUTCAR    = Deps::MakePath($CARDir, 'OUTCAR', 0);
#my $IBZKPT    = Deps::MakePath($CARDir, 'IBZKPT', 0);
my $EIGENVAL  = Deps::MakePath($CARDir, 'EIGENVAL', 0);
my $DOSCSV    = Deps::MakePath($CARDir, "ExactDOS.csv", 0);
my $ThermlCSV = Deps::MakePath($CARDir, "ThermalProperties.csv", 0);

print "\n";
print "Calculate EF(T), dEele(T), Sele(T), dEele(T) - TSele(T)\n";
print "\n";

print "Read files from [$CARDir]\n";
print "INCAR   : $INCAR\n";
print "POSCAR  : $POSCAR\n";
print "POTCAR  : $POTCAR\n";
print "OUTCAR  : $OUTCAR\n";
#print "IBZKPT  : $IBZKPT\n";
print "EIGENVAL: $EIGENVAL\n";
print "DOS CSV : $DOSCSV\n";
print "\n";

my $vasp = new VASP;

my $Crystal;
my ($a, $b, $c, $alpha, $beta, $gamma);
my ($V, $RV);
if(-f $INCAR and -f $POSCAR and -f $POTCAR) {
	$Crystal = $vasp->ReadStructureFromCARFiles($CARDir, 0);
	($a, $b, $c, $alpha, $beta, $gamma) = $Crystal->LatticeParameters();
	$V = $Crystal->Volume();
	$RV = $Crystal->ReciprocalVolume();
	print "Lattice parameters: $a, $b, $c A   $alpha, $beta, $gamma degrees\n";
	print "  V  = $V A^3\n";
	print "  RV = $RV A^-3\n";
	print "\n";
}

my $pOUT   = $vasp->ReadOUTCARtoHash($OUTCAR);
my $NELECT = $pOUT->{NELECT};
my $EF     = $pOUT->{EF};
#my $EF = $vasp->ReadFermiEnergy($OUTCAR);
print "from $OUTCAR:\n";
print "  NELECT = $NELECT\n";
print "  EF     = $EF eV\n";
print "\n";

#my ($nK, $pkArray) = $vasp->ReadIBZKPTToArray($IBZKPT);
#print "from $IBZKPT:\n";
#print "  nk = $nK\n";
#for(my $i = 0 ; $i < $nK ; $i++) {
#	my $p = $pkArray->[$i];
#	print "  k[$i] ($p->{kx}, $p->{ky}, $p->{kz}) w = $p->{w}\n";
#}
#print "\n";

print "Calculate exact D(E) / Ne(E)\n";
print "  Erange: $Emin - $Emax at $Estep step for $nE points\n";
print "  Note: Fermi level is adjusted to 0.0 eV hereafter.\n";

my ($Dunder, $Dover, $pD, $pNe) = &CalculateExactDOS($EIGENVAL, $EF, $Emin, $Emax, $Estep, $nE);
print "  DOS < $Emin = $Dunder\n";
print "  DOS > $Emax = $Dover\n";

my $Ne0 = &CalNe($Dunder, $pD, 0.0, 0.0);
print "  Ne up to  0.0 eV @   0 K = $Ne0\n";
if(0) {
	$Ne0 = &CalNe($Dunder, $pD, 0.0, -1.0);
	print "  Ne up to -1.0 eV @   0 K = $Ne0\n";
	$Ne0 = &CalNe($Dunder, $pD, 0.0,  1.0);
	print "  Ne up to  1.0 eV @   0 K = $Ne0\n";
	$Ne0 = &CalNe($Dunder, $pD, 100.0, 0.0);
	print "  Ne up to  0.0 eV @ 100 K = $Ne0\n";
	$Ne0 = &CalNe($Dunder, $pD, 300.0, 0.0);
	print "  Ne up to  0.0 eV @ 300 K = $Ne0\n";
	$Ne0 = &CalNe($Dunder, $pD, 600.0, 0.0);
	print "  Ne up to  0.0 eV @ 600 K = $Ne0\n";
}
print "\n";

print "NeSelect = $NeSelect\n";
if($NeSelect eq 'NELECT') {
	$Ne0 = $NELECT;
	print "  Ne0 = $Ne0 given by NELECT is used for following calculations.\n";
}
elsif($NeSelect eq 'GivenEF') {
#	$Ne0 = $Ne0;
	print "  Ne0 = $Ne0 at given EF = 0.0 eV is used for following calculations.\n";
}
elsif($NeSelect eq 'SpecifiedEF') {
	$Ne0 = &CalNe($Dunder, $pD, 0.0, $SpecifiedEF);
	print "  Ne0 = $Ne0 at specified EF = $SpecifiedEF eV is used for following calculations.\n";
}
print "\n";

print "  Save exact DOS(E) and Ne(E) to [$DOSCSV]\n";
my $out = JFile->new($DOSCSV, 'w') or die "Can not write to [$DOSCSV].\n";
$out->print("E(eV),D_i,Ne_i,fe($Temp K),fh($Temp K),Dunder=$Dunder,Dover=$Dover\n");
for(my $i = 0 ; $i < $nE ; $i++) {
	my $Ei = $Emin + $i * $Estep;
	my $fe = &FEe($Ei, 0.0, $Temp);
	my $fh = &FEh($Ei, 0.0, $Temp);
	$out->print("$Ei,$pD->[$i],$pNe->[$i],$fe,$fh\n");
}
$out->Close();
print "\n";

print "Calculate thermal properties\n";
my $Eele0;
print "  Save thermal properties to [$ThermlCSV]\n";
my $tcsv = JFile->new($ThermlCSV, 'w') or die "Can not write to [$ThermlCSV].\n";
$tcsv->print("T(K),EF(eV),dNe,Eele(eV),dEele(eV),Sele(eV/K),T*Sele(eV),dFele(eV)\n");

for(my $T = $Tmin ; $T <= $Tmax + $gEPS ; $T += $Tstep) {
print "  T = $T K\n";
	my $NeT = &CalNe($Dunder, $pD, $T, $EF);
#print " Ne(T=$T, EF=$EF) = $NeT\n";
#print "EF range: $EF+$dEFmin, $EF+$dEFmax\n";
	my ($EFT, $dNe) = &CalEF($Dunder, $pD, $T, $Ne0, $dEFmin, $dEFmax, $EFEPS);

	my $Eele = &CalEele($pD, $T, $EFT);
	$Eele0 = $Eele if(!defined $Eele0);
	my $dEele = $Eele - $Eele0;

	my $Sele  = &CalSele($pD, $T, $EFT);
	my $TSele = $T * $Sele;

	my $dFele = $dEele - $TSele;

	printf "    EF    = %12.8f eV (dNe = %g)\n", $T, $EFT, $dNe;
	printf "    Eele  = $Eele eV (E >= $Emin eV)\n";
	printf "    dEele = $dEele eV\n";
	printf "    Sele  = $Sele eV/K   TSele = $TSele eV\n";
	printf "    dFele = $dFele eV\n";
	$tcsv->print("$T,$EFT,$dNe,$Eele,$dEele,$Sele,$TSele,$dFele\n");
}
$tcsv->Close();

#print "log(10)=", log(10.0), "\n";

exit;


#===============================
# functions
#===============================
sub CalSele
{
	my ($pD, $T, $EFT) = @_;

	my $Sele = 0.0;
	for(my $ie = 0 ; $ie < $nE ; $ie++) {
		my $Ei = $Emin + $ie * $Estep;
		my $fe = &FEe($Ei, $EF, $T);
		my $fh = 1.0 - $fe;
		next if($fe < 1.0e-10 or $fh < 1.0e-10);

		$Sele += $pD->[$ie] * ($fe * log($fe) + $fh * log($fh));
	}
	return -$kB * $Sele / $e;
}

sub CalEele
{
	my ($pD, $T, $EFT) = @_;

	my $Eele = 0.0;
	for(my $ie = 0 ; $ie < $nE ; $ie++) {
		my $Ei = $Emin + $ie * $Estep;
		$Eele += $Ei * $pD->[$ie] * &FEe($Ei, $EF, $T);
	}
	return $Eele;
}

sub CalNe
{
	my ($Dunder, $pD, $T, $EF) = @_;

	my $Ne = $Dunder;
	for(my $ie = 0 ; $ie < $nE ; $ie++) {
		my $Ei = $Emin + $ie * $Estep;
		$Ne += $pD->[$ie] * &FEe($Ei, $EF, $T);
	}
	return $Ne;
}

sub CalEF
{
	my ($Dunder, $pD, $T, $NELECT, $EFmin, $EFmax, $EFEPS) = @_;

	my $dNe0 = &CalNe($Dunder, $pD, $T, $EFmin) - $NELECT;
	my $dNe1 = &CalNe($Dunder, $pD, $T, $EFmax) - $NELECT;
#print "  dNe($EFmin) = $dNe0  dNe($EFmax) = $dNe1\n";
		if($dNe0 * $dNe1 > 0.0) {
			print "Error in CalEF: dNe($EFmin)=$dNe0 and dNe($EFmax)=$dNe1 has the same signs.\n";
			print "  Reconsider dEFmin and dEFmax applicable for binary search.\n";
			print "\n";
			exit;
		}

	my ($EFnew, $dNe0, $dNe1, $dNenew);
	my $EF0 = $EFmin;
	my $EF1 = $EFmax;
	my $nMaxIter = 100;
	for(my $i = 0; $i < $nMaxIter ; $i++) {
		$dNe0   = &CalNe($Dunder, $pD, $T, $EF0) - $NELECT;
		$dNe1   = &CalNe($Dunder, $pD, $T, $EF1) - $NELECT;
#print "T=$T  EFrange = ($EF0, $EF1)  dNe = ($dNe0, $dNe1)\n";
		if(abs($dNe0) < $NeEPS) {
			return ($EF0, $dNe0);
		}
		if(abs($dNe1) < $NeEPS) {
			return ($EF1, $dNe1);
		}

#		if($EF1 - $EF0 < $EFEPS and abs($dNe0) < $NeEPS and abs($dNe1) < $NeEPS) {
		if($EF1 - $EF0 < $EFEPS) {
			$EFnew = ($EF0 + $EF1) / 2.0;
			$dNenew = &CalNe($Dunder, $pD, $T, $EFnew) - $NELECT;
#exit;
			return ($EFnew, $dNenew);
		}

		$EFnew  = ($EF0 + $EF1) / 2.0;
		$dNenew = &CalNe($Dunder, $pD, $T, $EFnew) - $NELECT;
		if($dNe0*$dNenew < 0.0) {
			$EF1 = $EFnew;
		}
		else {
			$EF0 = $EFnew;
		}
	}
	print "Error in CalEF: Not converged after $nMaxIter iterations\n";
exit;
	return (undef, undef);
}

sub CalculateExactDOS
{
	my ($EIGENVAL, $EF, $Emin, $Emax, $Estep, $nE) = @_;

	my (@D, @Ne);
	my ($Dunder, $Dover) = (0.0, 0.0);

	my $in = JFile->new($EIGENVAL, 'r') or die "Can not read [$EIGENVAL].\n";
	my $line;
	for(my $i = 0 ; $i < 5 ; $i++) {
		$in->ReadLine();
	}
#print "line: $line\n";
	my $SampleName = $line;
	$SampleName =~ s/^\s*//;
	$SampleName =~ s/\s*$//;

	$line = $in->ReadLine();
	my ($nn, $nK, $nLevels) = Utils::Split("\\s+", $line);
print "    nn,nK,nLevels=$nn, $nK, $nLevels\n";

	for(my $ie = 0 ; $ie < $nLevels ; $ie++) {
		$D[$ie] = 0.0;
	}

	for(my $i = 0 ; $i < $nK ; $i++) {
		$in->ReadLine();
		$line = $in->ReadLine();
		my ($kx, $ky, $kz, $w) = Utils::Split("\\s+", $line);
		for(my $ie = 0 ; $ie < $nLevels ; $ie++) {
			$line = $in->ReadLine();
			my ($idx, $Ei, $occ) = Utils::Split("\\s+", $line);
			$Ei -= $EF;
			my $iE = int(($Ei - $Emin) / $Estep);
#print "i=$i ie=$ie   Ei = $Ei ($Emin - $Emax)  iE=$iE (< $nE)\n";
			if($iE < 0.0) {
				$Dunder += 2.0 * $w;
			}
			elsif($nE <= $iE) {
				$Dover += 2.0 * $w;
			}
			else {
				$D[$iE] += 2.0 * $w;
#print "   D[$iE] = $D[$iE]\n";
			}
		}
	}
	$in->Close();
	
	$Ne[0] = $Dunder;
	for(my $ie = 0 ; $ie < $nE ; $ie++) {
		$Ne[$ie+1] += $Ne[$ie] + $D[$ie];
	}
	
	return ($Dunder, $Dover, \@D, \@Ne);
}

# 電子のFermi-Dirac分布
sub FEe
{
	my ($E, $EF, $T) = @_;

	if($T == 0.0) {
		return 0.0 if($E > $EF);
		return 0.5 if($E == $EF);
		return 1.0;
	}
	return 1.0 / ( 1.0 + exp(($E-$EF) / $T / $KToeV) );
}

# 正孔のFermi-Dirac分布
sub FEh
{
	my ($E, $EF, $T) = @_;

	if($T == 0.0) {
		return 0.0 if($E < $EF);
		return 0.5 if($E == $EF);
		return 1.0;
	}
	return 1.0 / ( 1.0 + exp(($EF-$E) / $T / $KToeV) );
}