#============================================================
# SemiconductorFromDFT
#============================================================
package SemiconductorFromDFT;
use Common;
@ISA = qw(Common);

#公開したいサブルーチン
@EXPORT = qw();

use strict;

use Deps;
use Utils;
use JFile;

use Sci qw($kB $JToeV $eVToJ $e);
use Sci::Optimize;
use Crystal::VASP;
use Math::Matrix;
use Math::MatrixReal;

#================================================
# 静的メンバー関数
#================================================

#==================================================================
# 変数取得関数
#==================================================================


#==================================================================
# コンストラクタ、デストラクト
#==================================================================
sub new
{
	my ($module) = @_;
	my $this = {};
	bless $this;
	return $this;
}

sub DESTROY
{
	my $this = shift;
}

#==================================================================
# 一般メンバー関数
#==================================================================
sub pDefectHash
{
	my ($this) = @_;
	return $this->{pDefectHash};
}
sub pDefect
{
	my ($this, $name, $q) = @_;

	my $p = $this->{pDefectHash};
	return undef if(!$p);
	return undef if(!$p->{$name});

	return $p->{$name} if(!defined $q);
	return $p->{$name}{$q};
}

sub pOrderedNames
{
	my ($this, $name) = @_;

	my $p = $this->{pDefectHash};
	if(!defined $name) {
		return $p->{pOrderedNames};
	}
	return $p->{$name}{pOrderedQ};
}

sub HfatEF
{
	my ($this, $name, $q, $T, $EF, $IgnoreS) = @_;
	$IgnoreS = 0 if(!defined $IgnoreS);

	my $p = $this->pDefect($name, $q);
	if(!$p) {
		print("Error in SemiconductorFromDFT::HfatEF: Can not find Hf data for [$name][q=$q].\n");
		exit;
	}
	
# Hf(EF) = Hf0 + q(EF - E0) - dS*T = (Hf0 - qE0) + qEF - dS*T
	my $dG = $p->{Hf00} + $p->{q} * $EF; # in eV
	$dG -= $p->{dS} * $kB * $T * $JToeV if(!$IgnoreS);

	return $dG;
}

sub CalPartitionFunction
{
	my ($this, $site, $T, $EF, $IgnoreS) = @_;

	my $p = $this->pDefectHash();
	my $Z = 1.0; # for site with no defect
	foreach my $name (keys %$p) {
			next if($name =~ /p/);
		my $p2 = $p->{$name};
		foreach my $q (keys %$p2) {
			next if($q =~ /p/);
			my $p3 = $p2->{$q};

			next if($p3->{Site} ne $site);
			my $Hf = $this->HfatEF($name, $q, $T, $EF, $IgnoreS);
			$Z += $this->BoltzmannProbability($Hf, $T);
		}
	}
	return $Z;
}

sub BoltzmannProbability
{
	my ($this, $dG, $T) = @_; # $dG in eV
	return exp(-$dG * $eVToJ / ($kB * $T));
}

sub DefectDensity
{
	my ($this, $N0, $dG, $T, $Z) = @_; # $dG in eV

	$Z = 1.0 if(!defined $Z);
	my $B = $this->BoltzmannProbability($dG, $T);
	return $N0 * $B / $Z; # same unit as $N0
}

sub AllFindTransitionEnergies
{
	my ($this, $App, $Tdefect, $Emin, $Emax, $IgnoreSForMinHf) = @_;

	my @inf;

	my $pNames = $this->pDefectHash()->{pOrderedNames};
	for(my $i = 0 ; $i < @$pNames ; $i++) {
		my $name = $pNames->[$i];
		my $pInf = $this->FindTransitionEnergies(undef, $name, $Tdefect, $Emin, $Emax, $IgnoreSForMinHf);
		$inf[$i] = {
			name => $name,
			$pInf => $pInf
			};
	}
	return \@inf;
}

sub FindTransitionEnergies
{
	my ($this, $App, $name, $Tdefect, $Emin, $Emax, $IgnoreSForMinHf) = @_;

	my $pDefect = $this->pDefectHash($name);
#print("  $name:\n");
	my $pQName  = $pDefect->{pOrderedQ};

	my @inf;
	my $idx = 0;

	my $EF = $Emin;
	my ($Minq, $MinHf) = $this->FindMinimumHf($App, $name, $EF, $Tdefect, $IgnoreSForMinHf);
#print "    EF=$EF: Min q=$Minq  Hf=$MinHf\n";

	$inf[$idx] = {
		E  => $EF,
		Hf => $MinHf,
		q  => $Minq,
		};
	$idx++;

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

#print "    Et=$Et => q=$Nextq  Hf=$Hft\n";
		$inf[$idx] = {
			E  => $Et,
			Hf => $Hft,
			qlower  => $Minq,
			qhigher => $Nextq
			};
		$idx++;

		$Minq = $Nextq;
		$EF   = $Et;

		last if($Et >= $Emax);
	}

	if($EF < $Emax) {
		($Minq, $MinHf) = $this->FindMinimumHf($App, $name, $Emax, $Tdefect, $IgnoreSForMinHf);
#print "    EF=$Emax: Min q=$Minq  Hf=$MinHf\n";
		$inf[$idx] = {
			E  => $Emax,
			Hf => $MinHf,
			q  => $Minq,
			};
	}
	
	return \@inf;
}

	
sub CalNeEF
{
	my ($this, $App, $Emin, $Emax, $nE, $T, $ND, $ED, $NA, $EA, $pEDOS, $pTDOS) = @_;

	my $Estep = ($Emax - $Emin) / ($nE - 1);

	my (@Ne, @Nh, @NDp, @NAm, @Ntot);
	for(my $i = 0 ; $i < $nE ; $i++) {
		my $EF    = $Emin + $i * $Estep;
		my $T1000 = 1000.0 / $T;

		$this->CalDistributionFunctions($App, $T, $EF);
		my ($Ne, $Nh, $NDp, $NAm, $TotalCharge) = $this->CalTotalCharge($App, $T, $EF, $ND, $ED, $NA, $EA, $pEDOS, $pTDOS);
		$Ne[$i]   = $Ne;
		$Nh[$i]   = $Nh;
		$NDp[$i]  = $NDp;
		$NAm[$i]  = $NAm;
		$Ntot[$i] = $TotalCharge;
	}
	
	return (\@Ne, \@Nh, \@NDp, \@NAm, \@Ntot);
}

sub CalHfNdEF
{
	my ($this, $App, $Emin, $Emax, $nE, $T, $pNe, $pNh, $IgnoreS) = @_;

	my $Estep = ($Emax - $Emin) / ($nE - 1);

	my @NeTotal;
	my $pNames = $this->pDefectHash()->{pOrderedNames};

	for(my $iE = 0 ; $iE < $nE ; $iE++) {
		$NeTotal[$iE] = $pNe->[$iE] - $pNh->[$iE];
	}

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

#$App->print("    q=$q:\n");
			my $site = $pDefectQ->{Site};
			my (@Hf, @Nd);
			for(my $iE = 0 ; $iE < $nE ; $iE++) {
				my $EF  = $Emin + $iE * $Estep;
				my $Hf  = $this->HfatEF($name, $q, $T, $EF, $IgnoreS);
				my $Z   = $this->CalPartitionFunction($site, $T, $EF, $IgnoreS);
				my $Nd  = $this->DefectDensity($pDefectQ->{N0}, $Hf, $T, $Z);
				$Hf[$iE] = $Hf;
				$Nd[$iE] = $Nd;
				$NeTotal[$iE] += $pDefectQ->{q} * $Nd[$iE];
			}
			$pDefectQ->{pHf} = \@Hf;
			$pDefectQ->{pNd} = \@Nd;
		}
	}
	
	return (\@NeTotal);
}

sub FindEFequilibrium
{
	my ($this, $App, $Emin, $Emax, $nE, $pNeTotal) = @_;

	my $Estep = ($Emax - $Emin) / ($nE - 1);

	my ($idx0, $idx1) = (0, $nE-1);
	my ($Ne0, $Ne1) = ($pNeTotal->[$idx0], $pNeTotal->[$idx1]);
	while(1) {
#print "idx=($idx0, $idx1) ($Ne0, $Ne1)\n";	
		if(abs($idx1 - $idx0) <= 1) {
			last;
		}

		my $idx2 = int(($idx0 + $idx1) / 2.0);
		my $Ne2  = $pNeTotal->[$idx2];
		
		if($Ne0 * $Ne2 == 0.0) {
			if($Ne2 == 0.0) {
				$idx0 = $idx1 = $idx2;
				last;
			}
			if($Ne0 == 0.0) {
				$idx1 = $idx0;
				last;
			}
		}
		elsif($Ne0 * $Ne2 < 0.0) {
			$idx1 = $idx2;
			($Ne0, $Ne1)   = ($pNeTotal->[$idx0], $pNeTotal->[$idx1]);
			next;
		}
		elsif($Ne1 * $Ne2 < 0.0) {
			$idx0 = $idx2;
			($Ne0, $Ne1)   = ($pNeTotal->[$idx0], $pNeTotal->[$idx1]);
			next;
		}
		else {
#$App->print("Error to find EF,e: Can not find Ne(total) = 0\n");
			return (undef);
		}
	}

	my $EF0 = $Emin + $idx0 * $Estep;
	my $EF1 = $Emin + $idx1 * $Estep;
	($Ne0, $Ne1) = ($pNeTotal->[$idx0], $pNeTotal->[$idx1]);
	my $EFe = $EF0 - $Ne0 * ($EF1 - $EF0) / ($Ne1 - $Ne0);
#$App->print("  Found: idx=($idx0, $idx1) EF=($EF0, $EF1) Ne=($Ne0, $Ne1)\n");

	return ($idx0, $idx1, $EF0, $EF1, $Ne0, $Ne1, $EFe);
}

sub FindMinimumHf
{
	my ($this, $App, $name, $EF, $T, $IgnoreS) = @_;

	my ($Minq, $MinHf);

	my $pDefect = $this->pDefect($name);
	foreach my $q (%$pDefect) {
		next if($q !~ /^[+-\d]+$/);
#print "q=$q\n";

#		my $pDefectQ = $pDefect->{$q};
		my $Hf  = $this->HfatEF($name, $q, $T, $EF, $IgnoreS);
		next if(!defined $Hf);

		if(!defined $Minq) {
			$Minq  = $q;
			$MinHf = $Hf;
		}
		elsif($MinHf > $Hf) {
			$Minq  = $q;
			$MinHf = $Hf;
		}
		
	}

	return ($Minq, $MinHf);
}

sub FindNextTransition
{
	my ($this, $App, $name, $q, $StartEF, $T, $IgnoreS) = @_;

	my $NextEt;
	my $Nextq;

	my $pDefect  = $this->pDefect($name);
	my $pDefectQ = $pDefect->{$q};
# Hf(EF) = Hf00 + q*EF
	my $Hf00 = $pDefectQ->{Hf00};

	foreach my $q2 (%$pDefect) {
		next if($q2 !~ /^[+-\d]+$/);
		next if($q2 >= $q);
#print "q=$q2\n";

		my $pDefectQ2 = $pDefect->{$q2};
		my $Hf002 = $pDefectQ2->{Hf00};

		my $Et = -($Hf002 - $Hf00) / ($q2 - $q);
#print "     q=$q2  Et=$Et\n";
		next if($Et < $StartEF);
		next if(defined $NextEt and $Et > $NextEt);

		$NextEt = $Et;
		$Nextq  = $q2;
#print "       =>$Nextq, $NextEt\n";
	}

	if(!defined $NextEt) {
		return (undef);
	}

	my $Hft  = $this->HfatEF($name, $Nextq, $T, $NextEt, $IgnoreS);

	return ($Nextq, $NextEt, $Hft);
}

sub ReadInputCSV
{
	my ($this, $InFile, $ReadToThis) = @_;
	$ReadToThis = 1 if(!defined $ReadToThis);

	my (%p, %d);
	my $in = JFile->new($InFile, 'r');
	if(!$in) {
#print "Error in SemiconductorFromDFT::ReadInputCSV: Can not read [$InFile]\n";
		return undef;
	}

	my $IsDefect = 0;
	while(1) {
		my $line = $in->ReadLine();
		last if(!defined $line);
		
		Utils::DelSpace($line);
		next if($line eq '');
		if($line =~ /^#DefectData/) {
			$IsDefect = 1;
			last;
		}
		next if($line =~ /^#/);

		my @a = split(/\s*,\s*/, $line);
		next if($a[0] eq '');
		$p{$a[0]} = $a[1];
	}

	my $HfKey    = $p{HfKey};
	my $idxHfKey = -1;
	my @DefectDataLabels;
	if($IsDefect) {
		my $line = $in->ReadLine();
		@DefectDataLabels = split(/\s*,\s*/, $line);
		for(my $i = 0 ; $i < @DefectDataLabels ; $i++) {
			if($DefectDataLabels[$i] eq $HfKey) {
				$idxHfKey = $i;
#print "found at i=$i\n";
#exit;
				last;
			}
		}
	}
	my @OrderedNames;
	my %NameAdded;
	while(1) {
		my $line = $in->ReadLine();
		last if(!defined $line);

		Utils::DelSpace($line);
		next if($line eq '');
		next if($line =~ /^#/);
		
		my @a = split(/\s*,\s*/, $line);
		next if($a[0] eq '');

#print "$a[$0]: $a[$1]\n";
		if(!$NameAdded{$a[0]}) {
			push(@OrderedNames, $a[0]);
			$NameAdded{$a[0]}++;
		}
		$d{$a[0]}{pOrderedQ} = [] if(!$d{$a[0]}{pOrderedQ});
		my $p = $d{$a[0]}{pOrderedQ};
		push(@$p, $a[1]);

		for(my $i = 0 ; $i < @a ; $i++) {
			$d{$a[0]}{$a[1]}{$DefectDataLabels[$i]} = $a[$i];
		}
		if($idxHfKey >= 0) {
			$d{$a[0]}{$a[1]}{Hf0} = $a[$idxHfKey];
		}
# Hf(EF) = Hf0 + q(EF - E0) - dS*T = (Hf0 - qE0) + qEF - dS*T
		$d{$a[0]}{$a[1]}{Hf00} = $d{$a[0]}{$a[1]}{Hf0} - $d{$a[0]}{$a[1]}{q} * $d{$a[0]}{$a[1]}{E0};
	}
	$d{pOrderedNames} = \@OrderedNames;

	$in->Close();

	if($ReadToThis) {
		$this->{pParameterHash}    = \%p;
		$this->{pDefectHash}       = \%d;
		$this->{pDefectDataLabels} = \@DefectDataLabels;
		return (\%p, \%d);
	}
	else {
		my %in;
		$this->{pParameterHash}    = $in{pParameterHash}    = \%p;
		$this->{pDefectHash}       = $in{pDefectHash}       = \%d;
		$this->{pDefectDataLabels} = $in{pDefectDataLabels} = \@DefectDataLabels;
		return (\%in, \%p, \%d);
	}
}

sub PrintParameters
{
	my ($this, $App, $pParameters, $pDefect) = @_;

	$App->print("Parameters\n");
	foreach my $key (sort keys %$pParameters) {
		$App->print("$key: $pParameters->{$key}\n");
	}
	$App->print("\n");
	$App->print("Defects\n");
	my $pDefect = $this->pDefectHash();
	if(!$pDefect) {
		$App->print("\nError in SemiconductorFromDFT::PrintParameters: Can not get pDefectHash\n");
		return;
	}
	my $pNames  = $pDefect->{pOrderedNames};
	if(!$pNames) {
		$App->print("\nError in SemiconductorFromDFT::PrintParameters: Can not get pOrderedNames\n");
		return;
	}
	for(my $i = 0 ; $i < @$pNames ; $i++) {
		my $name = $pNames->[$i];
		$App->print("  $name:\n");
		my $pDefect = $this->pDefect($name);
		my $pQName  = $pDefect->{pOrderedQ};
		for(my $j = 0 ; $j < @$pQName ; $j++) {
			my $q = $pQName->[$j];
			my $pDefectQ = $this->pDefect($name, $q);
#			my $pDefectQ = $p->{$name}{$q};
			next if(!$pDefectQ);

			$App->print("    q=$q:\n");
			foreach my $key (sort keys %$pDefectQ) {
				$App->print("      $key: $pDefectQ->{$key}\n");
			}
		}
	}
}

sub ReadDFTData
{
	my ($this, $App, $CARDir) = @_;

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

	$this->{nDOSData}   = $pDOSupArray->nData();
	$this->{DOSTitle}   = $pDOSupArray->Title();
	$this->{CellVolume} = $POSCARHash->{Volume};
	$this->{EFDOS}      = $vasp->{EFDOS};
	$VBM               -= $this->{EFDOS};
	$CBM               -= $this->{EFDOS};
	$this->{VBM}       = $VBM;
	$this->{CBM}       = $CBM;
#print("Title: $title\n");
#print("    X: $xlabel\n");
#print("nData in DOS(up): $nData\n");

	my $pEDOS = $this->{pEDOS} = $pDOSupArray->GetDataArray("x0");
	my $pTDOS;
#last if(not defined $pX0);
	my $nData = @$pEDOS;
	my $xlabel = $pDOSupArray->{"x0_Name"};
	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;
	}

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

	$this->{pTDOS} = $pTDOS;
	$this->{dEDOS} = $dEDOS;
	$this->{pNe}   = $pNe;

	return ($EF, $VBM, $CBM, $Eg, $pDOSupArray, $pDOSdnArray, $pEDOS, $pTDOS,
			$INCARHash, $KPOINTSHash, $VCRelaxKPOINTSHash, $POTCARHash, $POSCARHash, $SCFOUTCARHash, $VCROUTCARHash);
}

sub InitializeDOS
{
	my ($this, $App, $T, $EF) = @_;

	my $pEDOS = $this->{pEDOS};
	my $dEDOS = $pEDOS->[1] - $pEDOS->[0];
	my $pTDOS = $this->{pTDOS};
	my $pNe   = $this->{pNe};

	my $EF    = ($this->{CBM} + $this->{VBM}) / 2.0;
	my $idxEF = ($EF - $pEDOS->[0]) / $dEDOS;
	$idxEF = 0        if($idxEF < 0);
	$idxEF = $dEDOS-1 if($idxEF >= $dEDOS);

	my $Ne0K = Algorism::Interpolate($pEDOS, $pNe, $EF);

	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;
	}
	$this->{pNeFDe} = Algorism::IntegrateByTrapezoid($pEDOS, \@DOSFDe, $idxEF, undef);
	$this->{pNeFDh} = Algorism::IntegrateByTrapezoid($pEDOS, \@DOSFDh, 0, $idxEF);

	$this->{pDOSFDe} = \@DOSFDe;
	$this->{pDOSFDh} = \@DOSFDh;
	$this->{Ne0K}    = $Ne0K;
	
	return ($EF, $Ne0K);
}

sub CalDistributionFunctions
{
	my ($this, $App, $T, $EF) = @_;

	my $pEDOS = $this->{pEDOS};
	my $dEDOS = $pEDOS->[1] - $pEDOS->[0];
	my $pTDOS = $this->{pTDOS};
	my $pNe   = $this->{pNe};

	my $idxEF = ($EF - $pEDOS->[0]) / $dEDOS;
	$idxEF = 0        if($idxEF < 0);
	$idxEF = $dEDOS-1 if($idxEF >= $dEDOS);

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

#	$this->{pDOSFDe} = \@DOSFDe;
#	$this->{pDOSFDh} = \@DOSFDh;
	return ($this->{pDOSFDe}, $this->{pDOSFDh}, $this->{pNeFDe}, $this->{pNeFDh});
}

sub CalTotalCharge
{
	my ($this, $App, $T, $EF, $ND, $ED, $NA, $EA, $pEDOS, $pTDOS) = @_;

	$pEDOS = $this->{pEDOS} if(!defined $pEDOS);
	$pTDOS = $this->{pTDOS} if(!defined $pTDOS);

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

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

	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 CalSForEF
{
	my ($this, $App, $T, $ND, $ED, $NA, $EA, $pEDOS, $pTDOS, $EF) = @_;

	$this->CalDistributionFunctions($App, $T, $EF);
	my ($Ne, $Nh, $NDp, $NAm, $TotalCharge) = $this->CalTotalCharge($App, $T, $EF, $ND, $ED, $NA, $EA, $pEDOS, $pTDOS);
#print("EF = $EF eV  Total Charge=$TotalCharge cm3\n");

	return $TotalCharge;
}

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

	$pEDOS = $this->{pEDOS} if(!defined $pEDOS);
	$pTDOS = $this->{pTDOS} if(!defined $pTDOS);

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

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


1;