#===============================================
# TranSIESTA
#===============================================
package TranSIESTA;
use Exporter;
@ISA = qw(Exporter);

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

use strict;
use File::Basename;

use Sci::Science;
use Sci::EnergyBand;
use GraphData;
use MyTk::GraphFrameArray;

use Crystal::Crystal;
use Crystal::AtomType;
use Crystal::AtomSite;
use Crystal::CIF;
use Crystal::XCrySDen;

#===============================================
# パス（読み込みDB）
# Web関係変数
# CGI の仮想パス
# プログラム名
#===============================================
my $Program         = basename($0);
my $ProgramDir      = "d:\\Programs";
my $WebRootDir      = "d:\\MyWebs";
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 new
{
	my ($module) = @_;
	my $this = {};
	bless $this;
	return $this;
}

sub DESTROY
{
	my $this = shift;
}

#============================================================
#　継承クラスで定義しなおす関数
#============================================================
sub CheckFileType
{
	my ($path) = @_;
	my $in = new JFile($path, "r");
	return undef unless($in);
	my $line1 = $in->ReadLine();
	my $line2 = $in->ReadLine();
	$in->Close();
#	if($line1 =~ /^# ATK initiated at/) {
	if($line1 =~ /^# ATK / or $line2 =~ /^# ATK /) {
		return "ATK 2.0 Output";
	}
	return undef;
}

sub ReadForceConvergenceFromATK2Output
{
	my ($this, $path, $target) = @_;
	my ($drive, $dir, $filename, $ext, $lastdir, $filebody)
		= Deps::SplitFilePath($path);
	my $SampleName = $filebody;

print "  TranSIESTA::ReadForceConvergenceFromATK2Output: Reading [$path].\n";
	my $in = new JFile($path, "r");
	unless($in) {
		print "  Error: Can not ead [$path].\n";
		return undef;
	}

	$this->SetDataArray(new GraphDataArray);
	my $Data = new GraphData;
	$Data->SetTitle($SampleName);

	if(defined $target and $target =~ /ForceConvergence/i) {
	}
	else { #'EnergyConvergence'
		my @Iteration;
		my @Force;
		my $iter = 0;
		my $HitCount = 0;
		$in->rewind();
		while(!$in->eof()) {
			my $line = $in->SkipTo("^# Largest force component ");
			last if($in->eof());
			$HitCount++;
			next if($HitCount % 2);
			my ($f) = ($line =~ /=\s+([\d\.+\-Ee]+)/i);
			$Force[$iter] = $f;
			$Iteration[$iter] = $iter+1;
print "i=$iter: F(largest)=$Force[$iter]\n";
			$iter++;
		}
		$Data->{'x0'}      = \@Iteration;
		$Data->{'x0_Name'} = "Iteration";
		$Data->{'y0'}      = \@Force;
		$Data->{'y0_Name'} = "Force / eV/Ang";
	}

	$in->Close();

	$Data->CalMinMax();
	$this->{'DataArray'}->AddGraphData($Data);
	$this->{'DataArray'}->{'TargetData'} = $target;

	$in->Close();
print "  TranSIESTA::ReadForceConvergenceFromATK2Output: finished.\n";
	return $path;
}

sub ReadEnergyConvergenceFromATK2Output
{
	my ($this, $path, $target) = @_;
	my ($drive, $dir, $filename, $ext, $lastdir, $filebody)
		= Deps::SplitFilePath($path);
	my $SampleName = $filebody;

print "  TranSIESTA::ReadEnergyConvergenceFromATK2Output: Reading [$path].\n";
	my $in = new JFile($path, "r");
	unless($in) {
		print "  Error: Can not ead [$path].\n";
		return undef;
	}

	$this->SetDataArray(new GraphDataArray);
	my $Data = new GraphData;
	$Data->SetTitle($SampleName);

	if(defined $target and $target =~ /ForceConvergence/i) {
	}
	else { #'EnergyConvergence'
		my @Iteration;
		my @Energy;
		my $iter = 0;
		$in->rewind();
		while(!$in->eof()) {
			my $line = $in->SkipTo("^# sc ");
			last if($in->eof());
			my ($DETot) = ($line =~ /dEbs =\s+([\d\.+\-Ee]+)/i);
			$Energy[$iter] = $DETot;
			$Iteration[$iter] = $iter+1;
print "i=$iter: dE(total)=$Energy[$iter]\n";
			$iter++;
		}
		$Data->{'x0'}      = \@Iteration;
		$Data->{'x0_Name'} = "Iteration";
		$Data->{'y0'}      = \@Energy;
		$Data->{'y0_Name'} = "Energy / eV";
	}

	$in->Close();

	$Data->CalMinMax();
	$this->{'DataArray'}->AddGraphData($Data);
	$this->{'DataArray'}->{'TargetData'} = $target;


	$in->Close();
print "  TranSIESTA::ReadEnergyConvergenceFromATK2Output: finished.\n";
	return $path;
}

sub ReadBandFromATK2Output
{
	my ($this, $path) = @_;
	my ($drive, $dir, $filename, $ext, $lastdir, $filebody)
		= Deps::SplitFilePath($path);
	my $SampleName = $filebody;

print "  Read [$path].\n";
	my $in = new JFile($path, "r");
	unless($in) {
		print "  Error: Can not ead [$path].\n";
		return undef;
	}

	my $SimulationType = $in->SkipTo("SimulationType::Type ");
	$SimulationType =~ s/^.*\s(\S+)[\s\r\n]*?$/$1/;
	chomp($SimulationType);
print "SimulationType: $SimulationType\n";
	$in->rewind();
	my $line = $in->SkipTo("# BandLine : Number of segments ");
	my $nBand;
	$nBand = ($nBand =~ /\s(\d+)\s*?$/) if($line);
print "nBand: $nBand\n" if(defined $nBand);
	$in->rewind();
	$line = $in->SkipTo("# FermiEnergy : ");
	my ($EF, $unit);
	($EF, $unit) = ($EF =~ /:\s*([\d\.+-eE]+?)\s+(\w+)\s*$/) if($line);
print "EF: $EF $unit\n" if(defined $unit);

	my @BandArray;
	my @BandBoundaryDistance;
	my @BoundaryPositions;
	my $DistanceOffset = 0.0;
	my $nLevels = 0;
#$ib番目の計算したk線のバンドの組を読み込むループ
	for(my $ib = 0 ; $ib < $nBand ; $ib++) {
#次のバンドデータの先頭に位置づける
		my $nK = $in->SkipTo("# Number of points = ");
		($nK) = ($nK =~ /\s(\d+)\s*$/);
#print "iBand=$ib  nK=$nK\n";
		$line = $in->SkipTo("# -----");

		my ($d, $e, $ne, $kx, $ky, $kz);
		my $MaxD = 0.0;
		my $strK;
		my $iBand = 0;
#$ib番面に計算したバンドの、異なるエネルギー準位のバンドを読み込むループ
		while(!$in->eof()) {
			my $i;
			my $d0 = 0.0;
#あるバンドのk点の組を読み込む
			for($i = 0 ; $i < $nK ; $i++) {
				$line = $in->ReadLine();
				last if($in->eof());
				last if($line =~ /# -----/);
				Utils::DelSpace($line);
				last if($line eq '');

				$BandArray[$nLevels] = new EnergyBand 
						unless($BandArray[$nLevels]);
				my ($d, $e, $ne, $kx, $ky, $kz) = Utils::Split("\\s+", $line);
				if($i == 0) {
					$d0 = $d;
				}
				$d = $d - $d0;
#if($iBand == 0) {
#	print "$ib: $iBand: $i: $d\n";
#}
				$MaxD = $d if($MaxD < $d);
				$strK = "$kx $ky $kz";
				$BandArray[$nLevels]->AddByKDEN($kx, $ky, $kz, $d+$DistanceOffset, $e, $ne);

# バンド境界の、最初の情報を入れておく
				if($nLevels == 0 and $i == 0) {
					push(@BandBoundaryDistance, 0.0);
					push(@BoundaryPositions, $strK);
				}
			}
#次のK点群のバンドを読み込む
			last if($line =~ /# -----/);
			if($i >= $nK-1) {
				$nLevels++;
				$iBand++;
			}
		}
#$MaxD: $ib番目に計算したバンドでの、最大のΔk
		$DistanceOffset += $MaxD;
		push(@BandBoundaryDistance, $DistanceOffset);
		push(@BoundaryPositions, $strK);

#print "ib=$ib MaxD: $MaxD / Off: $DistanceOffset\n";
	}
	$in->Close();

	$this->SetDataArray(new GraphDataArray);
	my $Data = new GraphData;
	$Data->SetTitle($SampleName);

	my $band = $BandArray[0];
	for(my $i = 0 ; ; $i++) {
		$band = $BandArray[$i];
		last unless($band);
		$band->CalMinMax();
		$Data->SetXDataArray($i, $band->pDistance());
		$Data->SetYDataArray($i, $band->pEnergy());
		$Data->{"Ne$i"} = $band->pNe();
		my ($ymin, $ymax) = $band->GetYMinMax();
#		my $Ne = $band->Ne(0); #Ne($i);
#		my $s = "Band $i ($ymin - $ymax) Ne=$Ne";
		my $s = "Band $i ($ymin - $ymax)";
		$Data->SetXName(0, "k");
		$Data->SetYName($i, $s);
	}

	$Data->CalMinMax();
	$this->DataArray()->{'BandBoundaryDistances'} = \@BandBoundaryDistance;
	$this->DataArray()->{'BandBoundaryPositions'} = \@BoundaryPositions;
	$this->DataArray()->AddGraphData($Data);

	return $path;
}

sub ReadFiles
{
	my ($this, $filename, $target) = @_;
	$this->ClearAll();
	my $FileType = $this->{'FileType'} = CheckFileType($filename);
	return undef unless($FileType);
print("  TranSIESTA::ReadFiles: Read [$filename] for [$target]\n") if($target);
print("  TranSIESTA::ReadFiles: Read [$filename]\n") unless($target);

	$this->SetFileName($filename);
	if($FileType eq "ATK 2.0 Output") {
		if(not defined $target or $target eq 'EnergyConvergence') {
			return $this->ReadEnergyConvergenceFromATK2Output($filename);
		}
		elsif($target =~ /ForceConvergence/i) {
			return $this->ReadForceConvergenceFromATK2Output($filename);
		}
		elsif($target eq 'EnergyBand') {
			return $this->ReadBandFromATK2Output($filename);
		}
	}
	return undef;
}

#============================================================
#　一般関数
#============================================================

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

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

sub IsCrystalATKFile
{
	my ($this, $filename) = @_;
	my $in = new JFile($filename, "r");
	return undef unless($in);
	while(!$in->eof()) {
		my $line = $in->ReadLine();
		if($line =~ /^\s*UnitCell::/i) {
			$in->Close();
			return 1;
		}
	}
	$in->Close();
	return 0;
}

sub IsTwoProbeATKFile
{
	my ($this, $filename) = @_;
	my $in = new JFile($filename, "r");
	return undef unless($in);
	while(!$in->eof()) {
		my $line = $in->ReadLine();
		if($line =~ /^\s*System::Type\s+TwoProbe/i) {
			$in->Close();
			return 1;
		}
	}
	$in->Close();
	return 0;
}

sub ReadATKFile
{
	my ($this, $filename, $ElectrodeSpacing, $IsPrint) = @_;
#print "ReadATKFile: $filename\n";
	return $this->ReadCrystalATKFile($filename, $IsPrint) 
			if($this->IsCrystalATKFile($filename));
	return $this->ReadTwoProbeATKFile($filename, $ElectrodeSpacing, $IsPrint) 
			if($this->IsTwoProbeATKFile($filename));
	return -1;
}

sub ReadTwoProbeATKFile
{
	my ($this, $filename, $ElectrodeSpacing, $IsPrint) = @_;
	$ElectrodeSpacing = 0 unless(defined $ElectrodeSpacing);
	$IsPrint = 1 unless(defined $IsPrint);
	my $a0 = Sci::a0();
#print "$filename is TwoProbe file.\n";

	my $crystal = new Crystal();
	my $LCrystal = new Crystal();
	my $RCrystal = new Crystal();

	my ($drive, $directory, $fname, $ext, $lastdir, $filebody)
		= Deps::SplitFilePath($filename);
	$crystal->SetCrystalName($filebody);
	$crystal->SetSampleName($filebody);

	my $in = new JFile($filename, "r");
	return undef unless($in);

#電極ファイルの読み込み
	my $line = $in->SkipTo("TwoProbe::LeftElectrode::ArchiveFile");
	my ($LeftElectrodeFile) = ($line =~ /(\S+)\s*?$/);
	$LeftElectrodeFile = Deps::ReplaceExtension($LeftElectrodeFile, ".atk");
	$in->rewind();
	$line = $in->SkipTo("TwoProbe::RightElectrode::ArchiveFile");
	my ($RightElectrodeFile) = ($line =~ /(\S+)\s*?$/);
	$RightElectrodeFile = Deps::ReplaceExtension($RightElectrodeFile, ".atk");

print "LeftElectrodeFile: $LeftElectrodeFile\n";
	$LCrystal = $this->ReadElectrodeATKFile($LeftElectrodeFile, $IsPrint);
print "RightElectrodeFile: $RightElectrodeFile\n";
	$RCrystal = $this->ReadElectrodeATKFile($RightElectrodeFile, $IsPrint);
print "\n";

#中心系、電極系の原子情報
	my (@AtomName, @XC, @YC, @ZC);

#格子定数の決定：まず、zの最大・最小値を調べる
	$in->rewind();
	$line = $in->SkipTo("AtomList::Format");
	my ($unit) = ($line =~ /(\S+)\s*?$/);
print "unit=$unit\n";
	my $l0 = 1.0;
	$l0 = $a0*1.0e10 unless($unit =~ /ang/i);

#%block TwoProbe::CentralAtomList
	$in->ReadLine();
	my $ZMin = 1.0e10;
	my $ZMax = -1.0e10;
	while(!$in->eof()) {
#       Au      1.442   0.8328  8.245
		$line = $in->ReadLine();
#	%endblock TwoProbe::CentralAtomList
		last if($line =~ /^\s*%endblock /i);
		my ($name, $xc, $yc, $zc) = ($line =~ /^\s*(\S+)\s+(\S+)\s+(\S+)\s+(\S+)/);
		last unless(defined $zc);

		$ZMin = $zc*$l0 if($ZMin > $zc);
		$ZMax = $zc*$l0 if($ZMax < $zc);
		push(@AtomName, $name);
		push(@XC, $xc*$l0);
		push(@YC, $yc*$l0);
		push(@ZC, $zc*$l0);
	}
# LeftElectrodeとあわせる原子
	my ($XFirst, $YFirst, $ZFirst) = ($XC[0], $YC[0], $ZC[0]);
# RightElectrodeとあわせる原子
	my $n1 = @XC - 1;
	my ($XLast, $YLast, $ZLast) = ($XC[$n1], $YC[$n1], $ZC[$n1]);

print "ZMinMax: $ZMin - $ZMax\n";
	my ($La, $Lb, $Lc, $Lalpha, $Lbeta, $Lgamma) = $LCrystal->LatticeParameters();
	my ($Ra, $Rb, $Rc, $Ralpha, $Rbeta, $Rgamma) = $RCrystal->LatticeParameters();
#格子定数は、両電極の和と中心領域のZ範囲、$ElectrodeSpacingの和
	my $c = $ZMax - $ZMin + $Lc + $Rc + $ElectrodeSpacing;
	$crystal->SetLatticeParameters($Ra, $Rb, $c, $Ralpha, $Rbeta, $Rgamma);

# 電極リストの追加
	if($LCrystal) {
		my @Sites = $LCrystal->GetCAsymmetricAtomSiteList();
		my $nAtom = $LCrystal->nAsymmetricAtomSite();
		my $dx = $XLast - $LCrystal->{"OriginalX[0]"};
		my $dy = $YLast - $LCrystal->{"OriginalY[0]"};
		my $dz = $ZLast - $LCrystal->{"OriginalZ[0]"};
print "  Shift for LeftElectrode: ($dx, $dy, $dz)\n";
		for(my $i = 0 ; $i < $nAtom ; $i++) {
			my $site = $Sites[$i];
			push(@AtomName, $site->AtomName());
			push(@XC, $LCrystal->{"OriginalX[$i]"} + $dx);
			push(@YC, $LCrystal->{"OriginalY[$i]"} + $dy);
			push(@ZC, $LCrystal->{"OriginalZ[$i]"} + $dz);
			my $zmin = $LCrystal->{"OriginalZ[$i]"};
			$ZMin = $zmin if($ZMin > $zmin);
		}
	}
	if($RCrystal) {
		my @Sites = $RCrystal->GetCAsymmetricAtomSiteList();
		my $nAtom = $RCrystal->nAsymmetricAtomSite();
		my $n1 = $nAtom-1;
		my $dx = $XFirst - $RCrystal->{"OriginalX[$n1]"};
		my $dy = $YFirst - $RCrystal->{"OriginalY[$n1]"};
		my $dz = $ZFirst - $RCrystal->{"OriginalZ[$n1]"};
print "  Shift for RightElectrode: ($dx, $dy, $dz)\n";
		for(my $i = 0 ; $i < $nAtom ; $i++) {
			my $site = $Sites[$i];
			push(@AtomName, $site->AtomName());
			push(@XC, $RCrystal->{"OriginalX[$i]"} + $dx);
			push(@YC, $RCrystal->{"OriginalY[$i]"} + $dy);
			push(@ZC, $RCrystal->{"OriginalZ[$i]"} + $dz);
			my $zmin = $RCrystal->{"OriginalZ[$i]"};
			$ZMin = $zmin if($ZMin > $zmin);
		}
	}
#my $n = @AtomName;
#print "n=$n\n";
#for(my $i = 0 ; $i < $n ; $i++) {
#	print $AtomName[$i], $XC[$i], $YC[$i], $ZC[$i], "\n";
#}

# 原子リストを$crystalに設定する
	my %AtomCount;
	for(my $i = 0 ; $i < @AtomName ; $i++) {
		my $name = $AtomName[$i];
		my $xc = $XC[$i];
		my $yc = $YC[$i];
		my $zc = $ZC[$i] - $ZMin;

		$crystal->AddAtomType($name);
		$AtomCount{$name}++;
		my $Label = $name . $AtomCount{$name};
		my ($x, $y, $z) = $crystal->CartesianToFractional($xc, $yc, $zc);
		next if(defined $crystal->FindIdenticalAsymmetricSite($name, $x, $y, $z));
		$x = 0 if(abs($x) < 1.0e-6);
		$y = 0 if(abs($y) < 1.0e-6);
		$z = 0 if(abs($z) < 1.0e-6);
		$x = CrystalObject::RoundSymmetricPosition($x);
		$y = CrystalObject::RoundSymmetricPosition($y);
		$z = CrystalObject::RoundSymmetricPosition($z);
		$x = Utils::Round($x, 6);
		$y = Utils::Round($y, 6);
		$z = Utils::Round($z, 6);
		$crystal->AddAtomSite($Label, $name, $x, $y, $z, 1.0);
print "$name:$Label ($x,$y,$z)\n";
	}

	$in->Close();

	return ($crystal, $LCrystal, $RCrystal);
}

sub ReadElectrodeATKFile
{
	my ($this, $filename, $IsPrint) = @_;
	$IsPrint = 1 unless(defined $IsPrint);
	my $a0 = Sci::a0();

	my $crystal = new Crystal();
	my ($drive, $directory, $fname, $ext, $lastdir, $filebody)
		= Deps::SplitFilePath($filename);
	$crystal->SetCrystalName($filebody);
	$crystal->SetSampleName($filebody);

	my $in = new JFile($filename, "r");
	return undef unless($in);

	my $line = $in->SkipTo("%block electrode::UnitCell");
	my $l0 = 1.0;
#2.88499566724 0.0 0.0
	$line = $in->ReadLine();
	my ($l11, $l12, $l13) = ($line =~ /\s*(\S+)\s+(\S+)\s+(\S+)/);
#-1.44249783362 2.49847953764 0.0
	$line = $in->ReadLine();
	my ($l21, $l22, $l23) = ($line =~ /\s*(\S+)\s+(\S+)\s+(\S+)/);
#0.0 0.0 7.06676729488
	$line = $in->ReadLine();
	my ($l31, $l32, $l33) = ($line =~ /\s*(\S+)\s+(\S+)\s+(\S+)/);
#printf("%12s %7.4f %7.4f %7.4f\n", "LattVectr:", $l11*$l0, $l12*$l0, $l13*$l0);
#printf("%12s %7.4f %7.4f %7.4f\n", "",           $l21*$l0, $l22*$l0, $l23*$l0);
#printf("%12s %7.4f %7.4f %7.4f\n", "",           $l31*$l0, $l32*$l0, $l33*$l0);
	$crystal->SetLatticeVector($l11*$l0, $l12*$l0, $l13*$l0, 
				   $l21*$l0, $l22*$l0, $l23*$l0, 
				   $l31*$l0, $l32*$l0, $l33*$l0);
	my ($a, $b, $c, $alpha, $beta, $gamma)
			= $crystal->CalculateLatticeConstantFromVector();
print "Lattice Parameters: $a $b $c  $alpha $beta $gamma\n";

	$in->rewind();
	$line = $in->SkipTo("%block electrode::AtomList");
	$l0 = 1.0;
	my %AtomCount;
	my $count = 0;
	while(!$in->eof()) {
#Au 1.44249783362 0.832826512546 1.17779454915
		$line = $in->ReadLine();
#print "line: $line\n";
#%endblock electrode::AtomList
		last if($line =~ /^\s*%endblock /i);
		my ($name, $xc, $yc, $zc) = ($line =~ /^\s*(\S+)\s+(\S+)\s+(\S+)\s+(\S+)/);
		last unless(defined $zc);

		$crystal->AddAtomType($name);
		$AtomCount{$name}++;
		my $Label = $name . $AtomCount{$name};
		$crystal->{"OriginalX[$count]"} = $xc;
		$crystal->{"OriginalY[$count]"} = $yc;
		$crystal->{"OriginalZ[$count]"} = $zc;
		my ($x, $y, $z) = $crystal->CartesianToFractional($xc, $yc, $zc);
		$x = 0 if(abs($x) < 1.0e-6);
		$y = 0 if(abs($y) < 1.0e-6);
		$z = 0 if(abs($z) < 1.0e-6);
		$x = CrystalObject::RoundSymmetricPosition($x);
		$y = CrystalObject::RoundSymmetricPosition($y);
		$z = CrystalObject::RoundSymmetricPosition($z);
		$x = Utils::Round($x, 6);
		$y = Utils::Round($y, 6);
		$z = Utils::Round($z, 6);
		$crystal->AddAtomSite($Label, $name, $x, $y, $z, 1.0);
print "$name:$Label ($x,$y,$z)\n";
		$count++;
	}
	$crystal->SetLatticeParameters($a, $b, $c, $alpha, $beta, $gamma);

	$in->Close();
	$crystal->ExpandCoordinates();

	return $crystal;
}

sub ReadCrystalATKFile
{
	my ($this, $filename, $IsPrint) = @_;
	$IsPrint = 1 unless(defined $IsPrint);
	my $a0 = Sci::a0();

	my $crystal = new Crystal();
	my ($drive, $directory, $fname, $ext, $lastdir, $filebody)
		= Deps::SplitFilePath($filename);
	$crystal->SetCrystalName($filebody);
	$crystal->SetSampleName($filebody);

	my $in = new JFile($filename, "r");
	return undef unless($in);

	my $line = $in->SkipTo("UnitCell::");
#print "line: $line\n";
	my ($l0, $unit) = ($line =~ /LatticeConstant\s*([\d\.]+)\s+(\w+)/);
#print "l0=$l0 $unit\n";
	$l0 *= $a0*1.0e10 if($unit !~ /ang/i);
print "l0=$l0 Ang\n";

#%block Bulk::UnitCell
	$line = $in->ReadLine();

#  1.000000   0.000000   0.000000
	$line = $in->ReadLine();
	my ($l11, $l12, $l13) = ($line =~ /\s*(\S+)\s+(\S+)\s+(\S+)/);
# -0.500000   0.866025   0.000000
	$line = $in->ReadLine();
	my ($l21, $l22, $l23) = ($line =~ /\s*(\S+)\s+(\S+)\s+(\S+)/);
#  0.000000   0.000000   1.601998
	$line = $in->ReadLine();
	my ($l31, $l32, $l33) = ($line =~ /\s*(\S+)\s+(\S+)\s+(\S+)/);
#printf("%12s %7.4f %7.4f %7.4f\n", "LattVectr:", $l11*$l0, $l12*$l0, $l13*$l0);
#printf("%12s %7.4f %7.4f %7.4f\n", "",           $l21*$l0, $l22*$l0, $l23*$l0);
#printf("%12s %7.4f %7.4f %7.4f\n", "",           $l31*$l0, $l32*$l0, $l33*$l0);
	$crystal->SetLatticeVector($l11*$l0, $l12*$l0, $l13*$l0, 
				   $l21*$l0, $l22*$l0, $l23*$l0, 
				   $l31*$l0, $l32*$l0, $l33*$l0);
	my ($a, $b, $c, $alpha, $beta, $gamma)
			= $crystal->CalculateLatticeConstantFromVector();
print "Lattice Parameters: $a $b $c  $alpha $beta $gamma\n";

	$in->rewind();
	$line = $in->SkipTo("AtomList::Format");
	($unit) = ($line =~ /(\S+)\s*?$/);
#print "Unit for position: $unit\n";
	$l0 = 1.0;
	$l0 = $a0*1.0e10 if($unit !~ /ang/i);
#%block Bulk::AtomList
	$in->ReadLine();
	my %AtomCount;
	my $count = 0;
	while(!$in->eof()) {
#Zn    0.000000000    1.872173699    0.000000000 
		$line = $in->ReadLine();
#print "line: $line\n";
#%endblock Bulk::AtomList
		last if($line =~ /^%endblock /i);
		my ($name, $xc, $yc, $zc) = ($line =~ /^\s*(\S+)\s+(\S+)\s+(\S+)\s+(\S+)/);
		last unless(defined $zc);

		$crystal->AddAtomType($name);
		$AtomCount{$name}++;
		my $Label = $name . $AtomCount{$name};
		$crystal->{'OriginalX[$count]'} = $xc;
		$crystal->{'OriginalY[$count]'} = $yc;
		$crystal->{'OriginalZ[$count]'} = $zc;
		my ($x, $y, $z) = $crystal->CartesianToFractional($xc, $yc, $zc);
		$x = 0 if(abs($x) < 1.0e-6);
		$y = 0 if(abs($y) < 1.0e-6);
		$z = 0 if(abs($z) < 1.0e-6);
		$x = CrystalObject::RoundSymmetricPosition($x);
		$y = CrystalObject::RoundSymmetricPosition($y);
		$z = CrystalObject::RoundSymmetricPosition($z);
		$x = Utils::Round($x, 6);
		$y = Utils::Round($y, 6);
		$z = Utils::Round($z, 6);
		$crystal->AddAtomSite($Label, $name, $x, $y, $z, 1.0);
print "$name:$Label ($x,$y,$z)\n";
		$count++;
	}
	$crystal->SetLatticeParameters($a, $b, $c, $alpha, $beta, $gamma);

	$in->Close();
	$crystal->ExpandCoordinates();

	return $crystal;
}

sub SaveTranSIESTAATKFile
{
	my ($this, $Crystal, $filename, 
		$Function, $BasisSet, $Version, 
		$XCFunctional, $InitialSpin, $IsChooseRandomly) = @_;
	my $LF = "<br>\n";

#ファイル作製開始
	unless(open(OUT,">$filename")) {
		print "Can not write to $filename.$LF$LF";
		return;
	}
	if($Function =~ /relax/i) {
		print OUT "Simulation::Type Relaxation\n";
		print OUT "Simulation::ConstrainAtoms 0\n";
		print OUT "Relaxation::ForceTolerance 5.e-2 eV/Ang\n";
		print OUT "Relaxation::MaxSteps 100\n";
		print OUT "Relaxation::MaxDisplacement 0.2 Bohr\n";
	}
	else {
		print OUT "Simulation::Type SingleConfiguration\n";
	}
	print OUT "System::Type Bulk\n";
	print OUT "Method::Type NumOrb\n";
	print OUT "NumOrb::MeshCutoff 100 Ry\n";
	if($Version eq 'Ver11' ) {
		print OUT "NumOrb::BasisType $BasisSet\n";
	}
	elsif($Version eq 'Ver20') {
		print OUT "NumOrb::BasisSet::Size $BasisSet\n";
		print OUT "NumOrb::XCFunctional $XCFunctional\n";
		if($XCFunctional =~ /s/i) {
			print OUT "SCF::InitialSpinpolarization $InitialSpin Hbar\n";
		}
	}
	print OUT "\n";

#	my $lattice = $CIFContent{"_symmetry_cell_setting"};
	my $SPGName = $Crystal->SPGNameByOutputMode();
	my $iSPG    = $Crystal->iSPGByOutputMode();
	my ($a,$b,$c,$alpha,$beta,$gamma) 
		= $Crystal->LatticeParametersByOutputMode(0);

	print OUT "UnitCell::LatticeConstant $a Ang\n";
	print OUT "%block Bulk::UnitCell\n";
	my ($ax,$ay,$az,$bx,$by,$bz,$cx,$cy,$cz) 
		= $Crystal->LatticeVectorsByOutputMode(0);
	printf OUT "%10.6f %10.6f %10.6f\n", $ax/$a, $ay/$a, $az/$a;
	printf OUT "%10.6f %10.6f %10.6f\n", $bx/$a, $by/$a, $bz/$a;
	printf OUT "%10.6f %10.6f %10.6f\n", $cx/$a, $cy/$a, $cz/$a;
	print OUT "%endblock Bulk::UnitCell\n";
	print OUT "\n";

	my @ExpandedAtomSiteList = $Crystal->GetCExpandedAtomSiteListByOutputMode();
	my $nTotalExpandedAtomSite = @ExpandedAtomSiteList;
#	my $nAsymmetricAtomSite = &GetCIFValueByKey("nAsymmetricAtomSite");

	print OUT "AtomList::Format         Ang\n";
	print OUT "%block Bulk::AtomList\n";

#Occupancyが1のサイト
	my $count = 0;
	for(my $i = 0 ; $i < @ExpandedAtomSiteList ; $i++) {
		my $atom      = $ExpandedAtomSiteList[$i];
		my $atomname  = $atom->AtomNameOnly();
		my $charge    = $atom->Charge();
		my ($x,$y,$z) = $atom->Position(1);
		my $occupancy = $atom->Occupancy();
		next if($occupancy < 0.9999);

		$count++;
		my ($xc,$yc,$zc) = $Crystal->FractionalToCartesianByOutputMode($x,$y,$z);

		printf OUT "%-2s %14.9f %14.9f %14.9f \n",
			$atomname, $xc, $yc, $zc;
	}
#Occupancyが1未満のサイト
	if(@ExpandedAtomSiteList > $count and !$IsChooseRandomly) {
		print "<H3>***Warning***: One line should be removed from $filename.</H3>\n";
		print OUT "*** Shared atoms: Remove this line for calculation: OCC=";
		for(my $i = 0 ; $i < @ExpandedAtomSiteList ; $i++) {
			my $atom      = $ExpandedAtomSiteList[$i];
			my $occupancy = $atom->Occupancy();
			next if($occupancy >= 0.9999);
			printf OUT "%lg ", $occupancy;
		}
		print OUT "\n";
	}
	for(my $i = 0 ; $i < @ExpandedAtomSiteList ; $i++) {
		my $atom      = $ExpandedAtomSiteList[$i];
		my $atomname  = $atom->AtomNameOnly();
		my $charge    = $atom->Charge();
		my ($x,$y,$z) = $atom->Position(1);
		my $occupancy = $atom->Occupancy();
		next if($occupancy >= 0.9999);

		my ($xc,$yc,$zc) = $Crystal->FractionalToCartesianByOutputMode($x,$y,$z);

		printf OUT "%-2s %14.9f %14.9f %14.9f \n",
			$atomname, $xc, $yc, $zc;
	}

	print OUT "%endblock Bulk::AtomList\n";
	print OUT "\n";

#	&ListPartialOccupancyAtoms();
	my $an = 10;
	my $nx = int($an / $a) + 1;
	my $ny = int($an / $b) + 1;
	my $nz = int($an / $c) + 1;
	print OUT "Bulk::NumKPoints::A $nx\n";
	print OUT "Bulk::NumKPoints::B $ny\n";
	print OUT "Bulk::NumKPoints::C $nz\n";

#	if($Function =~ /relax/i) {
#		print OUT "Bulk::NumKPoints::A 1\n";
#		print OUT "Bulk::NumKPoints::B 1\n";
#		print OUT "Bulk::NumKPoints::C 1\n";
#	}
#	else {
#		print OUT "Bulk::NumKPoints::A 5\n";
#		print OUT "Bulk::NumKPoints::B 5\n";
#		print OUT "Bulk::NumKPoints::C 5\n";
#	}
	print OUT "\n";

	my $XC = new XCrySDen;
	my $KListDBDir = $XC->SetKListDBDir($this->KListDBDir);
#print "klist DB dir: $KListDBDir$LF";
	my $SPGName0 = $Crystal->SPGName();
	my $iSPG0    = $Crystal->iSPG();
	my ($KListFilePath, @KList) = $XC->ReadKList($iSPG0, $SPGName0);
#print "KList DB file: $KListFilePath$LF\n";

	if($Function !~ /relax/i) {
		if($Version eq 'Ver11' ) {
			print OUT "Analysis::BandLine::K0 (1.0,0.0,0.0) FracRecip \n";
			print OUT "Analysis::BandLine::K1 (0.0,0.0,0.0) FracRecip \n";
			print OUT "Analysis::BandLine::NumPoints0 21\n";
		}
		elsif($Version eq 'Ver20' ) {
#			print OUT "Analysis::BandLine::A::Start (1.0,0.0,0.0) FracRecip \n";
#			print OUT "Analysis::BandLine::A::End   (0.0,0.0,0.0) FracRecip \n";
#			print OUT "Analysis::BandLine::A::NumPoints 21\n";
#			print OUT "\n";
#			print OUT "Analysis::BandLine::B::Start (0.0,0.0,0.0) FracRecip \n";
#			print OUT "Analysis::BandLine::B::End   (0.0,0.0,1.0) FracRecip \n";
#			print OUT "Analysis::BandLine::B::NumPoints 21\n";
			my $nKList = @KList - 1;
			my $Char = 'A';
			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'};
				print OUT "Analysis::BandLine::${Char}::Start ($kx1,$ky1,$kz1) FracRecip \n";
				print OUT "Analysis::BandLine::${Char}::End   ($kx2,$ky2,$kz2) FracRecip \n";
				print OUT "Analysis::BandLine::${Char}::NumPoints $nPoint\n";
				print OUT "\n";
				$Char++;
			}
		}
		print OUT "\n";
	}

	close(OUT);

	return ($KListFilePath, $filename);
}

1;