#!/usr/bin/perl

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

use strict;
#use warnings;
use File::Path;
use File::Basename;
use File::Find;

#use Math::Matrix;
#use Math::MatrixReal;

use Deps;
use Utils;
use JFile;

use MyApplication;
use Sci qw($pi $torad $todeg);
#use Sci::Algorism;
#use Sci::GeneralFileFormat;

use Crystal::CIF;
use Crystal::Crystal;
use Crystal::SpaceGroup;
use Crystal::AtomType;
use Crystal::PWSCF;
#use Crystal::VESTA;
use Crystal::XCrySDen;

#===============================================
# デバッグ関係変数
#===============================================
#$PrintLevelが大きいほど、情報が詳しくなる
my $PrintLevel = 0;

#===============================================
# 文字コード関係変数
#===============================================
# sjis, euc, jis, noconv
my $PrintCharCode      = Deps::PrintCharCode();
my $OSCharCode         = Deps::OSCharCode();
my $FileSystemCharCode = Deps::FileSystemCharCode();

my $LF        = Deps::LF();
my $DirSep    = Deps::DirSep();
my $RegDirSep = Deps::RegDirSep();

my $PWSCFPPDir = ProgVars::QEPPDir();
print "QE PP Dir: $PWSCFPPDir\n";
my $PPSaveDir  = 'PP';

#===============================================
# Applicationオブジェクト作成
#===============================================
my $App = new MyApplication;
exit if($App->Initialize() < 0);

#$App->SetLF("<br>\n");
#$App->SetPrintCharCode("sjis");
#$App->SetDebug($Debug);
$App->SetDeleteHTMLFlag(1);

#===============================================
# スクリプト大域変数
#===============================================
my $Program             = basename($0);
my $InitialDirectory    = Deps::GetWorkingDirectory();
my %ParamHash;

#==========================================
# コマンドラインオプション読み込み
#==========================================
$App->AddArgument("--Action",              
	"--Action=[SearchPP|FindPP|MakeInputFiles|MakeCIFFilesFromOutput|MakeXSF|MakeHistoryCSV|UpdateInput"
			."|SummarizeStatus]",
	'');
$App->AddArgument("--BurstToP1",           "--BurstToP1=[0|1] (Def: 0)", 0);
$App->AddArgument("--PPType",              "--PPType=[US|USPP|NC|PAW|SL|1/r] (Def: PAW)", 'PAW');
$App->AddArgument("--Functional",          "--Functional=[LDA|SLA|PZ|PBE|PBESOL|LYP|B88|BLYP|P86|TPSS"
										               ."|PW|PSX|PSC|PBC|GGC|GGX|NOG|NOGX|NOGC] (Def: PBE)", 'PBE');
$App->AddArgument("--PPCopyMode",          "--PPCopyMode=[All|Selected|No] (Def: Selected)", 'Selected');
$App->AddArgument("--UseBravaisLattice",   "--UseBravaisLattice=[0|1] (Def: 0)", 0);
$App->AddArgument("--UseCELLDM",           "--UseCELLDM=[0|1] (Def:1)",          1);
$App->AddArgument("--Function",            "--Function=[scf|nscf|band|edensity|relax|vc-relax|md|phonon] (Def:scf)", "scf");
$App->AddArgument("--SpinPolarized",       "--SpinPolarized=[No|Yes] (Def:No)",  "No");
$App->AddArgument("--aKProduct",           "--aKProduct=val (Def:1.5)",          1.5);
$App->AddArgument("--ecutwfc",             "--ecutwfc=val (Def:x1.0,30.0)",      'x1.0,30.0');
$App->AddArgument("--ecutrho",             "--ecutrho=val (Def:x1.0,140.0)",     'x1.0,140.0');
$App->AddArgument("--nbnd",                "--nbnd=val (Def:x1.2)",              'x1.2');
$App->AddArgument("--KPoints",             "--KPoints=[File|XG|GX|..] (Def:File)", 'File');
$App->AddArgument("--nk",                  "--nk=val (Def:50)",              50);
$App->AddArgument("--nbnd",                "--nbnd=val (Def:x1.2)",              'x1.2');
$App->AddArgument("--press",               "--press=val in kBar (Def:0.001)",          0.001);
$App->AddArgument("--etot_conv_thr",       "--etot_conv_thr=val (Def:0.5e-4)", 0.5e-4);
$App->AddArgument("--forc_conv_thr",       "--forc_conv_thr=val (Def:0.5e-3)", 0.5e-3);
$App->AddArgument("--conv_thr",            "--conv_thr=val (Def:0.5e-6)",      0.001);
$App->AddArgument("--press_conv_thr",      "--press_conv_thr=val in kBar (Def:0.001)", 0.001);
$App->AddArgument("--DebugMode",           "--DebugMode: Set DebugMode",         '');

exit 1 if($App->ReadArgs(1, "sjis", 0) != 1);
my $Args = $App->Args();

my %ArgHash = $Args->GetArgHash();
foreach my $key (keys %ArgHash) {
#print "key: $key: $ArgHash{$key}\n";
	if($key =~ /^Param:(.*?)$/i) {
		$ParamHash{$1} = $ArgHash{$key};
#print "$1:  $ArgHash{$key}\n";
	}
}
$App->{pParamHash} = \%ParamHash;

#==========================================
# メイン関数スタート
#==========================================

my $Debug = $Args->GetGetArg("DebugMode");
$App->SetDebug($Debug);
my $Action = $Args->GetGetArg("Action");

my $ret = 0;
if($Action =~ /SearchPP/i) {
	&SearchPP();
}
elsif($Action =~ /FindPP/i) {
	&FindPP();
}
elsif($Action =~ /MakeInputFiles/i) {
	&MakeInputFiles();
}
elsif($Action =~ /MakeCIFFilesFromOutput/i) {
	&MakeCIFFilesFromOutput();
}
elsif($Action =~ /MakeXSF/i) {
	&MakeXSF();
}
elsif($Action =~ /MakeHistoryCSV/i) {
	&MakeHistoryCSV();
}
elsif($Action =~ /UpdateInput/i) {
	&UpdateInput();
}
elsif($Action =~ /SummarizeStatus/i) {
	&SummarizeStatus();
}
else {
	$App->print("Error: Invalid Action: [$Action]\n");
}

#Utils::EndHTML();

exit $ret;

#===============================================
# スクリプト終了
#===============================================

#==========================================
# &Subroutines
#==========================================
sub MakeSampleNameFromPath
{
	my ($path) = @_;

	my @filenames  = fileparse($path, "\.[^\.]+");
	my $SampleName = $filenames[0];
	$SampleName    =~ s/\..*$//;
	return $SampleName;
}

sub GetCrystalFromCIF
{
	my ($CIFFile) = @_;

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

	my $CrystalName = $Crystal->CrystalName();
	print "CrystalName: $CrystalName\n";
	my ($a,$b,$c,$alpha,$beta,$gamma) = $Crystal->LatticeParameters();
	printf "cell: %9.6f %9.6f %9.6f   %9.6f %9.6f %9.6f \n", $a, $b, $c, $alpha, $beta, $gamma;
	my ($a11, $a12, $a13, $a21, $a22, $a23, $a31, $a32, $a33) = $Crystal->LatticeVectors();
	printf "  Vectors: (%9.6f %9.6f %9.6f)\n", $a11, $a12, $a13;
	printf "           (%9.6f %9.6f %9.6f)\n", $a21, $a22, $a23;
	printf "           (%9.6f %9.6f %9.6f)\n", $a31, $a32, $a33;
	my $vol = $Crystal->Volume();
	print "  Volume: $vol A^3\n";
#対称要素を使い、非対称単位中の原子を展開する
#この際、Densityも計算される
	$Crystal->ExpandCoordinates();
	my $Density = $Crystal->Density();
	print "  Density: $Density g/cm^3\n";
	print "\n";

	my $SPG = $Crystal->GetCSpaceGroup();
	my $SPGName = $SPG->SPGName();
	my $iSPG    = $SPG->iSPG();
	my $LatticeSystem = $SPG->LatticeSystem();
	print "Space Group: $SPGName ($iSPG)\n";
	print "Lattice System: $LatticeSystem\n";
	my $nTranslation = $SPG->nTranslation();
	print "nTranslation: $nTranslation$LF";
	for(my $i = 0 ; $i < $nTranslation ; $i++) {
		my ($x,$y,$z) = $SPG->TranslationVector($i+1);
		print "   ($x, $y, $z)$LF";
	}
	my $nSymmetryOperation = $SPG->nSymmetryOperation();
	print "nSymmetryOperation: $nSymmetryOperation$LF";
	for(my $i = 0 ; $i < $nSymmetryOperation ; $i++) {
		my $symop = $SPG->SymmetryOperation($i+1);
		print "   $symop$LF";
	
		my ($x1,$y1,$z1,$t1,
			$x2,$y2,$z2,$t2,
			$x3,$y3,$z3,$t3)
			= $SPG->SymmetryOperationByMatrix($i+1);
		print "      $x1 $y1 $z1  $t1$LF";
		print "      $x2 $y2 $z2  $t2$LF";
		print "      $x3 $y3 $z3  $t3$LF";
	}
	print "\n";

	my @AtomTypeList = $Crystal->GetCAtomTypeList();
	my $nAtomType = @AtomTypeList;
	print "nAtomType: $nAtomType$LF";
	for(my $i = 0 ; $i < $nAtomType ; $i++) {
		my $label    = $AtomTypeList[$i]->Label();
		my $atomname = $AtomTypeList[$i]->AtomName();
		my $charge   = $AtomTypeList[$i]->Charge();
		my $AtomicNumber = $AtomTypeList[$i]->AtomicNumber();
		my $AtomicMass   = $AtomTypeList[$i]->AtomicMass();
		my $i1 = $i+1;
		print "   #$i1: $label : $atomname [$AtomicNumber] ($charge) "
			."($AtomicMass)$LF";
	}
	print "$LF";
#Crystalクラス中の、非対称単位中の原子 AsymmetricAtomSiteクラスのリストをとる
	my @AtomSiteList = $Crystal->GetCAsymmetricAtomSiteList();
	my $nAsymmetricAtomSite = @AtomSiteList;
	print "nAsymmetricAtomSite: $nAsymmetricAtomSite$LF";
	for(my $i = 0 ; $i < $nAsymmetricAtomSite ; $i++) {
		my $label     = $AtomSiteList[$i]->Label();
		my $type      = $AtomSiteList[$i]->AtomName();
		my ($x,$y,$z) = $AtomSiteList[$i]->Position(1);
		my $occupancy = $AtomSiteList[$i]->Occupancy();
		my $iAtomType = $Crystal->FindiAtomType($type);

		print "   $label ($type)[#$iAtomType]: ($x, $y, $z) [$occupancy]$LF";
	}

	my @ExpandedAtomSiteList   = $Crystal->GetCExpandedAtomSiteList();
	my $nTotalExpandedAtomSite = @ExpandedAtomSiteList;
	print "nTotalExpandedAtomSite: $nTotalExpandedAtomSite\n";
	for(my $i = 0 ; $i < $nTotalExpandedAtomSite ; $i++) {
		my $label     = $ExpandedAtomSiteList[$i]->Label();
		my $type      = $ExpandedAtomSiteList[$i]->AtomName();
		my ($x,$y,$z) = $ExpandedAtomSiteList[$i]->Position(1);
		my $occupancy = $ExpandedAtomSiteList[$i]->Occupancy();
		my $iAtomType = $Crystal->FindiAtomType($type);

		print "   $label ($type)[#$iAtomType]: ($x, $y, $z) [$occupancy]$LF";
	}

	return $Crystal;
}

sub FindPP
{
	my $InputFile        = $Args->GetGetArg(0);
	my $TargetFunctional = $Args->GetGetArg('Functional');
	$TargetFunctional    = 'PBE' if(!defined $TargetFunctional);
	my $TargetPPType     = $Args->GetGetArg('PPType');
	$TargetPPType        = 'PAW' if(!defined $TargetPPType);

	my $IsPrint = 1;

	$App->print("\n");
	$App->print("Find pseudo potentials...\n");
	$App->print("Input File           : $InputFile\n");
	$App->print("Pseudo potential type: $TargetPPType\n");
	$App->print("Functional           : $TargetFunctional\n");
	$App->print("\n");

	my $PWSCF = new PWSCF;
	my $Crystal;
	if($InputFile =~ /\.cif$/i) {
		my $CIF = new CIF();
		unless($CIF->Read($InputFile)) {
			$App->print("Error: Can not read [$InputFile].\n\n");
			return 0;
		}
		$Crystal = $CIF->GetCCrystal();
	}
	else {
		$Crystal = new Crystal;
		$PWSCF->ReadCrystalStructureFromInput($App, $InputFile, $Crystal);
	}
	my @AtomTypeList = $Crystal->GetCAtomTypeList();

	return $PWSCF->FindPPFiles(\@AtomTypeList, $TargetFunctional, $TargetPPType, $IsPrint);
}

sub SearchPP
{
	my $Functional      = $Args->GetGetArg('Functional');
	$Functional = 'PBE' if(!defined $Functional);
	my $PPType          = $Args->GetGetArg('PPType');
	$PPType = 'PAW' if(!defined $PPType);

	$App->print("\n");
	$App->print("Pseudo potential type: $PPType\n");
	$App->print("Functional           : $Functional\n");
	$App->print("\n");

	my $PWSCF = new PWSCF;
	my ($pPPHashArray, $pPPTypes, $pFunctionals) = $PWSCF->ReadAllPP($Debug);
}

sub SummarizeStatus
{
	my $InpFile         = $Args->GetGetArg(0);
	my $WoutFile        = $Args->GetGetArg(1);
	my $SummaryFile     = $Args->GetGetArg(2);

	$App->print("\n");
	$App->print("Input (.winp) File : $InpFile\n");
	$App->print("Input (.wout) File : $WoutFile\n");
	$App->print("Summary File       : $SummaryFile\n");
	$App->print("\n");
	
	my $PWSCF = new PWSCF;
print "Read from [$WoutFile]\n";
	my $pOut  = $PWSCF->ReadOutputToHash($WoutFile, 0);
print "Write summary to [$SummaryFile]\n";
	my $out = JFile->new($SummaryFile, 'w');
	if(!$out) {
		$App->print("\nError in SummarizeStatus: Can not write to [$SummaryFile].\n\n");
		exit;
	}

	$out->mprint("  EF   = $pOut->{EF} Ry\n");
	$out->mprint("  Etot = $pOut->{Etot} Ry\n");
	if($pOut->{RelaxationStatus} =~ /Error/i) {
		$out->mprint("  Relaxation $pOut->{RelaxationStatus}\n");
		return 1;
	}
	if($pOut->{SCFStatus} =~ /Error/i) {
		$out->mprint("  SCF $pOut->{SCFStatus}\n");
		return 2;
	}
	return 0;
}

sub MakeCIFFilesFromOutput
{
	my $InpFile         = $Args->GetGetArg(0);
	my $WoutFile        = $Args->GetGetArg(1);
	my $CIFInitialFile  = $Args->GetGetArg(2);
	my $CIFFinalFile    = $Args->GetGetArg(3);

	$App->print("\n");
	$App->print("Input (.winp) File : $InpFile\n");
	$App->print("Input (.wout) File : $WoutFile\n");
	$App->print("Initial CIF File   : $CIFInitialFile\n");
	$App->print("Final CIF File     : $CIFFinalFile\n");
	$App->print("\n");

	$App->print("\nDelete [$CIFInitialFile].");
	unlink($CIFInitialFile);
	$App->print("\nDelete [$CIFFinalFile].");
	unlink($CIFFinalFile);

	printf("\n");

	my $PWSCF      = new PWSCF;
	my $Crystal    = new Crystal;
	my $CIF        = new CIF;

	my $IniCrystal = $PWSCF->ReadCrystalStructureFromInput($App, $InpFile, $Crystal);
	$App->printf("Make CIFFile of initial structure [%s].\n", $CIFInitialFile);
	$CIF->CreateCIFFileFromCCrystal($IniCrystal, $CIFInitialFile, 0, 'unix');

	my $CIFHeader = $WoutFile;
	$CIFHeader =~ s/\.*$//;

	my $nStep     = $PWSCF->GetnFrameInOutput($App, $WoutFile);
	$App->print("nStep in [$WoutFile]: $nStep\n");

	$App->printf("\n");
	my $in  = JFile->new($WoutFile, "r") or die "$!: Can not read [$WoutFile].\n";
#	$PWSCF->ReadNextCrystalStructureFromOutput($App, $in, 0, $Crystal);
#	$App->printf("...Frame %d (D=%lg)\n", 0, $Crystal->Density());
#	$App->printf("Make CIFFile of initial structure [%s].\n", $CIFInitialFile);
#	$CIF->CreateCIFFileFromCCrystal($Crystal, $CIFInitialFile, 0, 'unix');

	for(my $i = 0 ; $i < $nStep ; $i++)
#	for(my $i = 1 ; $i < $nStep ; $i++)
	{
		last if(!$PWSCF->ReadNextCrystalStructureFromOutput($App, $in, $i, $Crystal));

		$App->printf("...Frame %d (D=%lg)\n", $i, $Crystal->Density());
	}
	$in->Close();

	$App->printf("Make CIFFile of final structure [%s].\n", $CIFFinalFile);
	$CIF->CreateCIFFileFromCCrystal($Crystal, $CIFFinalFile, 0, 'unix');
	
	$App->print("\n");
}

sub MakeHistoryCSV
{
	print "<b>Make History CSV File:</b>$LF";

	my $InpFile          = $Args->GetGetArg(0);
	my $WoutFile         = $Args->GetGetArg(1);
	my $TotECSVFile      = $Args->GetGetArg(2);
	my $StructureCSVFile = $Args->GetGetArg(3);

	my $nLastStep        = $Args->GetGetArg('nLastStep');
	my $nOutputInterval  = $Args->GetGetArg('nOutputInterval');
	   $nOutputInterval  = 1 if($nOutputInterval <= 0);
	my $iAvrFirstStep    = $Args->GetGetArg('iAvrFirstStep');
	   $iAvrFirstStep    = 1 if($iAvrFirstStep <= 0);
	my $iAvrLastStep     = $Args->GetGetArg('iAvrLastStep');

	$App->print("\n");
	$App->print("Input File        : $InpFile\n");
	$App->print("Output File       : $WoutFile\n");
	$App->print("Energy CSV File   : $TotECSVFile\n");
	$App->print("Structure CSV File: $StructureCSVFile\n");

	$App->print("nLastStep          : $nLastStep\n");
	$App->print("nOutputInterval    : $nOutputInterval\n");
	$App->print("iAvrFirstStep      : $iAvrFirstStep\n");
	$App->print("iAvrLastStep       : $iAvrLastStep\n");
	$App->print("\n");

	$App->print("\nDelete [$TotECSVFile].\n");
	unlink($TotECSVFile);
	$App->print("\nDelete [$StructureCSVFile].\n");
	unlink($StructureCSVFile);

	my $PWSCF   = new PWSCF;
	my $Crystal = new Crystal;
	$PWSCF->ReadCrystalStructureFromInput($App, $InpFile, $Crystal);

	my $nStep     = $PWSCF->GetnFrameInOutput($App, $WoutFile);
	$App->print("nStep in [$WoutFile]: $nStep\n");
	$nStep        = $nLastStep if($nLastStep > 0 && $nLastStep < $nStep);
	$App->print("nStep to be saved to [$StructureCSVFile]: $nStep\n");
	my $nOutput   = ($nStep-1) / $nOutputInterval + 1;
	$iAvrLastStep = $nStep if($iAvrLastStep <= 0 and $iAvrLastStep > $nStep);

	my @AtomTypeList         = $Crystal->GetCAtomTypeList();
	my $nAtomType            = @AtomTypeList;
	my @AtomSiteList         = $Crystal->GetCAsymmetricAtomSiteList();
	my $nAtomSite            = @AtomSiteList;
#	my @ExpandedAtomSiteList = $Crystal->GetCExpandedAtomSiteList();
#	my $n = @AtomSiteList;

	$App->printf("nAtomSites: %d  nAtomType: %d\n", $nAtomSite, $nAtomType);
	$App->print("nStep          : $nStep\n");
	$App->print("nOutputInterval: $nOutputInterval\n");
	$App->print("nOutput        : $nOutput\n");

	my $in   = JFile->new($WoutFile, "r")          or die "$!: Can not read [$WoutFile].\n";
	my $out  = JFile->new($StructureCSVFile, "wb") or die "$!: Can not write to [$StructureCSVFile].\n";
	my $out2 = JFile->new($TotECSVFile, "wb")      or die "$!: Can not write to [$TotECSVFile].\n";
	$out->print("iStep,TotE(Ry),a(A),b,c,alpha,beta,gamma,V,density(g/cm3)\n");
	$out2->print("iStep,iIter,TotE(Ry)\n");

	my $count = 0;
	my $AvrCount = 0;
	for(my $i = 0 ; $i < $nStep ; $i++)
	{
		my $pHash = $PWSCF->ReadNextCrystalStructureFromOutput($App, $in, $i, $Crystal);
		last if(!$pHash);

		my $pE = $pHash->{pEConvergence};
		for(my $j = 0 ; $j < @$pE ; $j++) {
			$out2->print("$i,$j,$pE->[$j]\n");
		}
	
		if($i % $nOutputInterval == 0) {
			$count++;
			my ($a, $b, $c, $alpha, $beta, $gamma) = $Crystal->LatticeParameters();
			my $V = $Crystal->Volume();
			my $D = $Crystal->Density();
			$out->print("$i,$pHash->{TotalEnergy},$a,$b,$c,$alpha,$beta,$gamma,$V,$D\n");
			$App->printf("...Frame %d (D=%lg)\n", $i, $Crystal->Density());
#			printf("Step %d / %d: Saved.\n", i, nStep);
		}
		else {
#			printf("Step %d / %d\n", i, nStep);
		}
	}

	$in->Close();
	$out->Close();
	$out2->Close();

	return 1;
}

sub MakeXSF
{
	print "<b>Make XSF File:</b>$LF";

	my $InpFile  = $Args->GetGetArg(0);
	my $WoutFile = $Args->GetGetArg(1);
	my $XSFFile  = $Args->GetGetArg(2);

	my $nLastStep       = $Args->GetGetArg('nLastStep');
	my $nOutputInterval = $Args->GetGetArg('nOutputInterval');
	   $nOutputInterval = 1 if($nOutputInterval <= 0);
	my $iAvrFirstStep   = $Args->GetGetArg('iAvrFirstStep');
	   $iAvrFirstStep   = 1 if($iAvrFirstStep <= 0);
	my $iAvrLastStep    = $Args->GetGetArg('iAvrLastStep');

	$App->print("\n");
	$App->print("Input File : $InpFile\n");
	$App->print("Output File: $WoutFile\n");
	$App->print("XSF File   : $XSFFile\n");
	$App->print("nLastStep          : $nLastStep\n");
	$App->print("nOutputInterval    : $nOutputInterval\n");
	$App->print("iAvrFirstStep      : $iAvrFirstStep\n");
	$App->print("iAvrLastStep       : $iAvrLastStep\n");
	$App->print("\n");

	$App->print("\nDelete [$XSFFile].\n");
	unlink($XSFFile);

	my $PWSCF   = new PWSCF;
	my $Crystal = new Crystal;
	$PWSCF->ReadCrystalStructureFromInput($App, $InpFile, $Crystal);

	my $nStep     = $PWSCF->GetnFrameInOutput($App, $WoutFile);
	$App->print("nStep in [$WoutFile]: $nStep\n");
	$nStep        = $nLastStep if($nLastStep > 0 && $nLastStep < $nStep);
	$App->print("nStep to be saved to [$XSFFile]: $nStep\n");
	my $nOutput   = ($nStep-1) / $nOutputInterval + 1;
	$iAvrLastStep = $nStep if($iAvrLastStep <= 0 and $iAvrLastStep > $nStep);

	my @AtomTypeList         = $Crystal->GetCAtomTypeList();
	my $nAtomType            = @AtomTypeList;
	my @AtomSiteList         = $Crystal->GetCAsymmetricAtomSiteList();
	my $nAtomSite            = @AtomSiteList;
#	my @ExpandedAtomSiteList = $Crystal->GetCExpandedAtomSiteList();
#	my $n = @AtomSiteList;

	$App->printf("nAtomSites: %d  nAtomType: %d\n", $nAtomSite, $nAtomType);
	$App->print("nStep          : $nStep\n");
	$App->print("nOutputInterval: $nOutputInterval\n");
	$App->print("nOutput        : $nOutput\n");

	my $xc = new XCrySDen;
	$xc->WriteXSFAnimationFileHeader($XSFFile, $nStep);

	my $in  = JFile->new($WoutFile, "r") or die "$!: Can not read [$WoutFile].\n";
	my $out = JFile->new($XSFFile, "wb") or die "$!: Can not write to [$XSFFile].\n";

	my $count = 0;
	my $AvrCount = 0;
	for(my $i = 0 ; $i < $nStep ; $i++)
	{
		last if(!$PWSCF->ReadNextCrystalStructureFromOutput($App, $in, $i, $Crystal));

		if($i % $nOutputInterval == 0) {
			$count++;
			$xc->WriteXSFFileFromCrystal($XSFFile, $count, $Crystal);
			$App->printf("...Frame %d (D=%lg)\n", $i, $Crystal->Density());
#			printf("Step %d / %d: Saved.\n", i, nStep);
		}
		else {
#			printf("Step %d / %d\n", i, nStep);
		}
	}

	$in->Close();
	$out->Close();

	return 1;
}

sub UpdateInput
{
	my $InpFile     = $Args->GetGetArg(0);
	my $WoutFile    = $Args->GetGetArg(1);
	my $NewInpFile  = $Args->GetGetArg(2);

	$App->print("\n");
	$App->print("Input (.winp) File : $InpFile\n");
	$App->print("Input (.wout) File : $WoutFile\n");
	$App->print("New Input File     : $NewInpFile\n");
	$App->print("\n");

	$App->print("\nDelete [$NewInpFile].");
	unlink($NewInpFile);
	printf("\n");

	my $PWSCF      = new PWSCF();
	my $Crystal    = new Crystal;
	$PWSCF->ReadCrystalStructureFromInput($App, $InpFile, $Crystal);

	my $nStep     = $PWSCF->GetnFrameInOutput($App, $WoutFile);
	$App->print("nStep in [$WoutFile]: $nStep\n");

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

		$App->printf("...Frame %d (D=%lg)\n", $i, $Crystal->Density());
	}
	$in->Close();

	$App->printf("Create update INP file: %s\n", $NewInpFile);
	$PWSCF->UpdateInpFile($App, $InpFile, $NewInpFile, $Crystal);

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

sub MakeInputFiles
{
	print "<b>Make PWSCF Input Files:</b>$LF";

	my $IsPrint = 1;

	my $CIFFile           = $Args->GetGetArg(0);
	if($CIFFile eq 'auto') {
		my @files = glob("*.cif");
		@files = glob("*.CIF") if(@files == 0);
		for(my $i = 0 ; $i < @files ; $i++) {
			next if($files[$i] =~ /-initial\.cif$/i or $files[$i] =~ /-final\.cif$/i);
			$CIFFile = $files[$i];
			last;
		}
		if($CIFFile eq 'auto') {
			for(my $i = 0 ; $i < @files ; $i++) {
				if($files[$i] =~ /-final\.cif$/i) {
					$CIFFile = $files[$i];
					last;
				}
			}
		}
		if($CIFFile eq 'auto') {
			$CIFFile = $files[0];
		}
	}

	my $InpFile           = $Args->GetGetArg(1);
	my $BurstToP1         = $Args->GetGetArg('BurstToP1');
# scf nscf relax vc-relax md phonon
	my @Function          = $Args->GetGetArg('Function');
	my $Function          = $Function[0];
	my $PPType            = $Args->GetGetArg('PPType');
	my $Functional        = $Args->GetGetArg('Functional');
	my $PPCopyMode        = $Args->GetGetArg('PPCopyMode');
	$PPCopyMode           = 0 if(!defined $PPCopyMode);
	my $UseBravaisLattice = $Args->GetGetArg('UseBravaisLattice');
	$UseBravaisLattice    = 0 if(!defined $UseBravaisLattice);
	my $UseCELLDM         = $Args->GetGetArg('UseCELLDM');
	$UseCELLDM            = 1 if(!defined $UseCELLDM);
	my @PP                = $Args->GetGetArg('PseudoPotential');
	my $PP                = $PP[0];
	my $SpinPolarized     = $Args->GetGetArg('SpinPolarized');
	$SpinPolarized        = 1 if(!defined $SpinPolarized);
	my $aKProduct         = $Args->GetGetArg('aKProduct');
	$aKProduct            = 1.5 if(!defined $aKProduct);
	my $ecutwfc           = $Args->GetGetArg('ecutwfc');
	$ecutwfc              = 'x1.0,30.0', if(!defined $ecutwfc);
	my $ecutrho           = $Args->GetGetArg('ecutrho');
	$ecutrho              = 'x1.0,140.0' if(!defined $ecutrho);
	my $nbnd              = $Args->GetGetArg('nbnd');
	$nbnd                 = 'x1.2' if(!defined $nbnd);
	my $nk                = $Args->GetGetArg('nk');
	$nk                   = 10 if(!defined $nk);
	my $nstep             = $Args->GetGetArg('nstep');
	$nstep                = 50 if(!defined $nstep);
	my $press             = $Args->GetGetArg('press');
	$press                = 0.001 if(!defined $press);
	my $etot_conv_thr    = $Args->GetGetArg('etot_conv_thr');
	$etot_conv_thr       = 0.5e-4 if(!defined $etot_conv_thr);
	my $forc_conv_thr    = $Args->GetGetArg('forc_conv_thr');
	$forc_conv_thr       = 0.5e-3 if(!defined $forc_conv_thr);
	my $conv_thr         = $Args->GetGetArg('conv_thr');
	$conv_thr            = 0.5e-6 if(!defined $conv_thr);
	my $press_conv_thr    = $Args->GetGetArg('press_conv_thr');
	$press_conv_thr       = 0.001 if(!defined $press_conv_thr);
	my $KPoints           = $Args->GetGetArg('KPoints');
	$KPoints              = 'File' if(!defined $KPoints);

	$App->print("\n");
	$App->print("CIF File              : $CIFFile\n");
	$App->print("Burst to P1           : $BurstToP1\n");
	$App->print("Input File to be saved: $InpFile\n");
	$App->print("Function              : $Function\n");
	$App->print("Pseudo potential type : $PPType\n");
	$App->print("Functional            : $Functional\n");
	$App->print("Pseudo Potential      : $PP\n");
	$App->print("Use Bravais lattice   : $UseBravaisLattice\n");
	$App->print("SpinPolarized         : $SpinPolarized\n");
	$App->print("aKProduct             : $aKProduct\n");
	$App->print("ecutwfc               : $ecutwfc Ry\n");
	$App->print("ecutrho               : $ecutrho Ry\n");
	$App->print("nbnd                  : $nbnd\n");
	$App->print("nk                    : $nk\n");
	$App->print("nstep                 : $nstep\n");
	$App->print("press                 : $press kBar\n");
	$App->print("etot_conv_thr         : $etot_conv_thr\n");
	$App->print("forc_conv_thr         : $forc_conv_thr\n");
	$App->print("conv_thr              : $conv_thr\n");
	$App->print("\n");

	$App->print("\nDelete [$InpFile].");
	unlink($InpFile);
	$App->print("\n");

	my $Crystal = &GetCrystalFromCIF($CIFFile);
	my $PWSCF   = new PWSCF();

	my $IsChooseRandomly = 0;
	my $SampleName = &MakeSampleNameFromPath($CIFFile);
	$PWSCF->SetSampleName($SampleName);
	my @Files = $PWSCF->SavePWSCFInpFile($Function, $PP, $Crystal, $InpFile, 
						$IsChooseRandomly, $UseBravaisLattice, $UseCELLDM, 
						$PPType, $Functional, $PPCopyMode, $PPSaveDir, $IsPrint,
						BurstToP1      => $BurstToP1,
						PrintForFindPP => 0,
						aKProdct       => $aKProduct,
						SpinPolarized  => $SpinPolarized,
						ecutwfc        => $ecutwfc,
						ecutrho        => $ecutrho,
						nbnd           => $nbnd,
						KPoints        => $KPoints,
						nk             => $nk,
						nstep          => $nstep,
						press          => $press,
						etot_conv_thr  => $etot_conv_thr,
						forc_conv_thr  => $forc_conv_thr,
						conv_thr       => $conv_thr,
						press_conv_thr => $press_conv_thr,
						);

	return;
}