#!/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 Sci qw($torad $eVToJ $RyToeV);
use Sci::ElasticConstant;
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 $ec = new ElasticConstant;

my $ConfigFile = 'configure.ini';
my $vars = new ProgVars;
my $pParams = {};
$pParams = $vars->ReadConfigure($ConfigFile, $pParams, 1, 1);

my $OS         = Utils::GetOSName();
my $SGROUPPath = ProgVars::SGROUPPath();
my $TEEPath    = ProgVars::TEEPath();

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

# eiJを計算するための格子定数の歪率
my $TableFile           = "table.txt";
my $SummaryCSV          = "summary.csv";
my $InitScriptPath = 'initialize.bat';
my $RunScriptPath  = 'CalculateElasticConstant.bat';

my $StrainedCIFTemplate = 'e{i}{j}-{sign}.cif';

#===============================================
# 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("--OriginalDir",       "--OriginalDir=path: Default original",       'original');
$App->AddArgument("--OriginalCIF",       "--OriginalCIF=path: Default original.cif",   'original.cif');
$App->AddArgument("--OriginalInput",     "--OriginalInput=path: Default original.cif", 'original.cif');
$App->AddArgument("--BurstToP1",         "--BurstToP1=[0|1]: Default: 1",               1);
$App->AddArgument("--UseBravaisLattice", "--UseBravaisLattice=[0|1]: Default: 1",       1);
$App->AddArgument("--Mode",              "--Mode=[Cabc|CXYZ|Vconst|Vconp|init]: Default Vconst", 'Vconst');
$App->AddArgument("--Index",             "--Index=[ijkl]: Default 1111",                  '1111');
$App->AddArgument("--nStep",             "--nStep=[integer >= 2]: Default 3",             3);
$App->AddArgument("--nRelaxationStep",   "--nRelaxationStep=[integer]: Default 100",      100);
$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("--e_tol",             "--e_tol=val     EPS for Etot convergence  : Default 1.0", 0.5e-4);
$App->AddArgument("--F_tol",             "--F_tol=val     EPS for force convergence : Default 1.0", 0.5e-3);
$App->AddArgument("--conv_thr",          "--conv_thr=val  EPS for SCF convergence   : Default 1.0", 0.5e-6);
$App->AddArgument("--pressure_tol",      "--pressure_tol=val in atm. EPS for P 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";

$App->SetDefault();

my $pvar = {};
%$pvar = $App->GetArgHash(1);
Utils::MergeHash($pvar, $pParams);
#foreach my $key (sort keys %$pvar) {
#print "  $key: $pvar->{$key}\n";
#}
#exit;

$ec->SetApplication($App);
$ec->SetpVar($pvar);

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

$App->AddVar('OS', $OS);
$App->AddVar('CIFFile', $App->GetGetArg(0));
$App->AddVar('SampleName', $SampleName);

#my $Action        = $App->GetGetArg('Action',        undef, 1, 1);
my $Engine        = $App->GetGetArg('Engine',        undef, 1, 1);
my $Mode          = $App->GetGetArg('Mode',          undef, 1, 1);

$ec->SetEngine($Engine);
$ec->SetMode($Mode);
$ec->SetTablePath($TableFile);
my $pEngine = $ec->SetpEngine();

my $Index         = $App->GetGetArg('Index',         undef, 1, 1);
my $OriginalDir   = $App->GetGetArg('OriginalDir',   undef, 1, 1);
my $OriginalCIF   = $App->GetGetArg('OriginalCIF',   undef, 1, 1);
my $OriginalInput = $App->GetGetArg('OriginalInput', undef, 1, 1);
my $BurstToP1     = $App->GetGetArg('BurstToP1', undef, 1, 1);
my $UseBravaisLattice = $App->GetGetArg('UseBravaisLattice', undef, 1, 1);
my $nStep         = $App->GetGetArg('nStep',         undef, 1, 1);
my $nRelaxationStep = $App->GetGetArg('nRelaxationStep',         undef, 1, 1);
my $nOrder        = $App->GetGetArg('nOrder',        undef, 1, 1);
my $pressure0     = $App->GetGetArg('pressure0',     undef,  1, 1);
my $e_tol         = $App->GetGetArg('e_tol',         undef,  1, 1);
my $F_tol         = $App->GetGetArg('F_tol',         undef,  1, 1);
my $conv_thr      = $App->GetGetArg('conv_thr',      undef,  1, 1);
my $pressure_tol  = $App->GetGetArg('pressure_tol',  undef,  1, 1);
my $fdiff         = $App->GetGetArg('fdiff',         undef,  1, 1);

$ec->print("\n");
$App->PrintArgs("%-30s: %s\n", 
	[qw(OS CIFFile Action Engine Mode Index OriginalCIF OriginalInput OriginalDir
	    BurstToP1 UseBravaisLattice Library
	    fdiff nStep nOrder nRelaxationStep pressure0
	    conv_thr e_tol F_tol pressure_tol), sort $App->ArgKeys()]);

my ($RunScriptTemplate) = $ec->GetFileConfig($TemplateSource);

#===============================================
# 実行
#===============================================
if($pvar->{Action} =~ /^MakeInputs$/i) {
	&MakeInputs();
}
elsif($pvar->{Action} =~ /^MakeRunScript$/i) {
	&MakeRunScript();
}
elsif($pvar->{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("Original CIF File: $OriginalCIF\n");
	my $CIF = $App->{'CIF'} = new CIF;
	unless($CIF->Read($OriginalCIF, 0)) {
		$App->print("Error: Can not read [$OriginalCIF]\n\n");
		exit 2;
	}
	my $Crystal = $CIF->GetCCrystal();

	if($Mode =~ /^Cabc/) {
		&CalculateElasticTensorCabc($Crystal, $Mode);
		return;
	}
	elsif($Mode =~ /^CXYZ/) {
		&CalculateElasticTensorCabc($Crystal, $Mode);
		return;
	}
	elsif($Mode =~ /^Vconst/) {
		&CalculateElasticTensorVconst($Crystal, $Mode);
		return;
	}
	elsif($Mode =~ /^Vconp/) {
		&CalculateElasticTensorVconp($Crystal, $Mode);
		return;
	}
	else {
		$App->InvalidParameter('Mode', 'CalculateElasticTensor', 1);
		exit;
	}
}

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

	my $table = $ec->OpenTable(undef, '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 = $ec->ReadTableLine();
		last if(!defined $p);

		my $li   = $p->{idx1};
		my $lj   = $p->{idx2};
		my $P    = $p->{P};
		my $file = $p->{file};
		my $dir  = $p->{dir};

		print "  data $idx in dir [$dir]: idx=$li  P=$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 $pOut     = $PWSCF->ReadOutputToHash($OutFile, 0);
		if($pOut->{RelaxationStatus} =~ /Error/i) {
			$App->print("\n");
			$App->print("***********************************************\n");
			$App->print("***Relaxation $pOut->{RelaxationStatus}\n");
			$App->print("***********************************************\n");
			$App->print("\n");
			next;
		}
		if($pOut->{SCFStatus} =~ /Error/i) {
			$App->print("\n");
			$App->print("***********************************************\n");
			$App->print("***SCF $pOut->{SCFStatus}\n");
			$App->print("***********************************************\n");
			$App->print("\n");
			next;
		}

		my $nStep    = $PWSCF->GetnFrameInOutput($App, $OutFile);
		my $Crystal2 = $PWSCF->ReadCrystalStructureFromInput($App, $InputFile, $Crystal, 0);
		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, 0));
		}
		$in->Close();

		my @latt = $Crystal2->LatticeParameters();
		my $V    = $Crystal2->Volume();
		my $Etot = $ec->ReadEtot($dir, $OutFile);
		$Emin = $Etot if($Emin > $Etot);
		
		$data[$idx] = {
			li    => $li,
			lj    => $lj,
			P     => $P,
			Pcal  => $pOut->{Pkbar},
			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++;
	}
	$ec->CloseTable();
	$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");
	if($ndata < $nOrder+1) {
		$App->print("  Warning: ndata=$ndata is smaller than nOrder+1=", $nOrder+1, "\n");
		$nOrder = $ndata - 1;
		$App->print("      nOrder is redured to $nOrder\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(set),P(cal),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(set),P(cal),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},$p->{Pcal},$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},$p->{Pcal},$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 = $ec->OpenTable(undef, '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 = $ec->ReadTableLine();
		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 = $ec->ReadEtot($dir);
		my $V = $vol * (1.0 + $fdV);
		$data[$idx] = {
			li   => $li,
			lj   => $lj,
			fdV  => $fdV,
			V    => $V,
			Etot => $Etot,
			file => $file,
			};
		$idx++;
	}
	$ec->CloseTable();

	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();

	$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 = $ec->OpenTable(undef, 'r');
	if(!defined $table) {
		$App->print("\nError in CalculateElasticTensorVconst: Can not read [$TableFile].\n\n");
		exit;
	}

	my @data;
	my $idx = 0;
	my $Emin = 1.0e100;
	my @idata;
	while(1) {
		my $p = $ec->ReadTableLine();
		last if(!defined $p);

		my $idx1 = $p->{idx1};
		my $idx2 = $p->{idx2};
		my $i1   = $p->{i1};
		my $i2   = $p->{i1};
		my $j1   = $p->{i2};
		my $j2   = $p->{j2};
		my $e    = $p->{e};
		my $e2   = $p->{e2};
		my $file = $p->{file};
		my $dir  = $p->{dir};

		my $Etot = $ec->ReadEtot($dir);
		$Emin = $Etot if($Emin > $Etot);
		$data[$idx] = {
			idx1 => $idx1 + 0,
			idx2 => $idx2 + 0,
			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,
			dir  => $dir,
			};
		$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();
#$Emin = 0.0;
	print "  Emin = $Emin eV\n";
	for(my $i = 0 ; $i < @data ; $i++) {
		$data[$i]{Etot} -= $Emin;
	}

# 対角成分の場合は専用ルーティンに飛ばす
	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;
	$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;
	$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;
	$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 =~ /^Cabc/) {
#		&MakeRunScriptCabc();
		$ec->MakeRunScriptCXYZ($CIFFile, $SampleName, $TableFile, $RunScriptTemplate, $RunScriptPath);
		return;
	}
	elsif($Mode =~ /^CXYZ/) {
		$ec->MakeRunScriptCXYZ($CIFFile, $SampleName, $TableFile, $RunScriptTemplate, $RunScriptPath);
		return;
	}
	elsif($Mode =~ /^Vconst/) {
		&MakeRunScriptVconst();
		return;
	}
	elsif($Mode =~ /^Vconp/) {
		&MakeRunScriptVconp();
		return;
	}
	elsif($Mode eq 'init') {
		&MakeRunScriptInit();
		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 = $ec->OpenTable(undef, 'r');
	if(!defined $table) {
		$App->print("\nError in MakeRunScriptCabc: Can not read [$TableFile].\n\n");
		exit;
	}

	my ($fStr, $rStr, $zStr);
	my ($ir, $if, $iz) = (0, 0, 0);
	while(1) {
		my $p = $ec->ReadTableLine();
		last if(!defined $p);

		my $i1 = $p->{i1};
		my $j1 = $p->{j1};
		my $i2 = $p->{i2};
		my $j2 = $p->{j2};
		my $e  = $p->{e};
		my $e2 = $p->{e2};
		my $file = $p->{file};
		my $dir  = $p->{dir};
		$e  += 0.0;
		$e2 += 0.0;
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 
					 . "\n"; # \@echo off\n";
#	$str =~ s/\{ArrayString\}/$ArrayString/i;
#	$str =~ s/\{SampleName\}/$SampleName/i;
#$v{$key} = Template->new()->ReplaceByHash($var, \%v, '{', '}');

	my %var = (
		Engine       => $Engine,
		Mode         => $Mode,
		Function     => $ec->RelaxMode(),
		SampleName   => $SampleName,
		CIFFile      => $CIFFile,
		UseBravaisLattice => $UseBravaisLattice,
		BurstToP1    => $BurstToP1,
		OriginalDir  => $OriginalDir,
		OriginalCIF  => $OriginalCIF,
		ArrayString  => $ArrayString,
		pressure     => $pressure0 * 0.001,
		nRelaxationStep => $nRelaxationStep,
		e_tol        => $e_tol,
		F_tol        => $F_tol,
		conv_thr     => $conv_thr,
		pressure_tol => $pressure_tol * 0.001,
		);
	$str = Template->new()->ReplaceByHash($str, \%var, '{', '}');

	$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 MakeRunScriptInit
{
	my $table = $ec->OpenTable(undef, 'w');
	if(!defined $table) {
		$App->print("\nError in MakeRunScriptInit: Can not write to [$TableFile].\n\n");
		exit;
	}
	$ec->WriteTableLine(
			idx1 => 0,
			P    => $pressure0 + 0.0,
			file => $OriginalCIF,
			dir  => $OriginalDir,
			);
	$ec->CloseTable();

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

#	$str =~ s/\{SampleName\}/$SampleName/i;

	my $ArrayString = "" #\@echo on\n"
					 . "set nr=0\n"
					 . "set nz=1\n"
					 . "set eArrayZ[0]=0\n"
					 . "set cifArrayZ[0]=original.cif\n"
					 . "set dirArrayZ[0]=original\n"
					 . "set nf=0\n"
					 . "\n"; # \@echo off\n";
#	$str =~ s/\{ArrayString\}/$ArrayString/i;
#	$str =~ s/\{SampleName\}/$SampleName/i;
#$v{$key} = Template->new()->ReplaceByHash($var, \%v, '{', '}');

	my %var = (
		Engine       => $Engine,
		Mode         => $Mode,
		Function     => $ec->RelaxMode(),
		SampleName   => $SampleName,
		CIFFile      => $CIFFile,
		UseBravaisLattice => $UseBravaisLattice,
		BurstToP1    => $BurstToP1,
		OriginalDir  => $OriginalDir,
		OriginalCIF  => $OriginalCIF,
		ArrayString  => $ArrayString,
		pressure     => $pressure0 * 0.001,
		nRelaxationStep => $nRelaxationStep,
		e_tol        => $e_tol,
		F_tol        => $F_tol,
		conv_thr     => $conv_thr,
		pressure_tol => $pressure_tol * 0.001,
		);
	$str = Template->new()->ReplaceByHash($str, \%var, '{', '}');

	$App->print("Save to [$InitScriptPath].\n");
	my $ret = JFile->SaveFile($str, $InitScriptPath, "w");
	if(!defined $ret) {
		$App->print("\nError in MakeRunScriptInit: Can not write to [$InitScriptPath].\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 = $ec->OpenTable(undef, '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 = $ec->ReadTableLine();
		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 
					 . "\n"; #\@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, '{', '}');

	my %var = (
		Engine       => $Engine,
		Mode         => $Mode,
		Function     => $ec->RelaxMode(),
		SampleName   => $SampleName,
		CIFFile      => $CIFFile,
		UseBravaisLattice => $UseBravaisLattice,
		BurstToP1    => $BurstToP1,
		OriginalDir  => $OriginalDir,
		OriginalCIF  => $OriginalCIF,
		ArrayString  => $ArrayString,
		pressure     => $pressure0 * 0.001,
		nRelaxationStep => $nRelaxationStep,
		e_tol        => $e_tol,
		F_tol        => $F_tol,
		conv_thr     => $conv_thr,
		pressure_tol => $pressure_tol * 0.001,
		);
	$str = Template->new()->ReplaceByHash($str, \%var, '{', '}');

	$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 = $ec->OpenTable(undef, '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 = $ec->ReadTableLine();
		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++;
		}
	}
	$ec->CloseTable();

	my $ArrayString = "" # \@echo on\n"
					 . "set nr=$ir\n" . $rStr 
					 . "set nz=$iz\n" . $zStr 
					 . "set nf=$if\n" . $fStr 
					 . "\n"; #\@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, '{', '}');

	my %var = (
		Engine        => $Engine,
		Mode          => $Mode,
		Function     => $ec->RelaxMode(),
		SampleName    => $SampleName,
		CIFFile       => $CIFFile,
		UseBravaisLattice => $UseBravaisLattice,
		BurstToP1     => $BurstToP1,
		OriginalDir   => $OriginalDir,
		OriginalCIF   => $OriginalCIF,
		OriginalInput => $OriginalInput,
		ArrayString   => $ArrayString,
		pressure      => $pressure0 * 0.001,
		nRelaxationStep => $nRelaxationStep,
		e_tol         => $e_tol,
		F_tol         => $F_tol,
		conv_thr      => $conv_thr,
		pressure_tol  => $pressure_tol * 0.001,
		);
	$str = Template->new()->ReplaceByHash($str, \%var, '{', '}');

	$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;
	}
}

#export "WIEN2kplCMD=$PerlPATH $TkPerlDir/WIEN2k/WIEN2k.pl"
#MakeSymmetrizedCIFCMD="$WIEN2kplCMD --Action=MakeSymmetrizedStructFile --TOL=0.001 *-final.cif out.struct"

sub MakeSymmetrizedStructFile
{
	my ($App, $CIFile, $StructFile, $TOL, $Prim) = @_;
	$TOL  = 0.0001 if(!defined $TOL or $TOL <= 0.0);
	$Prim = 0 if(!defined $Prim);

	$App->print("\n\n<b>Make Symmetrized WIEN2k Struct File:</b>\n");
	$App->print("Torrelance: $TOL\n");

#ファイル名を（ベース名, ディレクトリ名, 拡張子）に分解
	my ($drive, $dir, $filename, $ext, $lastdir, $filebody)
		= Deps::SplitFilePath($CIFFile);
	my $SampleName = $filebody;
	my $BasePath = "$dir";
	$BasePath = "$drive$dir" if($drive ne '');

	my $StructPath    = Deps::MakePath($BasePath, "${filebody}.struct");
	my $SymStructPath = Deps::MakePath($BasePath, "${filebody}-Symmetrized.struct");
	my $SymCIFPath    = Deps::MakePath($BasePath, "${filebody}-Symmetrized.cif");
	$App->print("Original(non-symmetric) CIF   : $CIFFile\n");
	$App->print("Original(non-symmetric) struct: $StructPath\n");
	$App->print("Symmetrized struct            : $SymStructPath\n");
	$App->print("Symmetrized CIF               : $SymCIFPath\n");

	my $CIF = new CIF();
	unless($CIF->Read($CIFFile)) {
		$App->print("Error: Can not read $CIFFile.\n\n");
		return 0;
	}
	my $Crystal = $CIF->GetCCrystal();
	$Crystal->ExpandCoordinates();

	$App->print("Save to $StructPath.\n");
	my $WIEN = new WIEN2k;
	$WIEN->SetSampleName($SampleName);
	$WIEN->SaveStructFile($Crystal, $StructPath, 0, 0, 0);

#	$App->print("instgen\n");
#	system("instgen");
	my $SGROUPCommand = "$SGROUPPath -set-TOL=$TOL -wi $StructPath";
	$SGROUPCommand    = "$SGROUPPath -prim -set-TOL=$TOL -wi $StructPath" if($Prim);
	$App->print("$SGROUPCommand\n");
	system("$SGROUPCommand | $TEEPath $StructPath.out");

	$SGROUPCommand = "$SGROUPPath -set-TOL=$TOL -wi $StructPath -wo $SymStructPath";
	$SGROUPCommand = "$SGROUPPath -prim -set-TOL=$TOL -wi $StructPath -wo $SymStructPath" if($Prim);
	$App->print("$SGROUPCommand\n");
	system("$SGROUPCommand");

	my $WIEN2k2 = new WIEN2k;
	my $SymCrystal = $WIEN2k2->ReadStructFile($SymStructPath);
	unless($SymCrystal) {
		$App->print("Error: Can not read [$SymStructPath].\n");
		return 0;
	}
	
	my $SPG = $SymCrystal->GetCSpaceGroup();
	my $SPGName = $SPG->SPGName();
	if($SPGName =~ /^R/i or $SPGName =~ /^P\s*\-?3/i) {
		$SPG->SetSPGName("$SPGName R");
		my ($a,$b,$c,$alpha,$beta,$gamma) = $SymCrystal->LatticeParameters();
		$App->print("  Hexagonal axis: a=$a   c=$c\n");
		my $ar = sqrt(3.0*$a*$a + $c*$c) / 3.0;
		my $alphah = 1.5 / sqrt(3.0 + ($c/$a)*($c/$a));
		$alphah = 2.0 * Sci::asin($alphah) * Sci::todeg();
		$App->print("  Rhombohedral axis: a=$ar   alpha=$alphah\n");
		$SymCrystal->SetLatticeParameters($ar,$ar,$ar,$alphah,$alphah,$alphah);
	}

	my $CIF2 = new CIF;
	if(!$CIF2->CreateCIFFileFromCCrystal($SymCrystal, $SymCIFPath, 0, 'unix')) {
		$App->print("Error: Can not create [$SymCIFPath].\n");
		return 0;
	}

	$App->print("\nMake Symmetrized Struct file Finished.\n\n");

	return;
}

#\Programs\Perl\Crystal\SPG\FindSPG.pl
#$App->AddArgument("--TolD",                "--TolD=val: Torelance for atomic position (A) [Def 0.01 A]", "0.01");
#$App->AddArgument("--xc0",                 "--xc0=val [Def 0.0]",                 "0.0");
#$App->AddArgument("--yc0",                 "--yc0=val [Def 0.0]",                 "0.0");
#$App->AddArgument("--zc0",                 "--zc0=val [Def 0.0]",                 "0.0");
#$App->AddArgument("--ScanAll",             "--ScanAll=[0|1] [Def 0]",             "0");
#$App->AddArgument("--FindSiteSymmetry",    "--FindSiteSymmetry=[0|1] [Def 0]",    "0");
#$App->AddArgument("--ShowCrystalInf",      "--ShowCrystalInf=[0|1] [Def 0]",      "0");
#$App->AddArgument("--ShowAllCoordination", "--ShowAllCoordination=[0|1] [Def 0]", "0");
#$App->AddArgument("--RCoordMax",           "--RCoordMax=val [Def 7.0]",           "7.0");
#$App->AddArgument("--ShowFoundStructure",  "--ShowFoundStructure=[0|1] [Def 0]",  "0");

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

	my $Crystal;

	my $RelaxedInput  = $ec->GetInputPath($OriginalDir, $OriginalCIF);
	my $RelaxedOutput = $ec->GetOutputPath($OriginalDir, $OriginalCIF);
	if($RelaxedInput and $RelaxedOutput) {
		$App->print("Read original structure from output File [$RelaxedOutput]\n");
		$Crystal = $ec->ReadFinalCrystalStructureFromOutput($App, $RelaxedInput, $RelaxedOutput, 0);
		if(!$Crystal) {
			$App->print("Error: Can not read [$RelaxedInput] or [$RelaxedOutput]\n\n");
#			exit 2;
		}
	}
	if(!$Crystal) {
		$App->print("Read original structure from 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;
		}
		$Crystal = $CIF->GetCCrystal();
	}
	$App->print("\n");

#$Crystal->SplitPartialSites();
	$Crystal->ExpandCoordinates();
	my $CrystalName = $Crystal->CrystalName();
	$App->print("CrystalName: $CrystalName\n");

	if($Mode =~ /^Cabc/) {
		$ec->MakeInputsCXYZ($Crystal);
#		&MakeInputsCabc($Crystal);
		return;
	}
	elsif($Mode =~ /^CXYZ/) {
		$ec->MakeInputsCXYZ($Crystal);
		return;
	}
	elsif($Mode =~ /^Vconst/){
		&MakeInputsVconst($Crystal);
		return;
	}
	elsif($Mode =~ /^Vconp/) {
		&MakeInputsVconp($Crystal);
		return;
	}
	elsif($Mode eq 'init') {
		&MakeInputsInit($Crystal);
		return;
	}
	else {
		$App->InvalidParameter('Mode', 'MakeInputs', 1);
		exit;
	}
}

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

	$ec->SaveCIF($Crystal, $OriginalCIF, 0, 'unix');

	my $table = $ec->OpenTable(undef, 'w');
	$ec->WriteTableLine(
			idx1 => 0,
			P    => $pressure0 + 0.0,
			e    => 0.0,
			file => $OriginalCIF,
			dir  => $OriginalDir,
			);
	$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 $dir = $StrainedCIFTemplate;
		if($ec->RelaxMode() eq 'scf') {
			$dir = "dV" . sprintf("%+05d", int($dV*10000));
		}
		else {
			$dir = "dVr" . sprintf("%+05d", int($dV*10000));
		}
		my $OutCIF = "$dir.cif";
		print "Save strained CIF to [$OutCIF]\n";
		$ec->SaveStrainedStructureV($Crystal, $i, 0, $k, $OutCIF);
		$ec->WriteTableLine(
				idx1 => $i,
				e    => $dV,
				file => $OutCIF,
				dir  => $dir,
				);
	}
	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 $dir = $StrainedCIFTemplate;
		if($ec->RelaxMode() eq 'scf') {
			$dir = "dV" . sprintf("%+05d", int($dV*10000));
		}
		else {
			$dir = "dVr" . sprintf("%+05d", int($dV*10000));
		}
		my $OutCIF = "$dir.cif";
		print "Save strained CIF to [$OutCIF]\n";
		$ec->SaveStrainedStructureV($Crystal, $i, 0, $k, $OutCIF);
		$ec->WriteTableLine(
				idx1 => $i,
				e    => $dV,
				file => $OutCIF,
				dir  => $dir,
				);
	}

	$ec->CloseTable();
}

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

	print "Save original CIF to [$OriginalCIF]\n";
	my $NewCrystal = $Crystal->Duplicate();
	$NewCrystal->BurstToP1();
	my @latt = $NewCrystal->LatticeParameters();
	my @lv   = $NewCrystal->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);
	$App->print("\n");

	$ec->SaveCIF($NewCrystal, $OriginalCIF, 0, 'unix');
	my $table = $ec->OpenTable(undef, 'w');
	$ec->WriteTableLine(
			idx1 => 0,
			idx2 => 0,
			i1   => $index[0],
			j1   => $index[1],
			i2   => $index[2],
			j2   => $index[3],
			P    => $pressure0 + 0.0,
			e    => 0.0,
			e2   => 0.0,
			file => $OriginalCIF,
			dir  => $OriginalDir,
			);

	$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 $dir;
		if($ec->RelaxMode() eq 'scf') {
			$dir    = "e${Index}_${i}_0_${e}_0";
		}
		else {
			$dir    = "er${Index}_${i}_0_${e}_0";
		}
		my $OutCIF = "$dir.cif";
		print "Save strained CIF to [$OutCIF]\n";
		if($index[0] == $index[2] and $index[1] == $index[3]) {
			$ec->SaveStrainedStructureCabc($Crystal, $OutCIF, 
				$e, $index[0]-1, $index[1]-1);
			$ec->WriteTableLine(
				idx1 => $i,
				idx2 => 0,
				i1   => $index[0],
				j1   => $index[1],
				i2   => $index[2],
				j2   => $index[3],
				P    => $pressure0 + 0.0,
				e    => $e,
				e2   => 0.0,
				file => $OutCIF,
				dir  => $dir,
				);
		}
		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 $dir;
				if($ec->RelaxMode() eq 'scf') {
					$dir = "e${Index}_${i}_${j}_${e}_$e2";
				}
				else {
					$dir = "er${Index}_${i}_${j}_${e}_$e2";
				}
				my $OutCIF = "$dir.cif";
				print "  Save strained CIF to [$OutCIF]\n";
				$ec->SaveStrainedStructureCabc($Crystal, $OutCIF, 
					$e, $index[0]-1, $index[1]-1, $e2, $index[2]-1, $index[3]-1);
				$ec->WriteTableLine(
					idx1 => $i,
					idx2 => $j,
					i1   => $index[0],
					j1   => $index[1],
					i2   => $index[2],
					j2   => $index[3],
					P    => $pressure0 + 0.0,
					e    => $e,
					e2   => $e2,
					file => $OutCIF,
					dir  => $dir,
					);
			}
		}
	}
	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 $dir;
			if($ec->RelaxMode() eq 'scf') {
				$dir = "e${Index}_${i}_${j}_${e}_${e2}";
			}
			else {
				$dir = "er${Index}_${i}_${j}_${e}_${e2}";
			}
			my $OutCIF = "$dir.cif";
			print "  Save strained CIF to [$OutCIF]\n";
			$ec->SaveStrainedStructureCabc($Crystal, $OutCIF, 
				$e, $index[0]-1, $index[1]-1, $e2, $index[2]-1, $index[3]-1);
			$ec->WriteTableLine(
				idx1 => $i,
				idx2 => $j,
				i1   => $index[0],
				j1   => $index[1],
				i2   => $index[2],
				j2   => $index[3],
				P    => $pressure0 + 0.0,
				e    => $e,
				e2   => $e2,
				file => $OutCIF,
				dir  => $dir,
				);
		}
	}
	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 $dir;
			if($ec->RelaxMode() eq 'scf') {
				$dir = "e${Index}_${i}_0_${e}_0";
			}
			else {
				$dir = "er${Index}_${i}_0_${e}_0";
			}
			my $OutCIF = "$dir.cif";
			print "  Save strained CIF to [$OutCIF]\n";
			$ec->SaveStrainedStructureCabc($Crystal, $OutCIF, 
				$e, $index[0]-1, $index[1]-1);
				$ec->WriteTableLine(
					idx1 => $i,
					idx2 => 0,
					i1   => $index[0],
					j1   => $index[1],
					i2   => $index[2],
					j2   => $index[3],
					P    => $pressure0 + 0.0,
					e    => $e,
					e2   => 0.0,
					file => $OutCIF,
					dir  => $dir,
					);
		}
		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 $dir;
				if($ec->RelaxMode() eq 'scf') {
					$dir = "e${Index}_${i}_${j}_${e}_$e2";
				}
				else {
					$dir = "er${Index}_${i}_${j}_${e}_$e2";
				}
				my $OutCIF = "$dir.cif";
				print "Save strained CIF to [$OutCIF]\n";
				$ec->SaveStrainedStructureCabc($Crystal, $OutCIF, 
					$e, $index[0]-1, $index[1]-1, $e2, $index[2]-1, $index[3]-1);
				$ec->WriteTableLine(
					idx1 => $i,
					idx2 => $j,
					i1   => $index[0],
					j1   => $index[1],
					i2   => $index[2],
					j2   => $index[3],
					P    => $pressure0 + 0.0,
					e    => $e,
					e2   => $e2,
					file => $OutCIF,
					dir  => $dir,
					);
			}
		}
	}

	$ec->CloseTable();
}

sub MakeInputsVconp
{
	my ($Crystal) = @_;
	
	my $table = $ec->OpenTable(undef, 'w');
	if(!$table) {
		die "$!: Can not write to [$TableFile].\n";
	}

	$ec->SaveCIF($Crystal, $OriginalCIF, 0, 'unix');
	$App->print("  0: Ptarget=0.0: $OriginalCIF in [original]\n");
	$pressure0 += 0.0;
	$ec->WriteTableLine(
			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");
		$ec->WriteTableLine(
				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");
		$ec->WriteTableLine(
				idx1 => $i,
				P    => $Ptarget,
				file => $OriginalCIF,
				dir  => $dir,
				);
	}

	$ec->CloseTable();
}