package TightBinding;
use Exporter;
@ISA = qw(Exporter);

#公開したいサブルーチン
@EXPORT = qw();

use strict;
use File::Basename;
#use warnings;
use Math::MatrixReal;
use Math::Complex;

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

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

#===============================================
# 変数取得関数
#===============================================

#===============================================
# コンストラクタ、デストラクタ
#===============================================
sub new
{
	my ($module) = @_;
	my $this = {};
	bless $this;
	return $this;
}

sub DESTROY
{
	my $this = shift;
}

#===============================================
# 一般関数
#===============================================
sub SetSampleName
{
	my ($this, $name) = @_;
	return $this->{'SampleName'} = $name;
}

sub SampleName
{
	my ($this) = @_;
	return $this->{'SampleName'};
}

sub ReadInput
{
	my ($this, $InFile) = @_;

	$this->{InFile} = $InFile;

	my %p;
	my $in = JFile->new($InFile, 'r') or die "$!: Can not read [$InFile].\n";
	while(!$in->eof()) {
		my $line = $in->ReadLine();
		next if($line =~ /^\s*#/);

		my ($key, $val) = ($line =~ /^\s*(\S+)\s*=\s*(\S.*)\s*$/);
		next if(!defined $key);
		$p{$key} = $val;
	}
	$in->Close();

	return $this->{pInput} = \%p;
}

sub PrintInputInf
{
	my ($this, $pInput) = @_;
	$pInput = $this->{pInput} if(!defined $pInput);

	print "# Control parameters\n";
	my @keylist = qw(CrystalName CIFFile UseConventionalCell KPOINTSFile BandFile DBFile
			 Rmax nRange 
			 UseZincBlendAnalyticalForm Diagonalization DiagonalizationEPS 
			 PrintLevel SortEigenValues);
	foreach my $key (@keylist) {
		printf "%-20s: $this->{pInput}{$key}\n", $key;
	}
	print "\n";
}

sub ReadDB
{
	my ($this, $DBFile) = @_;

	my (@lOrbName, @lmOrbName, %eta, %kF, %kd, %rc, %ri, %ro, %rd, %AtomOrbitals, %aii);
	$this->{DBFile} = $DBFile;

	my ($line, $key);
	my $in = JFile->new($DBFile, 'r') or die "$!: Can not read [$DBFile].\n";
	$in->rewind();
	$line = $in->SkipTo("lOrbName:");
	($key, @lOrbName) = Utils::Split("\\s+", $line);
#print "OrbName=", join(', ', @lOrbName, "\n");

	$in->rewind();
	$in->SkipTo("lmOrbName:");
	my $c = 0;
	while(1) {
		$line = $in->ReadLine();
		next if($line =~ /^\s*#/);

		my @b;
		($key, @b) = Utils::Split("\\s+", $line);
		last if(!defined $b[0]);

		$lmOrbName[$c++] = \@b;
	}
#print "OrbName=\n";
#for(my $i = 0 ; $i < @lmOrbName ; $i++) {
#	my $p = $lmOrbName[$i];
#	print "  ", join(', ', @$p), "\n";
#}

	$in->rewind();
	$in->SkipTo("eta:");
	while(1) {
		$line = $in->ReadLine();
		next if($line =~ /^\s*#/);

		Utils::DelSpace($line);
		my ($key, $val) = Utils::Split("\\s*=\\s*", $line);
		last if(!defined $key);

		$eta{$key} = $val;
	}
#foreach my $key (keys %eta) {
#print "eta: $key: $eta{$key}\n";
#}

	$in->rewind();
	$in->SkipTo("Atom specific");
	my $Atom = '';
	my $c    = 0;
	while(!$in->eof()) {
		$line = $in->ReadLine();
		next if($line =~ /^\s*#/);

		if($line =~ /^\s*([A-Z][a-z]?):/) {
			$Atom = $1;
#print "atom: $Atom\n";
			$c = 0;
		}

		if($line =~ /(\S.*?)\s*=\s*(\S.*?)\s*$/i) {
			my $key = $1;
			my $val = $2;
			if($key eq 'kF') {
				$kF{$Atom} = $val;
			}
			elsif($key eq 'kd') {
				$kd{$Atom} = $val;
			}
			elsif($key eq 'rc') {
				$rc{$Atom} = $val;
			}
			elsif($key eq 'ri') {
				$ri{$Atom} = $val;
			}
			elsif($key eq 'ro') {
				$ro{$Atom} = $val;
			}
			elsif($key eq 'rd') {
				$rd{$Atom} = $val;
			}
#print "  $key = $val\n";
		}
		elsif($line =~ /(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)/i) {
			my ($orb, $n, $l, $m, $E0) = ($1, $2, $3, $4, $5);
#print "  $orb: $n $l $m: $E0\n";
			next if(!defined $E0);
			$AtomOrbitals{$Atom}[$c] = [$n, $l, $m, $orb, $E0];
			$aii{$Atom}{$orb} = $E0;
			$c++;
		}
	}
	$in->Close();

if(0) {
foreach my $key (keys %kF) {
	print "  kF($key)=$kF{$key}\n";
}
foreach my $key (keys %rc) {
	print "  rc($key)=$kF{$key}\n";
}
foreach my $key (keys %ri) {
	print "  ri($key)=$ri{$key}\n";
}
foreach my $key (keys %aii) {
	my $p = $aii{$key};
	foreach my $key2 (keys %$p) {
		print "  aii($key)($key2)=$p->{$key2}\n";
	}
}
foreach my $key (keys %AtomOrbitals) {
	my $p = $AtomOrbitals{$key};
	for(my $i = 0 ; $i < @$p ; $i++) {
		print "  ($key): $p->[$i][0] $p->[$i][1] $p->[$i][2]\n";
	}
}
exit;
}

	$this->{pDB} = {
		plOrbName     => \@lOrbName,
		plmOrbName    => \@lmOrbName,
		peta          => \%eta,
		pkF           => \%kF,
		pkd           => \%kd,
		prc           => \%rc,
		pri           => \%ri,
		pro           => \%ro,
		prd           => \%rd,
		pAtomOrbitals => \%AtomOrbitals,
		paii          => \%aii,
		};
	
	return $this->{pDB};
}

sub PrintElectronicParameters
{
	my ($this, $pDB) = @_;
	$pDB = $this->{pDB} if(!defined $pDB);

	my @lOrbName     = @{$pDB->{plOrbName}};
	my @lmOrbName    = @{$pDB->{plmOrbName}};
	my %eta          = %{$pDB->{peta}};
	my %kF           = %{$pDB->{pkF}};
	my %kd           = %{$pDB->{pkd}};
	my %rc           = %{$pDB->{prc}};
	my %ri           = %{$pDB->{pri}};
	my %ro           = %{$pDB->{pro}};
	my %rd           = %{$pDB->{prd}};
	my %AtomOrbitals = %{$pDB->{pAtomOrbitals}};
	my %aii          = %{$pDB->{paii}};

	my $Crystal = $this->{pCrystal};
	my @AtomTypeList           = $Crystal->GetCAtomTypeList();
	my $nAtomType              = @AtomTypeList;
	my @ExpandedAtomSiteList   = $Crystal->GetCExpandedAtomSiteList();
	my $nTotalExpandedAtomSite = $Crystal->nTotalExpandedAtomSite();

	print "# Electronic parameters\n";
	for(my $i = 0 ; $i < $nAtomType ; $i++) {
		my $atomname = $AtomTypeList[$i]->AtomNameOnly();

if(!defined $kF{$atomname} and !defined $rc{$atomname} and !defined $ri{$atomname} and !defined $ro{$atomname}) {
	print "Error: Atomic data for [$atomname] is not given in [$this->{pInput}{DBFile}].\n";
	exit;
}
		print "  $i: $atomname\n";
		print "     rd=$rd{$atomname}\n" if(defined $rd{$atomname});
		print "     kd=$ri{$atomname}\n" if(defined $kd{$atomname});
		print "     kF=$kF{$atomname}\n" if(defined $kF{$atomname});
		print "     rc=$rc{$atomname}\n" if(defined $rc{$atomname});
		print "     ri=$ri{$atomname}\n" if(defined $ri{$atomname});
		print "     ro=$ro{$atomname}\n" if(defined $ro{$atomname});
		my $paii = $aii{$atomname};
		my $pao  = $AtomOrbitals{$atomname};

if(!defined $pao or !defined $pao) {
	print "Error: Atomic data for [$atomname] is not given in [$this->{pInput}{DBFile}].\n";
	exit;
}

		printf "     %2s %2s %2s [orbname]: E0 [eV]\n", 'n', 'l', 'm';
		for(my $i = 0 ; $i < @$pao ; $i++) {
			my $n        = $pao->[$i][0];
			my $l        = $pao->[$i][1];
			my $m        = $pao->[$i][2];
			my $orbname  = $pao->[$i][3];
			my $orbname2 = $n . $lmOrbName[$l][$l+$m];
			if($orbname ne $orbname2) {
				print "Error: Input orbit name [$atomname $orbname] is inconsistent "
					."with (n l m) = ($n $l $m) [$orbname2].\n";
				exit;
			}
			my $E0       = $paii->{$orbname};
			printf "     %2d %2d %2d [$orbname]: $E0\n", $n, $l, $m;
		}
	}
	print "\n";
}

sub ReadCIF
{
	my ($this, $CIFFile) = @_;
	
	$this->{CIFFile} = $CIFFile;

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

	my $SPG     = $Crystal->GetCSpaceGroup();
	my $SPGName = $SPG->SPGName();
	if($this->{pInput}{UseConventionalCell} or $SPGName !~ /^\s*[FIABC]/i) {
#		return $this->{pCrystal} = $Crystal;
	}
	else
	{
#Primitive Cellをつくる
		print "  Convert to Primitive Cell.\n";
		print "  SpaceGroup: $SPGName\n";

		my $T = new Math::MatrixReal(3, 3);
		if($SPGName =~ /^\s*F/i) {
			$T = Math::MatrixReal->new_from_cols( 
				[ [ 0.5, 0.5,   0], 
				  [ 0,   0.5, 0.5],
				  [ 0.5,   0, 0.5] ] );
			$T->transpose($T);
		}
		elsif($SPGName =~ /^\s*I/i) {
			$T = Math::MatrixReal->new_from_cols( 
				[ [-0.5, 0.5, 0.5], 
				  [ 0.5,-0.5, 0.5],
				  [ 0.5, 0.5,-0.5] ] );
			$T->transpose($T);
		}
		elsif($SPGName =~ /^\s*A/i) {
			$T = Math::MatrixReal->new_from_cols( 
				[ [ 1.0, 0.0, 0.0], 
				  [ 0.0, 0.5, 0.5],
				  [ 0.0,-0.5, 0.5] ] );
			$T->transpose($T);
		}
		elsif($SPGName =~ /^\s*B/i) {
			$T = Math::MatrixReal->new_from_cols( 
				[ [ 0.5, 0.0, 0.5], 
				  [ 0.0, 1.0, 0.0],
				  [-0.5, 0.0, 0.5] ] );
			$T->transpose($T);
		}
		elsif($SPGName =~ /^\s*C/i) {
			$T = Math::MatrixReal->new_from_cols( 
				[ [ 0.5, 0.5, 0.0], 
				  [-0.5, 0.5, 0.0],
				  [ 0.0, 0.0, 1.0] ] );
			$T->transpose($T);
		}

#実空間べクトル(ai)の変換行列
		print  "  Conversion matrix for real space vector: (a'i) = (Tij)(aj)\n";
		printf "  |%10.4f %10.4f %10.4f|\n", $T->element(1,1), $T->element(1,2), $T->element(1,3);
		printf "  |%10.4f %10.4f %10.4f|\n", $T->element(2,1), $T->element(2,2), $T->element(2,3);
		printf "  |%10.4f %10.4f %10.4f|\n", $T->element(3,1), $T->element(3,2), $T->element(3,3);
		print "\n";

		my ($a, $b, $c, $alpha, $beta, $gamma) = $Crystal->LatticeParameters();
		$Crystal = $Crystal->ConvertLattice($T, 0, 1);
	}

#my $CrystalName = $Crystal->CrystalName();
#my ($aa,$ba,$c,$alpha,$beta,$gamma) = $Crystal->LatticeParameters();
#my ($a11, $a12, $a13, $a21, $a22, $a23, $a31, $a32, $a33) = $Crystal->LatticeVectors();
#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();

	return $this->{pCrystal} = $Crystal;
}

sub PrintCrystalInf
{
	my ($this, $Crystal) = @_;
	$Crystal = $this->{pCrystal} if(!defined $Crystal);

	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 @AtomTypeList           = $Crystal->GetCAtomTypeList();
	my $nAtomType              = @AtomTypeList;
	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]->AtomNameOnly();
		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";
}

sub ReadKPOINTS
{
	my ($this, $InFile) = @_;

	my @KList;
	$this->{KPOINTSFile} = $InFile;

	my $in = JFile->new($InFile, 'r') or die "$!: Can not read [$InFile].\n";
	my $line = $in->ReadLine();
	my $n = $in->ReadLine() + 0;
	$line = $in->ReadLine();
	$line = $in->ReadLine();
	Utils::DelSpace($line);
	if($line !~ /R/i) {
		print "Error: KPOINTS should be given based on the reciprocal mode (mode=[$line]).\n";
		exit;
	}

	my ($kxp, $kyp, $kzp);
	my $c = 0;
	while(1) {
		last if($in->eof());
		$line = $in->ReadLine();
		my ($kx, $ky, $kz, $kname) = Utils::Split("\\s+", $line);
		next if(!defined $kz);

#print "l($kxp,$kyp,$kzp)/($kx,$ky,$kz):$line";
		if(!defined $kxp or ($kx != $kxp or $ky != $kyp or $kz != $kzp)) {
#print "c=$c\n";
			$KList[$c] = [$kx, $ky, $kz, $kname];
			$c++;
		}
		$kxp = $kx;
		$kyp = $ky;
		$kzp = $kz;
	}
#for(my $i = 0 ; $i < @KList ; $i++) {
#	my $p = $KList[$i];
#	print "$i: $p->[0], $p->[1], $p->[2], $p->[3]\n";
#}
	$in->Close();

	$this->{pKPOINTS} = {
		nK     => $n,
		pKList => \@KList,
		};

	return ($n, \@KList);
}

sub BuildOrbitals
{
	my ($this) = @_;

	my $Crystal = $this->{pCrystal};
	my @AtomTypeList           = $Crystal->GetCAtomTypeList();
	my $nAtomType              = @AtomTypeList;
	my @ExpandedAtomSiteList   = $Crystal->GetCExpandedAtomSiteList();
	my $nTotalExpandedAtomSite = $Crystal->nTotalExpandedAtomSite();

	my $pDB = $this->{pDB};
	my @lOrbName     = @{$pDB->{plOrbName}};
	my @lmOrbName    = @{$pDB->{plmOrbName}};
	my %AtomOrbitals = %{$pDB->{pAtomOrbitals}};

	my @Orbitals;
#print "n=$nTotalExpandedAtomSite\n";
	for(my $isite = 0 ; $isite < $nTotalExpandedAtomSite ; $isite++) {
		my $pAtom = $ExpandedAtomSiteList[$isite];
		my $atom  = $pAtom->AtomNameOnly();
		my $pAtom = $AtomOrbitals{$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 $io = @Orbitals - 1;
#print "$Orbitals[$io][0] $Orbitals[$io][1] $Orbitals[$io][2] $Orbitals[$io][3] \n";
		}
	}
#exit;

	return $this->{pOrbitals} = \@Orbitals;
}

sub BuildKList
{
	my ($this) = @_;

	my $Crystal  = $this->{pCrystal};
	my $pKPOINTS = $this->{pKPOINTS};
	my $n        = $pKPOINTS->{nK};
	my $pKPoints = $pKPOINTS->{pKList};

	my $BS = new BandStructure;
	my (@kl, @KList);

	for(my $i = 0 ; $i < @$pKPoints ; $i++) {
		$kl[$i] = {
			kx   => $pKPoints->[$i][0],
			ky   => $pKPoints->[$i][1],
			kz   => $pKPoints->[$i][2],
			name => $pKPoints->[$i][3],
		};
#print "k[$i]{name}=$kl[$i]{name}\n";
	}

	if($n > 2) {
		@$pKPoints = $BS->ReduceKList(\@kl);
#for(my $i = 0 ; $i < @$pKPoints ; $i++) {
#print "k[$i]{name}=$pKPoints->[$i]{name}\n";
#}
#exit;
		($n, $pKPoints) = $BS->DivideKList($Crystal, $pKPoints, $n);
	}

	for(my $i = 0 ; $i < @$pKPoints ; $i++) {
		$KList[$i] = [$pKPoints->[$i]{kx}, $pKPoints->[$i]{ky}, $pKPoints->[$i]{kz}, $pKPoints->[$i]{name}];
	}
#for(my $i = 0 ; $i < @$pKPoints ; $i++) {
#print "k[$i]{name}=$KList[$i][3]\n";
#}
#exit;

	my @k;
	my $ktot = 0.0;
	if($n <= 2) {
		for(my $i = 0 ; $i < @KList ; $i++) {
			my ($kx0, $ky0, $kz0, $name0) = @{$KList[$i]};
#printf "k[$i]: (%8.4f, %8.4f, %8.4f, %s)\n", $kx0, $ky0, $kz0, $name0;
			push(@k, [$kx0, $ky0, $kz0, 0.0, "$name0"]);
		}
	}
	else {
		for(my $i = 0 ; $i < @KList-1 ; $i++) {
			my ($kx0, $ky0, $kz0, $name0) = @{$KList[$i]};
			my ($kx1, $ky1, $kz1, $name1) = @{$KList[$i+1]};
#printf "k[$i]: (%8.4f, %8.4f, %8.4f, %s)-(%8.4f, %8.4f, %8.4f, %s)\n",
#		$kx0, $ky0, $kz0, $name0, $kx1, $ky1, $kz1, $name1;
			my $k2 = $Crystal->GetReciprocalDistanceFromK($kx0, $ky0, $kz0, $kx1, $ky1, $kz1);
			push(@k, [$kx0, $ky0, $kz0, $ktot, $name0]);
			$ktot += $k2;
			if($i == @KList-2) {
				push(@k, [$kx1, $ky1, $kz1, $ktot, $name1]);
			}
		}
	}
	for(my $i = 0 ; $i < @k ; $i++) {
printf "k[$i]: (%8.4f, %8.4f, %8.4f) k=%8.4f %s\n", $k[$i][0], $k[$i][1], $k[$i][2], $k[$i][3], $k[$i][4];
	}

	return $this->{pCalKList} = \@k;
}

sub PrintMatrix
{
	my ($this, $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 BuildFockMatrixForZincBlend
{
	my ($this, $kx, $ky, $kz) = @_;

	my $pInput  = $this->{pInput};

	my $Crystal = $this->{pCrystal};
	my @AtomTypeList           = $Crystal->GetCAtomTypeList();
	my $nAtomType              = @AtomTypeList;
	my @ExpandedAtomSiteList   = $Crystal->GetCExpandedAtomSiteList();
	my $nTotalExpandedAtomSite = $Crystal->nTotalExpandedAtomSite();

	my $pDB     = $this->{pDB};
	my @lOrbName     = @{$pDB->{plOrbName}};
	my @lmOrbName    = @{$pDB->{plmOrbName}};
	my %AtomOrbitals = %{$pDB->{pAtomOrbitals}};
	my %eta          = %{$pDB->{peta}};
	my %kF           = %{$pDB->{pkF}};
	my %kd           = %{$pDB->{pkd}};
	my %rc           = %{$pDB->{prc}};
	my %ri           = %{$pDB->{pri}};
	my %ro           = %{$pDB->{pro}};
	my %rd           = %{$pDB->{prd}};
	my %AtomOrbitals = %{$pDB->{pAtomOrbitals}};
	my %aii          = %{$pDB->{paii}};

	my $pOrbitals = $this->{pOrbitals};
	my @Orbitals  = @$pOrbitals;
	my $nOrb      = @Orbitals;

	my (@g, $ph111, $ph1mm, $phm1m, $phmm1, @hij);

#==============================================
# From analyatical equations
#==============================================
#Reciprocal W: 0.250000   0.500000   0.750000  W
#Caltesian  W: 0.5   0.0   1.0   W
	my ($sT, $T, $tRT, $RT, $tR) = Crystal->GetLatticeConversionMatrix('FCCPrim', 'ConvertedlToOriginal', '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 "T=\n";
	print $T;
	my $hkl = new Math::MatrixReal(3, 1);
	$hkl = Math::MatrixReal->new_from_cols([ [$kx, $ky, $kz] ]);
	my $chkl = $T * $hkl;
	my ($kxc, $kyc, $kzc) = ($chkl->element(1,1), $chkl->element(2,1), $chkl->element(3,1));
	print "k(Cartesian)=($kxc, $kyc, $kzc)\n";

	$ph111 =  $kxc * 0.25 + $kyc * 0.25 + $kzc * 0.25;
	$ph1mm =  $kxc * 0.25 - $kyc * 0.25 - $kzc * 0.25;
	$phm1m = -$kxc * 0.25 + $kyc * 0.25 - $kzc * 0.25;
	$phmm1 = -$kxc * 0.25 - $kyc * 0.25 + $kzc * 0.25;
	$g[0] = exp(i * $pi2 * $ph111) + exp(i * $pi2 * $ph1mm) + exp(i * $pi2 * $phm1m) + exp(i * $pi2 * $phmm1);
	$g[1] = exp(i * $pi2 * $ph111) + exp(i * $pi2 * $ph1mm) - exp(i * $pi2 * $phm1m) - exp(i * $pi2 * $phmm1);
	$g[2] = exp(i * $pi2 * $ph111) - exp(i * $pi2 * $ph1mm) + exp(i * $pi2 * $phm1m) - exp(i * $pi2 * $phmm1);
	$g[3] = exp(i * $pi2 * $ph111) - exp(i * $pi2 * $ph1mm) - exp(i * $pi2 * $phm1m) + exp(i * $pi2 * $phmm1);
for(my $i = 0 ; $i < @g ; $i++) {
print "g$i=$g[$i]\n";
}

	my $atomname0      = $ExpandedAtomSiteList[0]->AtomNameOnly();
	my ($x0, $y0, $z0) = $ExpandedAtomSiteList[0]->Position();
	my $atomname1      = $ExpandedAtomSiteList[1]->AtomNameOnly();
	my ($x1, $y1, $z1) = $ExpandedAtomSiteList[1]->Position();
	my $r = $Crystal->GetInterAtomicDistance($x0, $y0, $z0, $x1, $y1, $z1);

	my ($VssS, $VssP, $VssD) = $this->Getbij($r, $atomname0, 0, $atomname1, 0);
	my ($VspS, $VspP, $VspD) = $this->Getbij($r, $atomname0, 0, $atomname1, 1);
	my ($VppS, $VppP, $VppD) = $this->Getbij($r, $atomname0, 1, $atomname1, 1);

	my $Ess =  $VssS;
	my $Esp = -$VspS / sqrt(3.0);
	my $Exx =  $VppS / 3.0 + 2.0/3.0 * $VppP;
	my $Exy =  $VppS / 3.0 - $VppP / 3.0;
print "Ess=$VssS\n";
print "Esp=$Esp\n";
print "Exx=$Exx  Exy=$Exy\n";

	for(my $i = 0 ; $i < $nOrb ; $i++) {
		my $pOrbital = $Orbitals[$i];
		my $atomname = $pOrbital->[1];
		my $orbname  = $pOrbital->[2];
		my $E0   = $aii{$atomname}{$orbname};
print "$i: $E0, $atomname, $orbname\n";
		$hij[$i][$i] = $E0;
	}
	$hij[0][1] = 0.0;
	$hij[0][2] = 0.0;
	$hij[0][3] = 0.0;
	$hij[0][4] =  $Ess * $g[0];
#	$hij[0][5] =  $Esp * $g[2];
	$hij[0][5] = -$Esp * $g[2];
#	$hij[0][6] =  $Esp * $g[3];
	$hij[0][6] = -$Esp * $g[3];
#	$hij[0][7] =  $Esp * $g[1];
	$hij[0][7] = -$Esp * $g[1];
	$hij[1][2] =  0.0;
	$hij[1][3] =  0.0;
#	$hij[1][4] = -$Esp * $g[2];
	$hij[1][4] =  $Esp * $g[2];
	$hij[1][5] =  $Exx * $g[0];
	$hij[1][6] =  $Exy * $g[1];
	$hij[1][7] =  $Exy * $g[3];
	$hij[2][3] =  0.0;
#	$hij[2][4] = -$Esp * $g[3];
	$hij[2][4] =  $Esp * $g[3];
#	$hij[2][5] =  $Exy * $g[2]; #$g[2]
	$hij[2][5] =  $Exy * $g[1]; #$g[2]
	$hij[2][6] =  $Exx * $g[0];
#	$hij[2][7] =  $Exy * $g[1]; #$g[1]
	$hij[2][7] =  $Exy * $g[2]; #$g[1]
#	$hij[3][4] = -$Esp * $g[1];
	$hij[3][4] =  $Esp * $g[1];
	$hij[3][5] =  $Exy * $g[3];
#	$hij[3][6] =  $Exy * $g[1]; #$g[1];
	$hij[3][6] =  $Exy * $g[2]; #$g[1];
	$hij[3][7] =  $Exx * $g[0];
#  0   1   2   3  4   5   6   7
# s1 py1 pz1 px1 s2 py2 pz2 px2
	$hij[4][5] =  0.0;
	$hij[4][6] =  0.0;
	$hij[4][7] =  0.0;
	$hij[5][6] =  0.0;
	$hij[5][7] =  0.0;
	$hij[6][7] =  0.0;

	for(my $i = 0 ; $i < $nOrb ; $i++) {
		for(my $j = $i+1 ; $j < $nOrb ; $j++) {
#			$hij[$j][$i] = ~$hij[$i][$j];
			$hij[$j][$i] = Re($hij[$i][$j]) - i * Im($hij[$i][$j]);
#print "[$i,$j]: $hij[$j][$i] <=> $hij[$i][$j]\n";
		}
	}

	return @hij;
}

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

	my @s = qw(s1 py1 pz1 px1 s2 py2 pz2 px2);

	my $pInput  = $this->{pInput};
	my $Rmax    = $pInput->{Rmax};
	my $nrange  = $pInput->{nRange};

	my $Crystal = $this->{pCrystal};
	my @AtomTypeList           = $Crystal->GetCAtomTypeList();
	my $nAtomType              = @AtomTypeList;
	my @ExpandedAtomSiteList   = $Crystal->GetCExpandedAtomSiteList();
	my $nTotalExpandedAtomSite = $Crystal->nTotalExpandedAtomSite();

	my $pDB     = $this->{pDB};
	my @lOrbName     = @{$pDB->{plOrbName}};
	my @lmOrbName    = @{$pDB->{plmOrbName}};
	my %AtomOrbitals = %{$pDB->{pAtomOrbitals}};
	my %eta          = %{$pDB->{peta}};
	my %kF           = %{$pDB->{pkF}};
	my %kd           = %{$pDB->{pkd}};
	my %rc           = %{$pDB->{prc}};
	my %ri           = %{$pDB->{pri}};
	my %ro           = %{$pDB->{pro}};
	my %rd           = %{$pDB->{prd}};
	my %AtomOrbitals = %{$pDB->{pAtomOrbitals}};
	my %aii          = %{$pDB->{paii}};

	my $pOrbitals = $this->{pOrbitals};
	my @Orbitals  = @$pOrbitals;
	my $nOrb      = @Orbitals;

my (@g, @E, @hij);
if($pInput->{UseZincBlendAnalyticalForm}) {
	@hij = $this->BuildFockMatrixForZincBlend($kx, $ky, $kz);

	for(my $i = 0 ; $i < $nOrb ; $i++) {
		for(my $j = 0 ; $j < $nOrb ; $j++) {
#print "$i,$j $hij[$i][$j]\n";
			if(!defined $hij[$i][$j]) {
				print "hij[$i][$j] is not defined\n";
				exit;
			}
		}
	}
}
else {
	for(my $i = 0 ; $i < $nOrb ; $i++) {
	for(my $j = 0 ; $j < $nOrb ; $j++) {
		$hij[$i][$j] = 0;
	}
	}
}

#==============================================
# By summation
#==============================================
	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    = $AtomOrbitals{$name1};
my $IsPrintedHij = 0;
		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    = $AtomOrbitals{$name2};

my $IsPrinted = 0;
			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 $dx  = $x2r - $x1;
				my $dy  = $y2r - $y1;
				my $dz  = $z2r - $z1;
				my $r = $Crystal->GetInterAtomicDistance($x1, $y1, $z1, $x2r, $y2r, $z2r);
				next if(abs($r) > $Rmax);

				my ($VS, $VP, $VD);
				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($pInput->{PrintLevel} >= 1);
					}
					else {
						$HijR[$idx1][$idx2] = 0.0;
						$HijI[$idx1][$idx2] = 0.0;
print "Hij[$idx1][$idx2] = $HijR[$idx1][$idx2]\n" if($pInput->{PrintLevel} >= 1);
					}
				}
				else {
					($VS, $VP, $VD) = $this->Getbij($r, $name1, $l1, $name2, $l2);
if($pInput->{PrintLevel} >= 1 and !$IsPrinted) {
print "VS=$VS  VP=$VP  VD=$VD ($OrbName1-$OrbName2) r=$r\n";
$IsPrinted=1;
}
					$b = $this->CalVSPD($l1, $m1, $l2, $m2, $VS, $VP, $VD, $r, $dx, $dy, $dz);
#if($idx1 == 1 and $idx2 == 6) {
#print "$l1,$m1,$l2,$m2: ($ix,$iy,$iz): r=$r, b=$b (l=$lx,$ly,$lz) VS=$VS\n";
#exit;
#}
					my $phi  = $kx * $dx + $ky * $dy + $kz * $dz;
					my $ekrr = cos($pi2 * $phi);
					my $ekri = sin($pi2 * $phi);
					$HijR[$idx1][$idx2] += $ekrr * $b;
					$HijI[$idx1][$idx2] += $ekri * $b;

if($pInput->{UseZincBlendAnalyticalForm}) {
	my ($rx, $ry, $rz) = $Crystal->CalCartesianCoordinate($dx, $dy, $dz);
	my $lx  = $rx / $r;
	my $ly  = $ry / $r;
	my $lz  = $rz / $r;
	if($lx > 0.0) {
		$ly = $ly/$lx;
		$lz = $lz/$lx;
		$lx = 1;
		}
	else {
		$ly = $ly/abs($lx);
		$lz = $lz/abs($lx);
		$lx = -1;
	}
	printf "Hij[$idx1][$idx2]=(%8.4f,%8.4f)  b=%8.4f  lxyz=%2.0f%2.0f%2.0f  ekr=(%8.4f, %8.4f)\n", 
			$HijR[$idx1][$idx2], $HijI[$idx1][$idx2], $b, $lx, $ly, $lz, $ekrr, $ekri;
			my $err = abs($hij[$idx1][$idx2] - ($HijR[$idx1][$idx2] + i*$HijI[$idx1][$idx2]));
	if($err > 0.01) {
		print "  AF:H[$idx1][$idx2]($s[$idx1]-$s[$idx2])=$hij[$idx1][$idx2] ($HijR[$idx1][$idx2], $HijI[$idx1][$idx2]) err=$err\n";
#		print "     phi: 111=$ph111, 1mm=$ph1mm, m1m=$phm1m, mm1=$phmm1\n";
		print "     Error too large\n" if($err > 0.01);
		$IsPrintedHij = 1;
	}
}
else {
	my ($rx, $ry, $rz) = $Crystal->CalCartesianCoordinate($dx, $dy, $dz);
	my $lx  = $rx / $r;
	my $ly  = $ry / $r;
	my $lz  = $rz / $r;
	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($pInput->{PrintLevel} >= 1);
#printf "  Hij[%2d][%2d]=(%8.4f,%8.4f) = exp(2pi * i *%8.4f)) * %8.4f\n",
#		$idx1, $idx2, $HijR[$idx1][$idx2], $HijI[$idx1][$idx2], $phi, $b if($pInput->{PrintLevel} >= 1);
	printf "    [ekr=(%8.4f, %8.4f)  b=%8.4f (lxyz=%3.1f,%3.1f,%3.1f)]\n",
				$ekrr, $ekri,$b, $lx, $ly, $lz              if($pInput->{PrintLevel} >= 2);
#print "     x: ($x1,$y1,$z1)-($x2,$y2,$z2)/($x2r,$y2r,$z2r) = dx=($dx,$dy,$dz)\n";
	$IsPrintedHij = 1;
}

			}
			}
			}
}
#print "\n" if($IsPrintedHij);
		}
	}

	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];
		}
	}
#exit;

if($pInput->{UseZincBlendAnalyticalForm}) {
print "HijR($kx,$ky,$kz)=\n";
$this->PrintMatrix(\@HijR);
print "\n";
print "HijI($kx,$ky,$kz)=\n";
$this->PrintMatrix(\@HijI);
print "\n";

for(my $i = 0 ; $i < $nOrb ; $i++) {
for(my $j = 0 ; $j < $nOrb ; $j++) {
	$HijR[$i][$j] = Re($hij[$i][$j]);
	$HijI[$i][$j] = Im($hij[$i][$j]);
#print "$i,$j: ($HijR[$i][$j], $HijI[$i][$j])\n";
}
}
}

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

sub CalVSPD
{
	my ($this, $l1, $m1, $l2, $m2, $VS, $VP, $VD, $r, $dx, $dy, $dz) = @_;

	my $Crystal = $this->{pCrystal};

	my ($rx, $ry, $rz) = $Crystal->CalCartesianCoordinate($dx, $dy, $dz);
	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) { #sy
		$b =  $ly * $VS;
	}
	elsif($l1 == 1 and $m1 == -1 and $l2 == 0) { #ys
		$b = -$ly * $VS;
	}
	elsif($l1 == 0 and $l2 == 1 and $m2 == 0) { #sz
		$b =  $lz * $VS;
	}
	elsif($l1 == 1 and $m1 == 0 and $l2 == 0) { #zs
		$b = -$lz * $VS;
	}
	elsif($l1 == 0 and $l2 == 1 and $m2 == 1) { #sx
		$b =  $lx * $VS;
	}
	elsif($l1 == 1 and $m1 == 1 and $l2 == 0) { #xs
		$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 - $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 - $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 - $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 - $VP);
	}
	elsif(($l1 == 0 and $m1 ==  0 and $l2 == 2 and $m2 == -2) or
	      ($l1 == 2 and $m1 == -2 and $l2 == 0 and $m2 == 0)) { #sxy
		$b = sqrt(3.0)*$lx*$ly * $VS;
	}
	elsif(($l1 == 0 and $m1 ==  0 and $l2 == 2 and $m2 == -1) or
	      ($l1 == 2 and $m1 == -1 and $l2 == 0 and $m2 == 0)) { #syz
		$b = sqrt(3.0)*$ly*$lz * $VS;
	}
	elsif(($l1 == 0 and $m1 ==  0 and $l2 == 2 and $m2 ==  0) or
	      ($l1 == 2 and $m1 ==  0 and $l2 == 0 and $m2 == 0)) { #s3z2-r2
	$b = ($lz2 - 0.5*($lx2 + $ly2)) * $VS;
	}
	elsif(($l1 == 0 and $m1 ==  0 and $l2 == 2 and $m2 ==  1) or
	      ($l1 == 2 and $m1 ==  1 and $l2 == 0 and $m2 == 0)) { #szx
		$b = sqrt(3.0)*$lz*$lx * $VS;
	}
	elsif(($l1 == 0 and $m1 ==  0 and $l2 == 2 and $m2 ==  2) or
	      ($l1 == 2 and $m1 ==  2 and $l2 == 0 and $m2 == 0)) { #sx2-y2
		$b = 0.5 * sqrt(3.0) * ($lx2 - $ly2) * $VS;
	}
	elsif(($l1 == 1 and $m1 == -1 and $l2 == 2 and $m2 == -2)) { #y xy
		$b =  sqrt(3.0) * $ly2*$lx * $VS + $lx*($1.0-2.0*$ly2) * $VP;
	}
	elsif(($l1 == 2 and $m1 == -2 and $l2 == 1 and $m2 == -1)) { #xy y
		$b =  sqrt(3.0) * $ly2*$lx * $VS + $lx*($1.0-2.0*$ly2) * $VP;
		$b = -$b;
	}
	elsif(($l1 == 1 and $m1 == -1 and $l2 == 2 and $m2 == -1)) { #y yz
		$b =  sqrt(3.0) * $ly2*$lz * $VS + $lz*(1.0-2.0*$ly2) * $VP;
	}
	elsif(($l1 == 2 and $m1 == -1 and $l2 == 1 and $m2 == -1)) { #yz y
		$b =  sqrt(3.0) * $ly2*$lz * $VS + $lz*(1.0-2.0*$ly2) * $VP;
		$b = -$b;
	}
	elsif(($l1 == 1 and $m1 == -1 and $l2 == 2 and $m2 ==  0)) { #y 3z2-r2
		$b =  $ly*($lz2-0.5*($lx2+$ly2)) * $VS - sqrt(3.0) * $ly*$lz2 * $VP;
	}
	elsif(($l1 == 2 and $m1 ==  0 and $l2 == 1 and $m2 == -1)) { #3z2-r2 y
		$b =  $ly*($lz2-0.5*($lx2+$ly2)) * $VS - sqrt(3.0) * $ly*$lz2 * $VP;
		$b = -$b;
	}
	elsif(($l1 == 1 and $m1 == -1 and $l2 == 2 and $m2 ==  1)) { #y zx
		$b =  $lx*$ly*$lz * (sqrt(3.0) * $VS - 2.0 * $VP);
	}
	elsif(($l1 == 2 and $m1 ==  1 and $l2 == 1 and $m2 == -1)) { #zx y
		$b =  $lx*$ly*$lz * (sqrt(3.0) * $VS - 2.0 * $VP);
		$b = -$b;
	}
	elsif(($l1 == 1 and $m1 == -1 and $l2 == 2 and $m2 ==  2)) { #y x2-y2
		$b =  0.5*sqrt(3.0) * $ly*($lx2-$ly2) * $VS - $ly*(1.0+$lx2-$ly2) * $VP;
	}
	elsif(($l1 == 2 and $m1 ==  2 and $l2 == 1 and $m2 == -1)) { #x2-y2 y
		$b =  0.5*sqrt(3.0) * $ly*($lx2-$ly2) * $VS - $ly*(1.0+$lx2-$ly2) * $VP;
		$b = -$b;
	}
	elsif(($l1 == 1 and $m1 ==  0 and $l2 == 2 and $m2 == -2)) { #z xy
		$b =  $lx*$ly*$lz * (sqrt(3.0) * $VS - 2.0 * $VP);
	}
	elsif(($l1 == 2 and $m1 == -2 and $l2 == 1 and $m2 ==  0)) { #xy z
		$b =  $lx*$ly*$lz * (sqrt(3.0) * $VS - 2.0 * $VP);
		$b = -$b;
	}
	elsif(($l1 == 1 and $m1 ==  0 and $l2 == 2 and $m2 == -1)) { #z yz
		$b =  sqrt(3.0) * $lz2*$ly * $VS + $ly*(1.0-2.0*$lz2) * $VP;
	}
	elsif(($l1 == 2 and $m1 == -1 and $l2 == 1 and $m2 ==  0)) { #yz z
		$b = sqrt(3.0) * $lz2*$ly * $VS + $ly*(1.0-2.0*$lz2) * $VP;
		$b = -$b;
	}
	elsif(($l1 == 1 and $m1 ==  0 and $l2 == 2 and $m2 ==  0)) { #z 3z2-r2
		$b =  $lz*($lz2-0.5*($lx2+$ly2)) * $VS + sqrt(3.0)*$lz*($lx2+$ly2) * $VP;
	}
	elsif(($l1 == 2 and $m1 ==  0 and $l2 == 1 and $m2 ==  0)) { #3z2-r2 z
		$b =  $lz*($lz2-0.5*($lx2+$ly2)) * $VS + sqrt(3.0)*$lz*($lx2+$ly2) * $VP;
		$b = -$b;
	}
	elsif(($l1 == 1 and $m1 ==  0 and $l2 == 2 and $m2 ==  1)) { #z zx
		$b =  sqrt(3.0)*$lz2*$lx * $VS + $lx*(1.0-2.0*$lz) * $VP;
	}
	elsif(($l1 == 2 and $m1 ==  1 and $l2 == 1 and $m2 ==  0)) { #zx z
		$b =  sqrt(3.0)*$lz2*$lx * $VS + $lx*(1.0-2.0*$lz) * $VP;
		$b = -$b;
	}
	elsif(($l1 == 1 and $m1 ==  0 and $l2 == 2 and $m2 ==  2)) { #z x2-y2
		$b =  0.5*sqrt(3.0)*$lz*($lx2-$ly2) * $VS - $lz*($lx2-$ly2) * $VP;
	}
	elsif(($l1 == 2 and $m1 ==  2 and $l2 == 1 and $m2 ==  0)) { #x2-y2 z
		$b =  0.5*sqrt(3.0)*$lz*($lx2-$ly2) * $VS - $lz*($lx2-$ly2) * $VP;
		$b = -$b;
	}
	elsif(($l1 == 1 and $m1 ==  1 and $l2 == 2 and $m2 == -2)) { #x xy
		$b =  sqrt(3.0) * $lx2*$ly * $VS + $ly*(1.0-2.0*$lx2) * $VP;
	}
	elsif(($l1 == 2 and $m1 == -2 and $l2 == 1 and $m2 ==  1)) { #xy x
		$b =  sqrt(3.0) * $lx2*$ly * $VS + $ly*(1.0-2.0*$lx2) * $VP;
		$b = -$b;
	}
	elsif(($l1 == 1 and $m1 ==  1 and $l2 == 2 and $m2 == -1)) { #x yz
		$b =  $lx*$ly*$lz * (sqrt(3.0) * $VS - 2.0*$VP);
	}
	elsif(($l1 == 2 and $m1 == -1 and $l2 == 1 and $m2 ==  1)) { #yz x
		$b =  $lx*$ly*$lz * (sqrt(3.0) * $VS - 2.0*$VP);
		$b = -$b;
	}
	elsif(($l1 == 1 and $m1 ==  1 and $l2 == 2 and $m2 ==  0)) { #x 3z2-r2
		$b =  $lx*($lz2-0.5*($lx2+$ly2)) * $VS - sqrt(3.0) * $lx*$lz2 * $VP;
	}
	elsif(($l1 == 2 and $m1 ==  0 and $l2 == 1 and $m2 ==  1)) { #3z2-r2 x
		$b =  $lx*($lz2-0.5*($lx2+$ly2)) * $VS - sqrt(3.0) * $lx*$lz2 * $VP;
		$b = -$b;
	}
	elsif(($l1 == 1 and $m1 ==  1 and $l2 == 2 and $m2 ==  1)) { #x zx
		$b =  sqrt(3.0)*$lx2*$lz * $VS + $lz*(1.0-2.0*$lx2) * $VP;
	}
	elsif(($l1 == 2 and $m1 ==  1 and $l2 == 1 and $m2 ==  1)) { #zx x
		$b =  sqrt(3.0)*$lx2*$lz * $VS + $lz*(1.0-2.0*$lx2) * $VP;
		$b = -$b;
	}
	elsif(($l1 == 1 and $m1 ==  1 and $l2 == 2 and $m2 ==  2)) { #x x2-y2
		$b =  0.5*sqrt(3.0)*$lx*($lx2-$ly2) * $VS + $lx*(1.0-$lx2+$ly2) * $VP;
	}
	elsif(($l1 == 2 and $m1 ==  2 and $l2 == 1 and $m2 ==  1)) { #x x2-y2
		$b =  0.5*sqrt(3.0)*$lx*($lx2-$ly2) * $VS + $lx*(1.0-$lx2+$ly2) * $VP;
		$b = -$b;
	}
	
	return $b;
}

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

	my $pDB = $this->{pDB};
	my @lOrbName     = @{$pDB->{plOrbName}};
	my @lmOrbName    = @{$pDB->{plmOrbName}};
	my %AtomOrbitals = %{$pDB->{pAtomOrbitals}};
	my %eta          = %{$pDB->{peta}};
	my %kF           = %{$pDB->{pkF}};
	my %kd           = %{$pDB->{pkd}};
	my %rc           = %{$pDB->{prc}};
	my %ri           = %{$pDB->{pri}};
	my %ro           = %{$pDB->{pro}};
	my %rd           = %{$pDB->{prd}};
	my %AtomOrbitals = %{$pDB->{pAtomOrbitals}};
	my %aii          = %{$pDB->{paii}};

	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;
}
if(($l1 >= 1 or $l2 >= 1) and !defined $etaP) {
	print "Error: Can not find eta{$oname1${oname2}P}\n";
	exit;
}
if(($l1 >= 2 and $l2 >= 2) and !defined $etaD) {
	print "Error: Can not find eta{$oname1${oname2}D}\n";
	exit;
}
if($l1 >= 3 or $l2 >= 3) {
	print "Error: Orbital [l=$l1 or l=$l2] is not supported\n";
	exit;
}

# For sd, pd
	if($l1 < 2 and $l2 == 2) {
		my $c = $h2m * $rd{$name2}**1.5 / $d**3.5;
		return ($c*$etaS, $c*$etaP, $c*$etaD);
	}
	elsif($l1 == 2 and $l2 < 2) {
		my $c = $h2m * $rd{$name1}**1.5 / $d**3.5;
		return ($c*$etaS, $c*$etaP, $c*$etaD);
	}
# For dd
	elsif($l1 == 2 and $l2 == 2) {
		my $c = $h2m * (($rd{$name1}+$rd{$name2})/2.0)**3.0 / $d**5.0;
		return ($c*$etaS, $c*$etaP, $c*$etaD);
	}

# For ss, sp
	my $c = $h2m / $d / $d;
	return ($c*$etaS, $c*$etaP, $c*$etaD);
}

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

	my $pInput    = $this->{pInput};
	my $pOrbitals = $this->{pOrbitals};
	my @Orbitals  = @$pOrbitals;
	my $nOrb      = @Orbitals;

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

		$this->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 SaveMatrix
{
	my ($this, $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 PrintEigenValues
{
	my ($this, $pEV) = @_;

	my $pInput    = $this->{pInput};
	my $pOrbitals = $this->{pOrbitals};
	my @Orbitals  = @$pOrbitals;

	my $n = @$pEV;

	print "\n";
	print "E (eV): Cij";
	for(my $i = 0 ; $i < $n ; $i++) {
#push(@Orbitals, [$isite, $atom, $n.$lmOrbName[$l][$l+$m], $n, $l, $m]);
		printf "  %-18s", "$Orbitals[$i][0]:$Orbitals[$i][1]$Orbitals[$i][2]";
	}
	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,%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  = sqrt($vR*$vR + $vI*$vI);
			next if(abs($v) < $EPS);

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

			my $iSite = $Orbitals[$j][0];
			my $atom  = $Orbitals[$j][1];
			my $orb   = $Orbitals[$j][2];

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

1;