package ElasticConstant;
use Exporter;
@ISA = qw(Exporter);

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

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

use ProgVars;
use Deps;
use Utils;

use Sci qw($torad $eVToJ $RyToeV $a0 $torad $todeg asin acos);
#use Sci::Science;
use Sci::Vector;
use Sci::Optimize;

use Crystal::CIF;
use Crystal::Crystal;
use Crystal::PWSCF;
use Crystal::GULP;
use Crystal::VASP;
use Crystal::WIEN2k;
#use Crystal::MyUtility;


BEGIN {}

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


#===============================================
# Constructer, Destructer
#===============================================
sub new
{
	my ($module) = @_;
	my $this = {};
	bless $this;
	return $this;
}

sub DESTROY
{
	my $this = shift;
}

#===============================================
# Subroutines
#===============================================
sub App
{ 
	my ($this) = @_;
	if(!$this->{App}) {
		$this->{App} = new MyApplication();
	}
	return $this->{App};
}
sub SetApplication
{
	my ($this, $App) = @_;
	$App = new MyApplication() if(!$App);
	return $this->{App} = $App;
}
sub print
{
	shift->App()->print(@_);
}
sub print
{
	shift->App()->printf(@_);
}

sub Engine() { return shift->{Engine}; };
sub SetEngine
{
	my ($this, $Engine) = @_;
	return $this->{Engine} = $Engine;
}

sub Mode() 
{ 
	my ($this) = @_;

	my $Mode = $this->{Mode};
	if($Mode =~ /\+/) {
		$Mode =~ s/\+//g;
	}
	return $Mode;
}
sub RelaxMode() 
{ 
	my ($this) = @_;
	return $this->{RelaxMode} if($this->{RelaxMode});
	
}
sub SetMode
{
	my ($this, $Mode) = @_;

	$this->{Mode} = $Mode;
	if($Mode =~ /\+/) {
		$this->{Mode}      =~ s/\+//g;
		$this->{RelaxMode} = 'relax';
	}
	elsif($Mode =~ /^init/i or $Mode =~ /^Vconp/i) {
		$this->{RelaxMode} = 'vc-relax';
	}
	else {
		$this->{RelaxMode} = 'scf';
	}
	return $this->{Mode};
}

sub pEngine() { return shift->{pEngine}; };
sub SetpEngine
{
	my ($this, $Engine) = @_;
	$Engine = $this->{Engine} if(!defined $Engine);

	if($Engine =~ /^gulp$/i) {
		$this->{pEngine} = new GULP;
	}
	elsif($Engine =~ /^vasp$/i) {
		$this->{pEngine} = new VASP;
	}
	elsif($Engine =~ /^pwscf$/i) {
		$this->{pEngine} = new PWSCF;
	}
	return $this->{pEngine};
}

sub var
{
	my ($this, $key) = @_;
	my $val = $this->App()->var($key);
	return $val;
}

sub OpenSummaryCSV
{
	my ($this, $path, $mode) = @_;
	$path = $this->var('SummaryCSV') if(!defined $path);
	$mode = 'r' if(!defined $mode);

	$this->{pSummaryCSV} = JFile->new($path, $mode);
	if(!$this->{pSummaryCSV}) {
		if($mode =~ /r/) {
			print "Error in ElasticConstant::SummaryCSV: Can not read [$path].\n";
		}
		else {
			print "Error in ElasticConstant::SummaryCSV: Can not write to [$path].\n";
		}
		exit;
	}
	return $this->{pSummaryCSV};
}
sub CloseSummaryCSV
{
	my ($this) = @_;
	return 0 if(!$this->{pSummaryCSV});
	$this->{pSummaryCSV}->Close();
	$this->{pSummaryCSV} = undef;
	return 1;
}
sub WriteSummaryCSV
{
	my ($this, @args) = @_;
	my $out = $this->{pSummaryCSV};
	$out->print(@args);
}

sub ReadSummaryCSV
{
	my ($this) = @_;
	my $in = $this->{pSummaryCSV};
	return $in->ReadLine();
}

sub OpenTable
{
	my ($this, $path, $mode) = @_;
	$path = $this->var('TableFile') if(!defined $path);
	$mode = 'r' if(!defined $mode);

#print "Open table file [$path].\n";
	$this->{pTable} = JFile->new($path, $mode);
	if(!$this->{pTable}) {
		if($mode =~ /r/) {
			print "Error in ElasticConstant::OpenTable: Can not read [$path].\n";
		}
		else {
			print "Error in ElasticConstant::OpenTable: Can not write to [$path].\n";
		}
		exit;
	}
	return $this->{pTable};
}
sub CloseTable
{
	my ($this) = @_;
	return 0 if(!$this->{pTable});
	$this->{pTable}->Close();
	$this->{pTable} = undef;
	return 1;
}
sub WriteTableLine
{
	my ($this, %p) = @_;
	my $out = $this->{pTable};

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

	my @array = qw(dir file idx1 idx2 i1 j1 i2 j2 e e2);
	foreach my $key (@array) {
		$p{$key} = '' if(!defined $p{$key});
	}
	$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 ($this) = @_;
	my $in = $this->{pTable};
#print "in=$in\n";

	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 GetFileConfig
{
	my ($this, $TemplateSource, $Engine, $Mode) = @_;
	$Engine = $this->var('Engine') if(!defined $Engine);
	$Mode   = $this->('Mode')      if(!defined $Mode);
#print "T[$TemplateSource] E[$Engine] M[$Mode]\n";

	my $RunScriptTemplate;
	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-Vconst-template.bat");
			}	
			elsif($Mode eq 'init') {
				$RunScriptTemplate = Deps::MakePath($TemplateSource, "run-gulp-init-template.bat");
			}
			else {
				print("\nError in Main: 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");
			$RunScriptTemplate = Deps::MakePath($TemplateSource, "run-pwscf-Vconp-template.bat");
			return ($RunScriptTemplate);

			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-Vconst-template.bat");
			}
			elsif($Mode eq 'CXYZ') {
				$RunScriptTemplate = Deps::MakePath($TemplateSource, "run-pwscf-Vconst-template.bat");
			}
			elsif($Mode eq 'init') {
				$RunScriptTemplate = Deps::MakePath($TemplateSource, "run-pwscf-init-template.bat");
			}
			else {
				print("\nError in Main: Mode [$Mode] is not implemented for Engine [$Engine].\n\n");
				exit;
			}
		}
		else {
		}
	}
	elsif($Engine =~ /^vasp$/i) {
	}
	return ($RunScriptTemplate);
}

sub GetInputPath
{
	my ($this, $Dir, $FileName, $Engine, $Mode) = @_;
	$Engine = $this->{Engine};
	$Mode   = $this->Mode();

	if($Engine =~ /^pwscf$/i) {
		my ($drive, $directory, $filename, $ext1, $lastdir, $filebody) = Deps::SplitFilePath($FileName);
		return Deps::MakePath($Dir, "$filebody.pwin", 0);
	}
	return undef;
}

sub GetOutputPath
{
	my ($this, $Dir, $FileName, $Engine, $Mode) = @_;
	$Engine = $this->{Engine};
	$Mode   = $this->Mode();

	if($Engine =~ /^pwscf$/i) {
		my ($drive, $directory, $filename, $ext1, $lastdir, $filebody) = Deps::SplitFilePath($FileName);
		return Deps::MakePath($Dir, "$filebody.pwout", 0);
	}
	return undef;
}

sub ReadFinalCrystalStructureFromOutput
{
	my ($this, $App, $InpFile, $WoutFile, $Engine, $Mode, $PrintStructures) = @_;
	$Engine = $this->{Engine};
	$Mode   = $this->{Mode};
	$PrintStructures = 0 if(!defined $PrintStructures);

	if($Engine =~ /^pwscf$/i) {
		my $Crystal = $this->{pEngine}->ReadFinalCrystalStructureFromOutput(
				$App, $InpFile, $WoutFile, $PrintStructures);
		return $Crystal;
	}
	return undef;
}

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

	my ($Etot);
	my $dir = $file;
	$dir =~ s/\.cif$//i;
	if($Engine =~ /^gulp$/i) {
		my $OutFile = "$dir.out";
		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);
		print("  Read Etot from [$OutFile].\n");
		$Etot = $this->{pEngine}->ReadEtotIneV($OutFile);
#print "Etot=$Etot\n";
#exit if($OutFile =~ /P1/);
	}
	elsif($Engine =~ /^vasp$/i) {
		my $OutFile = "SCF/OUTCAR";
		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 {
		print("Error in ElasticConstant::ReadEtot: Invalid Engine [$Engine].\n");
		exit;
	}
	return $Etot;
}

# ei' = ei + epsij*ej
# ai = aij*ej
# ai' aij*e' = ai + aij*epsjk * ek
sub SaveStrainedStructureCXYZ
{
	my ($this, $Crystal, $OutCIF, $e, $li, $lj, $e2, $li2, $lj2) = @_;

	my $mode = $this->Mode();

my $li1  = $li  + 1;
my $lj1  = $lj  + 1;
my $li22 = $li2 + 1;
my $lj22 = $lj2 + 1;
if(defined $e2) {
print "  strained structure for $mode: e$li1$lj1 = $e  e$li22$lj22 = $e2\n";
}
else {
print "  strained structure for $mode: e$li1$lj1 = $e\n";
}
	my @latt = $Crystal->LatticeParameters();
	my $NewCrystal = $Crystal->Duplicate();
	$NewCrystal->BurstToP1();

	my @v;
	if($mode =~ /^CXYZ/i) {
		my @aij  = $Crystal->LatticeVectorsByVectors();
		my @a2ij = (Vector::zero(3), Vector::zero(3), Vector::zero(3));
		my @eij  = Vector::UnitVectors();
		my @e2ij = @eij; #Vector::UnitVectors();
#print "a1=", $aij[0], "\n";
#print "a2=", $aij[1], "\n";
#print "a3=", $aij[2], "\n";
#print "e1=", $eij[0], "\n";
#print "e2=", $eij[1], "\n";
#print "e3=", $eij[2], "\n";
#print "e1'=", $e2ij[0], "\n";
#print "e2'=", $e2ij[1], "\n";
#print "e3'=", $e2ij[2], "\n";
		$e2ij[$li] += 0.5 * $e * $eij[$lj];
		$e2ij[$lj] += 0.5 * $e * $eij[$li];
		if(defined $e2) {
			$e2ij[$li2] += 0.5 * $e2 * $eij[$lj2];
			$e2ij[$lj2] += 0.5 * $e2 * $eij[$li2];
		}
#print "e1'=", $e2ij[0], "\n";
#print "e2'=", $e2ij[1], "\n";
#print "e3'=", $e2ij[2], "\n";

		for(my $i = 0 ; $i < 3 ; $i++) {
#			$a2ij[$i] = Vector::zero(3);
			for(my $j = 0 ; $j < 3 ; $j++) {
				$a2ij[$i][$j] = 0.0;
				for(my $k = 0 ; $k < 3 ; $k++) {
					$a2ij[$i][$j] += $aij[$i][$k] * $e2ij[$k][$j];
				}
			}
		}
#print "a1'=", $a2ij[0], "\n";
#print "a2'=", $a2ij[1], "\n";
#print "a3'=", $a2ij[2], "\n";

		$NewCrystal->SetLatticeVectorsByVectors(@a2ij);
		$NewCrystal->CalMetrics();
	}
	elsif($mode =~ /^Cabc/i) {
		($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();

		$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();
	}
	else {
		$this->App()->print("Error in ElasticConstant::SaveStrainedStructureCXYZ Invalid mode [$mode]\n");
		exit;
	}

	($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->LatticeVectors();
	@latt = $NewCrystal->LatticeParameters();
	printf("  latt: %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f\n", @latt);
	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]);

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

	return ($e, $OutCIF);
}

sub MakeFileBody
{
	my ($this, $i, $j, $e1, $e2) = @_;

	my $Mode  = $this->Mode();
	my $Index = $this->var('Index');

	my $s = '';
	if($this->RelaxMode() eq 'scf') {
		$s = '';
	}
	else {
		$s = 'r';
	}

	if($Mode =~ /^CXY/i) {
		return "eXY$s${Index}_${i}_${j}_${e1}_$e2";
	}
	elsif($Mode =~ /^Cabc/i) {
		return "eab$s${Index}_${i}_${j}_${e1}_$e2";
	}
	elsif($Mode =~ /^Vconst/i) {
		return "dV$s" . sprintf("%+05d", int($e1*10000));
	}
	elsif($Mode =~ /^Vconp/i) {
		return "P$s" . sprintf("%g", $e1*0.001);
	}
	return undef
}

#\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
{
	my ($this) = @_;

	my $App  = $this->App();
	my $ec   = $this;

	my $CIFFile     = $this->var('CIFFile');
	my $Mode        = $this->var('Mode');
	my $Index       = $this->var('Index');
	my $fdiff       = $this->var('fdiff');
	my $nStep       = $this->var('nStep');
	my $OriginalCIF = $this->var('OriginalCIF');
	my $OriginalDir = $this->var('OriginalDir');
	my $pressure0   = $this->var('pressure0');

	$App->print("\n");
	$App->print("Make input files for calculating dielectric constant by finite method\n");
	$App->print("\n");
#$App->PrintArgs("%-30s: %s\n", [qw(TableFile)]);

	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 = 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);
		return;
	}
	elsif($Mode =~ /^CXYZ/) {
		$ec->MakeInputsCXYZ($Crystal);
		return;
	}
	elsif($Mode =~ /^Vconst/){
		$this->MakeInputsVconst($Crystal);
		return;
	}
	elsif($Mode =~ /^Vconp/) {
		$this->MakeInputsVconp($Crystal);
		return;
	}
	elsif($Mode eq 'init') {
		$this->MakeInputsInit($Crystal);
		return;
	}
	else {
		$App->InvalidParameter('Mode', 'MakeInputs', 1);
		exit;
	}
}

sub MakeInputsVconst
{
	my ($this, $Crystal) = @_;

	my $App = $this->App();
	my $ec   = $this;

	my $Mode        = $this->var('Mode');
	my $Index       = $this->var('Index');
	my $fdiff       = $this->var('fdiff');
	my $nStep       = $this->var('nStep');
	my $OriginalCIF = $this->var('OriginalCIF');
	my $OriginalDir = $this->var('OriginalDir');
	my $pressure0   = $this->var('pressure0');

	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, '2');
	exit if(!$table);

	$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;
		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;
		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 ($this, $Crystal) = @_;

	my $App = $this->App();
	my $ec   = $this;

	my $Mode        = $this->var('Mode');
	my $Index       = $this->var('Index');
	my $fdiff       = $this->var('fdiff');
	my $nStep       = $this->var('nStep');
	my $OriginalCIF = $this->var('OriginalCIF');
	my $OriginalDir = $this->var('OriginalDir');
	my $pressure0   = $this->var('pressure0');

	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');
	exit if(!$table);

	$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 ($this, $Crystal) = @_;

	my $App = $this->App();
	my $ec   = $this;

	my $Mode        = $this->var('Mode');
	my $Index       = $this->var('Index');
	my $fdiff       = $this->var('fdiff');
	my $nStep       = $this->var('nStep');
	my $OriginalCIF = $this->var('OriginalCIF');
	my $OriginalDir = $this->var('OriginalDir');
	my $pressure0   = $this->var('pressure0');

	my $table = $ec->OpenTable(undef, 'w');
	exit if(!$table);

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


sub MakeInputsCXYZ
{
	my ($this, $Crystal) = @_;

	my $App = $this->App();
	my $ec   = $this;

	my $Mode        = $this->var('Mode');
	my $Index       = $this->var('Index');
	my $fdiff       = $this->var('fdiff');
	my $nStep       = $this->var('nStep');
	my $OriginalCIF = $this->var('OriginalCIF');
	my $OriginalDir = $this->var('OriginalDir');
	my $pressure0   = $this->var('pressure0');

	my @index = ($Index =~ /^(\d)(\d)(\d)(\d)/i);
	$App->print("Calculate elastic tensor for CXYZ_{", 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 @v    = $NewCrystal->LatticeVectorsByMatrix();
	$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], );
	$App->print("\n");

	$this->SaveCIF($NewCrystal, $OriginalCIF, 0, 'unix');
	my $table = $ec->OpenTable(undef, 'w');
	exit if(!$table);

	$this->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("  index %3d: eXY$index[0]$index[1]=$e\n", $i);

		my $dir    = $this->MakeFileBody($i, 0, $e, 0.0);
		my $OutCIF = "$dir.cif";
		print "  Save strained CIF to [$OutCIF]\n";
		if($index[0] == $index[2] and $index[1] == $index[3]) {
			$this->SaveStrainedStructureCXYZ($Crystal, $OutCIF, 
				$e, $index[0]-1, $index[1]-1);
			$this->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("    index %3d: eXY$index[0]$index[1]=$e2\n", $j);
				my $dir    = $this->MakeFileBody($i, $j, $e, $e2);
				my $OutCIF = "$dir.cif";
				print "  Save strained CIF to [$OutCIF]\n";
				$this->SaveStrainedStructureCXYZ($Crystal, $OutCIF, 
					$e, $index[0]-1, $index[1]-1, $e2, $index[2]-1, $index[3]-1);
			$this->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("      index %3d: eXY$index[2]$index[3]=$e2\n", $j);
			my $dir    = $this->MakeFileBody($i, $j, $e, $e2);
			my $OutCIF = "$dir.cif";
			print "  Save strained CIF to [$OutCIF]\n";
			$this->SaveStrainedStructureCXYZ($Crystal, $OutCIF, 
				$e, $index[0]-1, $index[1]-1, $e2, $index[2]-1, $index[3]-1);
			$this->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("  index %3d: eXY$index[0]$index[1]=$e\n", $i);
		if($index[0] == $index[2] and $index[1] == $index[3]) {
			my $dir    = $this->MakeFileBody($i, 0, $e, 0.0);
			my $OutCIF = "$dir.cif";
			print "  Save strained CIF to [$OutCIF]\n";
			$this->SaveStrainedStructureCXYZ($Crystal, $OutCIF, 
				$e, $index[0]-1, $index[1]-1);
				$this->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("    index %3d: eXY$index[2]$index[3]=$e2\n", $j);
				my $dir    = $this->MakeFileBody($i, $j, $e, $e2);
				my $OutCIF = "$dir.cif";
				print "   Save strained CIF to [$OutCIF]\n";
				$this->SaveStrainedStructureCXYZ($Crystal, $OutCIF, 
					$e, $index[0]-1, $index[1]-1, $e2, $index[2]-1, $index[3]-1);
				$this->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,
					);
			}
		}
	}

	$this->CloseTable();
}

sub MakeRunScript
{
	my ($this) = @_;
	
	my $App             = $this->App();
	my $Mode            = $this->Mode();
	my $CIFFile         = $App->var('CIFFile');
	my $RunScriptTemplate = $App->var('RunScriptTemplate');
	my $RunScriptPath   = $App->var('RunScriptPath');
	my $InitScriptPath  = $App->var('InitScriptPath');
	my $SampleName      = $App->var('SampleName');
	my $TableFile       = $App->var('TableFile');
#$App->PrintArgs("%-30s: %s\n", [qw(TableFile RunScriptTemplate RunScriptPath)]);
#exit;

	$App->print("\n");
	if($Mode =~ /^init/i) {
		$App->print("Make run script [$InitScriptPath] from template [$RunScriptTemplate]\n");
	}
	else {
		$App->print("Make run script [$RunScriptPath] from template [$RunScriptTemplate] and table [$TableFile]\n");
	}
	$App->print("\n");

	if($Mode =~ /^Cabc/) {
		$this->MakeRunScriptCXYZ($CIFFile, $SampleName, $RunScriptTemplate, $RunScriptPath);
		return;
	}
	elsif($Mode =~ /^CXYZ/) {
		$this->MakeRunScriptCXYZ($CIFFile, $SampleName, $RunScriptTemplate, $RunScriptPath);
		return;
	}
	elsif($Mode =~ /^Vconst/) {
		$this->MakeRunScriptVconst($CIFFile, $SampleName, $RunScriptTemplate, $RunScriptPath);
		return;
	}
	elsif($Mode =~ /^Vconp/) {
		$this->MakeRunScriptVconp($CIFFile, $SampleName, $RunScriptTemplate, $RunScriptPath);
		return;
	}
	elsif($Mode eq 'init') {
		$this->MakeRunScriptInit($CIFFile, $SampleName, $RunScriptTemplate, $RunScriptPath);
		return;
	}
	else {
		$App->InvalidParameter('Mode', 'MakeRunScript', 1);
		exit;
	}
}

sub MakeRunScriptInit
{
	my ($this, $CIFFile, $SampleName, $RunScriptTemplate, $RunScriptPath) = @_;

	my $ec              = $this;
	my $App             = $this->App();
	my $Engine          = $this->var('Engine');
	my $Mode            = $this->Mode();
	my $pressure0       = $this->var('pressure0');
	my $BurstToP1       = $this->var('BurstToP1');
	my $UseBravaisLattice = $this->var('UseBravaisLattice');
	my $OriginalDir     = $this->var('OriginalDir');
	my $OriginalCIF     = $this->var('OriginalCIF');
	my $OriginalInput   = $this->var('OriginalInput');
	my $nRelaxationStep = $this->var('nRelaxationStep');
	my $e_tol           = $this->var('e_tol');
	my $F_tol           = $this->var('F_tol');
	my $conv_thr        = $this->var('conv_thr');
	my $pressure_tol    = $this->var('pressure_tol');
	my $InitScriptPath  = $App->var('InitScriptPath');

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

	my %var = (
		Engine       => $Engine,
		Mode         => $Mode,
		Function     => $ec->RelaxMode(),
		SampleName   => $SampleName,
		CIFFile      => $CIFFile,
		UseBravaisLattice => $UseBravaisLattice,
		BurstToP1    => $BurstToP1,
		OriginalDir  => $OriginalDir,
		OriginalCIF  => $OriginalCIF,
		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 ($this, $CIFFile, $SampleName, $RunScriptTemplate, $RunScriptPath) = @_;

	my $ec              = $this;
	my $App             = $this->App();
	my $Engine          = $this->var('Engine');
	my $Mode            = $this->Mode();
	my $pressure0       = $this->var('pressure0');
	my $BurstToP1       = $this->var('BurstToP1');
	my $UseBravaisLattice = $this->var('UseBravaisLattice');
	my $OriginalDir     = $this->var('OriginalDir');
	my $OriginalCIF     = $this->var('OriginalCIF');
	my $OriginalInput   = $this->var('OriginalInput');
	my $nRelaxationStep = $this->var('nRelaxationStep');
	my $e_tol           = $this->var('e_tol');
	my $F_tol           = $this->var('F_tol');
	my $conv_thr        = $this->var('conv_thr');
	my $pressure_tol    = $this->var('pressure_tol');

	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');
	exit if(!$table);

	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, '{', '}');

#print "		Function     =>", $ec->RelaxMode(), "\n";
	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 ($this, $CIFFile, $SampleName, $RunScriptTemplate, $RunScriptPath) = @_;

	my $ec              = $this;
	my $App             = $this->App();
	my $Engine          = $this->var('Engine');
	my $Mode            = $this->Mode();
	my $pressure0       = $this->var('pressure0');
	my $BurstToP1       = $this->var('BurstToP1');
	my $UseBravaisLattice = $this->var('UseBravaisLattice');
	my $OriginalDir     = $this->var('OriginalDir');
	my $OriginalCIF     = $this->var('OriginalCIF');
	my $OriginalInput   = $this->var('OriginalInput');
	my $nRelaxationStep = $this->var('nRelaxationStep');
	my $e_tol           = $this->var('e_tol');
	my $F_tol           = $this->var('F_tol');
	my $conv_thr        = $this->var('conv_thr');
	my $pressure_tol    = $this->var('pressure_tol');

	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');
	exit if(!$table);

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

sub MakeRunScriptCXYZ
{
	my ($this, $CIFFile, $SampleName, $RunScriptTemplate, $RunScriptPath) = @_;
	
	my $App  = $this->App();
	my $Mode = $this->Mode();
	my $Engine = $this->var('Engine');
	my $pressure0 = $this->var('pressure0');
	my $BurstToP1 = $this->var('BurstToP1');
	my $UseBravaisLattice = $this->var('UseBravaisLattice');
	my $OriginalDir = $this->var('OriginalDir');
	my $OriginalCIF = $this->var('OriginalCIF');
	my $nRelaxationStep = $this->var('nRelaxationStep');
	my $e_tol    = $this->var('e_tol');
	my $F_tol    = $this->var('F_tol');
	my $conv_thr = $this->var('conv_thr');
	my $pressure_tol = $this->var('pressure_tol');

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

	my $table = $this->OpenTable(undef, 'r');
	exit if(!$table);

	my ($fStr, $rStr, $zStr);
	my ($ir, $if, $iz) = (0, 0, 0);
	while(1) {
		my $p = $this->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     => $this->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 SaveStrainedStructureCabc
{
	my ($this, $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();
	printf("  latt: %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f\n", @latt);
	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]);

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

	return ($e, $OutCIF);
}

sub SaveStrainedStructureV
{
	my ($this, $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();
	printf("  latt: %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f\n", @latt);

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

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

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

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

sub CalculateElasticTensor
{
	my ($this) = @_;

	my $Mode            = $this->var('Mode');
	my $OriginalCIF     = $this->var('OriginalCIF');

	$this->print("\nCalculate elastic tensor for Mode [$Mode]\n");

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

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

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

	my $App             = $this->App();
	my $nOrder          = $this->var('nOrder');

	my $table = $ec->OpenTable(undef, 'r');
	exit if(!$table);

	$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 [", $this->var('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 = $this->OpenSummaryCSV(undef, 'w');
	exit if(!$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");
	$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");
	}
	$this->CloseSummaryCSV();

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

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

	my $App             = $this->App();
	my $nOrder          = $this->var('nOrder');

	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');
	exit if(!$table);

	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 = $this->OpenSummaryCSV(undef, 'w');
	exit if(!$csv);

	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");
	}
	$this->CloseSummaryCSV();

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

sub CalculateElasticTensorCabcDiagonal
{
	my ($this, $Crystal, $vol, $pdata) = @_;
	my $ec = $this;

	my $App             = $this->App();
	my $Index           = $this->var('Index');
	my $nOrder          = $this->var('nOrder');

	$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 = $this->OpenSummaryCSV(undef, 'w');
	exit if(!$csv);

	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");
	}
	$this->CloseSummaryCSV();

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


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

	my $App             = $this->App();
	my $Index           = $this->var('Index');
	my $nOrder          = $this->var('nOrder');

	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');
	exit if(!$table);

	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 = $this->CalculateElasticTensorCabcDiagonal($Crystal, $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 [", $this->var('SummaryCSV'), "]\n");

	my $csv = $this->OpenSummaryCSV(undef, 'w');
	exit if(!$csv);

	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");
	$this->CloseSummaryCSV();
	
	$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";
}

1;