#!/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 Crystal::GULP;
use Crystal::VASP;
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 $ConfigFile = 'configure.ini';
my $vars = new ProgVars;
my $pParams = {};
$pParams = $vars->ReadConfigure($ConfigFile, $pParams, 1, 1);

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|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 $RunScriptTemplate;
my $RunScriptPath = 'CalculateElasticConstant.bat';
if($Engine =~ /^gulp$/i) {
	if($OS eq 'MSWin32') {
#		@CopyList = ("GoGULP.bat", "runForDielectricConstant.bat");
		if($Mode eq 'Vconst') {
			$RunScriptTemplate = Deps::MakePath($TemplateSource, "run-gulp-Vconst-template.bat");
		}
		elsif($Mode eq 'Vconp') {
			$RunScriptTemplate = Deps::MakePath($TemplateSource, "run-gulp-Vconp-template.bat");
		}
		elsif($Mode eq 'Cabc') {
			$RunScriptTemplate = Deps::MakePath($TemplateSource, "run-gulp-Cabc-template.bat");
		}
		else {
			$App->print("\nError: Mode [$Mode] is not implemented for Engine [$Engine].\n\n");
			exit;
		}
	}
	else {
	}
}
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;
		}
	}
	else {
	}
}
elsif($Engine =~ /^vasp$/i) {
}

#===============================================
# 実行
#===============================================
if($Action =~ /^MakeInputs$/i) {
	&MakeInputs();
}
elsif($Action =~ /^MakeRunScript$/i) {
	&MakeRunScript();
}
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;
	my $Emin = 1.0e100;
	while(1) {
		my $p = &ReadTableLine($table);
		last if(!defined $p);

		my $li   = $p->{i1};
		my $P    = $p->{P};
		my $file = $p->{file};
		my $dir  = $p->{dir};

		print "  $idx: $dir: $li, $P, $file\n";
		$P *= 0.001;
		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, $dir, $OutFile);
		$Emin = $Etot if($Emin > $Etot);
		
		$data[$idx] = {
			li    => $li,
			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");
	$App->print("  Emin  : $Emin\n");
	
	my (@x, @y);
	for(my $i = 0 ; $i < $ndata ; $i++) {
#print "i=$i: data=$data[$i]\n";
		$x[$i] = $data[$i]->{V};
		$y[$i] = $data[$i]->{Etot} - $Emin;
	}

	my $iPrintLevel = 2;
#$nOrder=3;
	my $Opt = new Optimize;
#print "nx=", scalar @x, "\n";
#print "ny=", scalar @y, "\n";
#print "nOrder=$nOrder\n";
	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
		$App->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");
		$App->print("  BV(fit)=$BVfit GPa\n");
		$csv->print("$p->{li},$p->{lj},$p->{P},$V,$pl->[0],$pl->[1],$pl->[2],$pl->[3],$pl->[4],$pl->[5],"
			 	   , $p->{Etot} - $Emin, ",$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 $p = &ReadTableLine($table);
		last if(!defined $p);

		my $li   = $p->{i1};
		my $lj   = $p->{j1};
		my $fdV  = $p->{e};
		my $file = $p->{file};
		my $dir  = $p->{dir};

		my $Etot = &ReadEtot($Mode, $Engine, $dir);
		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");
		print("  BV(fit)=$BVfit GPa\n");
		print("  BV(cal)=$BVcal GPa\n") if(defined $BVcal);
		$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 CalculateElasticTensorCabcDiagonal
{
	my ($Crystal, $Mode, $vol, $pdata) = @_;

	$App->print("\n");
	$App->print("Calculate elastic tensor for diagonal element.\n");

	my @data  = @$pdata;
	my $ndata = @data;
	@data = sort { $a->{e} <=> $b->{e}; } @data;
	for(my $i = 0 ; $i < $ndata ; $i++) {
		my $p = $data[$i];
		print("file $p->{file}   e=$p->{e}  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]->{e};
		$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");

	my $csv = JFile->new($SummaryCSV, "w");
	if(!defined $csv) {
		$App->print("\nError in CalculateElasticTensorCabcDiagonal: Can not write to [$SummaryCSV].\n\n");
		exit;
	}

	my $key = "C$Index";
	$csv->print("i1,j1,i2,j2,e1,e2,Etot(eV),Etot(fit),d2E/de2(fit),d2E/de2(numdiff),${key}_fit(GPa),${key}_numdiff(GPa),file\n");
	my ($emin, $emax) = (1.0e100, -1.0e100);
#print "n=$ndata\n";
	for(my $i = 0 ; $i < $ndata ; $i++) {
		my $p = $data[$i];
		my $e = $p->{e};
		$emin = $e if($emin > $e);
		$emax = $e if($emax < $e);
		my $Ecal = $Opt->YCal($e);
		my ($d2Ede2cal, $d2Ede2fit, $Ccal, $Cfit);
		$d2Ede2fit = 2.0 * $pCi->[2];
		for(my $j = 3 ; $j <= $nOrder ; $j++) {
			$d2Ede2fit += $j * ($j-1) * $pCi->[$j] * $e**($j-2);
		}
		$Cfit = $d2Ede2fit * $eVToJ / ($vol * 1.0e-30) * 1.0e-9; # GPa
		if($i >= 1 and $i <= @data-2) {
			$d2Ede2cal = ($data[$i+1]->{Etot} + $data[$i-1]->{Etot} - 2.0*$data[$i]->{Etot}) / ($data[$i+1]->{e} - $data[$i]->{e})**2;
			$Ccal = $d2Ede2cal * $eVToJ / ($vol * 1.0e-30) * 1.0e-9; # GPa
		}
#print "i=$i\n";
		printf("$key for e1=%6.4f / e2=%6.4f: ${key}_fit=%10.4g GPa (%10.4g kBar)  ${key}_cal=%10.4g GPa (%10.4g kBar\n",
						$p->{e}, $p->{e2}, $Cfit, $Cfit * 10.0, $Ccal, $Ccal * 10.0);
#		print("$p->{i1},$p->{j1},$p->{i2},$p->{j2},$p->{ie1},$p->{ie2},$p->{e},$p->{Etot},$Ecal,$d2Ede2fit,$d2Ede2cal,$Cfit,$Ccal,$p->{file}\n");
		$csv->print("$p->{i1},$p->{j1},$p->{i2},$p->{j2},$p->{e},$p->{e2},$p->{Etot},$Ecal,$d2Ede2fit,$d2Ede2cal,$Cfit,$Ccal,$p->{file}\n");
	}
	$csv->print("\n");
	$csv->print("da/a,Etot(fit)\n");
	my $n = 101;
	my $estep = ($emax - $emin) / ($n - 1);
	for(my $i = 0 ; $i < $n ; $i++) {
		my $e = $emin + $i * $estep;
		my $Ecal = $Opt->YCal($e);
		$csv->print("$e,$Ecal\n");
	}
	$csv->Close();

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


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

	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]);
	}
	
	$App->print("Index: ", join(', ', @index), "\n");

	my ($a2, $b2, $c2, $alpha, $beta, $gamma) = $Crystal->LatticeParameters();
	$App->print("  Latt: $a2, $b2, $c2 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 CalculateElasticTensorCabc: Can not read [$TableFile].\n\n");
		exit;
	}
	my @data;
	my $idx = 0;
	my $Emin = 1.0e100;
	my @idata;
	while(1) {
		my $line = $table->ReadLine();
		last if(!defined $line);

		my ($i1, $j1, $i2, $j2, $e, $e2, $file) = Utils::Split("\\s+", $line);
		my $header = $file;
		$header =~ s/\.cif$//i;

		my $Etot = &ReadEtot($Mode, $Engine, $header);
		$Emin = $Etot if($Emin > $Etot);
		$data[$idx] = {
			i1   => $i1 + 0,
			j1   => $j1 + 0,
			i2   => $i2 + 0,
			j2   => $j2 + 0,
			e    => $e + 0.0,
			e2   => $e2 + 0.0,
			Etot => $Etot + 0.0,
			file => $file,
			};
		my ($idx1, $idx2) = ($file =~ /_([+\-\d]+)_([+\-\d]+)/);
		($idx1, $idx2) = (0, 0) if(!defined $idx1 and !defined $idx2);
		$idx1++;
		$idx2++;
		$idata[$idx1][$idx2] = $data[$idx] if($idx1 >= 0 and $idx2 >= 0);
#print "(idx)=($idx1, $idx2)  e=($e,$e2) Etot=$Etot  p=$idata[$idx1][$idx2]\n";
		$idx++;
	}
	$table->Close();

# 対角成分の場合は専用ルーティンに飛ばす
	if($index[0] == $index[2] and $index[1] == $index[3]) {
		my $ret = &CalculateElasticTensorCabcDiagonal($Crystal, $Mode, $vol, \@data);
		return $ret;
	}

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

	$App->print("\n");
	$App->print("Elastic tensor by least squres fitting\n");
	$App->print("  nData : $ndata\n");
#	$App->print("  nOrder: $nOrder\n");
	my (@x, @y, @f);
	for(my $i = 0 ; $i < $ndata ; $i++) {
		$x[$i] = $data[$i]->{e};
		$y[$i] = $data[$i]->{e2};
		$f[$i] = $data[$i]->{Etot} - $Emin;
	}
	my $Opt = new Optimize;
	my $Method = "LCFunction";
	print "  use $Method\n";
	my ($OptVars, $MinVal);
	my $m = 6;
	my $iPrintLevel = 2;
	my (@v);
	for(my $i = 0 ; $i < $ndata ; $i++) {
		$v[$i] = [$data[$i]->{e}, $data[$i]->{e2}];
		$f[$i] = $data[$i]->{Etot} - $Emin;
	}
	($OptVars, $MinVal) = $Opt->Optimize(
					"lsqLCFunctions", \@v, \@f, $m,
					sub { &LCFunc(@_); },
					$iPrintLevel
					);
	print "  Optimized at S = $MinVal:\n";
	$App->print("    a0  = $OptVars->[0]\n");
	$App->print("    ax  = $OptVars->[1]\n");
	$App->print("    ay  = $OptVars->[2]\n");
	$App->print("    axx = $OptVars->[3]\n");
	$App->print("    axy = $OptVars->[4]\n");
	$App->print("    ayy = $OptVars->[5]\n");

	$App->print("\n");
	$App->print("   Writing ftting results to [$SummaryCSV]\n");
	my $csv = JFile->new($SummaryCSV, "w");
	if(!defined $csv) {
		$App->print("\nError in CalculateElasticTensorCabc: Can not write to [$SummaryCSV].\n\n");
		exit;
	}

	my $CxxStr = "C$index[0]$index[1]$index[0]$index[1]";
	my $CyyStr = "C$index[2]$index[3]$index[2]$index[3]";
	my $CxyStr  = "C$index[0]$index[1]$index[2]$index[3]";
	$csv->print("e1,e2,Etot(eV),Etot(fit)\n");
	my ($a0, $ax, $ay, $axx, $axy, $ayy, $axxx, $axxy, $ayyy) = @$OptVars;
	for(my $i = 0 ; $i < $ndata ; $i++) {
		my $x = $data[$i]->{e};
		my $y = $data[$i]->{e2};
		my $Etot = $data[$i]->{Etot} - $Emin;
		my $ycal = &YCal($x, $y, @$OptVars);
		$App->printf("  $i: (e1, e2) = (%6.3f, %6.3f): Etot(DFT)=%16.12g  Etot(fit)=%16.12g eV\n",
				$x, $y, $Etot, $ycal);
		$csv->print("$x,$y,$Etot,$ycal\n");
	}
	$csv->print("\n");
	$App->print("\n");

	my $Cxx = 2.0 * $axx * $eVToJ / ($vol * 1.0e-30) * 1.0e-9; # GPa
	my $Cyy = 2.0 * $ayy * $eVToJ / ($vol * 1.0e-30) * 1.0e-9; # GPa
	my $Cxy =       $axy * $eVToJ / ($vol * 1.0e-30) * 1.0e-9; # GPa
	$csv->print("index,Cfit(GPa)\n");
	$csv->print("$CxxStr,$Cxx\n");
	$App->print("  $CxxStr = $Cxx GPa\n");
	$csv->print("$CyyStr,$Cyy\n");
	$App->print("  $CyyStr = $Cyy GPa\n");
	$csv->print("$CxyStr,$Cxy\n");
	$App->print("  $CxyStr = $Cxy GPa\n");
	$App->print("\n");
	
	$App->print("Elastic tensor by numerical differentiation\n");
	my @p = ($idata[0][1], $idata[1][1], $idata[2][1]);
	my $e = $p[0]->{e};
	my $d2Udex20 = ($p[2]->{Etot} + $p[0]->{Etot} - 2.0*$p[1]->{Etot}) / $e / $e;
	my $Cxx = $d2Udex20 * $eVToJ / ($vol * 1.0e-30) * 1.0e-9; # GPa

	@p = ($idata[1][0], $idata[1][1], $idata[1][2]);
	my $d2Udey20 = ($p[2]->{Etot} + $p[0]->{Etot} - 2.0*$p[1]->{Etot}) / $e / $e;
	my $Cyy = $d2Udey20 * $eVToJ / ($vol * 1.0e-30) * 1.0e-9; # GPa

	@p = ($idata[0][2], $idata[1][2], $idata[2][2]);
	my $dUdexp = ($p[2]->{Etot} - $p[0]->{Etot}) / 2.0 / $e;
	@p = ($idata[0][0], $idata[1][0], $idata[2][0]);
	my $dUdexm = ($p[2]->{Etot} - $p[0]->{Etot}) / 2.0 / $e;
	my $d2Udxdy = ($dUdexp - $dUdexm) / 2.0 / $e;
	my $Cxy = $d2Udxdy * $eVToJ / ($vol * 1.0e-30) * 1.0e-9; # GPa
	$csv->print("index,Cnum(GPa)\n");
	$csv->print("$CxxStr,$Cxx\n");
	$App->print("  $CxxStr = $Cxx GPa\n");
	$csv->print("$CyyStr,$Cyy\n");
	$App->print("  $CyyStr = $Cyy GPa\n");
	$csv->print("$CxyStr,$Cxy\n");
	$App->print("  $CxyStr = $Cxy GPa\n");
	$csv->Close();
	
	$App->print("\n");
}

sub LCFunc
{
	my ($idx, $v) = @_;
	if($idx == 0) {
		return 1.0;
	}
	elsif($idx == 1) {
		return $v->[0];
	}
	elsif($idx == 2) {
		return $v->[1];
	}
	elsif($idx == 3) {
		return $v->[0]*$v->[0];
	}
	elsif($idx == 4) {
		return $v->[0]*$v->[1];
	}
	elsif($idx == 5) {
		return $v->[1]*$v->[1];
	}

print "invalid idx=$idx\n";
exit;
}

sub YCal
{
	my ($x, $y, $a0, $ax, $ay, $axx, $axy, $ayy) = @_;
	return $a0 + $ax * $x +$ay * $y + $axx * $x*$x + $axy * $x*$y + $ayy * $y*$y;
}

sub LSQ2DFunc
{
	my ($px, $py, $pf, $pVars, $iPrintLevel) = @_;
#	my ($a0, $ax, $ay, $axx, $axy, $ayy) = @$pVars;
	my $nData = @$px;
	my $S = 0.0;
	for(my $i = 0 ; $i < $nData ; $i++) {
		my $ycal = &YCal($px->[$i], $py->[$i], @$pVars);
		my $d = $pf->[$i] - $ycal;
		$S +=  $d*$d;
	}
	$S = sqrt($S / $nData);
print("a=", join(', ', @$pVars), "  S=$S\n") if($iPrintLevel >= 2);
	return $S;
}

sub PDLDisplayFunc {
	my ($piddle) = @_;
#print "Display: piddle=$piddle\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 ($fStr, $rStr, $zStr);
	my ($ir, $if, $iz) = (0, 0, 0);
	while(1) {
		my $line = $table->ReadLine();
		last if(!defined $line);

		my ($i1, $j1, $i2, $j2, $e, $e2, $file) = Utils::Split("\\s+", $line);
		$e += 0.0;
		my $dir = $file;
		$dir =~ s/\.\w+?$//;
print "$i1: $file\n";		

		if($e < 0.0) {
			$App->print(" Negative strain: e=$e atm from [$file].\n");
			$rStr .= "set eArrayR[$ir]=$e\n";
			$rStr .= "set cifArrayR[$ir]=$file\n";
			$rStr .= "set dirArrayR[$ir]=$dir\n";
			$ir++;
		}
		elsif($e > 0.0) {
			$App->print(" Positive strain: e=$e atm from [$file].\n");
			$fStr .= "set eArrayF[$if]=$e\n";
			$fStr .= "set cifArrayF[$if]=$file\n";
			$fStr .= "set dirArrayF[$if]=$dir\n";
			$if++;
		}
		else {
			$App->print(" Zero strain: e=$e atm from [$file].\n");
			$zStr .= "set eArrayZ[$iz]=$e\n";
			$zStr .= "set cifArrayZ[$iz]=$file\n";
			$zStr .= "set dirArrayZ[$iz]=$dir\n";
			$iz++;
		}
	}
	$table->Close();

	my $ArrayString = "\@echo on\n"
					 . "set nr=$ir\n" . $rStr 
					 . "set nz=$iz\n" . $zStr 
					 . "set nf=$if\n" . $fStr 
					 . "\@echo off\n";
	$str =~ s/\{ArrayString\}/$ArrayString/i;
	$str =~ s/\{SampleName\}/$SampleName/i;
#$v{$key} = Template->new()->ReplaceByHash($var, \%v, '{', '}');

	$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 ($fStr, $rStr, $zStr);
	my ($ir, $if, $iz) = (0, 0, 0);
	while(1) {
		my $p = &ReadTableLine($table);
		last if(!defined $p);

		my $li  = $p->{i1};
		my $P   = $p->{P};    # in atm
		my $e   = $P * 0.001; # in kBar
		my $file = $p->{file};
		my $dir  = $p->{dir};
		if($P < 0.0) {
			$App->print(" Negative pressure: P=$P atm from [$file].\n");
			$rStr .= "set eArrayR[$ir]=$e\n";
			$rStr .= "set cifArrayR[$ir]=$file\n";
			$rStr .= "set dirArrayR[$ir]=$dir\n";
			$ir++;
		}
		elsif($P > 0.0) {
			$App->print(" Positive pressure: P=$P atm from [$file].\n");
			$fStr .= "set eArrayF[$if]=$e\n";
			$fStr .= "set cifArrayF[$if]=$file\n";
			$fStr .= "set dirArrayF[$if]=$dir\n";
			$if++;
		}
		else {
			$App->print(" Zero pressure: P=$P atm from [$file].\n");
			$zStr .= "set eArrayZ[$iz]=$e\n";
			$zStr .= "set cifArrayZ[$iz]=$file\n";
			$zStr .= "set dirArrayZ[$iz]=$dir\n";
			$iz++;
		}
	}
	$table->Close();

	my $ArrayString = "\@echo on\n"
					 . "set nr=$ir\n" . $rStr 
					 . "set nz=$iz\n" . $zStr 
					 . "set nf=$if\n" . $fStr 
					 . "\@echo off\n";
	$str =~ s/\{ArrayString\}/$ArrayString/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;
#$v{$key} = Template->new()->ReplaceByHash($var, \%v, '{', '}');

	$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 ($fStr, $rStr, $zStr);
	my ($ir, $if, $iz) = (0, 0, 0);
	while(1) {
		my $p = &ReadTableLine($table);
		last if(!defined $p);

		my $e    = $p->{e};
		my $file = $p->{file};
		my $dir  = $p->{dir};

		if($e < 0.0) {
			$App->print(" Shrunk structure: dV/V==$e from [$file].\n");
			$rStr .= "set eArrayR[$ir]=$e\n";
			$rStr .= "set cifArrayR[$ir]=$file\n";
			$rStr .= "set dirArrayR[$ir]=$dir\n";
			$ir++;
		}
		elsif($e > 0.0) {
			$App->print(" Expanded structure: dV/V==$e from [$file].\n");
			$fStr .= "set eArrayF[$if]=$e\n";
			$fStr .= "set cifArrayF[$if]=$file\n";
			$fStr .= "set dirArrayF[$if]=$dir\n";
			$if++;
		}
		else {
			$App->print(" Original structure: dV/V==$e from [$file].\n");
			$zStr .= "set eArrayZ[$iz]=$e\n";
			$zStr .= "set cifArrayZ[$iz]=$file\n";
			$zStr .= "set dirArrayZ[$iz]=$dir\n";
			$iz++;
		}
	}
	$table->Close();

	my $ArrayString = "\@echo on\n"
					 . "set nr=$ir\n" . $rStr 
					 . "set nz=$iz\n" . $zStr 
					 . "set nf=$if\n" . $fStr 
					 . "\@echo off\n";
	$str =~ s/\{ArrayString\}/$ArrayString/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;
#$v{$key} = Template->new()->ReplaceByHash($var, \%v, '{', '}');

	$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 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 %10.6f %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);

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

	$table->Close();
}

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

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

	if(defined $e2) {
		$v[$li2][0] += 0.5 * $e2 * $v[$lj2][0] * $latt[$li2] / $latt[$lj2];
		$v[$li2][1] += 0.5 * $e2 * $v[$lj2][1] * $latt[$li2] / $latt[$lj2];
		$v[$li2][2] += 0.5 * $e2 * $v[$lj2][2] * $latt[$li2] / $latt[$lj2];
		$v[$lj2][0] += 0.5 * $e2 * $v[$li2][0] * $latt[$lj2] / $latt[$li2];
		$v[$lj2][1] += 0.5 * $e2 * $v[$li2][1] * $latt[$lj2] / $latt[$li2];
		$v[$lj2][2] += 0.5 * $e2 * $v[$li2][2] * $latt[$lj2] / $latt[$li2];
	}

	$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 'Vconst') {
#		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');
	$App->print("  0: Ptarget=0.0: $OriginalCIF in [original]\n");
#	$table->printf("%2d %2d %10.6f %s\n", 0, 0, $pressure0, $OriginalCIF);
	$pressure0 += 0.0;
	&WriteTableLine($table,
			idx1 => 0,
			P    => $pressure0,
			file => $OriginalCIF,
			dir  => "P0", # "original",
			);

	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;
		$Ptarget += 0.0;
		my $dir = "P" . ($Ptarget * 0.001);
		$App->print("  $i: Ptarget=$Ptarget: $OriginalCIF in [$dir]\n");
#		$table->printf("%2d %2d %10.6f %s\n", $i, 0, $Ptarget, $OriginalCIF);
		&WriteTableLine($table,
				idx1 => $i,
				P    => $Ptarget,
				file => $OriginalCIF,
				dir  => $dir,
				);
	}
	for(my $i = 1 ; $i < $n + $nStep ; $i++) {
#print "i=$i\n";
		my $Ptarget = $pressure0 + $i * $fdiff;
		$Ptarget += 0.0;
		my $dir = "P" . ($Ptarget * 0.001);
		$App->print("  $i: Ptarget=$Ptarget: $OriginalCIF in [$dir]\n");
#		$table->printf("%2d %2d %10.6f %s\n", $i, 0, $Ptarget, $OriginalCIF);
		&WriteTableLine($table,
				idx1 => $i,
				P    => $Ptarget,
				file => $OriginalCIF,
				dir  => $dir,
				);
	}

	$table->Close();
}

sub WriteTableLine
{
	my ($out, %p) = @_;

	if(!defined $p{dir}) {
		$p{dir} = $p{file};
		$p{dir} =~ s/\.\w+?$//;
	}
	$p{e} = $p{P} if(!defined $p{e});

	$out->print("$p{dir}, $p{file}, $p{idx1}, $p{idx2}, "
			   ."$p{i1}, $p{j1}, $p{i2}, $p{j2}, $p{e}, $p{e2}\n");
}

sub ReadTableLine
{
	my ($in) = @_;
	
	my $line = $in->ReadLine();
	return undef if(!defined $line);
	my ($dir, $file, $idx1, $idx2, $i1, $j1, $i2, $j2, $e, $e2, @a) 
		= Utils::Split("\\s*,\\s*", $line, 0);
#print "$dir, $file, $idx1, $idx2, $i1, $j1, $i2, $j2, $e, $e2:", 
#		join(':', @a), "\n";
	my $p = {
		dir  => $dir,
		file => $file,
		i1   => $i1 + 0,
		j1   => $j1 + 0,
		i2   => $i2 + 0,
		j2   => $j2 + 0,
		idx1 => $idx1,
		idx2 => $idx2,
		P    => $e,
		e    => $e + 0.0,
		e2   => $e2 + 0.0,
		};
	return $p;
}

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);
	&WriteTableLine($table,
			idx1 => 0,
			e    => 0.0,
			file => $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";
		&SaveStrainedStructureV($Crystal, $i, 0, $k, $OutCIF);
#		$table->printf("%2d %2d %10.6f %s\n", $i, 0, $dV, $OutCIF);
		&WriteTableLine($table,
				idx1 => $i,
				e    => $dV,
				file => $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";
		&SaveStrainedStructureV($Crystal, $i, 0, $k, $OutCIF);
#		$table->printf("%2d %2d %10.6f %s\n", $i, 0, $dV, $OutCIF);
		&WriteTableLine($table,
				idx1 => $i,
				e    => $dV,
				file => $OutCIF,
				);
	}

	$table->Close();
}

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

	my ($Etot);
	my $dir = $file;
	$dir =~ s/\.cif$//i;
	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 SaveCIF
{
	my ($Crystal, $OutCIF, $IsChooseRandomly, $strCRLF) = @_;

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