#!/usr/bin/perl

use lib 'd:/Programs/Perl/lib';
#use lib 'd:/MyWebs/cgi-bin/lib';

use lib '/home/tkamiya/bin/lib';
use lib '/home/tkamiya/bin';

use strict;
#use warnings;

use Sci qw($hbar $me $e);# acos asin);
use JFile;
use Crystal::TightBinding;

#==============================================================
# General parameters
#==============================================================
my $CEk = $hbar * $hbar / 2.0 / $me * 1.0e20 / $e;

my $InputFile   = ($ARGV[0])? $ARGV[0] : 'TB.in';

#==============================================================
# Read files
#==============================================================
my $tb      = new TightBinding;
my $pParams = $tb->ReadInput($InputFile);
if($pParams->{Diagonalization} ne 'Math::MatrixReal' and !-f $pParams->{Diagonalization}) {
	print "Error: Can not find the diagonalization program [$pParams->{Diagonalization}].\n";
	exit;
}
if(!-f $pParams->{CIFFile}) {
	print "Error: Can not find the CIF file [$pParams->{CIFFile}].\n";
	exit;
}
if(!-f $pParams->{KPOINTSFile}) {
	print "Error: Can not find the KPOINTS file [$pParams->{KPOINTSFile}].\n";
	exit;
}

my ($n, $pKPoints) = $tb->ReadKPOINTS($pParams->{KPOINTSFile});
my $Crystal        = $tb->ReadCIF($pParams->{CIFFile});

#==============================================================
# Main routine
#==============================================================
print "#==========================================#\n";
print "# Free electron model                      #\n";
print "#==========================================#\n";
print "\n";

$tb->PrintInputInf();
my $krange = $tb->{pInput}{kRange};
my $me     = $tb->{pInput}{me};
my $UG     = $tb->{pInput}{UG};
$CEk = $CEk / $me;

print "# Crystal structure\n";
print "Crystal structure read from [$pParams->{CIFFile}]\n";
$tb->PrintCrystalInf();

#==============================================================
# Build Orbitals
#==============================================================
my @Orbitals;
my $c = 0;
print "# Build orbitals\n";
for(my $ikz = -$krange ; $ikz <= $krange ; $ikz++) {
for(my $iky = -$krange ; $iky <= $krange ; $iky++) {
for(my $ikx = -$krange ; $ikx <= $krange ; $ikx++) {
	$Orbitals[$c] = [ $ikx, $iky, $ikz, "$ikx $iky $ikz" ];
	$c++;
}
}
}
my $nOrb = $c;
print "nOrbitals: $nOrb\n";
for(my $io = 0 ; $io < $nOrb ; $io++) {
	printf "  %2d: ik=%2d %2d %2d\n", $io, $Orbitals[$io][0], $Orbitals[$io][1], $Orbitals[$io][2];
}
print "\n";
#exit;

#==============================================================
# Build k list
#==============================================================
my $pCalKList = $tb->BuildKList();
my @k         = @$pCalKList;
my $nCalK     = @k;

print  "# k points\n";
print  "nK(target)=$n\n";
printf "nK(adjusted)=$nCalK\n";
for(my $i = 0 ; $i < $nCalK ; $i++) {
	my ($kx0, $ky0, $kz0, $name0) = @{$k[$i]};
	printf "k[$i]: (%8.4f, %8.4f, %8.4f, %s)\n", $kx0, $ky0, $kz0, $name0;
}
print "\n";
#exit;

my $out = JFile->new($pParams->{BandFile}, 'w') or die "$!: Can not write to [$pParams->{BandFile}].\n";
$out->print("name,kx,ky,kz,k");
for(my $io = 0 ; $io < $nOrb ; $io++) {
#	$out->print(",E($Orbitals[$io][3])");
	$out->print(",E($io)");
}
$out->print("\n");

for(my $i = 0 ; $i < @k ; $i++) {
	my $p = $k[$i];
	my ($kx, $ky, $kz, $ktot, $name) = @$p;

	print "#================================================\n";
	print "# k = ($kx, $ky, $kz): $name\n";
	print "#================================================\n";
	print "\n";
	
	my (@HijR, @HijI);
	for(my $i1 = 0 ; $i1 < $nOrb ; $i1++) {
		for(my $i2 = 0 ; $i2 < $nOrb ; $i2++) {
			$HijR[$i1][$i2] = 0.0;
			$HijI[$i1][$i2] = 0.0;
		}
	}

	for(my $i1 = 0 ; $i1 < $nOrb ; $i1++) {
		my ($ikx1, $iky1, $ikz1, $name1) = @{$Orbitals[$i1]};
		my $k = $Crystal->GetReciprocalDistanceFromK($kx + $ikx1, $ky + $iky1, $kz + $ikz1, 0, 0, 0);
		$HijR[$i1][$i1] = $CEk * $k * $k;
		for(my $i2 = $i1+1 ; $i2 < $nOrb ; $i2++) {
			my ($ikx2, $iky2, $ikz2, $name2) = @{$Orbitals[$i2]};
			if(abs($ikx2-$ikx1) == 1 or abs($iky2-$iky1) == 1 or abs($ikz2-$ikz1) == 1) {
				$HijR[$i1][$i2] = $HijR[$i2][$i1] = $UG;
			}
		}
	}

	if($tb->{pInput}{PrintLevel} >= 0) {
		print "HijR($kx,$ky,$kz)=\n";
		$tb->PrintMatrix(\@HijR);
		print "\n";
#		print "HijI($kx,$ky,$kz)=\n";
#		$tb->PrintMatrix(\@HijI);
#		print "\n";
	}

	my $pEV = &Diagonalize(\@HijR, \@HijI);
	&PrintEigenValues($pEV);

	$out->print("$name,$kx,$ky,$kz,$ktot");
	for(my $io = 0 ; $io < $nOrb ; $io++) {
		$out->print(",$pEV->[$io]{E}");
	}
	$out->print("\n");
}
$out->Close();

exit;

sub Diagonalize
{
	my ($pHijR, $pHijI) = @_;

	my $pInput = $tb->{pInput};
	my $nOrb   = @$pHijR;

	my (@EV, @En, @VnR, @VnI);

	if($pInput->{Diagonalization} eq 'Math::MatrixReal') {
		for(my $i = 0 ; $i < $nOrb ; $i++) {
			for(my $j = $i ; $j < $nOrb ; $j++) {
				if(abs($pHijI->[$i][$j]) > 1.0e-3) {
					print "Error: Hij is not a real matrix: HijI[$i][$j] = $pHijI->[$i][$j]\n";
				}
			}
		}

		my $pHij = new Math::MatrixReal($nOrb, $nOrb);
		for(my $i = 0  ; $i < $nOrb ; $i++) {
			for(my $j = 0 ; $j < $nOrb ; $j++) {
				$pHij->assign($i+1, $j+1, $pHijR->[$j][$i]);
			}
		}

		my ($pEn, $pCnij) = $pHij->sym_diagonalize();

		for(my $i = 0  ; $i < $nOrb ; $i++) {
			$En[$i] = $pEn->element($i+1, 1);
			for(my $j = 0 ; $j < $nOrb ; $j++) {
				$VnR[$i][$j] = $pCnij->element($j+1, $i+1);
				$VnI[$i][$j] = 0.0;
			}
		}
	}
	else {
		my $MatrixIn  = 'matrix.in';
		my $MatrixOut = 'matrix.out';

		unlink($MatrixIn);
		unlink($MatrixOut);

		$tb->SaveMatrix($MatrixIn, $pHijR, $pHijI, $pInput->{DiagonalizationEPS});

		my $ret = system($pInput->{Diagonalization});
		if($ret != 0) {
			print "Error: Can not execute [$pInput->{Diagonalization}] with err=$ret.\n";
			exit;
		}

		my ($ii, $jj);
		my $in = JFile->new($MatrixOut, 'r') or die "$!: Can not read [$MatrixOut].\n";

		my $n = $in->ReadLine() + 0;
		$in->ReadLine();
		for(my $i = 0 ; $i < $n ; $i++) {
			($ii, $En[$i]) = Utils::Split("\\s+", $in->ReadLine());
if($i+1 != $ii) {
	print "Error in read [$MatrixOut]; index ", $i+1, " is different from the one in the file ii=$ii.\n";
	exit;
}
		}
		$in->ReadLine();
		for(my $i = 0 ; $i < $n ; $i++) {
			for(my $j = 0 ; $j < $n ; $j++) {
				($ii, $jj, $VnR[$j][$i], $VnI[$j][$i]) = Utils::Split("\\s+", $in->ReadLine());
if($i+1 != $ii or $j+1 != $jj) {
	print "Error in read [$MatrixOut]; index (", $i+1, ',', $j+1, ") is different from the one in the file ii=($ii,$jj).\n";
	exit;
}
#print "$i, $j ($VnR[$i][$j], $VnI[$i][$j])\n";
			}
		}
		$in->Close();
#exit;
	}

	for(my $i = 0 ; $i < $nOrb ; $i++) {
		$EV[$i] = {E => $En[$i], VR => $VnR[$i], VI => $VnI[$i]};
	}

	if($pInput->{SortEigenValues}) {
		@EV = sort { $a->{E} <=> $b->{E} } @EV;
	}

	return \@EV;
}


sub PrintEigenValues
{
	my ($pEV) = @_;

	my $pInput = $tb->{pInput};
	my $n      = @$pEV;

#	$Orbitals[$c] = [ $ikx, $iky, $ikz, "$ikx $iky $ikz" ];
	print "\n";
	print "E (eV): Cij";
	for(my $i = 0 ; $i < $n ; $i++) {
		printf "  %-10s", "$Orbitals[$i][3]";
	}
	print "\n";

	for(my $i = 0 ; $i < $n ; $i++) {
		my $p = $pEV->[$i];
		printf "  %8.4f:", $p->{E};
		for(my $j = 0 ; $j < $n ; $j++) {
			my $vR = $p->{VR}[$j];
#			my $vI = $p->{VI}[$j];
			$vR = 0.0 if(abs($vR) < 1.0e-6);
#			$vI = 0.0 if(abs($vI) < 1.0e-6);

			printf " %8.4f", $vR;
#			printf " (%8.4f,%8.4f)", $vR, $vI;
		}
		print "\n";
	}

	print "\n";
	print "Molecular orbitals:\n";
	my $EPS = ($pInput->{MOEPS} > 0.0)? $pInput->{MOEPS} : 0.01;
	for(my $i = 0 ; $i < $n ; $i++) {
		my $p = $pEV->[$i];
		printf "  E(%2i)=%8.4f: ", $i, $p->{E};
		for(my $j = 0 ; $j < $n ; $j++) {
			my $vR = $p->{VR}[$j];
#			my $vI = $p->{VI}[$j];
			my $v  = abs($vR);
#			my $v  = sqrt($vR*$vR + $vI*$vI);
			next if(abs($v) < $EPS);

			my $sv = sprintf("%6.2f", abs($v));
			Utils::DelSpace($sv);

			my $orb   = $Orbitals[$j][3];

			if($vR > 0.0) {
				printf "+%s%s", $sv, "[$orb]";
			}
			else {
				printf "-%s%s", $sv, "[$orb]";
			}
		}
		print "\n";
	}
	print "\n";
}