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

my $OutFile = "TB1DBand.csv";

my $Rmax   = 3.0;
my $nrange = 1;

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

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

# e0
my %Atom = (
#	'Si' => [ '2s', ],
#	'C'  => [ '2s', ],
	'Si' => [ '2s', '2p', '3s' ],
	'C'  => [ '2s', '2p' ],
	);

my @HijIndex;

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

my %bij = (
	'Si 2s-Si 2s' => 0.2,
	'Si 2s-Si 2p' => 0.2,
	'Si 2p-Si 2s' => 0.2,
	'Si 2s-Si 3s' => 0.2,
	'Si 3s-Si 2s' => 0.2,
	'Si 2p-Si 2p' => 0.2,
	'Si 2p-Si 3s' => 0.2,
	'Si 3s-Si 3s' => 0.2,

	'C 2s-C 2s'   => 0.2,
	'C 2s-C 2p'   => 0.2,
	'C 2p-C 2s'   => 0.2,
	'C 2p-C 2p'   => 0.2,

	'Si 2s-C 2s'  => 0.2,
	'C 2s-Si 2s'  => 0.2,
	'Si 2s-C 2p'  => 0.2,
	'C 2p-Si 2s'  => 0.2,
	'Si 2p-C 2s'  => 0.2,
	'C 2s-Si 2p'  => 0.2,
	'Si 2p-C 2p'  => 0.2,
	'C 2p-Si 2p'  => 0.2,
	'Si 3s-C 2s'  => 0.2,
	'C 2s-Si 3s'  => 0.2,
	'Si 3s-C 2p'  => 0.2,
	'C 2p-Si 3s'  => 0.2,
	);

#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, ''],
	);

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++) {
		push(@Orbitals, [$isite, $atom, $pAtom->[$iorb]]);
		$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;

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

	my @Hij;
	for(my $i = 0  ; $i < $nOrb ; $i++) {
		for(my $j = 0 ; $j < $nOrb ; $j++) {
			$Hij[$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 $OrbName1 = $pAtom1->[$io1];
					my $idx1     = $HijIndex[$is1][$io1];
					for(my $io2 = 0 ; $io2 < @$pAtom2 ; $io2++) {
						my $OrbName2 = $pAtom2->[$io2];
						my $idx2     = $HijIndex[$is2][$io2];
						next if($idx2 < $idx1);

						if($r == 0.0) {
							if($idx1 == $idx2) {
								$Hij[$idx1][$idx2] += $aii{"$name1 $OrbName1"};
print "Hij[$idx1][$idx2] = aii{$name1 $OrbName1}\n";
							}
						}
						else {
#							my $ekr = exp(2.0*$pi * $k * ($x1 - $x2r));
#							my $ekr = exp(2.0*$pi * i * $k * ($x1 - $x2r));
							my $ekr = cos(2.0*$pi * $k * ($x1 - $x2r));
							my $b   = $bij{"$name1 $OrbName1-$name2 $OrbName2"};
							if($OrbName1 =~ /s/ and $OrbName2 =~ /s/) {
#								$b = $b;
							}
							elsif($OrbName1 =~ /s/ and $OrbName2 =~ /p/ or
								  $OrbName1 =~ /p/ and $OrbName2 =~ /s/) {
								$b = ($x1 - $x2r) / $r * $b;
							}
							elsif($OrbName1 =~ /p/ and $OrbName2 =~ /p/) {
#								$b = $b;
							}
							
#print "ekr=$ekr  b=$b\n";
							$Hij[$idx1][$idx2] += $ekr * $b;
print "Hij[$idx1][$idx2] = exp(i * $k * ($x1 - $x2r)) * bij{$name1 $OrbName1-$name2 $OrbName2}\n";
#print "Hij[$idx1][$idx2] = $Hij[$idx1][$idx2]\n";
						}
print "Hij[$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, $Hij[$j][$i]);
		}
		for(my $j = $i ; $j < $nOrb ; $j++) {
			$pHij->assign($i+1, $j+1, $Hij[$i][$j]);
		}
	}

	return $pHij;
}