#!/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 Math::MatrixReal;
use Math::Complex;

use Sci qw($pi $a0 $torad $todeg);# acos asin);
use Crystal::CIF;
use Crystal::Crystal;
use Sci::MyMatrixReal;

#==============================================================
# General parameters
#==============================================================
my $OutFile  = "TB1DBand.csv";
my @lOrbName = ('s', 'p', 'd', 'f');

my $pi2 = 2.0 * $pi;
my $h2m = 7.62; # eV A^2
my $e2  = 14.40; # eV A

my %eta = (
	'ssS' => -1.40,
	'spS' =>  1.84,
	'ppS' =>  3.24,
	'ppP' => -0.81,
	'sdS' => -3.16,
	'pdS' => -2.95,
	'pdP' =>  1.36,
	'ddS' => -16.2,
	'ddP' =>  8.75,
	'ddD' =>  0.0,
	);

#==============================================================
# Crystal structure
#==============================================================
my $a    = 4.0;
my @Site = (
	[ 'Si', 0.0, 0.0, 0.0 ],
	[ 'C',  0.5, 0.0, 0.0 ],
	);

#==============================================================
# Electronic parameters for atoms
#==============================================================
# n, l, m
my %AtomOribitals = (
	'Si' => [ [2, 0, 0], [2, 1, 0] ],
	'C'  => [ [2, 0, 0], [2, 1, 0] ],
	);

# e0
my %aii = (
	'C'  => { '2s' => -17.52, '2p' => -8.97 },
	'Si' => { '2s' => -13.55, '2p' => -6.52 },
	);

my %kF = (
	'C'  => 2.76,
	'Si' => 1.81,
	);

my %rc = (
	'C'  => 0.37,
	'Si' => 0.56,
	);

my %ri = (
#	'C'  => ,
	'Si' => 0.38,
	);

#==============================================================
# Calculation condition
#==============================================================
my $Rmax   = 3.0;
my $nrange = 1;

my ($k0, $k1, $nk) = (0.0, 0.5, 11);
my $kstep = ($k1 - $k0) / ($nk - 1);

#FCC Cartesian
my @KList = (
	[0.0,  0.0,  0.0, 'GAMMA'],
	[0.0,  0.0,  1.0, 'X'],
	[0.5,  0.0,  1.0, 'W'],
	[0.5,  0.5,  0.5, 'L'],
	[0.75, 0.75, 0.0, 'K'],
	[0.5,  0.0,  0.0, ''],
	);

#==============================================================
# Main routine
#==============================================================
my @HijIndex;
my @Orbitals;
for(my $isite = 0 ; $isite < @Site ; $isite++) {
	my $atom  = $Site[$isite][0];
	my $pAtom = $AtomOribitals{$atom};
	for(my $iorb = 0 ; $iorb < @$pAtom ; $iorb++) {
		my $n = $pAtom->[$iorb][0];
		my $l = $pAtom->[$iorb][1];
		push(@Orbitals, [$isite, $atom, "$n$lOrbName[$l]"]);
		$HijIndex[$isite][$iorb] = scalar @Orbitals - 1;
	}
}
my $nOrb = @Orbitals;

print "nOrbitals: $nOrb\n";
print "Hij index:\n";
for(my $io = 0 ; $io < $nOrb ; $io++) {
	print "  $io: iSite=$Orbitals[$io][0]  $Orbitals[$io][1] $Orbitals[$io][2]\n";
}

my $out = JFile->new($OutFile, 'w') or die "$!: Can not write to [$OutFile].\n";
$out->print("k");
for(my $io = 0 ; $io < $nOrb ; $io++) {
	$out->print(",E$io");
}
$out->print("\n");

for(my $ik = 0 ; $ik < $nk ; $ik++) {
	my $k = $k0 + $ik * $kstep;

	my $pHij = &CalHij($k);
print "\n";
print "Hij($k)=\n";
print $pHij;

	my ($En, $Cnij) = $pHij->sym_diagonalize();
print "\n";
print "Diagonarized\n";
print "En=\n";
print $En;
print "Cnij=\n";
print $Cnij;

	$out->print("$k");
	for(my $io = 0 ; $io < $nOrb ; $io++) {
		$out->print(",", $En->element($io+1,1));
	}
	$out->print("\n");
}

$out->Close();

exit;

#================================================================
# Subroutines
#================================================================
sub CalHij
{
	my ($k) = @_;

	my (@HijR, @HijI);
	for(my $i = 0  ; $i < $nOrb ; $i++) {
		for(my $j = 0 ; $j < $nOrb ; $j++) {
			$HijR[$j][$i] = 0.0;
			$HijI[$j][$i] = 0.0;
		}
	}

	for(my $is1 = 0  ; $is1 < @Site ; $is1++) {
		my $site1  = $Site[$is1];
		my $name1  = $site1->[0];
		my $x1     = $site1->[1];
		my $pAtom1 = $AtomOribitals{$name1};
		for(my $is2 = 0  ; $is2 < @Site ; $is2++) {
			my $site2  = $Site[$is2];
			my $name2  = $site2->[0];
			my $x2     = $site2->[1];
			my $pAtom2 = $AtomOribitals{$name2};
			for(my $ix = -$nrange ; $ix <= $nrange ; $ix++) {
				my $x2r = $x2 + $ix;
				my $r = $a * ($x2r - $x1);
				next if(abs($r) > $Rmax);

				for(my $io1 = 0 ; $io1 < @$pAtom1 ; $io1++) {
					my $n1       = $pAtom1->[$io1][0];
					my $l1       = $pAtom1->[$io1][1];
					my $OrbName1 = "$n1$lOrbName[$l1]";
					my $idx1     = $HijIndex[$is1][$io1];
					for(my $io2 = 0 ; $io2 < @$pAtom2 ; $io2++) {
						my $n2       = $pAtom2->[$io2][0];
						my $l2       = $pAtom2->[$io2][1];
						my $OrbName2 = "$n2$lOrbName[$l2]";
						my $idx2     = $HijIndex[$is2][$io2];
						next if($idx2 < $idx1);

						if($r == 0.0) {
							if($idx1 == $idx2) {
								$HijR[$idx1][$idx2] += $aii{$name1}{$OrbName1};
print "Hij[$idx1][$idx2] = aii{$name1}{$OrbName1}\n";
							}
						}
						else {
							my $ekrr = cos($pi2 * $k * ($x1 - $x2r));
							my $ekri = sin($pi2 * $k * ($x1 - $x2r));
							my ($VS, $VP, $VD) = &Getbij($r, $name1, $l1, $name2, $l2);
# Slater-Koster conversion
							my $lx  = ($x1 - $x2r) / $r;
							my $lx2 = $lx * $lx;
							if($l1 == 0 and $l2 == 0) {
								$b = $VS;
							}
							elsif(($l1 == 0 and $l2 == 1) or ($l1 == 1 and $l2 == 0)) {
								$b = $lx * $VS;
							}
							elsif(($l1 == 1 and $l2 == 1)) {
								$b = $lx2 * $VS + (1.0 - $lx2) * $VP;
							}
print "  ekr=($ekrr, $ekri)  VS=$VS  b=$b\n";
							$HijR[$idx1][$idx2] += $ekrr * $b;
							$HijI[$idx1][$idx2] += $ekri * $b;
print "HijR[$idx1][$idx2] = exp(i * $k * ($x1 - $x2r)) * bij{$name1 $OrbName1-$name2 $OrbName2} = $HijR[$idx1][$idx2]\n";
						}
print "HijR[$idx1][$idx2] <= r = $r [$x1 - $x2r] ($is1:$name1, $io1:$OrbName1), ($is2:$name2, $io2:$OrbName2, ix=$ix)\n";
					}
				}
			}
		}
	}

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

	return $pHij;
}

sub Getbij
{
	my ($d, $name1, $l1, $name2, $l2) = @_;

	my $k      = $h2m / $d / $d;

	my $oname1 = $lOrbName[$l1];
	my $oname2 = $lOrbName[$l2];
	my $etaS = (defined $eta{"$oname1${oname2}S"})? $eta{"$oname1${oname2}S"} : $eta{"$oname2${oname1}S"};
	my $etaP = (defined $eta{"$oname1${oname2}P"})? $eta{"$oname1${oname2}P"} : $eta{"$oname2${oname1}P"};
	my $etaD = (defined $eta{"$oname1${oname2}D"})? $eta{"$oname1${oname2}D"} : $eta{"$oname2${oname1}D"};
	
if(!defined $etaS) {
	print "Error: Can not find eta{$oname1${oname2}S}\n";
	exit;
}

	return ($k*$etaS, $k*$etaP, $k*$etaD);
}