#!/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 $PrintLevel = 2; # from -2 to 2

my $OutFile    = "TB1DBand.csv";
my @lOrbName   = ('s', 'p', 'd', 'f');
my @lmOrbName  = (
	['s'],
	['py', 'pz', 'px'],							# for l=1, m=-1, 0, 1
	['dxy', 'dyz', 'd3z2-r2', 'dzx', 'dx2-y2'], # for l=2, m=-2,-1,0,1,2
	['fy(3x2-y2)', 'fxyz', 'fy(5z2-r2)', 'fz(5z2-3r2)',
	 'fx(5z2-r2)', 'fz(x2-y2)', 'fx(x2-3y2)'],	# for l=3, m=-3,-2,-1,0,1,2,3
	);

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,
	'spP' =>  0.0,
	'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 $CIFFile = 'Si-Primitive.cif';
#my $CIFFile = 'Si.cif';

#No P1 RhombHex HexRhomb HexOrtho FCCPrim BCCPrim A/B/CCenterPrim
#my $BravaisLattice = 'No';
my $BravaisLattice = 'FCCPrim';

#==============================================================
# 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.5,  0.0,  1.0, 'W'],
	[0.5,  0.5,  0.5, 'L'],
	[0.0,  0.0,  0.0, 'GAMMA'],
	[0.0,  0.0,  1.0, 'X'],
	[0.5,  0.0,  1.0, 'W'],
	[0.75, 0.0, 0.75, 'K'],
	);

my $InputMode      = 'Cartesian';
#my $InputMode     = 'Reciprocal';
my $Direction;
if($InputMode eq 'Cartesian') {
	$Direction = 'OriginalToConverted';
}
else {
	$Direction = 'ConvertedlToOriginal';
}
# 1st: No P1 RhombHex HexRhomb HexOrtho FCCPrim BCCPrim A/B/CCenterPrim
# 2nd: OriginalToConverted ConvertedlToOriginal
# 3rd: ai xyz a*i hkl
my ($sT, $T, $tRT, $RT, $tR) = Crystal->GetLatticeConversionMatrix($BravaisLattice, $Direction, 'ai', 1);
#1st argで指定された変換行列  : $sT
#2nd argで選択された変換行列  : $T
#実空間座標(x,y,z)の変換行列  : $tRT = transpose(T)
#逆空間ベクトル(a*i)の変換行列: $RT  = transpose(T)
#逆空間座標(h k l)の変換行列  : $tR  = $T
#print "Conv M specified by 1st arg: sT=\n", $sT;
#print "Conv M selected by 2nd arg :  T=\n", $T;
#print "Conv M for (x,y,z)         : RT=\n", $RT;
#print "Conv M for (a*_i)          :tRT=\n", $tRT;
#print "Conv M for (h k l)         : tT=\n", $tR;
my $hkl = new Math::MatrixReal(3, 1);

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

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

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

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

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

#==============================================================
# Main routine
#==============================================================
#==============================================================
# Read CIF
#==============================================================
print "\n";
print "Crystal structure read from [$CIFFile]\n";

my $cif = new CIF;
$cif->Read($CIFFile) or die "$!: Can not read [$CIFFile]\n";
my $Crystal = $cif->GetCCrystal();
$Crystal->ExpandCoordinates();

my $CrystalName = $Crystal->CrystalName();
print "  CrystalName: $CrystalName\n";
my ($aa,$ba,$c,$alpha,$beta,$gamma) = $Crystal->LatticeParameters();
printf "  cell: %9.6f %9.6f %9.6f   %9.6f %9.6f %9.6f \n", $aa, $ba, $c, $alpha, $beta, $gamma;
my ($a11, $a12, $a13, $a21, $a22, $a23, $a31, $a32, $a33) = $Crystal->LatticeVectors();
printf "  Vectors: (%9.6f %9.6f %9.6f)\n", $a11, $a12, $a13;
printf "           (%9.6f %9.6f %9.6f)\n", $a21, $a22, $a23;
printf "           (%9.6f %9.6f %9.6f)\n", $a31, $a32, $a33;
my ($g11, $g12, $g13, $g21, $g22, $g23, $g31, $g32, $g33) = $Crystal->Metrics();
my $vol = $Crystal->Volume();
my ($Ra,$Rb,$Rc,$Ralpha,$Rbeta,$Rgamma) = $Crystal->ReciprocalLatticeParameters();
my ($Ra11, $Ra12, $Ra13, $Ra21, $Ra22, $Ra23, $Ra31, $Ra32, $Ra33) = $Crystal->ReciprocalLatticeVectors();
my ($Rg11, $Rg12, $Rg13, $Rg21, $Rg22, $Rg23, $Rg31, $Rg32, $Rg33) = $Crystal->RMetrics();
my $Rvol = $Crystal->CalculateReciprocalVolume();
my $SPG           = $Crystal->GetCSpaceGroup();
my $SPGName       = $SPG->SPGName();
my $iSPG          = $SPG->iSPG();
my $LatticeSystem = $SPG->LatticeSystem();
my @AtomTypeList  = $Crystal->GetCAtomTypeList();
my $nAtomType     = @AtomTypeList;
my @AtomSiteList = $Crystal->GetCAsymmetricAtomSiteList();
my $nAsymmetricAtomSite = @AtomSiteList;
my @ExpandedAtomSiteList = $Crystal->GetCExpandedAtomSiteList();
my $nTotalExpandedAtomSite = $Crystal->nTotalExpandedAtomSite();

print "  nAtomType: $nAtomType\n";
for(my $i = 0 ; $i < $nAtomType ; $i++) {
	my $label        = $AtomTypeList[$i]->Label();
	my $atomname     = $AtomTypeList[$i]->AtomName();
	my $charge       = $AtomTypeList[$i]->Charge();
	my $AtomicNumber = $AtomTypeList[$i]->AtomicNumber();
	my $AtomicMass   = $AtomTypeList[$i]->AtomicMass();
	my $i1 = $i+1;
	print "    #$i1: $label : $atomname [$AtomicNumber] ($charge) ($AtomicMass)\n";
}

print "  nTotalExpandedAtomSite: $nTotalExpandedAtomSite\n";
for(my $i = 0 ; $i < $nTotalExpandedAtomSite ; $i++) {
	my $label     = $ExpandedAtomSiteList[$i]->Label();
	my $type      = $ExpandedAtomSiteList[$i]->AtomName();
	my ($x,$y,$z) = $ExpandedAtomSiteList[$i]->Position(1);
	my $occupancy = $ExpandedAtomSiteList[$i]->Occupancy();
	my $iAtomType = $Crystal->FindiAtomType($type);

	print "    $label ($type)[#$iAtomType]: ($x, $y, $z) [$occupancy]\n";
}
print "\n";

#==============================================================
# Analyze Hij dimension
#==============================================================
my @Orbitals;
for(my $isite = 0 ; $isite < $nTotalExpandedAtomSite ; $isite++) {
	my $pAtom = $ExpandedAtomSiteList[$isite];
	my $atom  = $pAtom->AtomNameOnly();
	my $pAtom = $AtomOribitals{$atom};
	for(my $iorb = 0 ; $iorb < @$pAtom ; $iorb++) {
		my $n = $pAtom->[$iorb][0];
		my $l = $pAtom->[$iorb][1];
		my $m = $pAtom->[$iorb][2];
		push(@Orbitals, [$isite, $atom, $n.$lmOrbName[$l][$l+$m], $n, $l, $m]);
	}
}
my $nOrb = @Orbitals;

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

#==============================================================
# Build k list
#==============================================================
my @k;
my ($kxp, $kyp, $kzp);
my $ktot = 0.0;
for(my $i = 0 ; $i < @KList-1 ; $i++) {
	my ($kx0, $ky0, $kz0, $name0) = @{$KList[$i]};
	my ($kx1, $ky1, $kz1, $name1) = @{$KList[$i+1]};
	my $dkx = ($kx1 - $kx0) / ($nk - 1);
	my $dky = ($ky1 - $ky0) / ($nk - 1);
	my $dkz = ($kz1 - $kz0) / ($nk - 1);
	if($i == 0) {
		($kxp, $kyp, $kzp) = ($kx0, $ky0, $kz0);
	}
	for(my $j = 0 ; $j < $nk ; $j++) {
		my $kx = $kx0 + $j * $dkx;
		my $ky = $ky0 + $j * $dky;
		my $kz = $kz0 + $j * $dkz;
		my $k2 = $Crystal->GetReciprocalDistanceFromK($kx, $ky, $kz, $kxp, $kyp, $kzp);
		$ktot += $k2;
		push(@k, [$kx, $ky, $kz, $ktot, "$name0-$name1"]);
		($kxp, $kyp, $kzp) = ($kx, $ky, $kz);
	}
}

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

for(my $i = 0 ; $i < @k ; $i++) {
	my $p = $k[$i];
	my ($kx, $ky, $kz, $ktot, $kname) = @$p;
print "\n";
print "k = ($kx, $ky, $kz): $kname\n";

	$hkl = Math::MatrixReal->new_from_cols([ [$kx, $ky, $kz] ]);
	my $chkl = $T * $hkl;
	($kx, $ky, $kz) = ($chkl->element(1,1), $chkl->element(2,1), $chkl->element(3,1));

	my ($pHijR, $pHijI) = &CalHij($kx, $ky, $kz);
if($PrintLevel >= 0) {
print "\n";
print "HijR($kx,$ky,$kz)=\n";
&PrintMatrix($pHijR);
print "HijI($kx,$ky,$kz)=\n";
&PrintMatrix($pHijR);
}

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 ($pEn, $pCnij) = &Diagonalize($pHijR, $pHijI);
&PrintEigenValues($pEn, $pCnij);

	my @Esorted = sort { $a <=> $b } @$pEn;
	my $p  = $k[$i];
	$out->print("$kname,$kx,$ky,$kz,$ktot");
	for(my $io = 0 ; $io < $nOrb ; $io++) {
		$out->print(",$Esorted[$io]");
	}
	$out->print("\n");
}

$out->Close();

exit;

#================================================================
# Subroutines
#================================================================
sub Diagonalize
{
	my ($pHijR, $pHijI) = @_;
	
	my $EPS = 1.0e-6;
	my $MatrixIn = 'matrix.in';
	&SaveMatrix($MatrixIn, $pHijR, $pHijI, $EPS);

	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();

	my (@En, @Cnij);
	for(my $i = 0  ; $i < $nOrb ; $i++) {
		$En[$i] = $pEn->element($i+1, 1);
		for(my $j = 0 ; $j < $nOrb ; $j++) {
			$Cnij[$i][$j] = $pHij->element($i+1, $j+1);
		}
	}
	
	return (\@En, \@Cnij);
}

sub SaveMatrix
{
	my ($MatrixIn, $pHijR, $pHijI, $EPS) = @_;
	
	my $n = @$pHijR;

	my $mtx = JFile->new($MatrixIn, 'w') or die "$!: Can not write to [$MatrixIn].\n";
	$mtx->print("$n $EPS\n");
	$mtx->print("i j HijR HijI\n");
	for(my $i = 0 ; $i < $n ; $i++) {
		for(my $j = 0 ; $j < $n ; $j++) {
			$mtx->print($i+1, " ", $j+1, "   $pHijR->[$i][$j]   $pHijI->[$i][$j]\n");
		}
	}
	$mtx->print("i j SijR SijI\n");
	for(my $i = 0 ; $i < $n ; $i++) {
		for(my $j = 0 ; $j < $n ; $j++) {
			my $S = ($i == $j)? 1 : 0;
			$mtx->print($i+1, " ", $j+1, "   $S   0.0\n");
		}
	}
	
	$mtx->Close();
}

sub PrintMatrix
{
	my ($pM) = @_;
	
	my $n = @$pM;
	for(my $i = 0 ; $i < $n ; $i++) {
		print "[";
		for(my $j = 0 ; $j < $n ; $j++) {
			my $v = $pM->[$i][$j];
			$v = 0.0 if(abs($v) < 1.0e-6);
			printf " %12.4g", $v;
		}
		print "]\n";
	}
}

sub PrintEigenValues
{
	my ($pEn, $pCnij) = @_;

	my $n = @$pEn;

	print "\n";
	print "E (eV): Cij";
	for(my $i = 0 ; $i < $n ; $i++) {
		printf " %12s", "$Orbitals[$i][0]:$Orbitals[$i][1]$Orbitals[$i][2]";
	}
	print "\n";

	for(my $i = 0 ; $i < $n ; $i++) {
		printf "  %8.4f:", $pEn->[$i];
		for(my $j = 0 ; $j < $n ; $j++) {
			my $v = $pCnij->[$i][$j];
			$v = 0.0 if(abs($v) < 1.0e-6);
			printf " %12.4g", $v;
		}
		print "\n";
	}
}

sub PrintMatrixReal
{
	my ($pM, $n) = @_;
	
	for(my $i = 0 ; $i < $n ; $i++) {
		print "[";
		for(my $j = 0 ; $j < $n ; $j++) {
			my $v = $pM->element($i+1, $j+1);
			$v = 0.0 if(abs($v) < 1.0e-6);
			printf " %12.4g", $v;
		}
		print "]\n";
	}
}

sub PrintVectorReal
{
	my ($pM, $n) = @_;
	
	print "[";
	for(my $i = 0 ; $i < $n ; $i++) {
		printf " %12.4g", $pM->element($i+1, 1);
	}
	print "]\n";
}

sub CalHij
{
	my ($kx, $ky, $kz) = @_;

	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 $idx1 = 0  ; $idx1 < $nOrb ; $idx1++) {
# [$isite, $atom, $n.$lmOrbName[$l][$l+$m], $n, $l, $m]);
		my $pOrbital1 = $Orbitals[$idx1];
		my $is1       = $pOrbital1->[0];
		my $name1     = $pOrbital1->[1];
		my $OrbName1  = $pOrbital1->[2];
		my $n1        = $pOrbital1->[3];
		my $l1        = $pOrbital1->[4];
		my $m1        = $pOrbital1->[5];
		my $pSite1    = $ExpandedAtomSiteList[$is1];
		my ($x1, $y1, $z1) = $pSite1->Position(1);
		my $pAtom1    = $AtomOribitals{$name1};
		for(my $idx2 = $idx1  ; $idx2 < $nOrb ; $idx2++) {
#print "Sites: $is1 - $is2\n";
			my $pOrbital2 = $Orbitals[$idx2];
			my $is2       = $pOrbital2->[0];
			my $name2     = $pOrbital2->[1];
			my $OrbName2  = $pOrbital2->[2];
			my $n2        = $pOrbital2->[3];
			my $l2        = $pOrbital2->[4];
			my $m2        = $pOrbital2->[5];
			my $pSite2    = $ExpandedAtomSiteList[$is2];
			my ($x2, $y2, $z2) = $pSite2->Position(1);
			my $pAtom2    = $AtomOribitals{$name2};
			for(my $iz = -$nrange ; $iz <= $nrange ; $iz++) {
			for(my $iy = -$nrange ; $iy <= $nrange ; $iy++) {
			for(my $ix = -$nrange ; $ix <= $nrange ; $ix++) {
				my $x2r = $x2 + $ix;
				my $y2r = $y2 + $iy;
				my $z2r = $z2 + $iz;
				my $r = $Crystal->GetInterAtomicDistance($x1, $y1, $z1, $x2r, $y2r, $z2r);
				next if(abs($r) > $Rmax);

				if($r == 0.0) {
					if($idx1 == $idx2) {
if(!defined $aii{$name1}{$OrbName1}) {
print "Error: aii is not defined for $name1 $OrbName1\n";
exit;
}

						$HijR[$idx1][$idx2] += $aii{$name1}{$OrbName1};
print "Hij[$idx1][$idx2] = aii{$name1}{$OrbName1} = $HijR[$idx1][$idx2]\n" if($PrintLevel >= 1);
					}
				}
				else {
					my ($VS, $VP, $VD) = &Getbij($r, $name1, $l1, $name2, $l2);
					my ($rx, $ry, $rz) = $Crystal->CalCartesianCoordinate($x1-$x2r, $y1-$y2r, $z1-$z2r);
					my $lx  = $rx / $r;
					my $ly  = $ry / $r;
					my $lz  = $rz / $r;
					my $lx2 = $lx * $lx;
					my $ly2 = $ly * $ly;
					my $lz2 = $lz * $lz;
# Slater-Koster conversion
#['py', 'pz', 'px']: m = (-1,0,1)
					if($l1 == 0 and $l2 == 0) { #ss
						$b = $VS;
					}
					elsif(($l1 == 0 and $l2 == 1 and $m2 == -1) or ($l1 == 1 and $m1 == -1 and $l2 == 0)) { #sy
						$b = $ly * $VS;
					}
					elsif(($l1 == 0 and $l2 == 1 and $m2 ==  0) or ($l1 == 1 and $m1 ==  0 and $l2 == 0)) { #sz
						$b = $lz * $VS;
					}
					elsif(($l1 == 0 and $l2 == 1 and $m2 ==  1) or ($l1 == 1 and $m1 ==  1 and $l2 == 0)) { #sx
						$b = $lx * $VS;
					}
					elsif(($l1 == 1 and $m1 == -1 and $l2 == 1 and $m2 == -1)) { #yy
						$b = $ly2 * $VS + (1.0 - $ly2) * $VP;
					}
					elsif(($l1 == 1 and $m1 ==  0 and $l2 == 1 and $m2 ==  0)) { #zz
						$b = $lz2 * $VS + (1.0 - $lz2) * $VP;
					}
					elsif(($l1 == 1 and $m1 ==  1 and $l2 == 1 and $m2 ==  1)) { #xx
						$b = $lx2 * $VS + (1.0 - $lx2) * $VP;
					}
					elsif(($l1 == 1 and $m1 == -1 and $l2 == 1 and $m2 ==  0) or 
						  ($l1 == 1 and $m1 ==  0 and $l2 == 1 and $m2 == -1)) { #yz
						$b = $ly*$lz * $VS - $ly*$lz * $VP;
					}
					elsif(($l1 == 1 and $m1 == -1 and $l2 == 1 and $m2 ==  1) or 
						  ($l1 == 1 and $m1 ==  1 and $l2 == 1 and $m2 == -1)) { #yx
						$b = $ly*$lx * $VS - $ly*$lx * $VP;
					}
					elsif(($l1 == 1 and $m1 ==  0 and $l2 == 1 and $m2 ==  1) or 
						  ($l1 == 1 and $m1 ==  1 and $l2 == 1 and $m2 ==  0)) { #zx
						$b = $lz*$lx * $VS - $lz*$lx * $VP;
					}

					my $phi  = $kx * ($x1 - $x2r) + $ky * ($y1 - $y2r) + $kz * ($z1 - $z2r);
					my $ekrr = cos($pi2 * $phi);
					my $ekri = sin($pi2 * $phi);
					$HijR[$idx1][$idx2] += $ekrr * $b;
					$HijI[$idx1][$idx2] += $ekri * $b;
printf "Hij[$idx1][$idx2]: r = %8.4f ($is1:$name1, $OrbName1)-($is2:$name2, $OrbName2, ixyz=$ix,$iy,$iz): VS=%8.4f VP=%8.4f\n",
		$r, $VS, $VP if($PrintLevel >= 1);
printf "  HijR[%2d][%2d]=%8.4f = exp(2pi * i *%8.4f)) * %8.4f\n",
		$idx1, $idx2, $HijR[$idx1][$idx2], $phi, $b if($PrintLevel >= 1);
printf "    [ekr=(%8.4f, %8.4f)  b=%8.4f (lxyz=%8.4f,%8.4f,%8.4f)]\n",
		$ekrr, $ekri,$b, $lx, $ly, $lz if($PrintLevel >= 2);
				}
			}
			}
			}
		}
	}

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

	return (\@HijR, \@HijI);
}

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"};
#print "eta=($etaS, $etaP, $etaD) for $oname1$oname2\n";
if(!defined $etaS) {
	print "Error: Can not find eta{$oname1${oname2}S}\n";
	exit;
}

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