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

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

use strict;
use File::Basename;

use Crystal::SpaceGroup;
use Crystal::DVXa;
use Crystal::XCrySDen;

#my $DirectorySeparator = "\\";

my $WebRootDir      = "d:\\MyWebs";
my $ProgramDir      = "d:\\Programs";
my $DVXaDir         = Deps::MakePath($ProgramDir, "DVXa97",   1);
my $NonRelDir       = Deps::MakePath($DVXaDir, "exec",   1);
my $NonRelPath      = Deps::MakePath($NonRelDir, "nonrel", 0);
my $KListDBDir      = Deps::MakePath($WebRootDir, "Research", 1);
$KListDBDir         = Deps::MakePath($KListDBDir, "klist", 1);

#===============================================
# 変数取得関数
#===============================================
sub KListDBDir {
	my $this = shift;
	return $this->{'KListDBDir'} if($this->{'KListDBDir'});
	return $this->{'KListDBDir'} = $KListDBDir;
}
sub SetKListDBDir { 
	my($this, $d) = @_; 
	$this->{'KListDBDir'} = $d;
	return $d;
}
sub NonRelDir {
	my $this = shift;
	return $this->{'NonRelDir'} if($this->{'NonRelDir'});
	return $this->{'NonRelDir'} = $NonRelDir;
}
sub SetNonRelDir { 
	my($this, $d) = @_; 
	$this->{'NonRelDir'} = $d;
	$this->SetNonRelFile();
	return $d;
}
sub SetNonRelFile {
	my ($this) = @_;
	my $file = "nonrel";
	$this->{'NonRelFile'} = $file;
	$NonRelDir = $this->NonRelDir();
	$this->{'NonRelPath'} 
		= $NonRelPath = Deps::MakePath($NonRelDir, $file, 0);
#print "dd: ", $this->NonRelDir(), "<br>\n";
#print "ff: ", $this->{'NonRelPath'}, "<br>\n";
	return $this->{'NonRelPath'};
}
sub NonRelPath
{
	my ($this) = @_;
	return $this->{'NonRelPath'};
}

#===============================================
# コンストラクタ、デストラクタ
#===============================================
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 SaveTBInputFile
{
	my ($this, $Crystal, $filename) = @_;
	my $LF = "<br>\n";

	my $DVXa = new DVXa;
	my $NonRelDir  = $DVXa->SetNonRelDir($this->NonRelDir());
	my $NonRelPath = $DVXa->NonRelPath();
#print "Nonrel path: $NonRelPath$LF";
	my $XC = new XCrySDen;
	my $KListDBDir = $XC->SetKListDBDir($this->KListDBDir);
#print "klist DB dir: $KListDBDir$LF";

#ファイル名を（ベース名, ディレクトリ名, 拡張子）に分解
	my @filenames = fileparse($filename, "\.[^\.]+");

	my @AtomTypeList = $Crystal->GetCAtomTypeList();

#ファイル作製開始
	unless(open(OUT,">$filename")) {
		print "Can not write to [$filename].$LF$LF";
		return;
	}

	my $SampleName = $this->SampleName();
	my ($a,$b,$c,$alpha,$beta,$gamma) 
		= $Crystal->LatticeParametersByOutputMode(0);
	my @AtomTypeList = $Crystal->GetCAtomTypeList();
	my $nAtomTypes = @AtomTypeList;
	my @ExpandedAtomSiteList 
		= $Crystal->GetCExpandedAtomSiteListByOutputMode();
	my $nAtoms = @ExpandedAtomSiteList;
#	my $SPGName = $Crystal->SPGNameByOutputMode();
#	my $iSPG    = $Crystal->iSPGByOutputMode();
	my $SPGName = $Crystal->SPGName();
	my $iSPG    = $Crystal->iSPG();
	my ($KListFilePath, @KList) = $XC->ReadKList($iSPG, $SPGName);
#print "KList DB file: $KListFilePath$LF\n";

	print OUT<<EOT1;
[Title]
$SampleName
[Atomic Basis functions]   n l m  n(electrons)  energy(eV)  ri  rd
EOT1
#print "nAtomTypes: $nAtomTypes<br>\n";
	print OUT "$nAtomTypes\n";
	for(my $i = 0 ; $i < $nAtomTypes ; $i++) {
		my $type   = $AtomTypeList[$i];
		my $name   = $type->AtomNameOnly();
		my $charge = $type->Charge();
#print "$i: $name ($charge)<br>\n";
		$charge = 0 unless(defined $charge);

#print "'$name'<br>\n";
		print OUT "'$name'  '$charge'\n";
		my ($nOrbitals, $AtomicNumber, @Orbitals) = $DVXa->ReadOrbitals($name);
		DVXa::ModifyIonOccupancy(\@Orbitals, $charge);
		print OUT "  $nOrbitals\n";
		for(my $j = 0 ; $j < $nOrbitals ; $j++) {
			my $pOrb = $Orbitals[$j];
			my $n = $pOrb->{'n'};
			my $l = $pOrb->{'l'};
			my $m = 0;
			my $Energy = $pOrb->{'Energy'};
			my $Ne = $pOrb->{'Ne'};
			my $ri = 0.0;
			my $rd = 0.0;
			if($l == 2.0) {
				$rd = 1.10;
			}
#n l m  n(electrons)  energy(eV)  ri  rd
			printf OUT "  %d %d %d  %8.4f  %12.4f %8.4f %8.4f\n",
				$n, $l, $m, $Ne, $Energy, $ri, $rd;
		}
#'Sr'
#  10
#  1 0 0   2. 1 0 0
#  2 0 0   2. 1 0 0
#  2 1 0   6. 1 0 0
#  3 0 0   2. 1 0 0
#  3 1 0   6. 1 0 0
#  3 2 0  10. 1 0 0
#  4 0 0   2. 1 0 0
#  4 1 0   6. 1 0 0
#  4 2 0   0. 1 0 0
#  5 0 0   0. 5.00 1.10 0
	}
	
	print OUT<<EOT2;

[Coefficients of Hamiltonian matrix] korb1, korb2, sigma(1) or pi(2) or delta(3), ehta
10
  1 1 1 -1.40
  1 2 1  1.84
  2 2 1  3.24
  2 2 2 -0.81
  1 3 1 -3.16
  2 3 1 -2.95
  2 3 2  1.36
  3 3 1 -16.2
  3 3 2  8.75
  3 3 3  0
 
[Crystal Structure] in A unit
EOT2

	print OUT "cell\n";
	printf OUT "%10.6f %10.6f %10.6f %12.6f %12.6f %12.6f\n",
		$a, $b, $c, $alpha, $beta, $gamma;

	print OUT "$nAtoms\n";
#Occupancyが1のサイト
	for(my $i = 0 ; $i < $nAtoms ; $i++) {
		my $atom      = $ExpandedAtomSiteList[$i];
		my $atomname  = $atom->AtomNameOnly();
		my ($x,$y,$z) = $atom->Position(1);
		my $occupancy = $atom->Occupancy();
#		my $mult      = $atom->Multiplicity();
		my $id        = $atom->Id();
		next if($occupancy < 0.9999);

		printf OUT "%2d %12.8f  %12.8f  %12.8f\n",
			$id, $x, $y, $z;
	}

#Occupancyが1未満のサイト
	for(my $i = 0 ; $i < $nAtoms ; $i++) {
		my $atom      = $ExpandedAtomSiteList[$i];
		my $atomname  = $atom->AtomNameOnly();
		my ($x,$y,$z) = $atom->Position(1);
		my $occupancy = $atom->Occupancy();
#		my $mult      = $atom->Multiplicity();
		my $id        = $atom->Id();
		next if($occupancy > 0.9999);

		printf OUT "%2d %12.8f  %12.8f  %12.8f\n",
			$id, $x, $y, $z;
	}

	print OUT "\n";
	print OUT "[Calculation points in Reciprocal Space]\n";
	my $nKList = @KList - 1;
	print OUT "$nKList\n";
	for(my $i = 1 ; $i < $nKList+1 ; $i++) {
		my $pk1 = $KList[$i-1];
		my $pk2 = $KList[$i];
		my $label1 = $pk1->{'Label'};
		my $kx1    = $pk1->{'kx'};
		my $ky1    = $pk1->{'ky'};
		my $kz1    = $pk1->{'kz'};
		my $label2 = $pk2->{'Label'};
		my $kx2    = $pk2->{'kx'};
		my $ky2    = $pk2->{'ky'};
		my $kz2    = $pk2->{'kz'};
		my $nPoint = $pk2->{'nPoint'};
#0 0 0     0 0 1        10
#0 0 0     0 1 1        10
		printf OUT "%10.4f %10.4f %10.4f   %10.4f %10.4f %10.4f   %3d   '%s' '%s'\n",
			$kx1, $ky1, $kz1, $kx2, $ky2, $kz2, $nPoint, $label1, $label2;
	}

	print OUT<<EOT3;

[Control] EPS for Eigen values, overlap integrals (ratio), overlap integrals (distance / A)
1E-6 1e-3 5.0
EOT3

	close OUT;

	return ($KListFilePath, $filename);
}


1;