#!/usr/bin/perl
##!d:/Perl/bin/perl

BEGIN {
my $BaseDir = $ENV{'TkPerlDir'};
#print "\nBaseDir: $BaseDir\n";
@INC = ("$BaseDir/lib", "$BaseDir/VNL", "c:/Programs/Perl/lib", "d:/Programs/Perl/lib", @INC);
}

use strict;
#use warnings;
#use CGI::Carp qw(fatalsToBrowser);
use Cwd;

use ProgVars;
use Deps;
use Utils;

use MyApplication;
use Crystal::CIF;
use Crystal::Crystal;
use Crystal::PWSCF;
use Sci qw($torad $eVToJ $RyToeV);
use Sci::Optimize;

BEGIN {
#	$| = 1;
	open (STDERR, ">&STDOUT");
}

#===============================================
# Definition
#===============================================

# eiJ = e_i^j = d a_i / d a_j
# U = 1/2 * Cijkl * eIJ * eKL = 1/2 * C_{ijkl} * e^{ij} * e^{IJ}
#   = 1/2 * CIjKl * eiJ * ekL
# CIjKl = d2U / deiJ / dekL
# U = U0 + a11*dx2 + a22*dy2 + 2a12*dx*dy
# (dx, dy) = ( h,  0) => Up0 =  a11
# (dx, dy) = (-h,  0) => Um0 =  a11
# (dx, dy) = ( 0,  h) => U0p =  a22
# (dx, dy) = ( 0, -h) => U0m =  a22
# (dx, dy) = ( h,  h) => Upp = 2a12
# (dx, dy) = (-h, -h) => Umm = 2a12

#===============================================
# 大域変数
#===============================================
my $OS = Utils::GetOSName();

# eiJを計算するための格子定数の歪率
my $OriginalCIF         = 'original.cif';
my $StrainedCIFTemplate = 'e{i}{j}-{sign}.cif';
my $TableFile           = "table.txt";
my $SummaryCSV          = "summary.csv";

my $ProgramDir     = ProgVars::ProgramDir();
my $TkPerlDir      = ProgVars::TkPerlDir();
my $ElasticDir     = Deps::MakePath($TkPerlDir,  ["ElasticConstant"], 0);
my $TemplateSource = Deps::MakePath($ElasticDir, ["Template"], 0);

my $BaseDir = cwd();
my $TemplateDir = Deps::MakePath($BaseDir, "Template", 0);

#===============================================
# Applicationクラス作製
#===============================================
my $App = new MyApplication;
exit if($App->Initialize() < 0);

# 実行プログラム名（デフォルトでは$0から取得）
my $Program = $App->Program();
# アプリケーション名
$App->SetAppName($Program);
# バージョン
my $Version = $App->SetVersion("Ver 0.1");

#===============================================
# コマンドラインオプション読み込み
#===============================================
$App->AddArgument("--Action",  "--Action=[MakeInputs|MakeRunScript|run|cal]: Default ''", '');
$App->AddArgument("--Mode",    "--Mode=[Cabc|CXYZ|Vconst|Vconp]: Default Vconst", 'Vconst');
$App->AddArgument("--Index",   "--Index=[ijkl]: Default 1111",              '1111');
$App->AddArgument("--nStep",   "--nStep=[integer >= 2]: Default 3",         3);
$App->AddArgument("--nOrder",  "--nOrder=[integer <= nStep-1]: Default 2",  2);
$App->AddArgument("--pressure0",    "--pressure0=val in atm. Origin of pressure: Default 1.0",           1.0);
$App->AddArgument("--pressure_tol", "--pressure_tol=val in atm. Torelance for convergence: Default 1.0", 1.0);
$App->AddArgument("--fdiff",   "--fdiff=da/a=dV/V=val: Default 0.01",     0.01);
$App->AddArgument("--fdiff",   "--fdiff=da/a=dV/V=val: Default 0.01",     0.01);
$App->AddArgument("--Engine",  "--Engine=[vasp|gulp]: Default ''",          '');
$App->AddArgument("--Library", "--Library=[lib path]: Default ''",          '');
$App->AddArgument("--help",    "--help    : Show this help");
#print "pIn=$App->{pIn}<br>\n";
exit 1 if($App->ReadArgs() != 1);
#print "pIn=$App->{pIn}<br>\n";

#===============================================
# 変数設定
#===============================================
my $Action       = $App->GetGetArg('Action', undef, 1, 1);
my $Mode         = $App->GetGetArg('Mode',   undef, 1, 1);
my $Index        = $App->GetGetArg('Index',  undef, 1, 1);
my $nStep        = $App->GetGetArg('nStep',  undef, 1, 1);
my $nOrder       = $App->GetGetArg('nOrder', undef, 1, 1);
my $pressure0    = $App->GetGetArg('pressure0',  undef,  1, 1);
my $pressure_tol = $App->GetGetArg('pressure_tol',  undef,  1, 1);
my $fdiff        = $App->GetGetArg('fdiff',  undef,  1, 1);
my $Engine       = $App->GetGetArg('Engine', undef, 1, 1);

my $CIFFile = $App->GetGetArg(0);
   $CIFFile = 'ZnO.cif' if($CIFFile eq '');
my $SampleName = $CIFFile;
$SampleName =~ s/\..*?$//;   

$App->print("\n");
$App->print("OS    : $OS\n");
$App->print("Action: $Action\n");
$App->print("Mode  : $Mode\n");
$App->print("Index : $Index\n");
$App->print("nStep : $nStep\n");
$App->print("nOrder: $nOrder = ");
$nOrder       = eval($nOrder);
$App->print("$nOrder\n");
$App->print("pressure0    : $pressure0 atm\n");
$App->print("pressure_tol : $pressure_tol atm\n");
$App->print("fdiff = da/a = dV/V: $fdiff\n");
$App->print("Engine: $Engine\n");

my (@CopyList, $RunScriptPath, $RunScriptTemplate);
if($Engine =~ /^gulp$/i) {
	@CopyList = ("GoGULP.bat", "runForDielectricConstant.bat");
	$RunScriptPath = "runForDielectricConstant.bat"
}
elsif($Engine =~ /^pwscf$/i) {
	if($OS eq 'MSWin32') {
#		@CopyList = ("GoPWSCF.bat", "runForDielectricConstant.bat");
		if($Mode eq 'Vconst') {
			$RunScriptTemplate = Deps::MakePath($TemplateSource, "run-pwscf-Vconst-template.bat");
		}
		elsif($Mode eq 'Vconp') {
			$RunScriptTemplate = Deps::MakePath($TemplateSource, "run-pwscf-Vconp-template.bat");
		}
		elsif($Mode eq 'Cabc') {
			$RunScriptTemplate = Deps::MakePath($TemplateSource, "run-pwscf-Cabc-template.bat");
		}
		else {
			$App->print("\nError: Mode [$Mode] is not implemented for Engine [$Engine].\n\n");
			exit;
		}
		$RunScriptPath = "CalculateElasticConstant.bat"
	}
	else {
		@CopyList = ("GoPWSCF.sh", "runForDielectricConstant.sh");
		$RunScriptPath = "runForDielectricConstant.sh"
	}
}
elsif($Engine =~ /^vasp$/i) {
	@CopyList = ("DoVASP.sh");
	$RunScriptPath = "bash DoVASP.sh"
}
if($Action =~ /^run$/i) {
	$App->print("Engine: $Engine\n") if($Action =~ /^run$/i);
	$App->print("Copy files [", join(', ', @CopyList), "] from [$TemplateDir].\n");
	$App->print("Run script: $RunScriptPath\n");
}

#===============================================
# 実行
#===============================================
if($Action =~ /^MakeInputs$/i) {
	&MakeInputs();
}
elsif($Action =~ /^MakeRunScript$/i) {
	&MakeRunScript();
}
elsif($Action =~ /^run$/i) {
	&run();
}
elsif($Action =~ /^cal$/i) {
	&CalculateElasticTensor();
}
else {
	$App->InvalidParameter('Action', 'Main', 1);
	exit 1;
}

exit;

#=========================================
# Subroutines
#=========================================
sub CalculateElasticTensor
{
	$App->print("\nCalculate elastic tensor for Mode {$Mode]\n");

	$App->print("CIF File: $CIFFile\n");
	my $CIF = $App->{'CIF'} = new CIF;
	unless($CIF->Read($CIFFile, 0)) {
		$App->print("Error: Can not read [$CIFFile]\n\n");
		exit 2;
	}
	my $Crystal = $CIF->GetCCrystal();

	if($Mode eq 'Cabc') {
		&CalculateElasticTensorCabc($Crystal, $Mode);
		return;
	}
	elsif($Mode eq 'CXYZ') {
print "\nError: Mode=$Mode is not implemented\n\n";	
		&CalculateElasticTensorCXYZ($Crystal, $Mode);
		return;
	}
	elsif($Mode eq 'Vconst') {
		&CalculateElasticTensorVconst($Crystal, $Mode);
		return;
	}
	elsif($Mode eq 'Vconp') {
		&CalculateElasticTensorVconp($Crystal, $Mode);
		return;
	}
	else {
		$App->InvalidParameter('Mode', 'CalculateElasticTensor', 1);
		exit;
	}
}

sub CalculateElasticTensorVconp
{
	my ($Crystal, $Mode) = @_;

	my $table = JFile->new($TableFile, "r");
	if(!defined $table) {
		$App->print("\nError in CalculateElasticTensorVconst: Can not read [$TableFile].\n\n");
		exit;
	}

	$App->print("\n");
	$App->print("Original cell\n");
	my ($a2, $b2, $c, $alpha, $beta, $gamma) = $Crystal->LatticeParameters();
	$App->print("  Latt: $a2, $b2, $c A   $alpha, $beta, $gamma\n");
	my $vol = $Crystal->Volume();
	$App->print("  Volume: $vol A^3\n");
	$App->print("\n");

	$App->print("Data in [$TableFile]\n");
	my @data;
	my $idx = 0;
	while(1) {
		my $line = $table->ReadLine();
		last if(!defined $line);

		my ($li, $lj, $P, $file) = Utils::Split("\\s+", $line);
		my $header = $file;
		$header =~ s/\.cif//i;
		print "  $idx: $header: $li, $lj, $P, $file\n";
		$P *= 0.001;
		my $dir = "P$P";
		my $InputFile = "$dir/original.pwin";
		my $OutFile   = "$dir/original.pwout";
		print "  Input file: $InputFile\n";

		my $PWSCF   = new PWSCF;
		my $nStep    = $PWSCF->GetnFrameInOutput($App, $OutFile);
		my $Crystal2 = $PWSCF->ReadCrystalStructureFromInput($App, $InputFile, $Crystal);
		my $in  = JFile->new($OutFile, "r") or die "$!: Can not read [$OutFile].\n";
		for(my $i = 0 ; $i < $nStep ; $i++)
		{
			last if(!$PWSCF->ReadNextCrystalStructureFromOutput($App, $in, $i, $Crystal));
		}
		$in->Close();

		my @latt = $Crystal2->LatticeParameters();
		my $V    = $Crystal2->Volume();
		my $Etot = &ReadEtot($Mode, $Engine, $header, $OutFile);

		$data[$idx] = {
			li    => $li,
			lj    => $lj,
			P     => $P,
			V     => $V,
			Etot  => $Etot,
			pLatt => \@latt,
			file  => $file,
			dir   => $dir,
			};
		$App->print("    V=$V  latt=($a, $b, $c, $alpha, $beta, $gamma)  Etot=$Etot eV\n");
		$idx++;
	}
	$table->Close();
	$App->print("\n");

	my $ndata = @data;
	@data = sort { $a->{V} <=> $b->{V}; } @data;
	for(my $i = 0 ; $i < $ndata ; $i++) {
		my $p = $data[$i];
		print("file $p->{file} in $p->{dir}   P=$p->{P}  V=$p->{V}  Etot=$p->{Etot} eV\n");
	}

	$App->print("\n");
	$App->print("Least squares fitting\n");
	$App->print("  nData : $ndata\n");
	$App->print("  nOrder: $nOrder\n");
	my (@x, @y);
	for(my $i = 0 ; $i < $ndata ; $i++) {
		$x[$i] = $data[$i]->{V};
		$y[$i] = $data[$i]->{Etot};
	}

	my $iPrintLevel = 0;
#$nOrder=3;
	my $Opt = new Optimize;
	my ($pCi, $S2) = $Opt->Optimize("mlsq", \@x, \@y, $nOrder, $iPrintLevel);
	$App->print("Optimized at (", join(',',@{$pCi}), ") with S2=$S2\n");
	$App->print("\n");

#	$App->print("dV/V=", ($data[1]->{V} - $data[0]->{V}) / $vol, "\n");
#	$App->print("dV/V=", ($data[2]->{V} - $data[1]->{V}) / $vol, "\n");
#	$App->print("dV/V=", ($data[3]->{V} - $data[2]->{V}) / $vol, "\n");

	my $csv = JFile->new($SummaryCSV, "w");
	if(!defined $csv) {
		$App->print("\nError in CalculateElasticTensorVconp: Can not write to [$SummaryCSV].\n\n");
		exit;
	}
	print("i,j,P,V(A^3),a(A),b(A),c(A),alpha,beta,gamma,Etot(eV),Etot(fit),d2E/dV2(fit),BVfit(GPa)\,file,dir\n");
	$csv->print("i,j,P,V(A^3),a(A),b(A),c(A),alpha,beta,gamma,Etot(eV),Etot(fit),d2E/dV2(fit),BVfit(GPa),file,dir\n");
	my ($Vmin, $Vmax) = (1.0e100, -1.0e100);
	for(my $i = 0 ; $i < @data ; $i++) {
		my $p = $data[$i];
		my $pl   = $p->{pLatt};
		my $V    = $p->{V};
		$Vmin = $V if($Vmin > $V);
		$Vmax = $V if($Vmax < $V);
		my $Ecal = $Opt->YCal($V);
		my ($d2EdV2fit, $BVfit);
		$d2EdV2fit = 2.0 * $pCi->[2];
		for(my $j = 3 ; $j <= $nOrder ; $j++) {
			$d2EdV2fit += $j * ($j-1) * $pCi->[$j] * $V**($j-2);
		}
		$BVfit = $d2EdV2fit * $eVToJ * $vol * 1.0e30 * 1.0e-9; # GPa
		print("$p->{li},$p->{lj},$p->{P},$V,$pl->[0],$pl->[1],$pl->[2],$pl->[3],$pl->[4],$pl->[5],"
			 ."$p->{Etot},$Ecal,$d2EdV2fit,$BVfit,$p->{file},$p->{dir}\n");
		$csv->print("$p->{li},$p->{lj},$p->{P},$V,$pl->[0],$pl->[1],$pl->[2],$pl->[3],$pl->[4],$pl->[5],"
			 	   ."$p->{Etot},$Ecal,$d2EdV2fit,$BVfit,$p->{file},$p->{dir}\n");
	}
	$csv->print("\n");
	$csv->print("V(A^3),Etot(fit)\n");
	my $n = 101;
	my $Vstep = ($Vmax - $Vmin) / ($n - 1);
	for(my $i = 0 ; $i < $n ; $i++) {
		my $V = $Vmin + $i * $Vstep;
		my $Ecal = $Opt->YCal($V);
		$csv->print("$V,$Ecal\n");
	}
	$csv->Close();

	$App->print("\n");
}

sub CalculateElasticTensorVconst
{
	my ($Crystal, $Mode) = @_;

	my ($a2, $b2, $c, $alpha, $beta, $gamma) = $Crystal->LatticeParameters();
	$App->print("  Latt: $a2, $b2, $c A   $alpha, $beta, $gamma\n");
	my $vol = $Crystal->Volume();
	$App->print("  Volume: $vol A^3\n");
	$App->print("\n");

	my $table = JFile->new($TableFile, "r");
	if(!defined $table) {
		$App->print("\nError in CalculateElasticTensorVconst: Can not read [$TableFile].\n\n");
		exit;
	}

	my @data;
	my $idx = 0;
	while(1) {
		my $line = $table->ReadLine();
		last if(!defined $line);

		my ($li, $lj, $fdV, $file) = Utils::Split("\\s+", $line);
		my $header = $file;
		$header =~ s/\.cif//i;

		my $Etot = &ReadEtot($Mode, $Engine, $header);
		my $V = $vol * (1.0 + $fdV);
		$data[$idx] = {
			li   => $li,
			lj   => $lj,
			fdV  => $fdV,
			V    => $V,
			Etot => $Etot,
			file => $file,
			};
#		print "  $file  dV=$fdV  Etot = $Etot eV\n";
#		$csv->print("$li,$lj,$fdV,$V,$Etot,$file\n");
		$idx++;
	}
	$table->Close();

	my $ndata = @data;
	@data = sort { $a->{V} <=> $b->{V}; } @data;
	for(my $i = 0 ; $i < $ndata ; $i++) {
		my $p = $data[$i];
		print("file $p->{file}   dV=$p->{fdV}  V=$p->{V}  Etot=$p->{Etot} eV\n");
	}

	$App->print("\n");
	$App->print("Least squares fitting\n");
	$App->print("  nData : $ndata\n");
	$App->print("  nOrder: $nOrder\n");
	my (@x, @y);
	for(my $i = 0 ; $i < $ndata ; $i++) {
		$x[$i] = $data[$i]->{V};
		$y[$i] = $data[$i]->{Etot};
	}
	my $iPrintLevel = 0;
#$nOrder=3;
	my $Opt = new Optimize;
	my ($pCi, $S2) = $Opt->Optimize("mlsq", \@x, \@y, $nOrder, $iPrintLevel);
	$App->print("Optimized at (", join(',',@{$pCi}), ") with S2=$S2\n");
	$App->print("\n");

	$App->print("dV/V=", ($data[1]->{V} - $data[0]->{V}) / $vol, "\n");
#	$App->print("dV/V=", ($data[2]->{V} - $data[1]->{V}) / $vol, "\n");
#	$App->print("dV/V=", ($data[3]->{V} - $data[2]->{V}) / $vol, "\n");

	my $csv = JFile->new($SummaryCSV, "w");
	if(!defined $csv) {
		$App->print("\nError in CalculateElasticTensorVconst: Can not write to [$SummaryCSV].\n\n");
		exit;
	}
	print("i,j,fdV,V(A^3),Etot(eV),Etot(fit),d2E/dV2(fit),d2E/dV2(numdiff),BVfit(GPa),BVnumdiff(GPa),file\n");
	$csv->print("i,j,fdV,V(A^3),Etot(eV),Etot(fit),d2E/dV2(fit),d2E/dV2(numdiff),BVfit(GPa),BVnumdiff(GPa),file\n");
	my ($Vmin, $Vmax) = (1.0e100, -1.0e100);
	for(my $i = 0 ; $i < @data ; $i++) {
		my $p = $data[$i];
		my $V    = $p->{V};
		$Vmin = $V if($Vmin > $V);
		$Vmax = $V if($Vmax < $V);
		my $Ecal = $Opt->YCal($V);
		my ($d2EdV2cal, $d2EdV2fit, $BVcal, $BVfit);
		$d2EdV2fit = 2.0 * $pCi->[2];
		for(my $j = 3 ; $j <= $nOrder ; $j++) {
			$d2EdV2fit += $j * ($j-1) * $pCi->[$j] * $V**($j-2);
		}
		$BVfit = $d2EdV2fit * $eVToJ * $vol * 1.0e30 * 1.0e-9; # GPa
		if($i >= 1 and $i <= @data-2) {
			$d2EdV2cal = ($data[$i+1]->{Etot} + $data[$i-1]->{Etot} - 2.0*$data[$i]->{Etot}) / ($data[$i+1]->{V} - $data[$i]->{V})**2;
			$BVcal = $d2EdV2cal * $eVToJ * $vol * 1.0e30 * 1.0e-9; # GPa
		}
		print("$p->{li},$p->{lj},$p->{fdV},$V,$p->{Etot},$Ecal,$d2EdV2fit,$d2EdV2cal,$BVfit,$BVcal,$p->{file}\n");
		$csv->print("$p->{li},$p->{lj},$p->{fdV},$V,$p->{Etot},$Ecal,$d2EdV2fit,$d2EdV2cal,$BVfit,$BVcal,$p->{file}\n");
	}
	$csv->print("\n");
	$csv->print("V(A^3),Etot(fit)\n");
	my $n = 101;
	my $Vstep = ($Vmax - $Vmin) / ($n - 1);
	for(my $i = 0 ; $i < $n ; $i++) {
		my $V = $Vmin + $i * $Vstep;
		my $Ecal = $Opt->YCal($V);
		$csv->print("$V,$Ecal\n");
	}
	$csv->Close();
	$csv->Close();

	$App->print("\n");
}

sub ReadEtot
{
	my ($Mode, $Engine, $file, $OutFile) = @_;

	my ($Etot);
	my $dir = $file;
	$dir =~ s/\..*?$//;
	if($Engine =~ /^gulp$/i) {
		my $OutFile = "$dir.out";
		$App->print("  Read Etot from [$OutFile].\n");
		my $in = JFile->new($OutFile, "r") or die "$!: Can not read [$OutFile].\n\n";
		my $line = $in->SkipTo("Total lattice enthalpy\\s*=");
		($Etot) = ($line =~ /=\s*([+\-\.\d]+)/);
		$in->Close();
	}
	elsif($Engine =~ /^pwscf$/i) {
		$OutFile = "$dir/$dir.pwout" if(!defined $OutFile);
		$App->print("  Read Etot from [$OutFile].\n");
		my $in = JFile->new($OutFile, "r") or die "$!: Can not read [$OutFile].\n\n";
		while(1) {
			my $line = $in->SkipTo("!\\s*total\\s+energy\\s*=");
			last if(!defined $line);
			($Etot) = ($line =~ /=\s*([+\-\.\d]+)/);
		}
		$in->Close();
#print "Etot=$Etot\n";
		$Etot *= $RyToeV;
#print "Etot=$Etot\n";
#exit if($OutFile =~ /P1/);
	}
	elsif($Engine =~ /^vasp$/i) {
		my $OutFile = "SCF/OUTCAR";
		$App->print("  Read Etot from [$OutFile].\n");
		my $in = JFile->new($OutFile, "r") or die "$!: Can not read [$OutFile].\n\n";
		while(1) {
			my $line = $in->SkipTo("free\\s+energy\\s+TOTEN\\s*=");
			last if(!defined $line);
			($Etot) = ($line =~ /=\s*([+\-\.\d]+)/);
		}
		$in->Close();
	}
	else {
		$App->InvalidParameter('Engine', 'ReadEtot', 1);
		exit;
	}
	return $Etot;
}

sub CalculateElasticTensorCabc
{
	my ($Crystal, $Mode) = @_;

	my ($a, $b, $c, $alpha, $beta, $gamma) = $Crystal->LatticeParameters();
	$App->print("  Latt: $a, $b, $c A   $alpha, $beta, $gamma\n");
	my $vol = $Crystal->Volume();
	$App->print("  Volume: $vol A^3\n");
	$App->print("\n");

	my $csv = JFile->new($SummaryCSV, "r") or die "$!: Can not read [$SummaryCSV].\n";
	my $line = $csv->ReadLine();
	my ($Etot0, @eiJ, @EtotpiJ, @EtotmiJ);
	while(1) {
		$line = $csv->ReadLine();
		last if(!defined $line);

		my ($li, $lj, $eiJ, $Etot, $file) = Utils::Split("[\\s,]+", $line);
		$App->printf(" %6s=%9.6f  Etot=%12.6g eV [$file]\n", "e$li$lj", $eiJ, $Etot);
		if($eiJ == 0.0 or $li < 0 or $lj < 0) {
			$Etot0 = $Etot;
		}
		elsif($eiJ > 0.0) {
			$EtotpiJ[$li][$lj] = $Etot;
			$eiJ[$li][$lj] = $eiJ;
		}
		else {
			$EtotmiJ[$li][$lj] = $Etot;
			$eiJ[$li][$lj] = $eiJ;
		}
	}
	
	$App->print("\n");
	$App->print("Elastic tensor\n");
	$App->print("  Etot0 = $Etot0 eV\n");
	my @CIjKl;
	for(my $i = 0 ; $i < 3 ; $i++) {
		for(my $j = $i ; $j < 3 ; $j++) {
			$CIjKl[$i][$j][$i][$j] 
				= ($EtotpiJ[$i][$j] + $EtotmiJ[$i][$j] - 2.0 * $Etot0) / ($fdiff*$fdiff);
			$CIjKl[$i][$j][$i][$j] *= $eVToJ / ($vol * 1.0e-30);
			$App->printf("  CIjKl[$i][$j][$i][$j] = %12.6g GPa = %12.6g kBar\n", 
				$CIjKl[$i][$j][$i][$j] * 1.0e-9, $CIjKl[$i][$j][$i][$j] * 1.0e-8);
		}
	}
	
	$App->print("\n");
}

sub MakeRunScript
{
	$App->print("\n");
	$App->print("Make run script [$RunScriptPath] from template [$RunScriptTemplate] and table [$TableFile]\n");
	$App->print("\n");

	if($Mode eq 'Cabc') {
		&MakeRunScriptCabc();
		return;
	}
	elsif($Mode eq 'CXYZ') {
print "\nError: Mode=$Mode is not implemented\n\n";	
		&MakeRunScriptCXYZ();
		return;
	}
	elsif($Mode eq 'Vconst') {
		&MakeRunScriptVconst();
		return;
	}
	elsif($Mode eq 'Vconp') {
		&MakeRunScriptVconp();
		return;
	}
	else {
		$App->InvalidParameter('Mode', 'MakeRunScript', 1);
		exit;
	}
}

sub MakeRunScriptCabc
{
	my $str   = JFile->ReadFile($RunScriptTemplate);
	if(!defined $str) {
		$App->print("\nError in MakeRunScriptVconp: Can not read [$RunScriptTemplate].\n\n");
		exit;
	}

	my $table = JFile->new($TableFile, "r");
	if(!defined $table) {
		$App->print("\nError in MakeRunScriptVconp: Can not read [$TableFile].\n\n");
		exit;
	}

	my (@forward, @reverse);
	while(1) {
		my $line = $table->ReadLine();
		last if(!defined $line);

		my ($i1, $j1, $i2, $j2, $ie1, $ie2, $e, $file) = Utils::Split("\\s+", $line);
		$e += 0.0;
		if($e < 0.0) {
			$App->print(" Negative strain: e=$e atm from [$file].\n");
			$file =~ s/\.cif//i;
			push(@reverse, $e);
		}
		elsif($e > 0.0) {
			$App->print(" Positive strain: e=$e atm from [$file].\n");
			$file =~ s/\.cif//i;
			push(@forward, $e); # atm to kBar
		}
		else {
			$App->print(" Zero strain: e=$e atm from [$file].\n");
		}
	}
	$table->Close();

	my $sr = join(' ', @reverse);
	my $sf = join(' ', @forward);
	$str =~ s/\{ListReverse\}/$sr/i;
	$str =~ s/\{ListForward\}/$sf/i;
	$str =~ s/\{SampleName\}/$SampleName/i;
	$str =~ s/\{CIFHeader\}/e${Index}_/i;

	$App->print("Save to [$RunScriptPath].\n");
	my $ret = JFile->SaveFile($str, $RunScriptPath, "w");
	if(!defined $ret) {
		$App->print("\nError in MakeRunScriptCabc: Can not write to [$RunScriptPath].\n\n");
		exit;
	}
}

sub MakeRunScriptVconp
{
	my $str   = JFile->ReadFile($RunScriptTemplate);
	if(!defined $str) {
		$App->print("\nError in MakeRunScriptVconp: Can not read [$RunScriptTemplate].\n\n");
		exit;
	}

	my $table = JFile->new($TableFile, "r");
	if(!defined $table) {
		$App->print("\nError in MakeRunScriptVconp: Can not read [$TableFile].\n\n");
		exit;
	}

	my (@forward, @reverse);
	while(1) {
		my $line = $table->ReadLine();
		last if(!defined $line);

		my ($li, $lj, $P, $file) = Utils::Split("\\s+", $line);
		if($P <= 0.0) {
			$App->print(" Negative pressure: P=$P atm from [$file].\n");
			$file =~ s/\.cif//i;
			push(@reverse, $P * 0.001);
		}
		elsif($P > 0.0) {
			$App->print(" Positive pressure: P=$P atm from [$file].\n");
			$file =~ s/\.cif//i;
			push(@forward, $P * 0.001); # atm to kBar
		}
		else {
#			$App->print(" Zero pressure: P=$P atm from [$file].\n");
		}
	}
	$table->Close();

	my $sr = join(' ', @reverse);
	my $sf = join(' ', @forward);
	$str =~ s/\{PListReverse\}/$sr/i;
	$str =~ s/\{PListForward\}/$sf/i;
	$str =~ s/\{SampleName\}/$SampleName/i;
	$pressure0 *= 0.001;
	$pressure_tol *= 0.001;
	$str =~ s/\{pressure\}/$pressure0/i;
	$str =~ s/\{pressure_tol\}/$pressure_tol/i;

	$App->print("Save to [$RunScriptPath].\n");
	my $ret = JFile->SaveFile($str, $RunScriptPath, "w");
	if(!defined $ret) {
		$App->print("\nError in MakeRunScriptVconp: Can not write to [$RunScriptPath].\n\n");
		exit;
	}
}

sub MakeRunScriptVconst
{
	my $str   = JFile->ReadFile($RunScriptTemplate);
	if(!defined $str) {
		$App->print("\nError in MakeRunScriptVconst: Can not read [$RunScriptTemplate].\n\n");
		exit;
	}

	my $table = JFile->new($TableFile, "r");
	if(!defined $table) {
		$App->print("\nError in MakeRunScriptVconst: Can not read [$TableFile].\n\n");
		exit;
	}

	my (@forward, @reverse);
	while(1) {
		my $line = $table->ReadLine();
		last if(!defined $line);
		
		my ($li, $lj, $eiJ, $file) = Utils::Split("\\s+", $line);
		if($eiJ < 0.0) {
			$App->print(" Shrunk structure: dV/V==$eiJ from [$file].\n");
			$file =~ s/\.cif//i;
			push(@reverse, "$file");
		}
		elsif($eiJ > 0.0) {
			$App->print(" Expanded structure: dV/V==$eiJ from [$file].\n");
			$file =~ s/\.cif//i;
			push(@forward, "$file");
		}
		else {
			$App->print(" Original structure: dV/V==$eiJ from [$file].\n");
		}
	}
	$table->Close();

	my $sr = join(' ', @reverse);
	my $sf = join(' ', @forward);
	$str =~ s/\{FileListReverse\}/$sr/i;
	$str =~ s/\{FileListForward\}/$sf/i;
	$str =~ s/\{SampleName\}/$SampleName/i;

	$App->print("Save to [$RunScriptPath].\n");
	my $ret = JFile->SaveFile($str, $RunScriptPath, "w");
	if(!defined $ret) {
		$App->print("\nError in MakeRunScriptVconst: Can not write to [$RunScriptPath].\n\n");
		exit;
	}
}

sub run
{
	$App->print("\n");
	$App->print("Calculate total energies\n");
	$App->print("Base directory: $BaseDir\n");
	$App->print("\n");

	my $table = JFile->new($TableFile, "r") or die "$!: Can not read [$TableFile].\n";
	my $out   = JFile->new($SummaryCSV, "w") or die "$!: Can not write to [$SummaryCSV].\n";
	$out->print("i,j,eiJ,Etot(eV),file\n");
	
	if($Mode eq 'Cabc') {
		&runCabs($table, $out, $BaseDir);
		return;
	}
	elsif($Mode eq 'CXYZ') {
print "\nError: Mode=$Mode is not implemented\n\n";	
		&runCXYZ($table, $out, $BaseDir);
		return;
	}
	elsif($Mode eq 'Vconst') {
		&runVconst($table, $out, $BaseDir);
		return;
	}
	elsif($Mode eq 'Vconp') {
		&runVconp($table, $out, $BaseDir);
		return;
	}
	else {
		$App->InvalidParameter('Mode', 'run', 1);
		exit;
	}
}

sub runVconst
{
	my ($table, $out, $BaseDir) = @_;
	
	while(1) {
		my $line = $table->ReadLine();
		last if(!defined $line);
		
		my ($li, $lj, $eiJ, $file) = Utils::Split("\\s+", $line);
		$App->print(" for e$li$lj=$eiJ from [$file].\n");
		
		my $dir = $file;
		$dir =~ s/\..*?$//;
		$App->print("  Create dir [$dir].\n");
		Utils::CreateDirectory($dir);
#my $WorkDir = cwd();
#$App->print("  Working dir: [$WorkDir]\n");
#system("ls");
		$App->print("  Copy [$file, ", join(', ', @CopyList), "] to [$dir].\n");
		if(!Utils::CopyFile($file, $dir)) {
			$App->print("\nError: Can not copy [$file] to [$dir]\n\n");
			exit;
		}
		foreach my $f (@CopyList) {
			my $path = "$TemplateDir/$f";
			if(!Utils::CopyFile($path, $dir)) {
				$App->print("\nError: Can not copy [$path] to [$dir]\n\n");
				exit;
			}
		}

		$App->print("  cd [$dir]\n");
		if(!chdir($dir)) {
			$App->print("\nError: Can not change dir to [$dir]\n\n");
			exit;
		}
		my $WorkDir = cwd();
		$App->print("  Working dir: [$WorkDir]\n");
		my $cmd;
		if($Engine =~ /^gulp$/i) {
			$cmd = "$RunScriptPath $dir";
		}
		elsif($Engine =~ /^vasp$/i) {
			$cmd = "$RunScriptPath";
		}
		$App->print("  run [$cmd]\n");
#		system("chmod +x $RunScriptPath");
		system($cmd);

		my $Etot = &ReadEtot($Mode, $Engine, $file);
		$App->print("  Etot = $Etot eV\n");
		$out->print("$li,$lj,$eiJ,$Etot,$file\n");

		$App->print("  cd [$BaseDir]\n");
		if(!chdir($BaseDir)) {
			$App->print("\nError: Can not change dir to [$dir]\n\n");
			exit;
		}

		$App->print("\n");
	}
	
	$out->Close();
	$table->Close();
}


sub runCabc
{
	my ($table, $out) = @_;
	
	while(1) {
		my $line = $table->ReadLine();
		last if(!defined $line);
		
		my ($li, $lj, $eiJ, $file) = Utils::Split("\\s+", $line);
		$App->print(" for e$li$lj=$eiJ from [$file].\n");
		
		my $dir = $file;
		$dir =~ s/\..*?$//;
		$App->print("  Create dir [$dir].\n");
		Utils::CreateDirectory($dir);
#my $WorkDir = cwd();
#$App->print("  Working dir: [$WorkDir]\n");
#system("ls");
		$App->print("  Copy [$file, ", join(', ', @CopyList), "] to [$dir].\n");
		if(!Utils::CopyFile($file, $dir)) {
			$App->print("\nError: Can not copy [$file] to [$dir]\n\n");
			exit;
		}
		foreach my $f (@CopyList) {
			my $path = "$TemplateDir/$f";
			if(!Utils::CopyFile($path, $dir)) {
				$App->print("\nError: Can not copy [$path] to [$dir]\n\n");
				exit;
			}
		}

		$App->print("  cd [$dir]\n");
		if(!chdir($dir)) {
			$App->print("\nError: Can not change dir to [$dir]\n\n");
			exit;
		}
		my $WorkDir = cwd();
		$App->print("  Working dir: [$WorkDir]\n");
		my $cmd;
		if($Engine =~ /^gulp$/i) {
			$cmd = "$RunScriptPath $dir";
		}
		elsif($Engine =~ /^vasp$/i) {
			$cmd = "$RunScriptPath";
		}
		$App->print("  run [$cmd]\n");
#		system("chmod +x $RunScriptPath");
		system($cmd);

		my $Etot = &ReadEtot($Mode, $Engine, $file);
		$App->print("  Etot = $Etot eV\n");
		$out->print("$li,$lj,$eiJ,$Etot,$file\n");

		$App->print("  cd [$BaseDir]\n");
		if(!chdir($BaseDir)) {
			$App->print("\nError: Can not change dir to [$dir]\n\n");
			exit;
		}

		$App->print("\n");
	}
	
	$out->Close();
	$table->Close();
}

sub MakeInputs
{
	$App->print("\n");
	$App->print("Make input files for calculating dielectric constant by finite method\n");
	$App->print("\n");

	$App->print("CIF File: $CIFFile\n");
	$App->print("\n");
	
#CIFファイル名を引数に与えて読み込む
	my $CIF = $App->{'CIF'} = new CIF;
	unless($CIF->Read($CIFFile, 0)) {
		$App->print("Error: Can not read [$CIFFile]\n\n");
		exit 2;
	}
	$App->print("\n");

	my $Crystal = $CIF->GetCCrystal();
#$Crystal->SplitPartialSites();
	$Crystal->ExpandCoordinates();
#pCrystal->m_ASFDCDBFile  = ASFDCPath;
#pCrystal->m_NonrelDBFile = NONRELPath;
#pCrystal->m_SPGDBFile    = SPGDBPath;
#pCrystal->SetFileName(InputFile);
#(pCrystal->SpaceGroup).SetNoSymmetry();
	my $CrystalName = $Crystal->CrystalName();
	$App->print("CrystalName: $CrystalName\n");

	if($Mode eq 'Cabc') {
		&MakeInputsCabc($Crystal);
		return;
	}
	elsif($Mode eq 'CXYZ') {
print "\nError: Mode=$Mode is not implemented\n\n";	
		&MakeInputsCXYZ($Crystal);
		return;
	}
	elsif($Mode eq 'Vconst') {
		&MakeInputsVconst($Crystal);
		return;
	}
	elsif($Mode eq 'Vconp') {
		&MakeInputsVconp($Crystal);
		return;
	}
	else {
		$App->InvalidParameter('Mode', 'MakeInputs', 1);
		exit;
	}
}

sub MakeInputsCabc
{
	my ($Crystal) = @_;

	my @index = ($Index =~ /^(\d)(\d)(\d)(\d)/i);
	$App->print("Calculate elastic tensor for Cabc_{", join(',', @index), "} in lattice coordinate.\n");
	if($index[0] > $index[1]) {
		($index[0], $index[1]) = ($index[1], $index[0]);
	}
	if($index[2] > $index[3]) {
		($index[2], $index[3]) = ($index[3], $index[2]);
	}
	if($index[0] != $index[2] or $index[1] != $index[3]) {
		$App->print("\nError in MakeInputsCabc: Non-diagonal component C", join('', @index), " is not supported.\n\n");
		exit;
	}
	
	$App->print("Index: ", join(', ', @index), "\n");
	$App->print("fdiff: $fdiff\n");
	$App->print("nStep: $nStep\n");

	my ($a,$b,$c,$alpha,$beta,$gamma) = $Crystal->LatticeParameters();
	$App->printf("cell: %9.6f %9.6f %9.6f   %9.6f %9.6f %9.6f \n", $a, $b, $c, $alpha, $beta, $gamma);

#===============================================
# 歪CIF出力
#===============================================
	my $table = JFile->new($TableFile, "w") or die "$!: Can not write to [$TableFile].\n";

	print "\n";
	print "Save original CIF to [$OriginalCIF]\n";
	my @latt = $Crystal->LatticeParameters();
	my @lv   = $Crystal->LatticeVectors();
	$App->printf("  latt: %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f\n", @latt);
	$App->printf("  vec : %9.6f %9.6f %9.6f\n"
				."        %9.6f %9.6f %9.6f\n"
				."        %9.6f %9.6f %9.6f\n", @lv);
	my $NewCrystal = $Crystal->Duplicate();
	$NewCrystal->BurstToP1();
	&SaveCIF($NewCrystal, $OriginalCIF, 0, 'unix');
		$table->printf("%2d %2d %2d %2d %2d %2d %10.6f %s\n", 
			$index[0], $index[1], $index[2], $index[3], 0, 0, 0.0, $OriginalCIF);
	$App->print("\n");

	print "\n";
	$App->print("Strained models\n");
	my $n = -int($nStep / 2);
	for(my $i = -1 ; $i >= $n ; $i--) {
		my $e = $i * $fdiff;
		$App->printf("  %3d: e$index[0]$index[1]=$e\n", $i);
#		$App->printf("  %3d: e$index[2]$index[3]=$e\n", $i);

		my $OutCIF = "e${Index}_$e.cif";
		print "Save strained CIF to [$OutCIF]\n";
		my ($eiJ) = &SaveStrainedStructureCabc(
			$Crystal, $e, $index[0], $index[1], $OutCIF);
		$table->printf("%2d %2d %2d %2d %2d %2d %10.6f %s\n", 
			$index[0], $index[1], $index[2], $index[3], $i, 0, $e, $OutCIF);
	}
	for(my $i = 1 ; $i < $n + $nStep ; $i++) {
		my $e = $i * $fdiff;
		$App->printf("  %3d: e$index[0]$index[1]=$e\n", $i);
#		$App->printf("  %3d: e$index[2]$index[3]=$e\n", $i);

		my $OutCIF = "e${Index}_$e.cif";
		print "Save strained CIF to [$OutCIF]\n";
		my ($eiJ) = &SaveStrainedStructureCabc(
			$Crystal, $e, $index[0], $index[1], $OutCIF);
		$table->printf("%2d %2d %2d %2d %2d %2d %10.6f %s\n", 
			$index[0], $index[1], $index[2], $index[3], $i, 0, $e, $OutCIF);
	}

	$table->Close();
}

sub SaveStrainedStructureCabc
{
	my ($Crystal, $e, $li, $lj, $OutCIF) = @_;

	my @v;
	my @latt = $Crystal->LatticeParameters();
	($v[0][0], $v[0][1], $v[0][2], $v[1][0], $v[1][1], $v[1][2], $v[2][0], $v[2][1], $v[2][2])
		= $Crystal->LatticeVectors();
	my $NewCrystal = $Crystal->Duplicate();
	$NewCrystal->BurstToP1();

	$v[$li][0] += 0.5 * $e * $v[$lj][0] * $latt[$li] / $latt[$lj];
	$v[$li][1] += 0.5 * $e * $v[$lj][1] * $latt[$li] / $latt[$lj];
	$v[$li][2] += 0.5 * $e * $v[$lj][2] * $latt[$li] / $latt[$lj];
	$v[$lj][0] += 0.5 * $e * $v[$li][0] * $latt[$lj] / $latt[$li];
	$v[$lj][1] += 0.5 * $e * $v[$li][1] * $latt[$lj] / $latt[$li];
	$v[$lj][2] += 0.5 * $e * $v[$li][2] * $latt[$lj] / $latt[$li];

	$NewCrystal->SetLatticeVectors(
		$v[0][0], $v[0][1], $v[0][2], $v[1][0], $v[1][1], $v[1][2], $v[2][0], $v[2][1], $v[2][2]
		);
	$NewCrystal->CalMetrics();

	@latt = $NewCrystal->LatticeParameters();
	$App->printf("  latt: %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f\n", @latt);
	$App->printf("  vec : %9.6f %9.6f %9.6f\n"
				."        %9.6f %9.6f %9.6f\n"
				."        %9.6f %9.6f %9.6f\n",
			$v[0][0], $v[0][1], $v[0][2], $v[1][0], $v[1][1], $v[1][2], $v[2][0], $v[2][1], $v[2][2]);

	&SaveCIF($NewCrystal, $OutCIF, 0, 'unix');
	$App->print("\n");

	return ($e, $OutCIF);
}

sub SaveStrainedStructure
{
	my ($Mode, $Crystal, $li, $lj, $fdiff, $OutCIF) = @_;

	if($Mode eq 'CXYZ') {
		return &SaveStrainedStructureCXYZ($Crystal, $li, $lj, $fdiff, $OutCIF);
	}
#	elsif($Mode eq 'Cabc') {
#		return &SaveStrainedStructureCabc($Crystal, $li, $lj, $fdiff, $OutCIF);
#	}
	elsif($Mode eq 'V') {
		return &SaveStrainedStructureV($Crystal, $li, $lj, $fdiff, $OutCIF);
	}
	else {
		$App->InvalidParameter('Mode', 'SaveStrainedStructure', 1);
		exit;
	}
}

sub SaveStrainedStructureV
{
	my ($Crystal, $li, $lj, $k, $OutCIF) = @_;

	my @latt = $Crystal->LatticeParameters();

	my $NewCrystal = $Crystal->Duplicate();
	for(my $i = 0 ; $i < 3 ; $i++) {
		$latt[$i] *= $k;
	}
	$NewCrystal->SetLatticeParameters(@latt);
	$NewCrystal->CalMetrics();

	@latt = $NewCrystal->LatticeParameters();
	$App->printf("  latt: %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f\n", @latt);

	&SaveCIF($NewCrystal, $OutCIF, 0, 'unix');
	$App->print("\n");

	return ($k*$k*$k, $OutCIF);
}

sub MakeInputsVconp
{
	my ($Crystal) = @_;
	
	my $OriginalCIF = "original.cif";
	my $table = JFile->new($TableFile, "w") or die "$!: Can not write to [$TableFile].\n";
	&SaveCIF($Crystal, $OriginalCIF, 0, 'unix');
	$table->printf("%2d %2d %10.6f %s\n", 0, 0, $pressure0, $OriginalCIF);
	$App->print("\n");

#	print "\n";
	my $n = -int($nStep / 2);
print "n=$n\n";
	for(my $i = -1 ; $i >= $n ; $i--) {
#print "i=$i\n";
		my $Ptarget = $pressure0 + $i * $fdiff;
		$App->print("  $i: Ptarget=$Ptarget\n");
		$table->printf("%2d %2d %10.6f %s\n", $i, 0, $Ptarget, $OriginalCIF);
	}
	for(my $i = 1 ; $i < $n + $nStep ; $i++) {
#print "i=$i\n";
		my $Ptarget = $pressure0 + $i * $fdiff;
		$App->print("  $i: Ptarget=$Ptarget\n");
		$table->printf("%2d %2d %10.6f %s\n", $i, 0, $Ptarget, $OriginalCIF);
	}

	$table->Close();
}

sub MakeInputsVconst
{
	my ($Crystal) = @_;
	
	my ($a,$b,$c,$alpha,$beta,$gamma) = $Crystal->LatticeParameters();
	$App->printf("cell: %9.6f %9.6f %9.6f   %9.6f %9.6f %9.6f \n", $a, $b, $c, $alpha, $beta, $gamma);

	my $OriginalCIF = "original.cif";
	my $table = JFile->new($TableFile, "w") or die "$!: Can not write to [$TableFile].\n";
	&SaveCIF($Crystal, $OriginalCIF, 0, 'unix');
	$table->printf("%2d %2d %10.6f %s\n", 0, 0, 0.0, $OriginalCIF);
	$App->print("\n");

	print "\n";
	$App->print("Strained models\n");
	my $n = -int($nStep / 2);
#print "n=$n\n";
	for(my $i = -1 ; $i >= $n ; $i--) {
#print "i=$i\n";
		my $fVtarget = 1.0 + $i * $fdiff;
		my $k = $fVtarget**(1.0/3.0);
		$App->print("  $i: fVtarget=$fVtarget  klatt=$k\n");
		next if($fVtarget <= 0.0);

		my $dV = $fVtarget - 1.0;
		my $OutCIF = $StrainedCIFTemplate;
		$OutCIF = "dV" . sprintf("%05d", int($dV*10000)) . ".cif";
		print "Save strained CIF to [$OutCIF]\n";
		&SaveStrainedStructure($Mode, $Crystal, $i, 0, $k, $OutCIF);
		$table->printf("%2d %2d %10.6f %s\n", $i, 0, $dV, $OutCIF);
	}
	for(my $i = 1 ; $i < $n + $nStep ; $i++) {
#print "i=$i\n";
		my $fVtarget = 1.0 + $i * $fdiff;
		next if($fVtarget <= 0.0);

		my $k = $fVtarget**(1.0/3.0);
		$App->print("  $i: fVtarget=$fVtarget  klatt=$k\n");

		my $dV = $fVtarget - 1.0;
		my $OutCIF = $StrainedCIFTemplate;
		$OutCIF = "dV" . sprintf("%05d", int($dV*10000)) . ".cif";
		print "Save strained CIF to [$OutCIF]\n";
		&SaveStrainedStructure($Mode, $Crystal, $i, 0, $k, $OutCIF);
		$table->printf("%2d %2d %10.6f %s\n", $i, 0, $dV, $OutCIF);
	}

	$table->Close();
}

sub MakeInputsCabc0
{
	my ($Crystal) = @_;

	my ($a,$b,$c,$alpha,$beta,$gamma) = $Crystal->LatticeParameters();
	$App->printf("cell: %9.6f %9.6f %9.6f   %9.6f %9.6f %9.6f \n", $a, $b, $c, $alpha, $beta, $gamma);
#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 $Density = $Crystal->Density();
#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 @AtomAsymSiteList = $Crystal->GetCAsymmetricAtomSiteList();
	my $nAtomAsymSites   = @AtomAsymSiteList;
	my @AtomSiteList     = $Crystal->GetCExpandedAtomSiteList();
	my $nAtomSites       = $Crystal->nTotalExpandedAtomSite();

#===============================================
# 歪CIF出力
#===============================================

	my $table = JFile->new($TableFile, "w") or die "$!: Can not write to [$TableFile].\n";

	print "\n";
	print "Save original CIF to [$OriginalCIF]\n";
	my @latt = $Crystal->LatticeParameters();
	my @lv   = $Crystal->LatticeVectors();
	$App->printf("  latt: %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f\n", @latt);
	$App->printf("  vec : %9.6f %9.6f %9.6f\n"
				."        %9.6f %9.6f %9.6f\n"
				."        %9.6f %9.6f %9.6f\n", @lv);
	my $NewCrystal = $Crystal->Duplicate();
	$NewCrystal->BurstToP1();
	&SaveCIF($NewCrystal, $OriginalCIF, 0, 'unix');
	$table->printf("%2d %2d %10.6f %s\n", -1, -1, 0.0, $OriginalCIF);
	$App->print("\n");

	for(my $i = 0 ; $i < 3 ; $i++) {
		for(my $j = $i ; $j < 3 ; $j++) {

			my $OutCIF = $StrainedCIFTemplate;
			$OutCIF =~ s/{i}/$i/g;
			$OutCIF =~ s/{j}/$j/g;
			$OutCIF =~ s/{sign}/p/g;
			print "Save strained CIF to [$OutCIF]\n";
			my ($eiJ) = &SaveStrainedStructure($Mode, $Crystal, $i, $j, +$fdiff, $OutCIF);
			$table->printf("%2d %2d %10.6f %s\n", $i, $j, $eiJ, $OutCIF);
		
			$OutCIF = $StrainedCIFTemplate;
			$OutCIF =~ s/{i}/$i/g;
			$OutCIF =~ s/{j}/$j/g;
			$OutCIF =~ s/{sign}/m/g;
			print "Save strained CIF to [$OutCIF]\n";
			my ($eiJ) = &SaveStrainedStructure($Mode, $Crystal, $i, $j, -$fdiff, $OutCIF);
			$table->printf("%2d %2d %10.6f %s\n", $i, $j, $eiJ, $OutCIF);
		}
	}

	$table->Close();
}

sub SaveCIF
{
	my ($Crystal, $OutCIF, $IsChooseRandomly, $strCRLF) = @_;

	my $IsPrint = 1;
	my $cif = new CIF;
	$cif->CreateCIFFileFromCCrystal($Crystal, $OutCIF, $IsChooseRandomly, $strCRLF, $IsPrint);
}