#!/usr/bin/perl
##!c:/Perl/bin/perl

BEGIN {
#use lib 'c:/Programs/Perl/lib';
#use lib '/home/tkamiya/bin/lib';
my $BaseDir = $ENV{'TkPerlDir'};
#print "\nBaseDir: $BaseDir\n";
@INC = ("$BaseDir/lib", "$BaseDir/VNL", "c:/Programs/Perl/lib", @INC);
}

use strict;
#use warnings;
#use CGI::Carp qw(fatalsToBrowser);

use Deps;
use Utils;
use Sci qw($kB $JToeV);
use Sci::Optimize;
use MyHTMLApplication;
use Crystal::VASP;

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

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

my $T  = 300.0;   # EF
my ($Tmin, $Tmax, $Tstep) = (220.0, 300.0, 20.0);
my $ND = 0.103e22; # cm-3
my $ED = 0.4075;    # from 0 eV of DOS
my $NA = 0.0e20; # cm-3
my $EA = 0.43;    # from 0 eV of DOS

my $nMaxIter = 300;
my $xEPS     = 1.0e-15;
my $SEPS     = 1.0e3;
my $IsPrint  = 1;
#my $EF0 = 0.1; # eV
#my ($Eint0, $Eint1) = (-5.0, 8.0); # eV

my $CARDir = "SnS-Pnma-PBE";
$CARDir = $ARGV[0] if(defined $ARGV[0]);

my $OutDOSCSV = 'dos.csv';
my $OutputCSV = 'out.csv';

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

$App->SetOutputMode("console");


print("\nCalculate Ne-T from DFT TDOS\n");
my $vasp = new VASP;
my  ($EF, $VBM, $CBM, $Eg, $pDOSupArray, $pDOSdnArray,
		$INCARHash, $KPOINTSHash, $VCRelaxKPOINTSHash, $POTCARHash, $POSCARHash, $SCFOUTCARHash, $VCROUTCARHash)
			= $vasp->ReadVASPResults($CARDir);
if($EF eq 'Error') {
	$App->print("Error in VASP::ReadVASPResults: Abort due to [$VBM].\n");
	exit;
}

my $nData = $pDOSupArray->nData();
my $title = $pDOSupArray->Title();
my $pEDOS = $pDOSupArray->GetDataArray("x0");
my $pTDOS;
#last if(not defined $pX0);
my $nData = @$pEDOS;
my $xlabel = $pDOSupArray->{"x0_Name"};
print("Title: $title\n");
print("    X: $xlabel\n");
print("nData in DOS(up): $nData\n");
for(my $id = 0 ; ; $id++) {
	my $pY0   = $pDOSupArray->GetDataArray("y$id");
	last if(not defined $pY0);
	my $label0 = $pDOSupArray->{"y${id}_Name"};
print("    Y: $label0\n");
	if($label0 eq 'DOS') { # DOS, Interstitial, S total, Sn total, S s, S p, Sn s, Sn p, Sn, d etc
		$pTDOS = $pY0;
		last;
	}
}

for(my $i = 0 ; $i < @$pEDOS ; $i++) {
	$pTDOS->[$i] /= $POSCARHash->{Volume} * 1.0e-24;
}

# $pEDOS is measured from EF
print("EF(DOS OUTCAR): $vasp->{EFDOS} eV\n");
$VBM -= $vasp->{EFDOS};
$CBM -= $vasp->{EFDOS};
print("EVBM after EF correction: $VBM eV\n");
print("ECBM after EF correction: $CBM eV\n");

my $dEDOS = $pEDOS->[1] - $pEDOS->[0];
my $pNe = Algorism::IntegrateByTrapezoid($pEDOS, $pTDOS);

my $Emid = ($CBM + $VBM) / 2.0;
my $idxEmid = ($Emid - $pEDOS->[0]) / $dEDOS;

my $EF   = $Emid;
my $Ne0K = Algorism::Interpolate($pEDOS, $pNe, $EF);
print("EF (mid gap): $EF eV\n");
print("Ne(EF) at 0K: $Ne0K cm-3\n");

my @DOSFDe;
my @DOSFDh;
for(my $i = 0 ; $i < @$pEDOS ; $i++) {
	my $E = $pEDOS->[$i];
	my $FDe = Sci::FDTeV($T, $EF, $E);
	my $FDh = 1.0 - $FDe;
	$DOSFDe[$i] = $pTDOS->[$i] * $FDe;
	$DOSFDh[$i] = $pTDOS->[$i] * $FDh;
}
my $pNeFDe = Algorism::IntegrateByTrapezoid($pEDOS, \@DOSFDe, $idxEmid, undef);
my $pNeFDh = Algorism::IntegrateByTrapezoid($pEDOS, \@DOSFDh, 0, $idxEmid);

$App->print("\n");
$App->print("Write DOSup to [$OutDOSCSV].\n");

my $out = JFile->new($OutDOSCSV, 'w') or die "$!: Can not write to [$OutDOSCSV].\n";
$out->print("E(eV), total DOS (cm-3eV-1), Ne(cm-3), DOS*FDe, DOS*FDh, Ne*FDe, Ne*FDh\n");
for(my $i = 0 ; $i < @$pEDOS ; $i++) {
	$out->print("$pEDOS->[$i],$pTDOS->[$i],$pNe->[$i],$DOSFDe[$i],$DOSFDh[$i],$pNeFDe->[$i],$pNeFDh->[$i]\n");
}
$out->Close();

my ($Ne, $Nh, $NDp, $NAm, $TotalCharge) = &CalTotalCharge($T, $EF, $pEDOS, $pTDOS);
print("Ne  = $Ne cm-3\n");
print("Nh  = $Nh cm-3\n");
print("ND+ = $NDp cm-3\n");
print("NA- = $NAm cm-3\n");
print("Total charge: $TotalCharge cm-3\n");

#=====================
# Cal EF
#=====================
print("\n");
print("Optimize EF for T = $T K\n");
print("  Criteria: dEF=$xEPS   d(Total Charge)=$SEPS\n");

my $out = JFile->new($OutputCSV, 'w') or die "$!: Can not write to [$OutDOSCSV].\n";
$out->print("T(K),1000/T(K-1),EF(eV),Ne(cm-3),log(Ne),Nh(cm-3),ND+(cm-3),NA-(cm-3),total charge(cm-3)\n");

for(my $t = $Tmin ; $t <= $Tmax + 0.1 ; $t += $Tstep) {
	my $T = $t;
	my $T1000 = 1000.0 / $T;
	my ($EF, $S) = &CalEF($T, $nMaxIter, $xEPS, $SEPS, $IsPrint, $VBM, $CBM, $pEDOS, $pTDOS);
	my ($Ne, $Nh, $NDp, $NAm, $TotalCharge) = &CalTotalCharge($T, $EF, $pEDOS, $pTDOS);
	my $logNe = log($Ne) / log(10.0);
	printf("T=%4.1f K  EF=%6.4f eV: Ne=%8.4e  Nh=%8.4e  Ntot=%8.4e cm-3\n", $T, $EF, $Ne, $Nh, $TotalCharge);
	$out->print("$T,$T1000,$EF,$Ne,$logNe,$Nh,$NDp,$NAm,$TotalCharge\n");
}

$out->Close();

exit;

#===============================================
# Subroutines
#===============================================
sub CalTotalCharge
{
	my ($T, $EF, $pEDOS, $pTDOS) = @_;

	my $nData = @$pEDOS;
	my $dEDOS = $pEDOS->[1] - $pEDOS->[0];

	my @DOSFDe;
	my @DOSFDh;
	for(my $i = 0 ; $i < @$pEDOS ; $i++) {
		my $E = $pEDOS->[$i];
		my $FDe = Sci::FDTeV($T, $EF, $E);
		my $FDh = 1.0 - $FDe;
		$DOSFDe[$i] = $pTDOS->[$i] * $FDe;
		$DOSFDh[$i] = $pTDOS->[$i] * $FDh;
	}
	my $pNeFDe = Algorism::IntegrateByTrapezoid($pEDOS, \@DOSFDe, $idxEmid, undef);
	my $pNeFDh = Algorism::IntegrateByTrapezoid($pEDOS, \@DOSFDh, 0, $idxEmid);

	my $Ne = $pNeFDe->[$nData-1];
	my $Nh = $pNeFDh->[$nData-1];
	my $NDp = $ND * (1.0 - Sci::FDTeV($T, $EF, $ED));
	my $NAm = $NA * Sci::FDTeV($T, $EF, $EA);
	my $TotalCharge = $Ne + $NAm - $Nh - $NDp;

#print("EF = $EF eV\n");
#print("$TotalCharge = $Ne + $NAm - $Nh - $NDp cm-3\n");

	return ($Ne, $Nh, $NDp, $NAm, $TotalCharge);
}

sub CalS
{
	my ($T, $pEDOS, $pTDOS, $EF) = @_;

	my ($Ne, $Nh, $NDp, $NAm, $TotalCharge) = &CalTotalCharge($T, $EF, $pEDOS, $pTDOS);
#print("EF = $EF eV  Total Charge=$TotalCharge cm3\n");

	return $TotalCharge;
}

# 二分法でフェルミ準位を計算
sub CalEF
{
	my ($T, $nMaxIter, $xEPS, $SEPS, $IsPrint, $Ev, $Ec, $pEDOS, $pTDOS) = @_;

	my $Emax = 30.0 * $kB * $T * $JToeV;

	my $opt = new Optimize;
	my ($pVars, $S) = $opt->BinarySearch(
					$Ev-1.0, $Ec+5.0, sub { &CalS($T, $pEDOS, $pTDOS, @_); },
					$nMaxIter, $xEPS, $SEPS, $IsPrint
				);
	my $EF = $pVars->[0];
	return ($EF, $S);
}


1;