#!/usr/bin/perl

use strict;

use lib 'd:/Programs/Perl/lib';
use lib "$ENV{TkPerlDir}/lib";
use Deps;
use Utils;
use JFile;
use Sci qw($e $me $pi $h $hbar);
use Sci::Algorism;
use Crystal::VASP;


my $nSmooth = 11;
my $InputFile  = "EIGENVAL.band";
$InputFile  = "EIGENVAL" if(!-f $InputFile);
my $OutputFile = "EffectiveMassBand.csv";
my $POSCAR     = "POSCAR";
my $KPOINTS    = "KPOINTS";

my $vasp = new VASP;
my $pKPOINTSHash = $vasp->ReadKPOINTStoHash($KPOINTS);
print "Lattice mode in KPOINTS: $pKPOINTSHash->{LatticeMode}\n";
if($pKPOINTSHash->{LatticeMode} !~ /^R/i) {
	print "Error: Lattice mode [$pKPOINTSHash->{LatticeMode}] is not supported.\n";
	exit;
}

my $Crystal = $vasp->ReadStructureFromCARFiles($POSCAR, 1) or die "Error: Can not read [$POSCAR].\n";
my ($a, $b, $c, $alpha, $beta, $gamma) = $Crystal->LatticeParameters();
print "Latt: $a, $b, $c angstrom    $alpha, $beta, $gamma\n";
my $Volume = $Crystal->Volume();
print "Volume: $Volume A^3\n";

my $EIGENVAL = $InputFile; #$vasp->GetEIGENVALFileName($POSCAR);
print "Read eigenvalues from [$EIGENVAL].\n";

my $OUTCAR = $vasp->GetOUTCARFileName($POSCAR);

my $out = new JFile;
$out->Open($OutputFile, "w") or die "Error: Can not write to [$OutputFile].\n";

$vasp->ReadEIGENVAL($EIGENVAL) or die "Error: Can not read [$EIGENVAL].\n";

my $EnergyBandupArray = $vasp->{EnergyBandupArray};
my $nEnergyBand = $EnergyBandupArray->nEnergyBand();
print "nEnergyBand = $nEnergyBand\n";
my $iVBM;
my $nk;
for(my $il = 0 ; $il < $nEnergyBand ; $il++) {
	my $EB = $EnergyBandupArray->GetEnergyBand($il);
	$nk = $EB->nData();
	my $ik = 0;
	my ($kx, $ky, $kz, $dk, $e, $n) = $EB->GetData($ik);
	if($e > 0.0) {
		$iVBM = $il-1;
		last;
	}
}
my $iCBM = $iVBM + 1;
print "iVBM=$iVBM\n";
#my $VBM3 = $EnergyBandupArray->GetEnergyBand($iVBM-4);
#my $VBM2 = $EnergyBandupArray->GetEnergyBand($iVBM-2);
my $VBM3 = $EnergyBandupArray->GetEnergyBand($iVBM-3);
my $VBM2 = $EnergyBandupArray->GetEnergyBand($iVBM-1);
my $VBM1 = $EnergyBandupArray->GetEnergyBand($iVBM);
my $CBM1 = $EnergyBandupArray->GetEnergyBand($iCBM);
my $CBM2 = $EnergyBandupArray->GetEnergyBand($iCBM+1);

my $pBoundaryDistances = $VBM1->{BoundaryDistances};
my $pBoundaryPositions = $VBM1->{BoundaryPositions};
for(my $i = 0 ; $i < @$pBoundaryDistances ; $i++) {
	print "Boundary k $i: $pBoundaryDistances->[$i]\n";
}
for(my $i = 0 ; $i < @$pBoundaryPositions ; $i++) {
	print "Boundary Pos $i: $pBoundaryPositions->[$i]\n";
}

$out->print("kx,ky,kz,k,EVBM-2(eV),EVBM-1(eV),EVBM(eV),ECBM(eV),ECBM+1(eV),"
		   ."d2EVBM(-2)/dE2,d2EVBM(-1)/dE2,d2EVBM/dE2,d2ECBM/dE2,d2ECBM(+1)/dE2,"
		   ."mh*(-2),mh*(-1),mh*,me*,me*(+1)");
my ($kx, $ky, $kz, $dk, $EVBM3, $EVBM2, $EVBM1, $ECBM1, $ECBM2, $n);
my $iBoundary = 1;
my $NextBoundary = $pBoundaryDistances->[$iBoundary];
my (@line, @d, @ev3, @ev2, @ev1, @ec1, @ec2);
my $ic = 0;
my ($d2EVBM3dk2, $d2EVBM2dk2, $d2EVBM1dk2, $d2ECBM1dk2, $d2ECBM2dk2);
my $Kk = 2.0 * $pi / 1.0e-10;
for(my $ik = 0 ; $ik < $nk ; $ik++) {
	($kx, $ky, $kz, $dk, $EVBM3, $n) = $VBM3->GetData($ik);
	($kx, $ky, $kz, $dk, $EVBM2, $n) = $VBM2->GetData($ik);
	($kx, $ky, $kz, $dk, $EVBM1, $n) = $VBM1->GetData($ik);
	($kx, $ky, $kz, $dk, $ECBM1, $n) = $CBM1->GetData($ik);
	($kx, $ky, $kz, $dk, $ECBM2, $n) = $CBM2->GetData($ik);
	$kx += 0.0;
	$ky += 0.0;
	$kz += 0.0;
	if($ik == 0 or $d[$ic-1] != $dk) {
		$line[$ic] = "$kx,$ky,$kz,$dk,$EVBM3,$EVBM2,$EVBM1,$ECBM1,$ECBM2";
		$d[$ic]    = $dk;
		$ev3[$ic]  = $EVBM3;
		$ev2[$ic]  = $EVBM2;
		$ev1[$ic]  = $EVBM1;
		$ec1[$ic]  = $ECBM1;
		$ec2[$ic]  = $ECBM2;
		$ic++;
	}

	if($dk > $NextBoundary - 1.0e-3) {
		$out->print("\n");
		for(my $i = 0 ; $i < $ic ; $i++) {
print "ic=$i: $line[$i]\n";
		}
		for(my $i = 0 ; $i < $ic ; $i++) {
			$d2EVBM3dk2 = &Cald2Edk2($i, \@d, \@ev3);
			$d2EVBM2dk2 = &Cald2Edk2($i, \@d, \@ev2);
			$d2EVBM1dk2 = &Cald2Edk2($i, \@d, \@ev1);
			$d2ECBM1dk2 = &Cald2Edk2($i, \@d, \@ec1);
			$d2ECBM2dk2 = &Cald2Edk2($i, \@d, \@ec2);
			my $mhs3 = ($d2EVBM3dk2 == 0)? 0 : $hbar * $hbar / ($d2EVBM3dk2 * $e / $Kk / $Kk) / $me;
			my $mhs2 = ($d2EVBM2dk2 == 0)? 0 : $hbar * $hbar / ($d2EVBM2dk2 * $e / $Kk / $Kk) / $me;
			my $mhs1 = ($d2EVBM1dk2 == 0)? 0 : $hbar * $hbar / ($d2EVBM1dk2 * $e / $Kk / $Kk) / $me;
			my $mes1 = ($d2ECBM1dk2 == 0)? 0 : $hbar * $hbar / ($d2ECBM1dk2 * $e / $Kk / $Kk) / $me;
			my $mes2 = ($d2ECBM2dk2 == 0)? 0 : $hbar * $hbar / ($d2ECBM2dk2 * $e / $Kk / $Kk) / $me;
print "ic=$i: $d[$i]  $ev1[$i]  $ec1[$i]  $d2EVBM1dk2  $d2ECBM1dk2\n";
print "      m=$mhs1, $mes1\n";
			$out->print("$line[$i],$d2EVBM3dk2,$d2EVBM2dk2,$d2EVBM1dk2,$d2ECBM1dk2,$d2ECBM2dk2,$mhs3,$mhs2,$mhs1,$mes1,$mes2\n");
		}

		@d   = ();
		@ev3 = ();
		@ev2 = ();
		@ev1 = ();
		@ec1 = ();
		@ec2 = ();
		$ic  = 0;
		$line[$ic] = "$kx,$ky,$kz,$dk,$EVBM3,$EVBM2,$EVBM1,$ECBM1,$ECBM2";
		$d[$ic]  = $dk;
		$ev3[$ic] = $EVBM3;
		$ev2[$ic] = $EVBM2;
		$ev1[$ic] = $EVBM1;
		$ec1[$ic] = $ECBM1;
		$ec2[$ic] = $ECBM2;
		$ic++;
		$NextBoundary = $pBoundaryDistances->[++$iBoundary];
	}
}
$out->Close();

exit;

sub Cald2Edk2
{
	my ($i, $pd, $pe) = @_;
	my $n = @$pe;
	if($i == 0) {
		return (($pe->[1] - $pe->[0])/($pd->[1] - $pd->[0]) - ($pe->[0] - $pe->[1])/($pd->[1] - $pd->[0])) / ($pd->[1] - $pd->[0]);
	}
	elsif($i == $n-1) {
		return (($pe->[$i-1] - $pe->[$i])/($pd->[$i] - $pd->[$i-1]) - ($pe->[$i] - $pe->[$i-1])/($d[$i] - $d[$i-1])) / ($pd->[$i] - $pd->[$i-1]);
	}
	return (($pe->[$i+1] - $pe->[$i])/($pd->[$i+1] - $pd->[$i]) - ($pe->[$i] - $pe->[$i-1])/($pd->[$i] - $pd->[$i-1])) / ($pd->[$i+1] - $pd->[$i-1]) * 2;
}