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

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

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

use ProgVars;
use Deps;
use Utils;

use Sci qw($torad $pi $eVToJ $RyToeV $a0 $torad $todeg asin acos $mn);
#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 printf
{
	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 GetiAtomInfFromCIF
{
	my ($this, $cifpath, $iAtom, $IsPrint) = @_;
	$IsPrint = 1 if(!defined $IsPrint);

	my $Crystal = CIF->new()->ReadCrystal($cifpath, 0);
	if(!$Crystal) {
		$this->print("Error: Can not read [$cifpath]\n\n") if($IsPrint);
		exit 2;
	}

	my $iAtom        = $this->var('iAtom');
	$iAtom = $this->GetiAtom($Crystal, $iAtom);
#	my @AtomSiteList = $Crystal->GetCAsymmetricAtomSiteList();
	my @AtomSiteList = $Crystal->GetCExpandedAtomSiteList();
	my $atom      = $AtomSiteList[$iAtom];
	my $label     = $atom->Label();
	my $name      = $atom->AtomNameOnly();
	my $charge    = $atom->Charge();
	my $Mass      = $Crystal->FindMass($name);
	my ($x,$y,$z) = $atom->Position();
	my $occ       = $atom->Occupancy();
	$this->print("Target atom: #$iAtom at $label: $name ($x, $y, $z)  "
			."M=$Mass  charge=$charge\n") if($IsPrint);

	return ($iAtom, $atom, $Mass);
}

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

	if($iAtom =~ /^\d+$/) {
	}
	elsif($iAtom =~ /^[A-Za-z]+$/) {
		$iAtom = $Crystal->FindiAtom(1, $iAtom);
	}
	elsif($iAtom =~ /\(.*\)/) {
		my ($x, $y, $z) = Utils::Split("\\s*[\\(\\),]\\s*", $iAtom);
		$iAtom = $Crystal->FindNearestSite($x, $y, $z);
#print "iAtom: pos=($x, $y, $z)\n";
	}
#print "iAtom: $iAtom\n";

	return $iAtom;
}

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{idx}, $p{idy}, $p{idz}, $p{e}, $p{ex}, $p{ey}, $p{ez}\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, $idx, $idy, $idz, $e, $ex, $ey, $ez, @a) 
		= Utils::Split("\\s*,\\s*", $line, 0);
print "$dir, $file, $idx, $idy, $idz, $e:", 
		join(':', @a), "\n";
	my $p = {
		dir  => $dir,
		file => $file,
		idx  => $idx + 0,
		idy  => $idy + 0,
		idz  => $idz + 0,
		e    => $e  + 0.0,
		ex   => $ex + 0.0,
		ey   => $ey + 0.0,
		ez   => $ez + 0.0,
		};
	return $p;
}

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

	my $RunScriptTemplate;
	if($Engine =~ /^gulp$/i) {
		if($OS eq 'MSWin32') {
			$RunScriptTemplate = Deps::MakePath($TemplateSource, 
						"run-gulp-vib-template.bat");
		}
		else {
		}
	}
	elsif($Engine =~ /^pwscf$/i) {
		if($OS eq 'MSWin32') {
#			@CopyList = ("GoPWSCF.bat", "runForDielectricConstant.bat");
			$RunScriptTemplate = Deps::MakePath($TemplateSource, 
						"run-pwscf-vib-template.bat");
			return ($RunScriptTemplate);
		}
		else {
		}
	}
	elsif($Engine =~ /^vasp$/i) {
		if($OS eq 'MSWin32') {
			$RunScriptTemplate = Deps::MakePath($TemplateSource, 
						"run-vasp-vib-template.bat");
			return ($RunScriptTemplate);
		}
		else {
			$RunScriptTemplate = Deps::MakePath($TemplateSource, 
						"run-vasp-vib-template.sh");
			my $DoVASPTemplate = Deps::MakePath($TemplateSource, 
						"run-vasp-vib-DoVASP-template.sh");
			return ($RunScriptTemplate, $DoVASPTemplate);
		}
	}
	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 ReadCrystalStructure
{
	my ($this, $App, $InpFile, $WoutFile, $Engine, $Mode, $PrintStructures) = @_;
	$Engine = $this->{Engine};
	$Mode   = $this->{Mode};
	$PrintStructures = 0 if(!defined $PrintStructures);

	my $Crystal;
	if($Engine =~ /^pwscf$/i) {
#		$Crystal = $this->{pEngine}->ReadFinalCrystalStructureFromOutput(
#				$App, $InpFile, $WoutFile, $PrintStructures);
		return $Crystal;
	}
	elsif($Engine =~ /^vasp$/i) {
		my ($drive, $CARDir, $filename, $ext, $lastdir, $filebody) 
			= Deps::SplitFilePath($InpFile);
		my $IsPrint = 0;
		$Crystal = $this->{pEngine}->ReadStructureFromCARFiles($CARDir, $PrintStructures);
		return $Crystal;
	}

	return undef;
}

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

	my ($Etot, $pForces);
	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 = Deps::MakePath($file, ["SCF", "OUTCAR"], 0);
		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->ReadLine();
			last if(!defined $line);

#			my $line = $in->SkipTo("free\\s+energy\\s+TOTEN\\s*=");
			if($line =~ /free\s+energy\s+TOTEN\s*=/) {
				($Etot) = ($line =~ /=\s*([+\-\.\d]+)/);
			}
			elsif($line =~ /^\s*FORCES\s+acting/) {
				$line = $in->SkipTo("TOTAL-FORCE");
				$in->ReadLine();
				my $c = 0;
				while(1) {
					$line = $in->ReadLine();
					last if(!defined $line);
					last if($line =~ /---/);
					my ($x, $y, $z, @f) = Utils::Split("\\s+", $line);
					$pForces->[$c] = \@f;
#print "f=", @f, "\n";
					$c++;
				}
			}
		}
		$in->Close();
	}
	else {
		print("Error in ElasticConstant::ReadEtot: Invalid Engine [$Engine].\n");
		exit;
	}
	return ($Etot, $pForces);
}

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 =~ /^Cabc/i) {
		return "eab$s${Index}_${i}_${j}_${e1}_$e2";
	}
	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 $CIFFile     = $this->var('CIFFile');
	my $Mode        = $this->var('Mode');
	my $iAtom       = $this->var('iAtom');
	my $xdiff       = $this->var('xdiff');
	my $nStep       = $this->var('nStep');
	my $OriginalCIF = $this->var('OriginalCIF');
	my $OriginalDir = $this->var('OriginalDir');

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

	$this->print("Read original structure from CIF File [$CIFFile]\n");

	my $Crystal = CIF->new()->ReadCrystal($OriginalCIF, 0);
	if(!$Crystal) {
		$this->print("Error: Can not read [$OriginalCIF]\n\n");
		exit 2;
	}

	$Crystal->BurstToP1();
	my $CrystalName = $Crystal->CrystalName();
	$this->print("CrystalName: $CrystalName\n");

	$iAtom = $this->GetiAtom($Crystal, $iAtom);
	my $atom = $Crystal->GetAtomSite(1, $iAtom);
	my $label = $atom->Label();
	my $name  = $atom->AtomName();
	my ($x, $y, $z) = $atom->Position();
	$this->print("Target atom: #$iAtom: $label: $name ($x, $y, $z)\n");

	if($Mode =~ /^vib/) {
		$this->MakeInputsVib($Crystal, $iAtom, $xdiff);
		return;
	}
	elsif($Mode eq 'init') {
		$this->MakeInputsInit($Crystal);
		return;
	}
	else {
		$this->InvalidParameter('Mode', 'MakeInputs', 1);
		exit;
	}
}

sub MakeInputsVib
{
	my ($this, $Crystal, $iAtom, $xdiff) = @_;

	my $Mode  = $this->var('Mode');
	my $nStep = $this->var('nStep');
	$nStep = 3 if($Mode eq 'vib');
	my $xdiff = $this->var('xdiff');

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

	$this->print("Calculate phonon vibration for Atom #$iAtom with dx = $xdiff\n");
	print "Mode=$Mode\n";
	print "nStep=$nStep\n";
	print "xdiff=$xdiff\n";

	my @latt = $Crystal->LatticeParameters();
	$this->print("Latt: $latt[0]  $latt[2]  $latt[3]  $latt[4]  $latt[5]\n");
	$this->print("xdiff: $xdiff\n");
	$this->print("nStep: $nStep\n");
	$this->print("OriginalCIF: $OriginalCIF\n");
	$this->print("OriginalDir: $OriginalDir\n");

	print "Save original CIF to [$OriginalCIF]\n";
	my $NewCrystal = $Crystal->Duplicate();
	$NewCrystal->BurstToP1();
	$this->SaveCIF($NewCrystal, $OriginalCIF, 0, 'unix');

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

	$this->WriteTableLine(
		idx => 0,
		idy => 0,
		idz => 0,
		e    => $xdiff,
		ex   => 0.0,
		ey   => 0.0,
		ez   => 0.0,
		file => $OriginalCIF,
		dir  => $OriginalDir,
		);

	$this->print("Displaced models\n");
	my $n = int($nStep / 2);
print "Displacement step: $nStep from -$n to $n\n";
#print "Mode: $Mode\n";
	if($Mode eq 'vib11') {
		for(my $ix = -$n ; $ix <= $n ; $ix++) {
			next if($ix == 0);

			my ($iy, $iz) = (0, 0);
			my $ex = $ix * $xdiff / $latt[0];
			my $ey = 0.0;
			my $ez = 0.0;
			$this->printf("  i=($ix,$iy,$iz): e=($ex, $ey, $ez)\n");

			my $dir = "exyz${ix}_${iy}_${iz}";
			my $OutCIF = "$dir.cif";
			print "Save strained CIF to [$OutCIF]\n";
			$this->SaveStrainedStructureVib($Crystal, $OutCIF, $iAtom, $ex, $ey, $ez);
			$this->WriteTableLine(	
				idx => $ix,
				idy => $iy,
				idz => $iz,
				e    => $xdiff,
				ex   => $ix * $xdiff,
				ey   => $iy * $xdiff,
				ez   => $iz * $xdiff,
				file => $OutCIF,
				dir  => $dir,
				);
		}
	}
	else {
		for(my $ix = -$n ; $ix <= $n ; $ix++) {
		for(my $iy = -$n ; $iy <= $n ; $iy++) {
		for(my $iz = -$n ; $iz <= $n ; $iz++) {
			next if($ix == 0 and $iy == 0 and $iz == 0);

			my $ex = $ix * $xdiff / $latt[0];
			my $ey = $iy * $xdiff / $latt[1];
			my $ez = $iz * $xdiff / $latt[2];
			$this->printf("  i=($ix,$iy,$iz): e=($ex, $ey, $ez)\n");

			my $dir = "exyz${ix}_${iy}_${iz}";
			my $OutCIF = "$dir.cif";
			print "Save strained CIF to [$OutCIF]\n";
			$this->SaveStrainedStructureVib($Crystal, $OutCIF, $iAtom, $ex, $ey, $ez);
			$this->WriteTableLine(	
				idx => $ix,
				idy => $iy,
				idz => $iz,
				e    => $ix * $xdiff,
				ex   => $ix * $xdiff,
				ey   => $iy * $xdiff,
				ez   => $iz * $xdiff,
				file => $OutCIF,
				dir  => $dir,
				);
		}
		}
		}
	}

	$this->CloseTable();
}

sub SaveStrainedStructureVib
{
	my ($this, $Crystal, $OutCIF, $iAtom, $ex, $ey, $ez) = @_;

	my $NewCrystal = $Crystal->Duplicate();
	$NewCrystal->BurstToP1();
#	my @AtomTypeList = $NewCrystal->GetCAtomTypeList();
#	my @AtomSiteList = $NewCrystal->GetCAsymmetricAtomSiteList();
	my @AtomSiteList = $NewCrystal->GetCExpandedAtomSiteList();
#	my $nAtomSite    = $NewCrystal->nTotalExpandedAtomSite();
	my $atom      = $AtomSiteList[$iAtom];
	my $label     = $atom->Label();
	my $atomname  = $atom->AtomNameOnly();
#	my $charge    = $atom->Charge();
	my ($x,$y,$z) = $atom->Position();
	my $occ       = $atom->Occupancy();
	$this->print("  $label: $atomname ($x, $y, $z) => ");
#	$atom->SetLabel($label);
#	$atom->SetAtomName($atomname);
	($x,$y,$z) = ($x+$ex, $y+$ey, $z+$ez);
	$this->print("($x, $y, $z)\n");
	$NewCrystal->SetAtomSite($iAtom, $label, $atomname, $x, $y, $z, $occ);
#	$atom->SetPosition($x, $y, $z);
#	$atom->SetOccupancy($occ);

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

	return ($ex, $ey, $ez, $OutCIF);
}

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

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

	if($Mode =~ /^vib/) {
		$this->MakeRunScriptVib($CIFFile, $SampleName, $RunScriptTemplate, $RunScriptPath);
		return;
	}
	elsif($Mode eq 'init') {
		$this->MakeRunScriptInit($CIFFile, $SampleName, $RunScriptTemplate, $RunScriptPath);
		return;
	}
	else {
		$this->InvalidParameter('Mode', 'MakeRunScript', 1);
		exit;
	}
}

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

	my $App             = $this->App();
	my $Engine          = $this->var('Engine');
	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) {
		$this->print("\nError in FrozenPhonon::MakeRunScriptInit: "
				."Can not read [$RunScriptTemplate].\n\n");
		exit;
	}

	my %var = (
		Engine       => $Engine,
		Function     => $this->RelaxMode(),
		SampleName   => $SampleName,
		CIFFile      => $CIFFile,
		UseBravaisLattice => $UseBravaisLattice,
		OriginalDir  => $OriginalDir,
		OriginalCIF  => $OriginalCIF,
		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, '{', '}');

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

sub MakeRunScriptVib
{
	my ($this, $CIFFile, $SampleName, $RunScriptTemplate, $RunScriptPath) = @_;
	
	my $App  = $this->App();
	my $Engine = $this->var('Engine');
	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 $DoVASPTemplate    = $App->var('DoVASPTemplate');
	my $DoVASPPath        = $App->var('DoVASPPath');
	$this->print("RunScriptTemplate: $RunScriptTemplate\n");
	$this->print("RunScriptPath    : $RunScriptPath\n");
	$this->print("DoVASPTemplate   : $DoVASPTemplate\n");
	$this->print("DoVASPPath       : $DoVASPPath\n");

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

	my ($CIFList) = ('"');
	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 $file = $p->{file};
		my $dir  = $p->{dir};
		$e  += 0.0;
print "$i1: $file\n";		

		$this->print(" Displacement: e=$e atm from [$file].\n");
		$CIFList .= "$file ";
	}
	$CIFList .= "\"\n";
	$table->Close();

	my %var = (
		Engine       => $Engine,
		Function     => $this->RelaxMode(),
		SampleName   => $SampleName,
		CIFFile      => $CIFFile,
		OriginalDir  => $OriginalDir,
		OriginalCIF  => $OriginalCIF,
		CIFList      => $CIFList,
		nRelaxationStep => $nRelaxationStep,
		e_tol        => $e_tol,
		F_tol        => $F_tol,
		conv_thr     => $conv_thr,
		pressure_tol => $pressure_tol * 0.001,
		);

	my $str   = JFile->ReadFile($RunScriptTemplate);
	if(!defined $str) {
		$this->print("\nError in FrozenPhonon::MakeRunScriptVib: "
			."Can not read [$RunScriptTemplate].\n\n");
		exit;
	}
	$str = Template->new()->ReplaceByHash($str, \%var, '{', '}');
	$this->print("Save to [$RunScriptPath].\n");
	my $ret = JFile->SaveFile($str, $RunScriptPath, "w");
	if(!defined $ret) {
		$this->print("\nError in MakeRunScriptVib: Can not write to [$RunScriptPath].\n\n");
		exit;
	}

	$str   = JFile->ReadFile($DoVASPTemplate);
	if(!defined $str) {
		$this->print("\nError in FrozenPhonon::MakeRunScriptVib: "
			."Can not read [$DoVASPTemplate].\n\n");
		exit;
	}
	$str = Template->new()->ReplaceByHash($str, \%var, '{', '}');
	$this->print("Save to [$DoVASPPath].\n");
	my $ret = JFile->SaveFile($str, $DoVASPPath, "w");
	if(!defined $ret) {
		$this->print("\nError in MakeRunScriptVib: Can not write to [$DoVASPPath].\n\n");
		exit;
	}
	
	system("chmod +x *.sh");
}

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

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

sub GetDataArray
{
	my ($this, $iAtom, $Engine, $Mode) = @_;

	my $App    = $this->App();
	my $e      = $this->var('xdiff');
	my $nStep  = $this->var('nStep');

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

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

		my $idx = $p->{idx};
		my $idy = $p->{idy};
		my $idz = $p->{idz};
		my $e    = $p->{e};
		my $ex   = $p->{ex};
		my $ey   = $p->{ey};
		my $ez   = $p->{ez};
		my $file = $p->{file};
		my $dir  = $p->{dir};
#$this->print("(id)=($idx, $idy, $idz)  file=$file  dir=$dir\n");
		my $PrintStructures = 1;
		my $Input  = Deps::MakePath($dir, ["SCF", "INCAR"]);
		my $Output = Deps::MakePath($dir, ["SCF", "OUTCAR"]);
		if(!-f $Input or !-f $Output) {
			$this->print("Warning: Result files [$Input][$Output] are missing. Skip.\n\n");
			next;
		}
		
		my $cry = $this->ReadCrystalStructure($App, $Input, $Output, 
				$Engine, $Mode, $PrintStructures);
		my @sites = $cry->GetCExpandedAtomSiteList();
		my $atom2 = $sites[$iAtom];
		my ($x2, $y2, $z2) = $atom2->Position();
		my ($Etot, $pForces) = $this->ReadEtot($dir);
		my $pf = $pForces->[$iAtom];

		$Emin = $Etot if($Emin > $Etot);
		$data[$c] = {
			idx   => $idx + 0,
			idy   => $idy + 0,
			idz   => $idz + 0,
			e     => $e  + 0.0,
			ex    => $ex + 0.0,
			ey    => $ey + 0.0,
			ez    => $ez + 0.0,
			pos   => [$x2, $y2, $z2],
			Etot  => $Etot + 0.0,
			Force => [$pf->[0], $pf->[1], $pf->[2]],
			file  => $file,
			dir   => $dir,
			};
$this->print("  (id)=($idx, $idy, $idz)  e=$e: ($x2, $y2, $z2) Etot=$Etot eV  f=($pf->[0], $pf->[1], $pf->[2]) file=$file\n");
#print "c=$c e=$e  $data[$c]{e}\n";
		$c++;
	}
	$table->Close();

#print "e=$e\n";
#$Emin = 0.0;
	print "  Emin = $Emin eV\n";
	my @idata;
	my $off = int($nStep / 2);
	for(my $i = 0 ; $i < @data ; $i++) {
		$data[$i]{Etot} -= $Emin;
		my $idx = $data[$i]{idx} + $off;
		my $idy = $data[$i]{idy} + $off;
		my $idz = $data[$i]{idz} + $off;
		$idata[$idx][$idy][$idz] = $data[$i];
	}
	$this->print("\n");

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

	return ($Emin, $off, \@data, \@idata);
}

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

	my $App            = $this->App();
	my $Mode           = $this->var('Mode');
	my $Engine         = $this->var('Engine');
	my $nStep          = $this->var('nStep');
	my $nOrder         = $this->var('nOrder');
	my $e              = $this->var('xdiff');
	my $OriginalCIF    = $this->var('OriginalCIF');

	my $Crystal = CIF->new()->ReadCrystal($OriginalCIF, 0);
	if(!$Crystal) {
		$this->print("Error: Can not read [$OriginalCIF]\n\n");
		exit 2;
	}

	$this->print("\nCalculate atomic vibration in 1D [$Mode]\n");
	$this->print("Original CIF File: $OriginalCIF\n");

	my ($iAtom, $pAtom, $Mass) = $this->GetiAtomInfFromCIF($OriginalCIF, $this->var('iAtom'), $App);

	my ($Emin, $off, $pdata, $pidata) = $this->GetDataArray($iAtom, $Engine, $Mode);
	my @data  = @$pdata;
	my @idata = @$pidata;
	my $ndata = @data;

	$this->print("\n");
	$this->print("Secnd derivative by least squres fitting\n");
	$this->print("  nData : $ndata\n");
	$this->print("  nOrder: $nOrder\n");

	my $Opt = new Optimize;
	my $Method = "LCFunction";
	print "  use $Method\n";
	my ($OptVars, $MinVal);
	my $m = $nOrder;
	my $iPrintLevel = 2;
	my (@v, @f);
	for(my $i = 0 ; $i < $ndata ; $i++) {
		$v[$i] = $e * $data[$i]{idx};
		$f[$i] = $data[$i]{Etot};
print "i=$i: e=$v[$i]}  f=$f[$i]\n";
	}
	my ($OptVars, $MinVal) = $Opt->Optimize("mlsq", \@v, \@f, $m, $iPrintLevel);
	print "  Optimized at S = $MinVal:\n";
	$this->print("\n");
	$this->print("E(dx,dy,dz) = a[0] + a[1] * dx + a[2] * dx*dx +...\n");
	for(my $i = 0 ; $i < $m ; $i++) {
		$this->print("    a[$i]  = $OptVars->[$i]\n");
	}

	$this->print("\n");
	$this->print("   Writing fitting results to [", $this->var('SummaryCSV'), "]\n");

	my $csv = $this->OpenSummaryCSV(undef, 'w');
	$csv->print("x,y,z,a[2](fit),a[2](cal),Etot,Efit\n");
	my (@dx, @Etot);
	for(my $i = 0 ; $i < $ndata ; $i++) {
		my $p = $data[$i];
		$dx[$i]   = $p->{ex};
		$Etot[$i] = $p->{Etot};
	}
	for(my $i = 0 ; $i < $ndata ; $i++) {
		my $p = $data[$i];
		my $x = $dx[$i];
		my $y = 0.0;
		my $z = 0.0;
		my $Etot = $Etot[$i]; #$p->{Etot};
		my $Efit = $Opt->mlsqYCal($x, $OptVars);
		my $W = 2.0 * $OptVars->[2];
		for(my $j = 3 ; $j < $m ; $j++) {
			$W += $m * ($m-1) * $OptVars->[$j] * $x**($j-2);
		}
		my $Wcal = '';
		if($i > 0 and $i <= @data-2) {
			$Wcal = Algorism::Differentiate2ndByIndex($i, \@Etot, $e);
		}
		$csv->print("$x,$y,$z,$W,$Wcal,$Etot,$Efit\n");
	}
	exit if(!$csv);

	$this->CloseSummaryCSV();

	my $f00 = 2.0 * $pi * sqrt(2.0 * abs($OptVars->[2]) * $eVToJ / 1.0e-20 / $Mass / $mn);
	$f00 = -$f00 if($OptVars->[2] < 0.0);
	$this->print("  Eigen frequencies (M = $Mass mn)\n");
	$this->printf("    f00 = %14.6g Hz (Wi = $OptVars->[2])\n", $f00);

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

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

	my $Mode           = $this->var('Mode');
	if($Mode eq 'vib11') {
		return $this->Calculate1D();
	}

	my $App            = $this->App();
	my $Engine         = $this->var('Engine');
	my $e              = $this->var('xdiff');
	my $OriginalCIF    = $this->var('OriginalCIF');

	my $Crystal = CIF->new()->ReadCrystal($OriginalCIF, 0);
	if(!$Crystal) {
		$this->print("Error: Can not read [$OriginalCIF]\n\n");
		exit 2;
	}

	$this->print("\nCalculate atomic vibration\n");
	$this->print("Original CIF File: $OriginalCIF\n");

	my ($iAtom, $pAtom, $Mass) = $this->GetiAtomInfFromCIF($OriginalCIF, $this->var('iAtom'), $App);

	my ($Emin, $off, $pdata, $pidata) = $this->GetDataArray($iAtom, $Engine, $Mode);
	my @data  = @$pdata;
	my @idata = @$pidata;
	my $ndata = @data;

	$this->print("\n");
	$this->print("Hessian matrix by least squres fitting\n");
	$this->print("  nData : $ndata\n");
#	$this->print("  nOrder: $nOrder\n");

	my $Opt = new Optimize;
	my $Method = "LCFunction";
	print "  use $Method\n";
	my ($OptVars, $MinVal);
#	my $m = 7;
	my $m = 10;
	my $iPrintLevel = 2;
	my (@v, @f);
	for(my $i = 0 ; $i < $ndata ; $i++) {
		$v[$i] = [$e * $data[$i]{idx}, $e * $data[$i]{idy}, $e * $data[$i]{idz}];
		$f[$i] = $data[$i]->{Etot};
#		$f[$i] = $data[$i]->{Etot} - $Emin;
print "i=$i: ", join(',', @{$v[$i]}), "  f=$f[$i]\n";
	}
	($OptVars, $MinVal) = $Opt->Optimize(
		$Method, \@v, \@f, $m,
		sub { &LCFunc(@_); },
		$iPrintLevel
		);
	print "  Optimized at S = $MinVal:\n";
	$this->print("\n");
	$this->print("E(dx,dy,dz) = a0 + ax * dx + ay * dy + az * dz\n");
	$this->print("                 + axx * dx*dx + ayy * dy*dy + azz * dz*dz\n");
	$this->print("                 + 2 * (ayz * dy*dz + azx * dz*dx + axy * dx*dy)\n");

	my ($a0, $ax, $ay, $az, $axx, $ayy, $azz, $ayz, $azx, $axy) = @$OptVars;
	my (@f, @a1num, @a2num);
	$f[0] = $idata[$off][$off][$off]{Force}[0];
	$f[1] = $idata[$off][$off][$off]{Force}[1];
	$f[2] = $idata[$off][$off][$off]{Force}[2];
	$a1num[0] = ($idata[$off+1][$off][$off]{Etot} - $idata[$off-1][$off][$off]{Etot}) / 2.0 / $e;
	$a1num[1] = ($idata[$off][$off+1][$off]{Etot} - $idata[$off][$off-1][$off]{Etot}) / 2.0 / $e;
	$a1num[2] = ($idata[$off][$off][$off+1]{Etot} - $idata[$off][$off][$off-1]{Etot}) / 2.0 / $e;
	$a2num[0] = 0.5 * ($idata[$off+1][$off][$off]{Etot} +$idata[$off-1][$off][$off]{Etot}
			- 2.0 * $idata[$off][$off][$off]{Etot}) / $e / $e;
	$a2num[1] = 0.5 * ($idata[$off][$off+1][$off]{Etot} +$idata[$off][$off-1][$off]{Etot}
			- 2.0 * $idata[$off][$off][$off]{Etot}) / $e / $e;
	$a2num[2] = 0.5 * ($idata[$off][$off][$off+1]{Etot} +$idata[$off][$off][$off-1]{Etot}
			- 2.0 * $idata[$off][$off][$off]{Etot}) / $e / $e;

	my @label = qw(a0 ax ay az axx ayy azz ayz azx axy);
	$this->print("coeff\tfit\t(num diff)\t(H-F force)\n");
	for(my $i = 0 ; $i < @label ; $i++) {
		if($i >= 1 and $i <= 3) {
			$this->printf("    %3s = %14.6f (%14.6f) (%14.6f)\n", 
				$label[$i], $OptVars->[$i], $a1num[$i-1], $f[$i-1]);
		}
		elsif($i >= 4 and $i <= 6) {
			$this->printf("    %3s = %14.6f (%14.6f)\n", 
				$label[$i], $OptVars->[$i], $a2num[$i-4]);
		}
		else {
			$this->printf("    %3s = %14.6f\n", $label[$i], $OptVars->[$i]);
		}
	}

	$this->print("\n");
	$this->print("   Writing fitting results to [", $this->var('SummaryCSV'), "]\n");

	my $csv = $this->OpenSummaryCSV(undef, 'w');
	$csv->print("x,y,z,Etot,Efit\n");
	for(my $i = 0 ; $i < @data ; $i++) {
		my $p = $data[$i];
		my $x = $p->{idx} * $e;
		my $y = $p->{idy} * $e;
		my $z = $p->{idz} * $e;
		my $e   = $p->{e};
		my $Etot = $p->{Etot};
		my $Efit = &YCal($x, $y, $z, @$OptVars);
		$csv->print("$x,$y,$z,$Etot, $Efit\n");
	}
	exit if(!$csv);

	$this->CloseSummaryCSV();

	$this->print("\n");
	$this->print("Diagonalize Hessian matrix\n");
	my $pHij = new Math::MatrixReal(3, 3);
	$pHij->assign(1, 1, $axx);
	$pHij->assign(2, 2, $ayy);
	$pHij->assign(3, 3, $azz);
	$pHij->assign(1, 2, $axy);
	$pHij->assign(2, 1, $axy);
	$pHij->assign(2, 3, $ayz);
	$pHij->assign(3, 2, $ayz);
	$pHij->assign(3, 1, $azx);
	$pHij->assign(1, 3, $azx);
	my ($pWi, $pVij) = $pHij->sym_diagonalize();
	my (@Wi, @Vij);
	for(my $i = 0  ; $i < 3 ; $i++) {
		$Wi[$i] = $pWi->element($i+1, 1);
		for(my $j = 0 ; $j < 3 ; $j++) {
			$Vij[$i][$j] = $pVij->element($i+1, $j+1);
		}
	}
	$this->print("  Eigen values Wi:\n");
	print $pWi;
	$this->print("  Eigen vectors Vij: dE(dx',dy',dz') = Wi[0] * dx'*dx' + Wi[0] * dy'*dy' + Wi[2] * dz'*dz'\n");
	print $pVij;
	
	my $f00 = 2.0 * $pi * sqrt(2.0 * abs($Wi[0]) * $eVToJ / 1.0e-20 / $Mass / $mn);
	$f00 = -$f00 if($Wi[0] < 0.0);
	my $f11 = 2.0 * $pi * sqrt(2.0 * abs($Wi[1]) * $eVToJ / 1.0e-20 / $Mass / $mn);
	$f11 = -$f11 if($Wi[1] < 0.0);
	my $f22 = 2.0 * $pi * sqrt(2.0 * abs($Wi[2]) * $eVToJ / 1.0e-20 / $Mass / $mn);
	$f22 = -$f22 if($Wi[2] < 0.0);
	$this->print("  Eigen frequencies (M = $Mass mn)\n");
	$this->printf("    f00 = %14.6g Hz (Wi = $Wi[0])\n", $f00);
	$this->printf("    f11 = %14.6g Hz (Wi = $Wi[1])\n", $f11);
	$this->printf("    f22 = %14.6g Hz (Wi = $Wi[2])\n", $f22);

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

sub LCFunc
{
	my ($idx, $v) = @_;
#print "idx=$idx  v=", join(',', @$v), "\n";
	if($idx == 0) {
		return 1.0;
	}
	elsif($idx == 1) {
		return $v->[0];
	}
	elsif($idx == 2) {
		return $v->[1];
	}
	elsif($idx == 3) {
		return $v->[2];
	}
	elsif($idx == 4) {
		return 1.0 * $v->[0]*$v->[0];
	}
	elsif($idx == 5) {
		return 1.0 * $v->[1]*$v->[1];
	}
	elsif($idx == 6) {
		return 1.0 * $v->[2]*$v->[2];
	}
	elsif($idx == 7) {
		return 2.0 * $v->[1]*$v->[2];
	}
	elsif($idx == 8) {
		return 2.0 * $v->[2]*$v->[0];
	}
	elsif($idx == 9) {
		return 2.0 * $v->[0]*$v->[1];
	}

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

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

1;