#!/usr/bin/perl
##!c:/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", "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;
use Sci::SemiconductorFromDFT;

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

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

# パラメータハッシュ
my ($pIn, $pDefect);

my $InFile = "in.csv";
$InFile = $ARGV[0] if(defined $ARGV[0]);

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

$App->SetOutputMode("console");

my $s = new SemiconductorFromDFT;

#===============================================
# ファイル読み込み
#===============================================
($pIn, $pDefect) = $s->ReadInputCSV($InFile);
$s->PrintParameters($App, $pIn, $pDefect);

$App->print("Action [$pIn->{Action}]\n");
if($pIn->{Action} eq 'Ne-T') {
	&CalNeT();
}
elsif($pIn->{Action} eq 'Hf-EF') {
	&CalHfEF();
}
elsif($pIn->{Action} eq 'PlotHf-EF') {
	&PlotHfEF();
}
else {
	$App->print("Error: No action found for [$pIn->{Action}]\n");
}

exit;


#===============================================
# Subroutines
#===============================================
sub PlotHfEF
{
	$App->print("Make Hf-EF data for plotting most stable charge states\n");

	my $Estep = ($pIn->{Emax} - $pIn->{Emin}) / ($pIn->{nE} - 1);
	my $pNames = $s->pDefectHash()->{pOrderedNames};

	for(my $i = 0 ; $i < @$pNames ; $i++) {
		my $name = $pNames->[$i];
$App->print("  $name:\n");
		my $pDefect = $s->pDefect($name);
		my $pQName  = $pDefect->{pOrderedQ};
$App->print("    q=");
		for(my $j = 0 ; $j < @$pQName ; $j++) {
			my $q = $pQName->[$j];
			my $pDefectQ = $s->pDefect($name, $q);
			next if(!$pDefectQ);
$App->print("$q, ");
		}
$App->print("\n");
	}

	my $out = JFile->new($pIn->{OutMinHfCSV}, 'w') or die "$!: Can not write to [$pIn->{OutMinHfCSV}].\n";

	for(my $i = 0 ; $i < @$pNames ; $i++) {
		my $name = $pNames->[$i];
$App->print("  $name:\n");
		$out->print("$name\n");
		$out->print("EF(eV),Hf(eV),q\n");
		my $pDefect = $s->pDefect($name);
		my $pQName  = $pDefect->{pOrderedQ};

		my $EF = $pIn->{Emin};
		my ($Minq, $MinHf) = $s->FindMinimumHf($App, $name, $EF, $pIn->{T}, $pIn->{IgnoreSForMinHf});
print "    EF=$EF: Min q=$Minq  Hf=$MinHf\n";
		$out->print("$EF,$MinHf,$Minq\n");

		while(1) {
			my ($Nextq, $Et, $Hft) = $s->FindNextTransition($App, $name, $Minq, $EF, $pIn->{T}, $pIn->{IgnoreSForMinHf});
			last if(!defined $Nextq);

print "    Et=$Et => q=$Nextq  Hf=$Hft\n";
			$out->print("$Et,$Hft,$Minq => $Nextq\n");
			$Minq = $Nextq;
			$EF   = $Et;

			last if($Et >= $pIn->{Emax});
		}

		if($EF < $pIn->{Emax}) {
			($Minq, $MinHf) = $s->FindMinimumHf($App, $name, $pIn->{Emax}, $pIn->{T}, $pIn->{IgnoreSForMinHf});
print "    EF=$pIn->{Emax}: Min q=$Minq  Hf=$MinHf\n";
			$out->print("$pIn->{Emax},$MinHf,$Minq\n");
		}

		$out->print("\n");
	}

	$out->Close();
}

sub CalHfEF
{
	$App->print("Read DFT Data from [$pIn->{CARDir}]\n");
	my  ($EF, $VBM, $CBM, $Eg, $pDOSupArray, $pDOSdnArray, $pEDOS, $pTDOS,
		 $INCARHash, $KPOINTSHash, $VCRelaxKPOINTSHash, $POTCARHash, $POSCARHash, $SCFOUTCARHash, $VCROUTCARHash)
				= $s->ReadDFTData($App, $pIn->{CARDir});

	$App->print("Title: $s->{DOSTitle}\n");
#	$App->print("    X: $xlabel\n");
	$App->print("nData in DOS(up): $s->{nDOSData}\n");

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

	my ($Emid, $Ne0K) = $s->InitializeDOS($App, $pIn->{T});
	$App->print("ED: $pIn->{ED} eV\n");
	$App->print("ND: $pIn->{ND} cm-3\n");
	$App->print("EA: $pIn->{EA} eV\n");
	$App->print("NA: $pIn->{NA} cm-3\n");

	$App->print("E(midgap)        : $s->{Emid} eV\n");
	$App->print("Ne(Emidgap) at 0K: $s->{Ne0K} cm-3\n");

#=====================
# Save to OutDOSCSV
#=====================
	$App->print("\n");
	$App->print("Write DOSup to [$pIn->{OutDOSCSV}].\n");

	my $out = JFile->new($pIn->{OutDOSCSV}, 'w') or die "$!: Can not write to [$pIn->{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],$s->{pNe}->[$i],$s->{pDOSFDe}->[$i],$s->{pDOSFDh}->[$i],$s->{pNeFDe}->[$i],$s->{pNeFDh}->[$i]\n");
	}
	$out->Close();

#===============================================
# Non-selfconsistent計算
#===============================================
	$App->print("\nNon-selfconsistent calculation for EF = midgap\n");

	$s->CalDistributionFunctions($App, $pIn->{T}, $Emid);
	my ($Ne, $Nh, $NDp, $NAm, $TotalCharge) 
		= $s->CalTotalCharge($App, $pIn->{T}, $Emid, $pIn->{ND}, $pIn->{ED}, $pIn->{NA}, $pIn->{EA}, $pEDOS, $pTDOS);

	$out->print("Ne(Emid) at 0K,$Ne0K,cm-3\n");
	$App->print("Ne  = $Ne cm-3\n");
	$App->print("Nh  = $Nh cm-3\n");
	$App->print("ND+ = $NDp cm-3\n");
	$App->print("NA- = $NAm cm-3\n");
	$App->print("Total charge: $TotalCharge cm-3\n");

	my $Estep = ($pIn->{Emax} - $pIn->{Emin}) / ($pIn->{nE} - 1);

#===============================================
# Calculate Ne-EF
#===============================================
	$App->print("\nCalculate Ne-EF from DFT TDOS\n");
	my ($pNe, $pNh, $pNDp, $pNAm, $pNtot) 
		= $s->CalNeEF($App, $pIn->{Emin}, $pIn->{Emax}, $pIn->{nE}, $pIn->{T}, 0.0, $pIn->{ED}, 0.0, $pIn->{EA}, $pEDOS, $pTDOS);

#===============================================
# Calculate Hf-Nd-EF
#===============================================
	$App->print("\nCalculate Hf-Defect density-EF from DFT TDOS\n");
	my $pNeTotal = $s->CalHfNdEF($App, $pIn->{Emin}, $pIn->{Emax}, $pIn->{nE}, $pIn->{T}, $pNe, $pNh, $pIn->{IgnoreS});

#===============================================
# Find EF @ equilibrium by Bisection method
#===============================================
	$App->print("\nFind EF @ equilibrium\n");

	my ($idx0, $idx1, $EF0, $EF1, $Ne0, $Ne1, $EFe) = $s->FindEFequilibrium($App, $pIn->{EMin}, $pIn->{Emax}, $pIn->{nE}, $pNeTotal);
	if(!defined $idx0) {
		$App->print("Error to find EF,e: Can not find Ne(total) = 0\n");
	}
	else {
		$App->print("  Found: idx=($idx0, $idx1) EF=($EF0, $EF1) Ne=($Ne0, $Ne1)\n");
		$App->print("  EF,e = $EFe eV\n");
	}

#===============================================
# Save Hf-EF
#===============================================
	$App->print("\nSave Hf-EF data\n");

	my $pNames = $s->pDefectHash()->{pOrderedNames};
	my $out = JFile->new($pIn->{OutHfCSV}, 'w') or die "$!: Can not write to [$pIn->{OutHfCSV}].\n";
$out->print("EF(eV)");
	for(my $i = 0 ; $i < @$pNames ; $i++) {
		my $name = $pNames->[$i];
#$App->print("  $name:\n");
		my $pDefect = $s->pDefect($name);
		my $pQName  = $pDefect->{pOrderedQ};
		for(my $j = 0 ; $j < @$pQName ; $j++) {
			my $q = $pQName->[$j];
			my $pDefectQ = $s->pDefect($name, $q);
			next if(!$pDefectQ);

#$App->print("    q=$q:\n");
			$out->print(",$name($q)");
		}
	}
	$out->print("\n");

	for(my $i = 0 ; $i < $pIn->{nE} ; $i++) {
		my $EF    = $pIn->{Emin} + $i * $Estep;
		$out->print("$EF");
		for(my $in = 0 ; $in < @$pNames ; $in++) {
			my $name    = $pNames->[$in];
			my $pDefect = $s->pDefect($name);
			my $pQName  = $pDefect->{pOrderedQ};
			for(my $iq = 0 ; $iq < @$pQName ; $iq++) {
				my $q = $pQName->[$iq];
				my $pDefectQ = $s->pDefect($name, $q);
				next if(!$pDefectQ);

				my $pHf = $pDefectQ->{pHf};
				$out->print(",$pHf->[$i]");
			}
		}
		$out->print("\n");
	}
	$out->Close();

#===============================================
# Save Ndefect-EF
#===============================================
	$App->print("\nSave Defect density-EF data\n");

	my $out = JFile->new($pIn->{OutDefectDensityCSV}, 'w') or die "$!: Can not write to [$pIn->{OutDefectDensityCSV}].\n";

	$out->print("EF(eV),Ne(total)(cm-3),log(Ne(total),Ne(cm-3),Nh(cm-3)");
	for(my $i = 0 ; $i < @$pNames ; $i++) {
		my $name = $pNames->[$i];
#$App->print("  $name:\n");
		my $pDefect = $s->pDefect($name);
		my $pQName  = $pDefect->{pOrderedQ};
		for(my $j = 0 ; $j < @$pQName ; $j++) {
			my $q = $pQName->[$j];
			my $pDefectQ = $s->pDefect($name, $q);
			next if(!$pDefectQ);

#$App->print("    q=$q:\n");
			$out->print(",$name($q)");
		}
	}
	$out->print("\n");

	for(my $i = 0 ; $i < $pIn->{nE} ; $i++) {
		my $EF    = $pIn->{Emin} + $i * $Estep;
		my $abs = abs($pNeTotal->[$i]);
		my $logNeTotal = 0;
		if($pNeTotal->[$i] > 0) {
			$logNeTotal = log($pNeTotal->[$i]) / log(10.0);
		}
		elsif($pNeTotal->[$i] < 0) {
			$logNeTotal = -log(-$pNeTotal->[$i]) / log(10.0);
		}
		$out->print("$EF,$pNeTotal->[$i],$logNeTotal,$pNe->[$i],$pNh->[$i]");
		for(my $in = 0 ; $in < @$pNames ; $in++) {
			my $name    = $pNames->[$in];
			my $pDefect = $s->pDefect($name);
			my $pQName  = $pDefect->{pOrderedQ};
			for(my $iq = 0 ; $iq < @$pQName ; $iq++) {
				my $q = $pQName->[$iq];
				my $pDefectQ = $s->pDefect($name, $q);
				next if(!$pDefectQ);

				my $pNd = $pDefectQ->{pNd};
				$out->print(",$pNd->[$i]");
			}
		}
		$out->print("\n");
	}

	$out->Close();

#===============================================
# Save Output & Ne-EF
#===============================================
	$App->print("\nSave Output, (Ne,Nh,ND+,NA-)-EF data\n");

	my $out = JFile->new($pIn->{OutputCSV}, 'w') or die "$!: Can not write to [$pIn->{OutDOSCSV}].\n";

	$out->print("EF(DOS OUTCAR),$s->{EFDOS},eV\n");
	$out->print("EVBM after EF correction,$VBM,eV\n");
	$out->print("ECBM after EF correction,$CBM,eV\n");
	$out->print("EF(equilibrium),$EFe,eV\n");
	$out->print("\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 $i = 0 ; $i < $pIn->{nE} ; $i++) {
		my $EF    = $pIn->{Emin} + $i * $Estep;
		my $T1000 = 1000.0 / $pIn->{T};
		my $logNe = log($pNe->[$i]) / log(10.0);
		printf("T=%4.1f K  EF=%6.4f eV: Ne=%8.4e  Nh=%8.4e  Ntot=%8.4e cm-3\n", $pIn->{T}, $EF, $pNe->[$i], $pNh->[$i], $pNtot->[$i]);
		$out->print("$pIn->{T},$T1000,$EF,$Ne,$logNe,$Nh,$NDp,$NAm,$TotalCharge\n");
	}
	$out->Close();
}

sub CalNeT
{
	$App->print("Read DFT Data from [$pIn->{CARDir}]\n");
	my  ($EF, $VBM, $CBM, $Eg, $pDOSupArray, $pDOSdnArray, $pEDOS, $pTDOS,
		 $INCARHash, $KPOINTSHash, $VCRelaxKPOINTSHash, $POTCARHash, $POSCARHash, $SCFOUTCARHash, $VCROUTCARHash)
				= $s->ReadDFTData($App, $pIn->{CARDir});

	$App->print("Title: $s->{DOSTitle}\n");
#	$App->print("    X: $xlabel\n");
	$App->print("nData in DOS(up): $s->{nDOSData}\n");

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

	my $Emid = ($VBM + $CBM) / 2.0;
	my ($Emid, $Ne0K) = $s->InitializeDOS($App, $pIn->{T}, $Emid);
	my $pNe     = $s->{pNe};
#	my $pDOSFDe = $s->{pDOSFDe};
#	my $pDOSFDh = $s->{pDOSFDh};
#	my $pNeFDe  = $s->{pNeFDe};
#	my $pNeFDh  = $s->{pNeFDh};

	$App->print("ED: $pIn->{ED} eV\n");
	$App->print("ND: $pIn->{ND} cm-3\n");
	$App->print("EA: $pIn->{EA} eV\n");
	$App->print("NA: $pIn->{NA} cm-3\n");

	$App->print("E(midgap)        : $s->{Emid} eV\n");
	$App->print("Ne(Emidgap) at 0K: $s->{Ne0K} cm-3\n");

#=====================
# Save to OutDOSCSV
#=====================
	$App->print("\n");
	$App->print("Write DOSup to [$pIn->{OutDOSCSV}].\n");

	my $out = JFile->new($pIn->{OutDOSCSV}, 'w') or die "$!: Can not write to [$pIn->{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],$s->{pDOSFDe}->[$i],$s->{pDOSFDh}->[$i],$s->{pNeFDe}->[$i],$s->{pNeFDh}->[$i]\n");
	}
	$out->Close();

#===============================================
# Non-selfconsistent計算
#===============================================
	$App->print("\nNon-selfconsistent calculation for EF = midgap\n");

	$s->CalDistributionFunctions($App, $pIn->{T}, $Emid);
	my ($Ne, $Nh, $NDp, $NAm, $TotalCharge) 
		= $s->CalTotalCharge($App, $pIn->{T}, $Emid, $pIn->{ND}, $pIn->{ED}, $pIn->{NA}, $pIn->{EA}, $pEDOS, $pTDOS);

	$App->print("Ne  = $Ne cm-3\n");
	$App->print("Nh  = $Nh cm-3\n");
	$App->print("ND+ = $NDp cm-3\n");
	$App->print("NA- = $NAm cm-3\n");
	$App->print("Total charge: $TotalCharge cm-3\n");

#===============================================
# EF-T計算
#===============================================
	$App->print("\nCalculate Ne-T from DFT TDOS\n");
	$App->print("Optimize EF for T = $pIn->{T} K\n");
	$App->print("  Criteria: dEF=$pIn->{xEPS}   d(Total Charge)=$pIn->{SEPS}\n");

	my $out = JFile->new($pIn->{OutputCSV}, 'w') or die "$!: Can not write to [$pIn->{OutDOSCSV}].\n";

	$out->print("EF(DOS OUTCAR),$s->{EFDOS},eV\n");
	$out->print("EVBM after EF correction,$VBM,eV\n");
	$out->print("ECBM after EF correction,$CBM,eV\n");
	$out->print("Ne(Emid) at 0K,$Ne0K,cm-3\n");
	$out->print("Ne,$Ne,cm-3\n");
	$out->print("Nh,$Nh,cm-3\n");
	$out->print("ND+,$NDp,cm-3\n");
	$out->print("NA-,$NAm,cm-3\n");
	$out->print("Total charge,$TotalCharge,cm-3\n");
	$out->print("\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 = $pIn->{Tmin} ; $t <= $pIn->{Tmax} + 0.1 ; $t += $pIn->{Tstep}) {
		my $T = $t;
		my $T1000 = 1000.0 / $T;
		my ($EF, $S) 
			= $s->CalEF($App, $T, $pIn->{nMaxIter}, $pIn->{xEPS}, $pIn->{SEPS}, $pIn->{IsPrint}, $VBM, $CBM, 
						$pIn->{ND}, $pIn->{ED}, $pIn->{NA}, $pIn->{EA}, $pEDOS, $pTDOS);
		my ($Ne, $Nh, $NDp, $NAm, $TotalCharge) 
			= $s->CalTotalCharge($App, $T, $EF, $pIn->{ND}, $pIn->{ED}, $pIn->{NA}, $pIn->{EA}, $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();

	return;
}

1;