#!/usr/bin/perl

#require "Util.pl";
#require "MyCGI.pl";
use Cwd;

use lib './lib';
use lib 'd:/tkProg/tklib/Perl/lib';
use lib 'd:/Programs/Perl/lib';

use strict;
#use warnings;
#use File::Path;
use File::Basename;
#use File::Find;

use ProgVars;
use Crystal::AtomType;
use Utils;

#===============================================
# デバッグ関係変数
#===============================================
my $Debug = 0;
#$PrintLevelが大きいほど、情報が詳しくなる
my $PrintLevel = 0;

#===============================================
# スクリプト大域変数
#===============================================
my $DirectorySeparator = "\\";

my $InputFromSTDIN = 1;

my $Action = "SCFCalcDVXa";
my $CalcAtomOnly = 0;

my $BaseDir = "";

my $SpinPolarizedCalculation = 0;
my $ConvergenceThreshold = 1.0e-4;
my $MaxLoop = 5;

my $MinO1sEnergy = -480;
my $MaxO1sEnergy = -420;
my $O1sEnergy    = -521.6;

#===============================================
# 文字コード関係変数
#===============================================
# sjis, euc, jis, noconv
my $FileSystemCharCode = "sjis";
my $MySQLCharCode      = "sjis";
my $WebCharSet         = "sjis";
my $WebCharCode        = "sjis";

#===============================================
# パス（読み込みDB）
# Web関係変数
# CGI の仮想パス
# プログラム名
#===============================================
my $Program   = basename($0);
my $ProgramDir = dirname($0);

my $EditorPath = "notepad.exe %1";

my $DBDir  = ProgVars::DBDir();
my $BinDir = ProgVars::BinDir();

my $DVXaDir   = ProgVars::DVXaDir(); #"d:\\Programs\\Quantum\\DVXa97";

my $DVXaExecDir     = &MakePath($DVXaDir, "exec", $DirectorySeparator, 0);

my $DVXaSpinPath    = &MakePath($DVXaExecDir, "gspin-P4.exe", $DirectorySeparator, 0);
if($SpinPolarizedCalculation) {
	$DVXaSpinPath   = &MakePath($DVXaExecDir, "hspin-P4.exe", $DirectorySeparator, 0);
}

my $MakeF05Path     = &MakePath($DVXaExecDir, "makef05.exe", $DirectorySeparator, 0);

my $NonrelInputPath  = &MakePath($DBDir, "nonrel", $DirectorySeparator, 0);
my $RadiiInputPath   = &MakePath($DBDir, "radii",  $DirectorySeparator, 0);
my $F27InputPath     = &MakePath($DVXaDir, "blank", $DirectorySeparator, 0);
my $SpaceGroupDBPath = &MakePath($DBDir, "SPGRA", $DirectorySeparator, 0);

my $PeriodicTableDir     = ProgVars::PeriodicTableDir();
my $IonRadius1DBPath     = &MakePath($PeriodicTableDir, "IonRadius1.html", $DirectorySeparator, 0);
my $IonRadius2DBPath     = &MakePath($PeriodicTableDir, "IonRadius2.html", $DirectorySeparator, 0);
my $PeriodicTable1DBPath = &MakePath($PeriodicTableDir, "PeriodicTable1.html",
								$DirectorySeparator, 0);
my $PeriodicTable2DBDir = &MakePath($PeriodicTableDir, "PeriodicTable2",
								$DirectorySeparator, 0);

#==========================================
# メイン関数スタート
#==========================================

print "=================================================================\n";
print "=================================================================\n";
print "    $Program\n";
print "        Copyright by T. Kamiya, 2005\n";
print "=================================================================\n";
print "=================================================================\n";
print "\n";

my $InitialDirectory = cwd(); 
$InitialDirectory =~ s/\//$DirectorySeparator/g;
print "Initial Directory: $InitialDirectory\n";
print "\n";

#==========================================
# 引数を取得する
# parseInput関数はutil.plにある
#==========================================

for(my $i = 0 ; $i < @ARGV ; $i++) {
	if($ARGV[$i] =~ /^--CalcAtomOnly=(.*)$/i) {
		$CalcAtomOnly = $1;
	}
	elsif($ARGV[$i] =~ /^--InputFromSTDIN=(.*)$/i) {
		$InputFromSTDIN = $1;
	}
	elsif($ARGV[$i] =~ /^--SpinPolarizedCalculation=(.*)$/i) {
		$SpinPolarizedCalculation = $1;
	}
	elsif($ARGV[$i] =~ /^--Action=(.*)$/i) {
		$Action = $1;
	}
	elsif($ARGV[$i] =~ /^--BaseDir=(.*)$/i) {
		$BaseDir = $1;
	}
	elsif($ARGV[$i] =~ /^--Debug=(.*)$/i) {
		$Debug = $1;
	}
	elsif($ARGV[$i] =~ /^--help$/i) {
		print "usage: $Program [Options]\n";
		print "  Options:\n";
		print "    --Action=[SCFCalcDVXa|NormalizeF08e|\n";
		print "              ShowSpaceGroup|ShowIonInf|\n";
		print "              MakeLD01File|ShowCIFInf]\n";
		print "    --CalcAtomOnly=[0|1]\n";
		print "    --BaseDir=Dir\n";
		print "    --SpinPolarizedCalculation=[0|1]\n";
		print "    --InputFromSTDIN=[0|1]\n";
		print "    --Debug=[0|1]\n";
		exit;
	}
	else {
		print "Invalid argument: $ARGV[$i]\n";
		print "  Use \"$Program --help\" for help.\n";
		exit;
	}
}

print "[Parameters]\n";
print "  DVXaDir       : $DVXaDir\n";
print "  Action        : $Action\n";
print "  BaseDir       : $BaseDir\n";
print "  InputFromSTDIN: $InputFromSTDIN\n";
print "  CalcAtomOnly  : $CalcAtomOnly\n";
print "  Debug         : $Debug\n";
print "  SpinPolarizedCalculation: $SpinPolarizedCalculation\n";
print "\n";

if($Action =~ /^SCFCalcDVXa$/i) {
	&ExecuteSCFCalcDVXa();
}
elsif($Action =~ /^NormalizeF08e$/i) {
	&ExecuteNormalizeF08e();
}
elsif($Action =~ /^ShowSpaceGroup$/i) {
	&ShowSpaceGroup();
}
elsif($Action =~ /^ShowIonInf$/i) {
	&ShowIonInformation();
}
elsif($Action =~ /^MakeLD01File$/i) {
	&MakeLD01File();
}
elsif($Action =~ /^ShowCIFInf$/i) {
	&ShowCIFInf();
}
else {
	print "\n";
	print "Invalid Action: $Action\n";
	exit;
}

exit;

#===============================================
# スクリプト終了
#===============================================

#==========================================
# &Subroutines
#==========================================

sub ShowCIFInf
{
	my $CIFFile = "CaF2.cif";
	print "CIF File[$CIFFile]=";
	my $tmp = <STDIN>;
	&DelSpace($tmp);
	$CIFFile = $tmp if($tmp ne '');
	print "CIF File: $CIFFile\n";

	my %Content;
	open(IN, "<$CIFFile") or die "$!: $CIFFile\n";

#最初の一行
	my $header = <IN>;

	while(<IN>) {
		my $line = $_;
		&DelSpace($line);
# 空行はスキップ
		next if($line =~ /^$/);
# #Endを見つけたら終わる
		last if($line =~ /^#End/i);

#loop_の無いデータを読み込み
		if($line =~ /^_/) {
			my ($key, $content) = ($line =~ /^(\S+?)\s+(.*)$/);
			if($key ne '') {
				$Content{$key} = $content;
				print "$key: $content\n";
				next;
			}
		}
#loop_のデータを読み込み
		elsif($line =~ /^loop_$/i) {
			my @Params;
			my $Count = 0;
			my $line2 = '';
#loop_内のパラメータ名を取得
			while(<IN>) {
				$line2 = $_;
				&DelSpace($line2);
				last unless($line2 =~ /^_/);

				$Params[$Count] = $line2;
				$Count++;
			}
foreach my $s (@Params) {
	print "Key: $s\n";
}
#loop_内のパラメータを取得
			my $LoopContent = $line2;
			while(<IN>) {
				$line2 = $_;
				&DelSpace($line2);
#次のloop_を見つけたら、LoopContentを解析
				last if($line2 =~ /^(_|loop_|#End)/i);

				$LoopContent = "$LoopContent $line2";
			}

			my $pt = 0;
			my $c = 0;
			while($LoopContent =~ /(\'.*?\'|\S+)/g) {
				my $s = $1;
				my $i = int($c / $Count) + 1;
				$c++;
				print "  $i: $Params[$pt]: $s\n";
				$pt++;
				$pt = 0 if($pt >= $Count);
			}
		}
	}
	close(IN);
}

sub MakeLD01File
{
	print "SPG # (ex. 190)=";
	my $spgnum = <STDIN>;
	&DelSpace($spgnum);
	print "SPG Set #=";
	my $spgset = <STDIN>;
	&DelSpace($spgset);

	open(IN, "<$SpaceGroupDBPath") or die "$!: $SpaceGroupDBPath\n";
	while(<IN>) {
		my $line = $_;
		$line =~ /^\s*?(\d+?)\s+?(\d+?)\s+?(\d+?)\s+?(\d+?)\s+?(\d+?)\s+?(\S.*)$/;
		my $NSPGR  = $1;
		my $NSET   = $2;
		my $LAUEG  = $3;
		my $NCENTR = $4;
		my $NSYM   = $5;
		my $SPGR   = $6;
		my $BravaisCell = substr($SPGR, 0, 1);
		my $nTranslate = 1;
		$nTranslate = 2 if($BravaisCell eq 'I');
		$nTranslate = 2 if($BravaisCell eq 'A');
		$nTranslate = 2 if($BravaisCell eq 'B');
		$nTranslate = 2 if($BravaisCell eq 'C');
		$nTranslate = 4 if($BravaisCell eq 'F');
		
		if($spgnum == $NSPGR and $spgset == $NSET) {
			print "Space group #$NSPGR (Set $NSET): $SPGR\n";
			print "    Center Inversion       : $NCENTR\n";
			print "    # of translation       : $nTranslate\n";
			my $TotalSymOperation = $NSYM * (1+$NCENTR) * $nTranslate;
			print "    # of symmetry operation: $TotalSymOperation ($NSYM)\n";
			print "\n";
			
			my $F01Path = "F01.DAT";
			print "Write to $F01Path\n";
			open(F01, ">$F01Path") or die "$!: $F01Path";
			$SPGR =~ s/\s//g;
			print F01 "$SPGR\n";
			print F01 "$NCENTR $nTranslate $NSYM\n";

			print F01 "0.0 0.0 0.0\n";
			if($BravaisCell eq 'I') {
				print F01 "0.5 0.5 0.5\n";
			}
			if($BravaisCell eq 'A') {
				print F01 "0.0 0.5 0.5\n";
			}
			if($BravaisCell eq 'B') {
				print F01 "0.5 0.0 0.5\n";
			}
			if($BravaisCell eq 'C') {
				print F01 "0.5 0.5 0.0\n";
			}
			if($BravaisCell eq 'F') {
				print F01 "0.5 0.5 0.0\n";
				print F01 "0.5 0.0 0.5\n";
				print F01 "0.0 0.5 0.5\n";
			}
			
			<IN>; # skip diffraction condition

#			print "    Symmetry operation:\n";
			for(my $i = 0 ; $i < $NSYM ; $i++) {
				my $SymOp = <IN>;
				&DelSpace($SymOp);
				printf "      %2d: %s\n", $i+1, $SymOp;
				my @SymOpXYZ = split(/,\s*/, $SymOp);
				my @XMatrix = &SymOpToMatrix($SymOpXYZ[0]);
				my @YMatrix = &SymOpToMatrix($SymOpXYZ[1]);
				my @ZMatrix = &SymOpToMatrix($SymOpXYZ[2]);
				printf F01 "%7.4f %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f\n",
					$XMatrix[0], $XMatrix[1], $XMatrix[2],
					$YMatrix[0], $YMatrix[1], $YMatrix[2],
					$ZMatrix[0], $ZMatrix[1], $ZMatrix[2];
				printf F01 "%7.4f %7.4f %7.4f\n",
					$XMatrix[3], $YMatrix[3], $ZMatrix[3];
			}
			print "\n";

			close(OUT);
			last;
		}
#SPGがマッチしないとき、残りの行を読み飛ばす
		my $line2 = <IN>;
		for(my $i = 0 ; $i < $NSYM ; $i++) {
			my $SymOp = <IN>;
		}
	}
	close(IN);

}

sub ShowIonInformation
{
	print "Ion name/atomic number=";
	my $IonName = <STDIN>;
	chomp($IonName);
	my ($AtomicNumber, $IonName, $IonCharge)
		= &GetAtomInformation($IonName);
	my $IonRadius = &GetCationRadius($IonName, 
			$IonCharge, 0);
	print "Name         : $IonName\n";
	print "Atomic number: $AtomicNumber\n";
	print "Charge       : $IonCharge\n" 
		if(0 <= $IonCharge and $IonCharge < 99);
	print "IonRadius1.html\n";
	print "  Ion radius   : $IonRadius A\n"
			if($IonRadius > 0.0);

	my (@rad) = &GetIonRadius2($IonName, $IonCharge, 0);
	print "IoranRadius2.html\n";
	foreach my $s (@rad) {
		if($s =~ /^[\d\.]+?\s/) {
			print "  Ion Radius: $s A\n";
		}
		elsif($s =~ /[\+-]$/) {
			print "  Ion Charge: $s\n";
		}
		else {
			print "  nCoord    : $s\n";
		}
	}
	
	my ($AtmNum, $JName, $EName, $NameOrigin) =
		AtomType::GetIonInfFromPeriodicTable1($IonName);
	print "PeriodicTable1.html\n";
	print "  Atomic Number : $AtmNum\n";
	print "  Japanese Name : $JName\n";
	print "  English Name  : $EName\n";
	print "  Origin of Name: $NameOrigin\n";

	my ($FileName, $AtomJName, $AtomEName,
		$AtomicMass, $FoundYear, $AtomicRadius,
		$MeltingPoint, $BoilingPoint, $Density,
		$HeatCapacity, $IonizationPotential, $ElectronAffinity) =
			&GetIonInfFromPeriodicTable2($AtomicNumber);
	print "PeriodicTable2 ($FileName)\n";
	print "  AtomJName   : $AtomJName\n" if($AtomJName ne '');
	print "  AtomEName   : $AtomEName\n" if($AtomEName ne '');
	print "  AtomicMass  : $AtomicMass\n" if($AtomicMass ne '');
	print "  FoundYear   : $FoundYear\n" if($FoundYear ne '');
	print "  AtomicRadius: $AtomicRadius A\n" if($AtomicRadius ne '');
	print "  MeltingPoint: $MeltingPoint C\n" if($MeltingPoint ne '');
	print "  BoilingPoint: $BoilingPoint C\n" if($BoilingPoint ne '');
	print "  Density     : $Density g/cm3\n" if($Density ne '');
	print "  HeatCapacity: $HeatCapacity cal/g/C\n" if($HeatCapacity ne '');
	print "  IonizationEnergy: $IonizationPotential eV\n" 
				if($IonizationPotential ne '');
	print "  ElectronAffinity: $ElectronAffinity eV\n" 
				if($ElectronAffinity ne '');
	print "\n";
}

sub GetIonRadius2
{
	my ($name, $charge, $nCoord) = (@_);

	$name = uc $name;
	my @radius;
	
	open(IN, "<$IonRadius2DBPath") or die "$!: $IonRadius2DBPath\n";
	while(<IN>) {
		my $line = $_;
		$line =~ m|\<P\>\<B\>(.*?)\</B\>\</P\>|i;
		my $AtomName = $1;
		next unless(uc $AtomName eq $name);

		my $rline = '';
		while(<IN>) {
			my $line2 = $_;
			last if($line2 =~ m|\<TD\s+align=middle.*?\>|i);
			chomp($line2);
			$rline .= $line2;
		}
#print "rline: $rline\n";
#		$rline =~ m|\<P\>(.*?)\</TD\>|i;
#		my @a = split(/(\<\/?B\>|\<\/?FONT.*?\>|\<BR\>|\<\/?P\>|\<I\>)/i, $1);;
		my @a = split(/(\<\/?B\>|\<\/?FONT.*?\>|\<BR\>|\<\/?P\>|\<I\>|\<P.*?\>)/i, $rline);;
#		my @a = split(/\<.*?\>/, $rline);;
		my $IsPauling = 0;
		my $IsReliable = 0;
#print "a: @a\n";
		foreach my $s (@a) {
			$IsPauling = 1 if($s =~ /\<I\>/i);
			$IsReliable = 1 if($s =~ /\<B\>/i);
			$s =~ s/^\s+//;
			$s =~ s/\s+$//;
			next if($s eq '' or $s =~ /[\<\>]/);
			if($s =~ /^[\d\.]+$/) { # to A
				$s = $s * 0.01;
				$s = "$s (less reliable)"
					unless($IsReliable);
				$s = "$s (deduced by Pauling(1960) or Ahrens(1952))"
					if($IsPauling);
			}
			push(@radius, $s);
			$IsPauling = 0;
			$IsReliable = 0;
		}
		last;
	}
	close(IN);
	return @radius;
}


sub GetIonInfFromPeriodicTable2
{
	my ($atnum) = (@_);

	my $AtomJName  = '';
	my $AtomEName  = '';
	my $AtomicMass = '';
	my $FoundYear  = '';
	my $AtomicRadius = '';
	my $MeltingPoint = '';
	my $BoilingPoint = '';
	my $Density      = '';
	my $HeatCapacity = '';
	my $IonizationPotential = '';
	my $ElectronAffinity    = '';
#原子番号 94 
#原子量 239 
#発見年 (年) 1940 
#原子半径 (A) 1.64 
#融点 (oC) 639.5 
#沸点 (oC) 3235 
#密度 (g/cm3) 　 
#比熱 (cal/g oC) 　 
#イオン化エネルギー (eV) 5.8 
#電子親和力 (eV) 

	my $fname = sprintf "A%03dj.html", int($atnum+0.1);
	my $filepath = &MakePath($PeriodicTable2DBDir, $fname,
					$DirectorySeparator, 0);
	unless(open(IN, "<$filepath")) {
		print "$!: $filepath\n";
		return ($filepath);
	}
	while(<IN>) {
		my $line = $_;
#print "line: $line\n";
		if($line =~ m|<title>(.*?)<\/title>|i) {
			$AtomEName = $1;
#print "AtomEName: $AtomEName\n";
		}
		elsif($line =~ m|\<font SIZE=\"\+3\">(.*?)(\(.*?\))?\<\/font\>|i) {
			$AtomJName = $1;
#print "AtomJName: $AtomJName\n";
		}
		elsif($line =~ m|原子番号|i) {
			my $line2 = <IN>;
			$line2 =~ /\>(\d+)\<\/font/;
#			$AtomicMass = $1;
		}
		elsif($line =~ m|原子量|i) {
			my $line2 = <IN>;
			$line2 =~ /\>(\d*)\<\/font/;
			$AtomicMass = $1;
#print "AtomicMass: $AtomicMass\n";
		}
		elsif($line =~ m|発見年.*\>(\d*?)\<\/font|i) {
			$FoundYear = $1;
#print "FoundYear: $FoundYear\n";
		}
		elsif($line =~ m|原子半径.*\>([\d\.]*)\<\/font|i) {
			$AtomicRadius = $1;
		}
		elsif($line =~ m|融点.*\>([-\+\d\.]*)\<\/font|i) {
			$MeltingPoint = $1;
		}
		elsif($line =~ m|沸点.*\>([-\+\d\.]*)\<\/font|i) {
			$BoilingPoint = $1;
		}
		elsif($line =~ m|密度.*\>([\d\.]*)\<\/font|i) {
			$Density = $1;
		}
		elsif($line =~ m|比熱.*\>([\d\.]*)\<\/font|i) {
			$HeatCapacity = $1;
		}
		elsif($line =~ m|イオン化.*\>([-\+\d\.]*)\<\/font|i) {
			$IonizationPotential = $1;
		}
		elsif($line =~ m|電子.*\>([-\+\d\.]*)\<\/font|i) {
			$ElectronAffinity = $1;
		}
	}
	close(IN);
	return ($filepath, $AtomJName, $AtomEName,
		$AtomicMass, $FoundYear, $AtomicRadius,
		$MeltingPoint, $BoilingPoint, $Density,
		$HeatCapacity, $IonizationPotential, $ElectronAffinity);
}

sub GetIonInfFromPeriodicTable1
{
	my ($name) = (@_);

	$name = uc $name;
#print "name: $name\n";
	my $AtmNum = -1;
	my $JName = "";
	my $EName = "";
	my $NameOrigin = "";
#    <TD>番号</TD>
#    <TD>記号</TD>
#    <TD>元素名</TD>
#    <TD>英語名</TD>
#    <TD>名前の由来</TD></TR>

	open(IN, "<$PeriodicTable1DBPath") or die "$!: $PeriodicTable1DBPath\n";

	while(<IN>) {
		my $line = $_;
#print "line: $line\n";
		next if($line !~ m|<TD><A name=.*?>(\d+)<\/A><\/TD>|i);
		$AtmNum = $1;
#print "AtmNum: $AtmNum\n";
		my $line    = <IN>;
		$line =~ m|<TD.*?>(.*?)<\/TD>|i;
		my $atname   = $1;
#print "atname: [$atname] ($name)\n";
		next if($name ne uc $atname and $AtmNum != int($name+0.1));

		$line    = <IN>;
		$line =~ m|<TD.*?>(.*?)<\/TD>|i;
		$JName = $1;
#print "JName: $JName\n";
		$line    = <IN>;
		$line =~ m|<TD.*?>(.*?)<\/TD>|i;
		$EName = $1;
		$line    = <IN>;
		$line =~ m|<TD.*?>(.*?)<\/TD>|i;
		$NameOrigin = $1;
		last;
	}
	close(IN);

	return ($AtmNum, $JName, $EName, $NameOrigin);
}

sub ShowSpaceGroup
{
	my @LaueGroup = (
		"Triclinic (1)",	# 1
		"Monoclinic, a (2/m)",	# 2
		"Monoclinic, b (2/m)",	# 3
		"Monoclinic, c (2/m)",	# 4
		"Orthorhombic (mmm)",	# 5
		"Tetragonal (4/m)",	# 6
		"Tetragonal (4/mmm)",	# 7
		"Trigonal, Rhombohedral (3)",	# 8
		"Trigonal, Hexagonal (3)",	# 9
		"Trigonal, Rhombohedral (3m)",	# 10
		"Trigonal, Hexagonal (3m)",	# 11
		"Hexagonal (6/m)",	# 12
		"Hexagonal (6/mmm)",	# 13
		"Cubic (m3)",		# 14
		"Cubic (m3m)"		# 15
	);
	my @CentroSymmetry = (
		"No CentroSymmetry",
		"CentroSymmetric"
	);
	my @ICONDTypeStr = (
		"00L", "0K0", "0KL", "H00", "H0L",
		"HK0", "HKL", "0KK", "HH0", "HHL",
		"H0H", "HKK", "HKH", "HHH"
	);
	my @ICONDCondStr = (
		"None",
		"H=2N",          "K=2N",           "L=2N",     "K+L=2N",    "H+L=2N",
		"H+K=2N",        "H,K,L ODD/EVEN", "K+L=4N",   "H+L=4N",    "H+K=4N",
		"2H+L=2N",       "2H+L=4N",        "H+K+L=2N", "-H+K+L=3N", "H-K+L=3N",
		"H=4N",          "K=4N",           "L=3N",     "L=4N",      "L=6N",
		"|H|>=|K|>=|L|", "2H+K=2N",        "2H+K=4N",  "H+2K=2N",   "H+2K=4N",
		"H=2N,K=2N",     "K=2N,L=2N",      "H=2N,L=2N"
	);

	print "SPG # (ex. 190) or Name (ex. P 4 n m)=";
	my $spgname = <STDIN>;
	&DelSpace($spgname);
	my @SPGElement = split(/\s+/, $spgname);

	open(IN, "<$SpaceGroupDBPath") or die "$!: $SpaceGroupDBPath\n";
	while(<IN>) {
		my $line = $_;
		$line =~ /^\s*?(\d+?)\s+?(\d+?)\s+?(\d+?)\s+?(\d+?)\s+?(\d+?)\s+?(\S.*)$/;
		my $NSPGR  = $1;
		my $NSET   = $2;
		my $LAUEG  = $3;
		my $NCENTR = $4;
		my $NSYM   = $5;
		my $SPGR   = $6;
		my $BravaisCell = substr($SPGR, 0, 1);
		my $nTranslate = 1;
		$nTranslate = 2 if($BravaisCell eq 'I');
		$nTranslate = 2 if($BravaisCell eq 'A');
		$nTranslate = 2 if($BravaisCell eq 'B');
		$nTranslate = 2 if($BravaisCell eq 'C');
		$nTranslate = 4 if($BravaisCell eq 'F');

		my $Identical = 1;
		if($spgname =~ /^\d+$/) {
			$Identical = 0 unless($spgname == $NSPGR);
		}
		elsif(@SPGElement > 0) {
			foreach my $s (@SPGElement) {
				$s = &RegExpQuote($s);
				unless($SPGR =~ /$s/) {
					$Identical = 0;
					last;
				}
			}
		}
		
		if($Identical) {
			print "Space group #$NSPGR (Set $NSET): $SPGR\n";
			$LAUEG--;
			print "    Laue Group             : @LaueGroup[$LAUEG]\n";
			print "    Bravais Cell           : $BravaisCell\n";
			print "    Center Inversion       : @CentroSymmetry[$NCENTR]\n";
			print "    # of translation       : $nTranslate\n";
			my $TotalSymOperation = $NSYM * (1+$NCENTR) * $nTranslate;
			print "    # of symmetry operation: $TotalSymOperation ($NSYM)\n";
			
			my $ICONDLine = <IN>;
			my @ICOND = unpack("a2a2a2a2a2a2a2a2a2a2a2a2a2a2", $ICONDLine);
			print "    Diffraction condition:\n";
			for(my $i = 0 ; $i < 14 ; $i++) {
				printf "      %3s: %s\n",
					$ICONDTypeStr[$i], $ICONDCondStr[$ICOND[$i]]
						if($ICOND[$i] > 0);
			}

			print "    Symmetry operation:\n";
			for(my $i = 0 ; $i < $NSYM ; $i++) {
				my $SymOp = <IN>;
				&DelSpace($SymOp);
				printf "      %2d: %s\n", $i+1, $SymOp;
				my @SymOpXYZ = split(/,\s*/, $SymOp);
				my @XMatrix = &SymOpToMatrix($SymOpXYZ[0]);
				my @YMatrix = &SymOpToMatrix($SymOpXYZ[1]);
				my @ZMatrix = &SymOpToMatrix($SymOpXYZ[2]);
				printf "         |x'|=|%6.3f %6.3f %6.3f||x| |%6.3f|\n",
					$XMatrix[0], $XMatrix[1], $XMatrix[2], $XMatrix[3];
				printf "         |y'|=|%6.3f %6.3f %6.3f||y|+|%6.3f|\n",
					$YMatrix[0], $YMatrix[1], $YMatrix[2], $YMatrix[3];
				printf "         |z'|=|%6.3f %6.3f %6.3f||z| |%6.3f|\n",
					$ZMatrix[0], $ZMatrix[1], $ZMatrix[2], $ZMatrix[3];
			}
			print "\n";
			next;
		}
		my $line2 = <IN>;
		for(my $i = 0 ; $i < $NSYM ; $i++) {
			my $SymOp = <IN>;
		}
	}
	close(IN);
}

sub SymOpToMatrix
{
	my ($symopstr) = (@_);

if($Debug) {
	print "so: $symopstr\n";
}
	my ($x, $y, $z, $t);
	if($symopstr =~ /([-\+\/\.\d]*)x/) {
		$x = $1;
		$x = 1 if($1 eq '' or $1 eq '+');
		$x = -1 if($1 eq '-');
	}
	if($symopstr =~ /([-\+\/\.\d]*)y/) {
		$y = $1;
		$y = 1 if($1 eq '' or $1 eq '+');
		$y = -1 if($1 eq '-');
	}
	if($symopstr =~ /([-\+\/\.\d]*)z/) {
		$z = $1;
		$z = 1 if($1 eq '' or $1 eq '+');
		$z = -1 if($1 eq '-');
	}
	if($symopstr =~ /([^xyz]*)$/) {
		$t = $1;
	}
	$x = eval($x);
	$y = eval($y);
	$z = eval($z);
	$t = eval($t);
	return ($x, $y, $z, $t);
}

sub ExecuteNormalizeF08e
{
	my $F08ePath = &MakePath($BaseDir, 'f08e',
				$DirectorySeparator, 0);
	my $F08eNormlizedPath = 
		&MakePath($BaseDir, "f08e.Normalized",
				$DirectorySeparator, 0);
	my $command = $EditorPath;

	$command =~ s/\%1/$F08ePath/g;
	print "Execute $command...\n";
	system($command);

	print "O 1s min energy[$MinO1sEnergy eV]=";
	my $tmp = <STDIN>;
	&DelSpace($tmp);
	$MinO1sEnergy = $tmp if($tmp ne '');
	print "O 1s max energy[$MaxO1sEnergy eV]=";
	$tmp = <STDIN>;
	&DelSpace($tmp);
	$MaxO1sEnergy = $tmp if($tmp ne '');
	print "O 1s energy[$O1sEnergy eV]=";
	my $tmp = <STDIN>;
	&DelSpace($tmp);
	$O1sEnergy = $tmp if($tmp ne '');

	print "MinO1sEnergy: $MinO1sEnergy eV\n";
	print "MaxO1sEnergy: $MaxO1sEnergy eV\n";
	print "O1sEnergy   : $O1sEnergy eV\n";
	print "\n";

	open(IN, "<$F08ePath") or die "$!: $F08ePath\n";
if($Debug) {
	print "Search key line.\n";
}
	while(<IN>) {
if($Debug) {
	print "  line: $_\n";
}
		last if($_ =~ /\(RY\)\s+\(HR\)\s+\(EV\)/);
	}
if($Debug) {
	print "Key line found.\n";
}
	<IN>; # blank line
	my $FilePosition = tell(IN);
	my $nO1s = 0;
	my $AverageO1sEnergy = 0;
if($Debug) {
	print "Search O 1s.\n";
}
	while(<IN>) {
		my $line = $_;
		&DelSpace($line);
		last if($line =~ /^$/);

		my ($sn, $nLevel, $Spin, $RY, $HR, $EV, $nElectron);
		if($line =~ /(UP|DOWN)/) {
			($sn, $nLevel, $Spin, $RY, $HR, $EV, $nElectron)
				= split(/\s+/, $line);
		}
		else {
			($sn, $nLevel, $RY, $HR, $EV, $nElectron)
				= split(/\s+/, $line);
			$Spin = '    ';
		}
if($Debug) {
	print "line: $line\n";
	print "  $sn: $EV\n";
}
		if($MinO1sEnergy <= $EV and $EV <= $MaxO1sEnergy) {
			$nO1s++;
			$AverageO1sEnergy += $EV;
		}
	}
	if($nO1s <= 0) {
		print "\n";
		print "Error: O1s level was not found.\n";
		return;
	}
	$AverageO1sEnergy /= $nO1s;
	my $EnergyShift = $O1sEnergy - $AverageO1sEnergy;
	print "Average O1s Energy: $AverageO1sEnergy eV\n";
	print "Energy shift      : $EnergyShift\n";
	print "\n";

	seek(IN, $FilePosition, 0);

	open(OUT, ">$F08eNormlizedPath") or die "$!: $F08eNormlizedPath\n";
	print OUT "*** M.O. Eigenvalue (eV)\n";
	print OUT "  Average O1s Energy: $AverageO1sEnergy eV\n";
	print OUT "  Energy shift      : $EnergyShift\n";
	print OUT "\n";
	while(<IN>) {
		my $line = $_;
		&DelSpace($line);
		last if($line =~ /^$/);
		
		my ($sn, $nLevel, $Spin, $RY, $HR, $EV, $nElectron);
		if($line =~ /(UP|DOWN)/) {
			($sn, $nLevel, $Spin, $RY, $HR, $EV, $nElectron)
				= split(/\s+/, $line);
		}
		else {
			($sn, $nLevel, $RY, $HR, $EV, $nElectron)
				= split(/\s+/, $line);
			$Spin = '    ';
		}

		printf OUT " %4d %4d %-4s %12.5f %12.5f\n",
			$sn, $nLevel, $Spin, $EV+$EnergyShift, $nElectron;
	}
	close(OUT);
	close(IN);
	
	print "Write to $F08eNormlizedPath\n";
	print "\n";
}

sub ExecuteSCFCalcDVXa
{
	my $AnionName = "O";
	my $AnionAtomicNumber = &GetAtomicNumber($AnionName);
	my $AnionCharge = -2;
	my $AnionRadius = 1.60; # A

	my $CationName = "Al";
	my $CationAtomicNumber = 13;
	my $CationCharge = +3;
	my $CationRadius = 0.675;
	my $CoordinationNumber = 6;

	print "Cation name/atomic number=";
	$CationName = <STDIN> if($InputFromSTDIN);
	chomp($CationName);
	($CationAtomicNumber, $CationName, $CationCharge)
		= &GetAtomInformation($CationName);

	print "Cation name [$CationName]=";
	my $tmp = <STDIN>;
	&DelSpace($tmp);
	$CationName = $tmp if($tmp ne '');

	print "Cation atomic number [$CationAtomicNumber]=";
	my $tmp = <STDIN>;
	&DelSpace($tmp);
	$CationAtomicNumber = $tmp if($tmp ne '');

	print "Cation charge [$CationCharge]=";
	$tmp = <STDIN>;
	&DelSpace($tmp);
	$CationCharge = $tmp if($tmp ne '');

	print "   Cation name: $CationName\n";
	print "   Cation atomic number: $CationAtomicNumber\n";
	print "   Cation charge: $CationCharge\n";

	unless($CalcAtomOnly) {
		print "Coordination number [$CoordinationNumber]=";
		$tmp = <STDIN>;
		chomp($tmp);
		$CoordinationNumber = $tmp if($tmp ne '');

		$CationRadius = &GetCationRadius($CationName, 
					$CationCharge, $CoordinationNumber);
		print "Cation radius [$CationRadius A]=";
		$tmp = <STDIN>;
		chomp($tmp);
		$CationRadius = $tmp if($tmp ne '');

		print "   CoordinationNumber: $CoordinationNumber\n";
		print "   Cation radius: $CationRadius A\n";
	}

#=========================================
# ディレクトリィの作製
#=========================================
	my $TotalCharge = $AnionCharge * $CoordinationNumber + $CationCharge;
	if($TotalCharge == 0) {
		$TotalCharge = '';
	}
	elsif($TotalCharge > 0) {
		$TotalCharge = "$TotalCharge+";
	}
	else {
		$TotalCharge = -$TotalCharge;
		$TotalCharge = "$TotalCharge-";
	}
	my $CompoundName = "($CationName$AnionName$CoordinationNumber)$TotalCharge";
	$CompoundName = "$CationName$TotalCharge" if($CoordinationNumber == 0);
	$CompoundName = "($CationName$AnionName)$TotalCharge" if($CoordinationNumber == 1);
	if($CalcAtomOnly) {
		$CompoundName = "$CationName$CationCharge";
	}
	$BaseDir = &MakePath($InitialDirectory, $CompoundName, 
			$DirectorySeparator, 0);

	print "\n";
	print "Create $BaseDir...\n";
	&CreateDirecotry($BaseDir);

#=========================================
# nonrelをコピー
#=========================================
	print "\n";
	print "Copy nonrel...\n";
	my $NonrelOutputPath = &MakePath($BaseDir, "nonrel", 
				$DirectorySeparator, 0);
	print "  $NonrelInputPath to $NonrelOutputPath\n";
	&CopyFile($NonrelInputPath, $NonrelOutputPath);

#=========================================
# radiiをコピー
#=========================================
	print "\n";
	print "Copy radii...\n";
	my $RadiiOutputPath = &MakePath($BaseDir, "radii", 
				$DirectorySeparator, 0);
	print "  $RadiiInputPath to $RadiiOutputPath\n";
	&CopyFile($RadiiInputPath, $RadiiOutputPath);

	print "\n";
	print "Chage Directory to $BaseDir\n";
	chdir($BaseDir);

#=========================================
# f01を作製
#=========================================
	my $F01Path = "f01";
	print "Make f01: $F01Path\n";
	open(F01, ">$F01Path");
	printf F01 "| Z ||NEQ||   X    ||   Y    ||   Z    |\n";
#	printf F01 "    8    3  -1.43895  -1.43895  -1.43895";
	my $Distance = $AnionRadius + $CationRadius;
	my ($x, $y, $z);
	$CalcAtomOnly = 1 if($CoordinationNumber == 0);
	if($CalcAtomOnly) {
		$x = 0.0; $y = 0.0; $z = 0.0;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $CationAtomicNumber, 1, $x, $y, $z;
	}
	elsif($CoordinationNumber == 1) {
		$x =  0.0; $y = 0.0; $z = -$Distance/2.0;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $CationAtomicNumber, 1, $x, $y, $z;
		$x =  0.0; $y = 0.0; $z =  $Distance/2.0;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;
	}
	elsif($CoordinationNumber == 2) {
		$x =  0.0; $y = 0.0; $z = 0.0;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $CationAtomicNumber, 1, $x, $y, $z;
		$x =  0.0; $y = 0.0; $z =  $Distance;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;
		$x =  0.0; $y = 0.0; $z = -$Distance;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;
	}
	elsif($CoordinationNumber == 3) {
		$x =  0.0; $y = 0.0; $z = 0.0;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $CationAtomicNumber, 1, $x, $y, $z;
		$x =  $Distance; $y = 0.0; $z = 0.0;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;
		my $kx = $Distance/2.0;
		my $ky = $Distance*sqrt(3)/2.0;
		$x = -$kx; $y =  $ky; $z =  0.0;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;
		$x = -$kx; $y = -$ky; $z =  0.0;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;
	}
	elsif($CoordinationNumber == 4) {
# M : ( 0, 0, 0)
# O1: ( 0, 0, R)           R: M-O間距離
# O2: (2x, 0, z)           幾何学よりz=-R/3
# O3: (-x,  sqrt(3)*x, z)  幾何学からx,y座標は決まる
# O4: (-x, -sqrt(3)*x, z)  幾何学からx,y座標は決まる
		$x =  0.0; $y = 0.0; $z = 0.0;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $CationAtomicNumber, 1, $x, $y, $z;
		$x = 0.0; $y = 0.0; $z = $Distance;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;
		my $kz = -$Distance / 3.0;
		my $kx = sqrt(2.0) / 3.0 * $Distance;
		$x = 2.0*$kx; $y = 0.0; $z = $kz;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;
		$x = -$kx; $y =  sqrt(3.0)*$kx; $z = $kz;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;
		$x = -$kx; $y = -sqrt(3.0)*$kx; $z = $kz;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;
	}
	elsif($CoordinationNumber == 6) {
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $CationAtomicNumber, 1, 0.0, 0.0, 0.0;
		$x =  $Distance; $y = 0.0; $z = 0.0;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;
		$x = -$Distance; $y = 0.0; $z = 0.0;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;
		$x = 0.0; $y =  $Distance; $z = 0.0;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;
		$x = 0.0; $y = -$Distance; $z = 0.0;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;
		$x = 0.0; $y = 0.0; $z =  $Distance;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;
		$x = 0.0; $y = 0.0; $z = -$Distance;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;
	}
	elsif($CoordinationNumber == 8) {
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $CationAtomicNumber, 1, 0.0, 0.0, 0.0;
		my $kx = $Distance/sqrt(3);
		$x =  $kx; $y =  $kx; $z =  $kx;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;
		$x =  $kx; $y =  $kx; $z = -$kx;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;
		$x =  $kx; $y = -$kx; $z =  $kx;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;
		$x =  $kx; $y = -$kx; $z = -$kx;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;
		$x = -$kx; $y =  $kx; $z =  $kx;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;
		$x = -$kx; $y =  $kx; $z = -$kx;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;
		$x = -$kx; $y = -$kx; $z =  $kx;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;
		$x = -$kx; $y = -$kx; $z = -$kx;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;
	}
	elsif($CoordinationNumber == 12) {
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $CationAtomicNumber, 1, 0.0, 0.0, 0.0;
		my $kx = $Distance/sqrt(2);
		$x =  $kx; $y =  $kx; $z =  0.0;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;
		$x =  $kx; $y = -$kx; $z =  0.0;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;
		$x = -$kx; $y =  $kx; $z =  0.0;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;
		$x = -$kx; $y = -$kx; $z =  0.0;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;

		$x =  $kx; $y = 0.0; $z =  $kx;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;
		$x =  $kx; $y = 0.0; $z = -$kx;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;
		$x = -$kx; $y = 0.0; $z =  $kx;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;
		$x = -$kx; $y = 0.0; $z = -$kx;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;

		$x = 0.0; $y =  $kx; $z =  $kx;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;
		$x = 0.0; $y =  $kx; $z = -$kx;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;
		$x = 0.0; $y = -$kx; $z =  $kx;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;
		$x = 0.0; $y = -$kx; $z = -$kx;
		printf F01 " %4d %4d %9.5f %9.5f %9.5f\n",
			 $AnionAtomicNumber, 2, $x, $y, $z;
	}
	else {
		print "Error: CoordinationNumber ($CoordinationNumber) is"
			." not supported.\n";
		return;
	}
	printf F01 "---------------------------------------------\n";
	printf F01 "|NEQ||  CHG   ||U/D||   RD   ||   VD   |    1 	\n";
#	printf F01 "    1   2.00000\n";
	if($CalcAtomOnly) {
		printf F01 " %4d %9.5f\n", 1, $CationCharge;
	}
	else {
		printf F01 " %4d %9.5f\n", 1, $CationCharge;
		printf F01 " %4d %9.5f\n", 2, $AnionCharge;
	}
	printf F01 "---------------------------------------------\n";
	printf F01 "    0     Unit     (0:angstrom  1:atomic)\n";
	printf F01 "    %d     Spin     (0:non-spin  1:spin  )\n",
			$SpinPolarizedCalculation;
	printf F01 "    0     M.P.     (0:No        1:Yes   )\n";
	printf F01 "    0     Sample Point ( <100000, =0 autoset )\n";

	close(F01);

#=========================================
# makef05.exeを実行してf05を作製
#=========================================
	print "\n";
	print "Make f05: Executing $MakeF05Path ...\n";
	system($MakeF05Path);

	print "Copy F27...\n";
	print "Chage Directory to $InitialDirectory\n";
	my $F27OutputPath = &MakePath($BaseDir, "f27", $DirectorySeparator, 0);
	print "  $F27InputPath to $F27OutputPath\n";
	&CopyFile($F27InputPath, $F27OutputPath);

#=========================================
# 収束計算のループ
#=========================================
	my $Count = 0;
	while(1) {
		$Count++;
#=========================================
# spin.exeの実行
#=========================================
		if($Count < 2) {
			print "\n";
			print "Execute $DVXaSpinPath...\n";
			print "Chage Directory to $BaseDir\n";
			chdir($BaseDir);
		}
		system($DVXaSpinPath);

#=========================================
# f06zから、収束データを取得
#=========================================
		my $F06zPath = "f06z";
		if($Count < 2) {
			print "\n";
			print "Read from $F06zPath\n";
		}
		open(F06zPath, "<$F06zPath");
		my $content = '';
		while(<F06zPath>) {
			$content = $_;
		}
		my ($n1, $a, $b, $c, $d, $convergence, $n2) =
			($content =~ /(.+),(.+),(.+),(.+),(.+),(.+),(.+)/g);
		printf "   Convergence[$Count]: $n1($n2) $convergence\n";

#=========================================
# ファイルの更新
#=========================================
		CopyFile("f27", "f27.$Count.back");
		CopyFile("f17", "f27");
		CopyFile("f05", "f05.$Count.back");
		CopyFile("f26", "f05");

		if($convergence < $ConvergenceThreshold) {
			print "Calculation converged.\n";
			print "   ($convergence < $ConvergenceThreshold)\n";
			last;
		}
		if($Count >= $MaxLoop) {
			print "\n";
			print "Calculation did not converge.\n";
			print "   ($convergence > $ConvergenceThreshold)\n";
			print "Exit Loop because \$Count reaches "
				."maximum number ($MaxLoop).\n";
			last;
		}
	}
}

#====================================================================
#====================================================================
#====================================================================
#                  Subroutines
#====================================================================
#====================================================================
#====================================================================

sub GetCationRadius
{
	my ($name, $charge, $nCoordination) = (@_);

if($Debug) {
	print "Enter GetCationRadius\n";
}
	$name = uc $name;
	$charge = '1+' if($charge eq '+');
	$charge = '1-' if($charge eq '-');
	$charge = 0.0 + $charge if($charge =~ /\+/);
	$charge = 0.0 - $charge if($charge =~ /\-/);
	$charge = '' if($charge == 0.0);
#print "c: $charge\n";
	my $r = -1.0;
	open(IN, "<$IonRadius1DBPath") or die "$!: $IonRadius1DBPath\n";
	while(<IN>) {
		my $line = $_;
		next if($line !~ m|<TR>|i);
		my $atmline    = <IN>;
		my $atmnumline = <IN>;
		my $radiusline = <IN>;
		$atmline =~ m|<TD.*?>(.*?)<SUP>(.*?)<\/SUP><\/TD>|i;
		my $atname   = uc $1;
		my $atcharge = $2;
		$atcharge = 0.0 + $atcharge if($atcharge =~ /\+/);
		$atcharge = 0.0 - $atcharge if($atcharge =~ /\-/);
		$atmnumline =~ m|<TD.*?>(.*?)<\/TD>|i;
		my $atomicnumber   = $1;
		$radiusline =~ m|<TD.*?>(.*?)<\/TD>|i;
		my $radius   = $1;
if($Debug) {
	print "atmline: $atmline";
	print "   $atname ($atcharge) $radius\n";
}
		return $radius if($name eq $atname and $charge eq $atcharge);
	}
	close(IN);
	return $r;
}

sub GetAtomInformation
{
	my ($name) = (@_);

	my ($AtomName, $AtomicNumber, $Charge);
	chomp($name);
	$name =~ s/^\s*//;
	$name = uc $name;

#print "##GetAtomInformation: $name\n";

	if($name =~ /^\d+$/) {
		$AtomicNumber = $name;
		$AtomName = &GetAtomName($AtomicNumber);
		$Charge = 999;
	}
	else {
		my @n = ($name =~ /^(.*?)\s*([\d\+-]+)?$/);
#print "n: $n  : $1  : $2\n";
		if(@n == 0) {
			$AtomName = '';
			$Charge = 999;
		}
		elsif(@n == 1) {
			$AtomName = $1;
			$Charge = 999;
		}
		else {
			$AtomName = $1;
			$Charge = $2;
		}
		$AtomicNumber = &GetAtomicNumber($AtomName);
	}
	$AtomName = ucfirst lc $AtomName;
	return ($AtomicNumber, $AtomName, $Charge);
}

sub GetAtomName
{
	my ($AtomNumber) = (@_);

#print "##GetAtomName: $AtomNumber\n";

	my $name = '';
	open(IN, "<$NonrelInputPath") or die "$!: $NonrelInputPath\n";
# H  HFS ATOM CALC.
#    1  300    1
#   1.00000   0.97804  20.00000  32.00000   0.00000
#  1.0  0.0  0.0   -0.30       1.0
	while(<IN>) {
		my $line = $_;       # H  HFS ATOM CALC.
		my $nextline = <IN>; #    1  300    1
		$nextline = <IN>;    #   1.00000   0.97804  20.00000  32.00000   0.00000
#print "line    : $line";
#print "nextline: $nextline";
		my ($n) = ($nextline =~ /^\s*([\d\.]+?)\s/);
#print "  n: $n  an: $AtomNumber\n";
		if($n == $AtomNumber) {
			($name) = ($line =~ /^\s*(\w+?)\s/);
			last;
		}
# 次の空行までスキップ
		while(<IN>) {
			my $line2 = $_;
			chomp($line2);
			$line2 =~ s/^\s*//;
			if($line2 eq '') {
				last;
			}
		}
	}
	close(IN);
	$name = ucfirst lc $name;
	return $name;
}

sub GetAtomicNumber 
{
	my ($name) = (@_);

	$name = uc $name;
	my $an = -1;
	open(IN, "<$NonrelInputPath") or die "$!: $NonrelInputPath\n";
	while(<IN>) {
# H  HFS ATOM CALC.
#    1  300    1
#   1.00000   0.97804  20.00000  32.00000   0.00000
#  1.0  0.0  0.0   -0.30       1.0
  		if($_ =~ /^\s*$name\s/) { # H  HFS ATOM CALC.
			my $line = <IN>;  #    1  300    1
			$line = <IN>;     #   1.00000   0.97804  20.00000
			($an) = ($line =~ /^\s*(.*?)\s/);
			last;
		}
	}
	close(IN);
	return $an;
}

sub CopyFile
{
	my ($infile, $outfile) = (@_);
	
	open(IN, "<$infile") or die "$!: $infile\n";
	open(OUT, ">$outfile") or die "$!: $outfile\n";
	my @content = <IN>;
	print OUT @content;
	close(OUT);
	close(IN);
}

sub CreateDirecotry
{
	my ($dir) = @_;
	$dir =~ s/\\/\//g;

if($Debug) {
	print "CreateDirectory: $dir\n";
}
#	my $IsHeadSep = 0;
#	$IsHeadSep = 1 if($dir =~ /^\//);

	my $IsRelative = 1;
	$IsRelative = 0 if($dir =~ /^[A-Za-z]:\// or $dir =~ /^\//);
#print "IsRelative: $IsRelative\n";
	my @eachpath = split(/\//, $dir);
	
	my $i = 0;
	my $path = "";
	if($eachpath[$i] =~ /^[A-Za-z]:/) {
		$path = $eachpath[$i];
		$i++;
	}
	$path .= $DirectorySeparator unless($IsRelative);

	for( ; $i < @eachpath ; $i++) {
		$path .= $eachpath[$i] . $DirectorySeparator;
if($Debug) {
	print "path: $path\n";
}
		next if(-d $path);
		if(-f $path) {
			print "Can not create [$path]: $path is a file.\n";
			return;
		}
if($Debug) {
	print "    Create $path [$DirectorySeparator]\n";
}
		mkdir($path);
	}
}

sub DelSpace
{
	$_[0] =~ s/^\s*//;
	$_[0] =~ s/\s*$//;
	return $_[0];
}

sub MakePath {
	my ($dir, $fname, $separator, $DoTerminate) = (@_);

	my $RegSep = &RegExpQuote($separator);

#	my $IsHeadSep = 0;
#	$IsHeadSep = 1 if($dir =~ /^$RegSep/);
	$dir   =~ s/$RegSep$//;
	$fname =~ s/^$RegSep//;
	$fname =~ s/$RegSep$//;

	if($dir eq '') {
		$dir = $fname;
	}
	else {
		$dir = "$dir$separator$fname";
	}
	
	if($DoTerminate) {
		unless($dir =~ /$RegSep$/) {
			$dir .= $separator;
		}
	}
	return $dir;
}

sub RegExpQuote
{
	my($str) = (@_);

	$str =~ s/\\/\\\\/g;
	$str =~ s/\//\\\//g;
	return $str;
}