#!/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');

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

# n, l, m
my %Atom = (
	'Si' => [ [2, 0, 0], [2, 1, 0], [3, 0, 0] ],
	'C'  => [ [2, 0, 0], [2, 1, 0] ],
	);

# e0
my %aii = (
	'Si' => { '2s' =>  0.2, '2p' => 0.25, '3s' => 0.3 },
	'C'  => { '2s' =>  0.1, '2p' => 0.15 },
	);

# beta_ij
my %bij = (
	'Si 2s' => { 'Si 2s' => 0.2, 'Si 2p' => 0.2, 'Si 3s' => 0.2, 'C 2s'  => 0.2, 'C 2p'  => 0.2, },
	'Si 2p' => { 'Si 2p' => 0.2, 'Si 3s' => 0.2, 'C 2s'  => 0.2, 'C 2p'  => 0.2, },
	'Si 3s' => { 'Si 3s' => 0.2, 'C 2s'  => 0.2, 'C 2p'  => 0.2, },
	'C 2s'  => { 'C 2s'  => 0.2, 'C 2p'  => 0.2, },
	'C 2p'  => { 'C 2p'  => 0.2, },
	);

#==============================================================
# 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 = $Atom{$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 = $Atom{$name1};
		for(my $is2 = 0  ; $is2 < @Site ; $is2++) {
			my $site2  = $Site[$is2];
			my $name2  = $site2->[0];
			my $x2     = $site2->[1];
			my $pAtom2 = $Atom{$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 = exp(2.0*$pi * $k * ($x1 - $x2r));
#							my $ekrr = exp(2.0*$pi * i * $k * ($x1 - $x2r));
							my $ekrr = cos(2.0*$pi * $k * ($x1 - $x2r));
							my $ekri = sin(2.0*$pi * $k * ($x1 - $x2r));
							my $b   = &Getbij($name1, $OrbName1, $name2, $OrbName2);
							if($l1 == 0 and $l2 == 0) {
#								$b = $b;
							}
							elsif(($l1 == 0 and $l2 == 1) or ($l1 == 1 and $l2 == 0)) {
								$b = ($x1 - $x2r) / $r * $b;
							}
							elsif(($l1 == 1 and $l2 == 1)) {
#								$b = $b;
							}
							
#print "ekr=$ekr  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}\n";
#print "Hij[$idx1][$idx2] = $Hij[$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 ($name1, $OrbName1, $name2, $OrbName2) = @_;

	return $bij{"$name1 $OrbName1"}{"$name2 $OrbName2"} if(defined $bij{"$name1 $OrbName1"}{"$name2 $OrbName2"});
	return $bij{"$name2 $OrbName2"}{"$name1 $OrbName1"} if(defined $bij{"$name2 $OrbName2"}{"$name1 $OrbName1"});
	
print "Error: Can not find bij{$name2 $OrbName2}{$name1 $OrbName1}\n";
exit;
}