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

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

use strict;
use File::Basename;

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

#use Crystal::MyUtility;
use Crystal::Crystal;
use Crystal::AtomType;
use Crystal::AtomSite;
use Crystal::CIF;
use Crystal::VASP;
use Crystal::BandStructure;

#my $ATKPPDir   = "D:\\Programs\\Perl\\VNL\\pseudopotentials\\atk2.0";
my $PWSCFDir   = ProgVars::QEDir(); #Deps::MakePath($TkPerlDir, "QE", 0);
my $PWSCFPPDir = ProgVars::QEPPDir();
my $PWSCFPPURL = "http://www.pwscf.org/pseudo/1.3/html/{Atom}.html";

my @PPPriority = (
	'pbe-sp-kjpaw',
	'pbe-sp',
	'pbe-kjpaw',
	'pbe-(n-|nsp-)?kjpaw',
	'pbe-.*kjpaw',
	'pbe(-n)?-rrkjus',
	'pbe-.*rrkjus',
	'pbe-.*van_ak',
	'pbe-.*van_bm',
	'pbe-.*van',
	'pbe-mt_fhi',
	'pbe-nsp',
	'pbe-',
	'rel-pbesol(-n)?-kjpaw',
	'rel-pbesol(-n)?-rrkjus',
	'pw91',
	'pz-(-n)?kjpaw',
	'pz-(-n)?rrkjus',
	'pz-mt',
	'pz-rrkjus',
	'pz-van',
	'blyp-(sp-?)van',
	'blyp-mt',
	'tpss-mt',
	);

BEGIN {}

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

sub DESTROY
{
	my $this = shift;
}

sub SetSampleName
{
	my ($this, $name) = @_;
	return $this->{'SampleName'} = $name;
}

sub SampleName
{
	my ($this) = @_;
	return $this->{'SampleName'};
}

sub GetiPWSCFBravaisLattice
{
	my ($this, $strCrystalBravaisLattice) = @_;
	$strCrystalBravaisLattice = uc $strCrystalBravaisLattice;
	my %iBravaisLattice = (
			"FREELATTICE"      =>  0, # Free Lattice
			"P(CUBIC)"         =>  1, # Cubic P (sc)
			"FC(CUBIC)"        =>  2, # Cubic F (fcc)
			"BC(CUBIC)"        =>  3, # Cubic I (bcc)
			"P(HEXAGONAL)"     =>  4, # Hexagonal and Trigonal P
			"P(TRIGONAL)"      =>  4, # Hexagonal and Trigonal P
			"R(RHOMBOHEDRAL)"  =>  5, # Trigonal R
			"R(TRIGONAL)"      =>  5, # Trigonal R
			"P(TETRAGONAL)"    =>  6, # Tetragonal P (st)
			"BC(TETRAGONAL)"   =>  7, # Tetragonal I (bct)
			"P(ORTHORHOMBIC)"  =>  8, # Orthorhombic P
#			"A(ORTHORHOMBIC)"  => -9, # Orthorhombic base-centered(bco)
			"C(ORTHORHOMBIC)"  =>  9, # Orthorhombic base-centered(bco)
			"FC(ORTHORHOMBIC)" => 10, # Orthorhombic face-centered
			"BC(ORTHORHOMBIC)" => 11, # Orthorhombic body-centered
			"P(MONOCLINIC)"    => 12, # Monoclinick P
			"C(MONOCLINIC)"    => 13, # Monoclinick base-centered
			"P(TRICLINIC)"     => 14, # Triclinick P
		);
#foreach my $key (keys %iBravaisLattice) {
#	print "[$key]: $iBravaisLattice{$key}<br>\n";
#}

#print "[$strCrystalBravaisLattice] => $iBravaisLattice{$strCrystalBravaisLattice}<br>\n";

	my $ibrav = abs($iBravaisLattice{$strCrystalBravaisLattice});
	if($ibrav >= 12) {
		$ibrav = 14;
	}
	elsif($ibrav == 2 or $ibrav == 3) {
		$ibrav = 1;
	}
	elsif($ibrav == 5) {
		$ibrav = 4;
	}
	elsif($ibrav == 7) {
		$ibrav = 6;
	}
	elsif($ibrav >= 9 and $ibrav <= 11) {
		$ibrav = 8;
	}
	return $ibrav;
}

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

	print "\nCalculate total number of valence electrons\n";
	
#	my $pPPFiles = $this->FindPPFiles(\@AtomTypeList, $Functional, $PPType, $IsPrint);
#	my @AtomTypeList         = $Crystal->GetCAtomTypeList();
#	my $nAtomType            = @AtomTypeList;
	my @AtomSiteList         = $Crystal->GetCAsymmetricAtomSiteList();
	my @nMultiplicityExpandedAtomSiteList 
		= $Crystal->GetCnMultiplicityExpandedAtomSiteList();

	my %nVal;
	my $nTotEVal = 0;
	for(my $i = 0 ; $i < @AtomSiteList ; $i++) {
		my $atomname = $AtomSiteList[$i]->AtomNameOnly();
		my $PPFile = $pPPFiles->{$atomname};
		my $pHash = $this->ReadPPtoHash(Deps::MakePath($PWSCFPPDir, $PPFile, 0), 0);
#print "pHash=$pHash\n";
print "  $i: $atomname: $pHash->{z_valence} in [$PPFile]\n";
		$nTotEVal += $nMultiplicityExpandedAtomSiteList[$i] * $pHash->{z_valence};
#$hash{wfc_cutoff} = $1;
#$hash{rho_cutoff} = $1;
	}
print "  Total number of valence electrons: $nTotEVal\n";

	return $nTotEVal;
}

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

	print "\nCalculate Ecut for base functions (wfc_cutoff) and electron density (rho_cutoff)\n";

	my @AtomSiteList         = $Crystal->GetCAsymmetricAtomSiteList();

	my ($MaxEcut, $MaxRhoEcut) = (0.0, 0.0);
	for(my $i = 0 ; $i < @AtomSiteList ; $i++) {
		my $atomname = $AtomSiteList[$i]->AtomNameOnly();
		my $PPFile = $pPPFiles->{$atomname};
		my $pHash = $this->ReadPPtoHash(Deps::MakePath($PWSCFPPDir, $PPFile, 0), 0);
print "  $i: $atomname: wfc_Ecut = $pHash->{wfc_cutoff}, rho_Ecut = $pHash->{rho_cutoff} in [$PPFile]\n";
		$MaxEcut    = $pHash->{wfc_cutoff} if($MaxEcut    < $pHash->{wfc_cutoff});
		$MaxRhoEcut = $pHash->{rho_cutoff} if($MaxRhoEcut < $pHash->{rho_cutoff});
	}

	return ($MaxEcut, $MaxRhoEcut);
}

sub ReadPPtoHash
{
	my ($this, $fname, $DebugMode) = @_;
	$DebugMode = 0 if(!defined $DebugMode);

	my ($drive, $directory, $filename, $ext1, $lastdir, $filebody) 
		= Deps::SplitFilePath($fname);
print "Read PP from [$filename]\n" if($DebugMode);

	my %hash = (path => $fname, filename => $filename);
	my $in = JFile->new($fname, 'r');
	if(!$in) {
		print "Can not read [$fname].\n" if($DebugMode);
		return 0;
	}

	my $line = $in->SkipTo("Generated ");
	($hash{method}) = ($line =~ /using\s+(\S.*),/);
	($hash{method}) = ($line =~ /by\s+(\S.*),/)          if($hash{method} eq '');
	($hash{method}) = ($line =~ /using\s+(\S+)/)         if($hash{method} eq '');
	($hash{method}) = ($line =~ /by\s+(\S+)/)            if($hash{method} eq '');
	($hash{method}) = ($line =~ /UPF file from\s+(\S+)/) if($hash{method} eq '');
	Utils::DelQuote($hash{method});

	$line = $in->SkipTo("Pseudopotential type:");
	($hash{PPType}) = ($line =~ /type:\s*(\S.+)\s*$/);
	Utils::DelQuote($hash{PPType});
	if($hash{PPType} eq '') {
		$hash{PPType} = 'PAW' if($hash{method} =~ /paw/i);
	}
	$hash{PPType} = 'USPP' if($hash{PPType} eq 'US');

	$line = $in->SkipTo("Element:");
	($hash{Element}) = ($line =~ /Element:\s*(\S.*)\s*$/);
	if($hash{Element} eq '') {
		$in->rewind();
		for(my $i = 0 ; $i < 20 ; $i++) {
			$line = $in->ReadLine();
			if($line =~ /^\s*([A-Z][a-z]?)\s/) {
				$hash{Element} = $1;
#print "e:$hash{Element}\n";
#exit;
				last;
			}
		}
	}
	Utils::DelQuote($hash{Element});
	$hash{Element} = ucfirst(lc($hash{Element}));

	$line = $in->SkipTo("Functional:");
	($hash{Functional}) = ($line =~ /Functional:\s*(\S.+)\s*$/);
	if($hash{Functional} eq '') {
		$in->rewind();
		for(my $i = 0 ; $i < 20 ; $i++) {
			$line = $in->ReadLine();
			if($line =~ /^\s*GGA-([A-Za-z0-9]+)\s/) {
				$hash{Functional} = $1;
#print "e:$hash{Element}\n";
#exit;
				last;
			}
		}
	}
	Utils::DelQuote($hash{Functional});

	while(1) {
		$line = $in->ReadLine();
		last if(!$line);
		last if($line =~ /\<PP_MESH/i);
		last if($line =~ /\<PP_R type=/i);

		if($line =~ /wfc_cutoff=\s*"?(.*?)"?\s*$/i) {
			$hash{wfc_cutoff} = $1;
			Utils::DelQuote($hash{wfc_cutoff});
		}
		elsif($line =~ /rho_cutoff=\s*"?(.*?)"?\s*$/i) {
			$hash{rho_cutoff} = $1;
			Utils::DelQuote($hash{rho_cutoff});
		}
		elsif($line =~ /z_valence=\s*"?(.*?)"?\s*$/i) {
			$hash{z_valence} = $1;
			Utils::DelQuote($hash{z_valence});
		}
	}
#print "  [$hash{Element}] by [$hash{method}]  PP=[$hash{PPType}]  Func=[$hash{Functional}]\n";
#print "    z_valence=$hash{z_valence}  wfc_cutoff=$hash{wfc_cutoff}  rho_cutoff=$hash{rho_cutoff}\n";
#exit;

	if($DebugMode) {
print "  [$hash{Element}] by [$hash{method}]  PP=[$hash{PPType}]  Func=[$hash{Functional}]\n";
		if($hash{Element} eq '' or $hash{method} eq '' 
		   or $hash{PPType} eq '' or $hash{Functional} eq '') {
			$in->rewind();
			for(my $j = 0 ; $j < 7 ; $j++) {
				$line = $in->ReadLine();
				exit if($line =~ /^Valence/);
print "  $line";
			}
exit;
		}
#exit if($files[$i] =~ /Ga\.pbe.*nsp/);
	}

	$in->Close();

	return \%hash;
}

sub ReadAllPP
{
	my ($this, $Debug) = @_;
	$Debug = 0 if(!defined $Debug);

	print "Read all pseudo potentials from [$PWSCFPPDir]\n";

	my (@PPHashArray, %PPTypeHash, %FunctionalHash, %ElementHash, %PPHash);
	my @files = sort glob(Deps::MakePath($PWSCFPPDir, "*.UPF", 0));
	for(my $i = 0 ; $i < @files ; $i++) {
#last if($i >= 90);
print "  Read [$files[$i]]\n" if($Debug);
		$PPHashArray[$i] = $this->ReadPPtoHash($files[$i], $Debug);

		my $Element    = $PPHashArray[$i]->{Element};
		my @PPType     = Utils::Split("\\s+", $PPHashArray[$i]->{PPType});
		my @Functional = Utils::Split("\\s+", $PPHashArray[$i]->{Functional});
		$ElementHash{$Element}++;
		for(my $i = 0 ; $i < @PPType ; $i++) {
			$PPTypeHash{$PPType[$i]}++;
		}
		for(my $i = 0 ; $i < @Functional ; $i++) {
			$FunctionalHash{$Functional[$i]}++;
		}
	}

	my @ElementArray    = sort keys %ElementHash;
	my @PPTypeArray     = sort keys %PPTypeHash;
	my @FunctionalArray = sort keys %FunctionalHash;
if($Debug) {
	print "\n";
	print "PPTypes: ";
	for(my $i = 0 ; $i < @PPTypeArray ; $i++) {
		print "$PPTypeArray[$i] ";
	}
	print "\n";

	print "Functionals: ";
	for(my $i = 0 ; $i < @FunctionalArray ; $i++) {
		print "$FunctionalArray[$i] ";
	}
	print "\n";
}

	for(my $l = 0 ; $l < @PPHashArray ; $l++) {
		my $p           =  $PPHashArray[$l];
		my $e           = $p->{Element};
		my @PPTypes     = Utils::Split("\\s+", $p->{PPType});
		my @Functionals = Utils::Split("\\s+", $p->{Functional});

		my %ph;
		for(my $i = 0 ; $i < @PPTypes ; $i++) {
			my $pptype = $PPTypes[$i];
			next if($ph{$pptype});
			$ph{$pptype}++;
			my %fh;
			for(my $j = 0 ; $j <@Functionals ; $j++) {
				my $f = $Functionals[$j];
				next if($fh{$f});
				$fh{$f}++;
				$PPHash{$e}{$pptype}{$f} = [] if(!defined $PPHash{$e}{$pptype}{$f});
				my $p = $PPHash{$e}{$pptype}{$f};
#print "Add $PPHashArray[$l] to [$e]$[pptype][$f]\n" if($e eq 'Ba');
				push(@$p, $PPHashArray[$l]);
			}
		}
	}

	return (\@PPHashArray, \@PPTypeArray, \@FunctionalArray, \%PPHash);
}

sub FindPPFiles
{
	my ($this, $pAtomTypeList, $TargetFunctional, $TargetPPType, $IsPrint) = @_;
	$TargetFunctional = 'PBE' if(!defined $TargetFunctional);
	$TargetPPType     = 'PAW' if(!defined $TargetPPType);
	$IsPrint          = 1 if(!defined $IsPrint);

	my @AtomTypeList = @$pAtomTypeList;

	my ($pPPHashArray, $pPPTypes, $pFunctionals, $pPPHash) = $this->ReadAllPP($IsPrint);

	print "\n";
	for(my $i = 0 ; $i < @AtomTypeList ; $i++) {
		my $atomname = $AtomTypeList[$i]->AtomName();
		my $p = $pPPHash->{$atomname};
		if(!defined $p) {
			print "Error: No pseudo potential is available for [$atomname].\n";
			last;
		}
#print "$i [$atomname][$PPType][$Functional]: p=$p ";
		print "Available pseudo potentials for [$atomname]:\n";
		foreach my $PPType (sort keys %$p) {
			my $p2 = $p->{$PPType};
print "  PPType=$PPType Functionals=";
#			if(!defined $p2) {
#				print "Error: No pseudo potential type is available for [$atomname].\n";
#				last;
#			}
#print "pp2=$p2 ";
#print "hash=", %$p2, "\n";
			foreach my $Functional (sort keys %$p2) {
print "$Functional  ";
				my $p3 = $p2->{$Functional};
#				if(!defined $p3) {
#					print "Error: No pseudo potential for functional [$Functional] is available.\n";
#					last;
#				}
#print "pp3=$p3 [$PPType][$Functional]";
			}
print "\n";
		}
	}

print "PPT:$TargetPPType:";
	my @PPTypes = Utils::Split(",", $TargetPPType);
print join(', ', @PPTypes), "\n";
	for(my $i = 0 ; $i < @PPTypes ; $i++) {
		my $pptype = $PPTypes[$i];
#print "$i: $ppt\n";
		print "\n";
		print "Available pseudo potentials for [$pptype][$TargetFunctional]\n";
		my ($pPPToBeUsed, $atomname, $ppt2) 
			= $this->FindPPFilesByPPType($pPPHash, $pAtomTypeList, $TargetFunctional, $pptype, $IsPrint);
		if($pPPToBeUsed == 0) {
			print "  Error: No pseudo potential is available for [$atomname][$ppt2][$TargetFunctional].\n";
			next;
		}
		else {
			my $MatchedIndex = $atomname;
			print " ***All PPs [$pptype][$TargetFunctional] are available for "
			     ."[PPPriority=$PPPriority[$MatchedIndex]]\n";# if($IsPrint);
			foreach my $atomname (sort keys %$pPPToBeUsed) {
				printf "  %2s: $pPPToBeUsed->{$atomname}\n", $atomname;
			}
			return $pPPToBeUsed;
		}
	}
}

sub FindPPFilesByPPType
{
	my ($this, $pPPHash, $pAtomTypeList, $TargetFunctional, $TargetPPType, $IsPrint) = @_;
	my @AtomTypeList = @$pAtomTypeList;

#print "\n";
#print "Enter PWSCF::FindPPFilesByPPType\n";
#print "Available pseudo potentials for [$TargetPPType][$TargetFunctional]\n";
	for(my $i = 0 ; $i < @AtomTypeList ; $i++) {
		my $atomname = $AtomTypeList[$i]->AtomName();
		my $p = $pPPHash->{$atomname};
		if(!defined $p) {
			return (0, $atomname);
		}
		my $p2 = $p->{$TargetPPType};
		if(!defined $p2) {
			print "  Error in PWSCF::FindPPFilesByPPType: "
			     ."No pseudo potential is available for [$atomname][$TargetPPType].\n";
			return (0, $atomname, $TargetPPType);
		}
		my $p3 = $p2->{$TargetFunctional};
		if(!defined $p3) {
			print "  Error in PWSCF::FindPPFilesByPPType: "
			     ."No pseudo potential is available for [$atomname][$TargetPPType][$TargetFunctional].\n";
			return (0, $atomname, $TargetPPType, $TargetFunctional);
#			exit;
		}
#print "  [$atomname]: ", join(', ', @$p3), "\n";
		for(my $i = 0 ; $i < @$p3 ; $i++) {
			my $p4 = $p3->[$i];
			print "     $p4->{filename}\n";
			$p->{PPCandidates} = [] if(!defined $p->{PPCandidates});
			my $ppc = $p->{PPCandidates};
			push(@$ppc, $p4->{filename});
			$p->{PPCandidates} = $ppc;
		}
	}

#	print "\n";
#	print "PP files priority:\n";
	my %PPToBeUsed;
	my $MatchedIndex = -1;
	for(my $i = 0 ; $i < @PPPriority ; $i++) {
		print "  $i: $PPPriority[$i]\n";
		my $AllPassed = 1;
		%PPToBeUsed = ();
		for(my $j = 0 ; $j < @AtomTypeList ; $j++) {
			my $atomname = $AtomTypeList[$j]->AtomName();
			my $p = $pPPHash->{$atomname};
			my $ppc = $p->{PPCandidates};
			my $AtomPassed = 0;
			for(my $k = 0 ; $k < @$ppc ; $k++) {
				my $fname = $ppc->[$k];
				if($fname =~ /$PPPriority[$i]/i) {
					print "     check $PPPriority[$i] vs $fname: Passed\n";
					$PPToBeUsed{$atomname} = $fname;
					$AtomPassed = 1;
					last;
				}
			}
			$AllPassed = 0 if(!$AtomPassed);
		}

		if($AllPassed) {
			$MatchedIndex = $i;
#			print " ***All PPs are available for $PPPriority[$i]\n" if($IsPrint);
			last;
		}
	}

#	print "\n";
#	print "PP files to be used:\n";
#	foreach my $atomname (sort keys %PPToBeUsed) {
#		printf "  %2s: $PPToBeUsed{$atomname}\n", $atomname;
#	}

	return (\%PPToBeUsed, $MatchedIndex);
}

sub MakeKPOINTSForBand
{
	my ($this, $Crystal, $klistfile, $KPoints, $nk) = @_;

#klistファイル読み込み
print "PWSCF::MakeKPOINTSForBand: *Read k points from [$klistfile]\n";
	my @klist;
	if($klistfile) {
		@klist = BandStructure->ReadKListFile($klistfile);
	}
	else {
		for(my $i = 0 ; $i < length($KPoints) ; $i++) {
			my $c0 = substr($KPoints, $i, 1);
			if($c0 =~ /G/i) {
				push(@klist, {name => 'GAMMA', kx => 0.0, ky => 0.0, kz => 0.0});
			}
			elsif($c0 =~ /X/i) {
				push(@klist, {name => 'X',     kx => 0.5, ky => 0.0, kz => 0.0});
			}
			elsif($c0 =~ /Y/i) {
				push(@klist, {name => 'Y',     kx => 0.0, ky => 0.5, kz => 0.0});
			}
			elsif($c0 =~ /Z/i) {
				push(@klist, {name => 'Z',     kx => 0.0, ky => 0.0, kz => 0.5});
			}
			elsif($c0 =~ /R/i) {
				push(@klist, {name => 'R',     kx => 0.5, ky => 0.5, kz => 0.5});
			}
			elsif($c0 =~ /M/i) {
				push(@klist, {name => 'M',     kx => 0.5, ky => 0.5, kz => 0.0});
			}
		}
	}
	if(@klist == 0) {
		print "Can not read $klistfile.\n";
		print "Use default klist.\n";
		push(@klist, {name => 'R',     kx => 0.5, ky => 0.5, kz => 0.5});
		push(@klist, {name => 'GAMMA', kx => 0.0, ky => 0.0, kz => 0.0});
		push(@klist, {name => 'X',     kx => 0.5, ky => 0.0, kz => 0.0});
		push(@klist, {name => 'M',     kx => 0.5, ky => 0.5, kz => 0.0});
		push(@klist, {name => 'Z',     kx => 0.0, ky => 0.0, kz => 0.5});
		push(@klist, {name => 'GAMMA', kx => 0.0, ky => 0.0, kz => 0.0});
	}

	return \@klist;
}

sub SavePWSCFInpFile
{
	my ($this, $Function, $PseudoPotential, $Crystal, $filename,
			$IsChooseRandomly, $PWSCFUseBravaisLattice, $PWSCFUseCELLDM,
			$PPType, $Functional, $PPCopyMode, $SavePPDir, $IsPrint, %args) = @_;
	$PPType                  = 'PAW'      if($PPType eq '');
	$Functional              = 'PBE'      if($Functional eq '');
	$PPCopyMode              = 'Selected' if($PPCopyMode eq '');
	if($PPCopyMode eq 'No') {
		$SavePPDir = $PWSCFPPDir;
	}
	elsif($SavePPDir eq '') {
		$SavePPDir              = './'       
	}
	
	$args{BurstToP1}        =   0    if(!defined $args{BurstToP1});
	$PWSCFUseBravaisLattice =   0    if(!defined $PWSCFUseBravaisLattice);
	$PWSCFUseCELLDM         =   1    if(!defined $PWSCFUseCELLDM);
	$args{PrintForFindPP}   =   1    if(!defined $args{PrintForFindPP});
	$args{aKProduct}        =   1.5  if(!defined $args{aKProduct});
	$args{ecutwfc}          =  30.0  if(!defined $args{ecutwfc});
	$args{ecutrho}          = 140.0  if(!defined $args{ecutrho});
	$args{SpinPolarized}    = 1      if(!defined $args{SpinPolarized});
	$args{KPoints}          = 50     if(!defined $args{KPoints});
	$args{press}            = 0.001  if(!defined $args{press});
	$args{etot_conv_thr}    = 0.5e-4 if(!defined $args{etot_conv_thr});
	$args{forc_conv_thr}    = 0.5e-3 if(!defined $args{forc_conv_thr});
	$args{conv_thr}         = 0.5e-6 if(!defined $args{conv_thr});
	$args{press_conv_thr}   = 1.0    if(!defined $args{press_conv_thr});
	$IsPrint                = 0      if(!defined $IsPrint);
	my $LF = "<br>\n";
	if($args{BurstToP1}) {
		$Crystal->BurstToP1();
	}

	my $kListFile;
	if($Function =~ /bands/i) {
		print "\n";

		my @klists = glob("*.KPOINTS");
		$kListFile = (@klists > 0)? $klists[0] : '';
		if($kListFile eq '') {
			@klists = glob("*.klist");
			$kListFile = (@klists > 0)? $klists[0] : '';
		}
		if($kListFile) {
			print("\nkList file: $kListFile\n");
		}
		else {
			print("\nLlist file is not found. Use default.\n");
		}
	}

	if($PPCopyMode ne 'No' and !-e $SavePPDir) {
		mkdir($SavePPDir);
	}

	my $vasp = new VASP;
	
	my @AtomTypeList         = $Crystal->GetCAtomTypeList();
	my $nAtomType            = @AtomTypeList;
	my @AtomSiteList         = $Crystal->GetCAsymmetricAtomSiteList();
	my $nAtomSite            = @AtomSiteList;
	my @ExpandedAtomSiteList = $Crystal->GetCExpandedAtomSiteList();
	my $nTotalExpandedAtomSite = @ExpandedAtomSiteList;
#print "nTEA=$nTotalExpandedAtomSite/$nAtomType/$nAtomSite\n";
#exit;

	my $pPPFiles = $this->FindPPFiles(\@AtomTypeList, $Functional, $PPType, $args{PrintForFindPP});
	my $nEVal = $this->GetnEVal($Crystal, $pPPFiles);
	if($args{nbnd} eq '') {
		$args{nbnd} = int($nEVal / 2 - 1.0e-6) + 1;
	}
	else {
		my @a = (Utils::Split("\s*,\s*", $args{nbnd}));
		for(my $i = 0 ; $i < @a ; $i++) {
			if($args{nbnd} =~ /^x(.*)$/) {
				my $v = $1 * (int($nEVal / 2 - 1.0e-6) + 1);
				$args{nbnd} = $v if($args{nbnd} < $v);
			}
			else {
				my $v = $1;
				$args{nbnd} = $v if($args{nbnd} < $v);
			}
		}
	}
	$args{nbnd} = (int($nEVal / 2 - 1.0e-6) + 1) if($args{nbnd} < (int($nEVal / 2 - 1.0e-6) + 1));
	$args{nbnd} = int($args{nbnd} - 1e-6) + 1;
print "# of valence electrons: $nEVal from PP files. nbnd = $args{nbnd} will be used.\n";

	my ($MaxEcut, $MaxRhoEcut) = $this->GetMaxEcut($Crystal, $pPPFiles);
	if($args{ecutwfc} eq '') {
		$args{ecutwfc} = $MaxEcut;
	}
	else {
		my @a = (Utils::Split("\s*,\s*", $args{ecutwfc}));
		for(my $i = 0 ; $i < @a ; $i++) {
			if($a[$i] =~ /^x(.*)$/) {
				my $v = $1 * $MaxEcut;
				$args{ecutwfc} = $v if($args{ecutwfc} < $v);
#print "$v = $1 * $MaxEcut, $a[$i]\n";
			}
			else {
				$args{ecutwfc} = $a[$i] if($args{ecutwfc} < $a[$i]);
			}
		}
	}
	$args{ecutwfc} = $MaxRhoEcut if($args{ecutwfc} < $MaxRhoEcut);

print "Ecut for basis functions: $MaxEcut Ry from PP files, $args{ecutwfc} Ry will be used.\n";
	if($args{ecutrho} eq '') {
		$args{ecutrho} = $MaxRhoEcut;
	}
	else {
		my @a = (Utils::Split("\s*,\s*", $args{ecutrho}));
		for(my $i = 0 ; $i < @a ; $i++) {
			if($a[$i] =~ /^x(.*)$/) {
				my $v = $1 * $MaxRhoEcut;
				$args{ecutrho} = $v if($args{ecutrho} < $v);
#print "$v = $1 * $MaxRhoEcut, $a[$i]\n";
			}
			else {
				$args{ecutrho} = $a[$i] if($args{ecutrho} < $a[$i]);
			}
		}
	}
	$args{ecutrho} = $MaxRhoEcut if($args{ecutrho} < $MaxRhoEcut);
print "Ecut for rho            : $MaxRhoEcut Ry from PP files, $args{ecutrho} Ry will be used.\n";
#exit;
	print "\n";

	my @Files;
	my @filenames = fileparse($filename, "\.[^\.]+");
	my $OutputDir = $filenames[1];

#ファイル作製開始
	unless(open(OUT,">$filename")) {
		print "Can not write to $filename.$LF$LF";
		return;
	}
	binmode(OUT);

	my ($dt, $nstep) = (20, 100);
	if($Function =~ /md/i) {
		$dt    = 20;
		$nstep = 100;
	}
	elsif($Function =~ /relax/i) {
		$dt    = 10;
		$nstep = 50;
	}
	$nstep = $args{nstep} if(defined $args{nstep});

	print OUT " &CONTROL\n";
	my $SampleName = $this->SampleName();
	print OUT "    title         = '$SampleName'\n";
	print OUT "    prefix        = 'wm'\n";
#	print OUT "    prefix        = '$SampleName'\n";
	print OUT "    outdir        = '${SampleName}_qe_data',\n";
#	print OUT "    outdir        = './',\n";
	print OUT "    pseudo_dir    = '$SavePPDir',\n";
	print OUT "    verbosity     = 'high',\n";
# scf nscf relax vc-relax md phonon
	if($Function =~ /^dos$/i) {
		print OUT "    calculation   = 'nscf',\n";
	}
	else {
		print OUT "    calculation   = '$Function',\n";
	}
	if($Function =~ /^scf$/i or $Function =~ /relax/i) {
		print OUT "    restart_mode  = 'from_scratch',\n";
	}
	else {
		print OUT "    restart_mode  = 'restart',\n";
	}
	print OUT "    wf_collect    = .True.,\n";
	print OUT "    tprnfor       = .True.,\n";
	print OUT "    iprint        = 1,\n";
	print OUT "    etot_conv_thr = $args{etot_conv_thr},\n";
	print OUT "    forc_conv_thr = $args{forc_conv_thr},\n";
	print OUT "    dt            = $dt,\n";
	print OUT "    nstep         = $nstep,\n";
	print OUT "    tstress       = .True.,\n";
#	print OUT "    max_seconds = 6000,\n";
	print OUT "    disk_io       = 'high',\n";
	print OUT " /\n";

	print OUT " &SYSTEM\n";
#	my $SPGName = $Crystal->SPGNameByOutputMode();
#	my $iSPG    = $Crystal->iSPGByOutputMode();

	my ($a,$b,$c,$alpha,$beta,$gamma) = $Crystal->LatticeParameters();
#	my ($ax,$ay,$az,$bx,$by,$bz,$cx,$cy,$cz) = $Crystal->LatticeVectors();
	my $cosG = Sci::dcos($gamma);
	my $cosB = Sci::dcos($beta);
	my $cosA = Sci::dcos($alpha);
	my $ka = 20.0;
	my ($nkx, $nky, $nkz) = $vasp->AnalyzeaKProduct($Crystal, $args{aKProduct});
#	my ($nkx, $nky, $nkz);
#	$nkx = int($ka / $a + 1.001);
#	$nky = int($ka / $b + 1.001);
#	$nkz = int($ka / $c + 1.001);

	my $BravaisLattice       = 'P';
	my $iPWSCFBravaisLattice = 14;
	if($PWSCFUseBravaisLattice) {
#		@ExpandedAtomSiteList   = $Crystal->GetCBravaisLatticeAsymmetricAtomSiteList();
		$nTotalExpandedAtomSite = @ExpandedAtomSiteList;
		$BravaisLattice         = $Crystal->GetBravaisLattice();
		$iPWSCFBravaisLattice   = $this->GetiPWSCFBravaisLattice($BravaisLattice);
		$iPWSCFBravaisLattice   = 14 if(!defined $iPWSCFBravaisLattice);

		print "<H3>Use Bravais Lattice [$BravaisLattice: Reduced to ibrav = $iPWSCFBravaisLattice]</H3>\n";
#		print "   Atomic sites are reduced to $nTotalExpandedAtomSite.<br>\n";
		print "\n";
	}
#print "nAS=$nTotalExpandedAtomSite\n"+
#exit;
#	my $nAsymmetricAtomSite = &GetCIFValueByKey("nAsymmetricAtomSite");

	print OUT "    ibrav = $iPWSCFBravaisLattice,\n";
# celldm()で指定するときは、a.u.。 A,B,Cの時はオングストローム
	if($PWSCFUseCELLDM) {
		my $k = $a0 * 1.0e10;
#print "latt=$a,$b,$c\n";
		$a /= $k;
		$b /= $k;
		$c /= $k;
		if($iPWSCFBravaisLattice == 5) {	# Trigonal R
			my $aR = $a;
			if(abs($cosA) < 1.0e-6) {
				print "This is a hexagonal unit cell [a=$a, c=$c bohr].<br>\n";
				$aR = sqrt(3.0 * $a*$a + $c*$c) / 3.0;
				print "  Convert to a rhombohedral lattice [aR=$aR bohr, ";
				printf OUT "    celldm(1)   = %10.6f,\n", $aR;
			}
		}
		else {
			printf OUT "    celldm(1)   = %10.6f,\n", $a;
		}
		if(1 <= $iPWSCFBravaisLattice and $iPWSCFBravaisLattice <= 3) {
		}
		elsif($iPWSCFBravaisLattice == 4 or $iPWSCFBravaisLattice == 6 or $iPWSCFBravaisLattice == 7) {
						# Hexagonal and Trigonal P, Tetragonal P, I
			printf OUT "    celldm(3)   = %10.6f,\n", $c / $a;
		}
		elsif($iPWSCFBravaisLattice == 5) {	# Trigonal R
			my $cosAR = $cosA;
			if(abs($cosA) < 1.0e-6) {
				my $ca = $c / $a;
				my $sina2 = 1.5 / sqrt(3.0 + $ca*$ca);
				my $alpha = asin($sina2) * 2.0;
				my $alphadegree = $alpha * $todeg;
				$cosAR = cos($alpha);
				print "alpha=$alphadegree].<br>\n";
				print "\n";
			}
			printf OUT "    celldm(4)   = %10.6f,\n", $cosAR;
		}
		elsif(8 <= $iPWSCFBravaisLattice and $iPWSCFBravaisLattice <= 11) {	# Orthorhombic P, base/face/body-centered
			printf OUT "    celldm(2)   = %10.6f,\n", $b / $a;
			printf OUT "    celldm(3)   = %10.6f,\n", $c / $a;
		}
		elsif($iPWSCFBravaisLattice == 12 or $iPWSCFBravaisLattice == 13) {	# Monoclinic P, base-centered
			printf OUT "    celldm(2)   = %10.6f,\n", $b / $a;
			printf OUT "    celldm(3)   = %10.6f,\n", $c / $a;
			printf OUT "    celldm(4)   = %10.6f,\n", $cosB;
		}
		else {
			printf OUT "    celldm(2)   = %10.6f,\n", $b / $a;
			printf OUT "    celldm(3)   = %10.6f,\n", $c / $a;
			printf OUT "    celldm(4)   = %10.6f,\n", $cosA;
			printf OUT "    celldm(5)   = %10.6f,\n", $cosB;
			printf OUT "    celldm(6)  = %10.6f,\n", $cosG;
		}
	}
	else {
		printf OUT "      A = %10.6f,\n", $a;
		printf OUT "      B = %10.6f,\n", $b;
		printf OUT "      C = %10.6f,\n", $c;
		printf OUT "      cosAB = %10.6f,\n", $cosG;
		printf OUT "      cosAC = %10.6f,\n", $cosB;
		printf OUT "      cosBC = %10.6f,\n", $cosA;
	}

	print OUT "    nat         = $nTotalExpandedAtomSite,\n";
	print OUT "    nspin       = $args{SpinPolarized},\n";
	print OUT "    ntyp        = $nAtomType,\n";
	print OUT "    ecutwfc     = $args{ecutwfc},\n";
	print OUT "    ecutrho     = $args{ecutrho},\n";
	print OUT "    nosym       = .False.,\n";
	print OUT "    noinv       = .False.,\n";
	print OUT "    tot_charge  = 0,\n";
	print OUT "    nbnd        = $args{nbnd},\n";
#	print OUT "    starting_magnetization(1)=0.5,\n";
	print OUT "    occupations = 'smearing',\n";
	print OUT "    smearing    = 'mp',\n";
#	print OUT "    smearing    = 'methfessel-paxton',\n";
	print OUT "    degauss     = 0.02,\n";
	print OUT " /\n";

	print OUT " &ELECTRONS\n";
	print OUT "    conv_thr         = $args{conv_thr},\n";
	print OUT "    mixing_beta      = 0.7,\n";
	print OUT "    mixing_mode      = 'plain',\n";
	print OUT "    electron_maxstep = 100,\n";
	print OUT "    diagonalization  = 'cg',\n";
	print OUT " /\n";
	print OUT "#ion_dynamics: damp for relax, verlet for md\n";
	print OUT " &IONS\n";
	if($Function =~ /relax/i) {
		print OUT "    ion_dynamics  = 'damp',\n";
	}
	elsif($Function =~ /md/i) {
		print OUT "   ion_dynamics='verlet',\n";
		print OUT "   ion_temperature='rescaling',\n";
		print OUT "   tempw=600.D0,\n";
		print OUT "   refold_pos=.TRUE.,\n";
	}
	else {
		print OUT "    ion_dynamics='damp',\n";
	}
	print OUT " /\n";
	print OUT "#cell_dofree: all, none, bfgs, cg, sd, damp-w\n";
	if($Function =~ /md/i or $Function =~ /vc.*relax/i) {
		print OUT " &CELL\n";
		print OUT "    cell_dofree   = 'all',\n";
#		print OUT "    cell_dynamics = 'none',\n";
		print OUT "    cell_dynamics = 'bfgs',\n";
#		print OUT "    cell_dynamics = 'sd',\n";
#		print OUT "    cell_dynamics = 'damp-w',\n";
		print OUT "    press          = $args{press},\n";
		print OUT "    press_conv_thr = $args{press_conv_thr},\n";
		print OUT "    wmass          = 0.20,\n";
		print OUT "    cell_factor    = 2.0,\n";
		print OUT " /\n";
	}
	if($Function =~ /phonon/i) {
		print OUT " &phonon\n";
		print OUT "    xqq(1) = .00, xqq(2) = .00, xqq(3) = 1.00\n";
		print OUT " /\n";
	}
	print OUT "ATOMIC_SPECIES\n";
	for(my $i = 0 ; $i < $nAtomType ; $i++) {
		my $atom = $AtomTypeList[$i];
		my $name = $atom->AtomNameOnly();
		my ($mass, $charge, $ai, $bi, $ci, $rad) = AtomType::GetBusingParameter($name);

		my $PPURL = $PWSCFPPURL;
		$PPURL =~ s/{Atom}/$name/g;
		push(@Files, $PPURL);

		if($PPCopyMode eq 'All' or $PPCopyMode eq 'Selected') {
			my $SourcePPPath   = Deps::MakePath($PWSCFPPDir, $pPPFiles->{$name});
			my $DestPPPath   = Deps::MakePath($SavePPDir, $pPPFiles->{$name});
			unless(Utils::CopyFile($SourcePPPath, $DestPPPath)) {
				print "Error in PWSCF::SavePWSCFInpFile: "
				    . "Can not copy [$SourcePPPath] to [$DestPPPath]$LF";
			}
			push(@Files, $DestPPPath);
		}
		printf OUT "   %4s %10.5f  %s \n", uc($name), $mass, $pPPFiles->{$name};
	}
	print OUT "\n";
	print OUT "ATOMIC_POSITIONS crystal \n";
#	if($PWSCFUseBravaisLattice) {
#		print OUT "ATOMIC_POSITIONS alat \n";
#	}
#	else {
#		print OUT "ATOMIC_POSITIONS crystal \n";
#	}
#Occupancyが1のサイト
	for(my $i = 0 ; $i < @ExpandedAtomSiteList ; $i++) {
		my $atom      = $ExpandedAtomSiteList[$i];
		my $atomname  = $atom->AtomNameOnly();
		my $charge    = $atom->Charge();
		my ($x,$y,$z) = $atom->Position(1);
		my $occupancy = $atom->Occupancy();
		next if($occupancy < 0.9999);

#		my ($xc,$yc,$zc) = $Crystal->FractionalToCartesian($x,$y,$z);
		printf OUT "   %4s  %14.8f %14.8f %14.8f\n",
			uc($atomname), $x, $y, $z;
#		printf OUT "   %4s  %14.8f %14.8f %14.8f   %2d %2d %2d\n",
#			uc($atomname), $x, $y, $z, 1, 1, 1;
	}
#Occupancyが1未満のサイト
	for(my $i = 0 ; $i < @ExpandedAtomSiteList ; $i++) {
		my $atom      = $ExpandedAtomSiteList[$i];
		my $atomname  = $atom->AtomNameOnly();
		my $charge    = $atom->Charge();
		my ($x,$y,$z) = $atom->Position(1);
		my $occupancy = $atom->Occupancy();
		next if($occupancy >= 0.9999);

#		my ($xc,$yc,$zc) = $Crystal->FractionalToCartesian($x,$y,$z);
		printf OUT "   %4s  %14.8f %14.8f %14.8f\n",
			uc($atomname), $x, $y, $z;
	}
	print OUT "\n";
	if($Function =~ /bands/i) {
		my $pKList = $this->MakeKPOINTSForBand($Crystal, $kListFile, $args{KPoints}, $args{nk});
		print OUT "K_POINTS {tpiba_b} \n";
#		print OUT "K_POINTS {crystal} \n";
		my $nk = @$pKList;
		print OUT "$nk\n";
		for(my $i = 0 ; $i < @$pKList ; $i++) {
			my $p = $pKList->[$i];
#			print OUT "  $i: $p->{name}: ($p->{kx}, $p->{ky}, $p->{kz})\n";
			printf OUT "  %9.6f %9.6f %9.6f   %d\n", $p->{kx}, $p->{ky}, $p->{kz}, $args{nk};
		}
#		print OUT "K_POINTS {crystal_b} \n";
#		print OUT "11\n";
#		print OUT "X  20\n";
#		print OUT "gG 20\n";
#		print OUT "Y   0\n";
#		print OUT "L  20\n";
#		print OUT "Z   0\n";
#		print OUT "N  20\n";
#		print OUT "gG 20\n";
#		print OUT "M   0\n";
#		print OUT "R  20\n";
#		print OUT "gG  0\n";
	}
	else {
		print OUT "K_POINTS {automatic} \n";
		print OUT "  $nkx $nky $nkz   1 1 1 \n";
	}

	close(OUT);

	my $TemplateDir = Deps::MakePath($PWSCFDir, 'Template', 0);
	if($Function =~ /band/i) {
		my $Source = Deps::MakePath($TemplateDir, 'bands.in', 0);
		my $Target = 'bands.in';
		print "Copy [$Source] to [$Target]\n";
		if(Utils::CopyFile($Source, $Target) <= 0) {
			print "\nError in PWSCF::SaveInputFiles: Can not copy [$Source] to [$Target].\n\n\n";
			exit;
		}
	}
	elsif($Function =~ /do/i) {
		my $Source = Deps::MakePath($TemplateDir, 'dos.in', 0);
		my $Target = 'dos.in';
		print "Copy [$Source] to [$Target]\n";
		if(Utils::CopyFile($Source, $Target) <= 0) {
			print "\nError in PWSCF::SaveInputFiles: Can not copy [$Source] to [$Target].\n\n\n";
			exit;
		}
		$Source = Deps::MakePath($TemplateDir, 'projwfc.in', 0);
		$Target = 'projwfc.in';
		print "Copy [$Source] to [$Target]\n";
		if(Utils::CopyFile($Source, $Target) <= 0) {
			print "\nError in PWSCF::SaveInputFiles: Can not copy [$Source] to [$Target].\n\n\n";
			exit;
		}
	}

	return ($filename, @Files);
}

sub GetInpFileName
{
	my ($this, $App, $path) = @_;

	my @filenames  = fileparse($path, "\.[^\.]+");
	my $infile = $filenames[0] . ".winp";
	   $infile = $filenames[0] . ".inp" if(!-e $infile);
	return $infile;
};

sub GetnFrameInOutput
{
	my ($this, $App, $InpFile) = @_;

	my $in  = JFile->new($InpFile, "r") or die "$!: Can not read [$InpFile].\n";

	my $nStep = 0;
	while(1) {
		my $line = $in->ReadLine();
		last if(!defined $line);

		$nStep++ if($line =~ /ATOMIC_POSITIONS/);
	}
	$in->Close();

	return $nStep;
}

sub ReadNextCrystalStructureFromOutput
{
	my ($this, $App, $in, , $iStep, $Crystal, $PrintStructures) = @_;
	$PrintStructures = 1 if(!defined $PrintStructures);

print "\nRead Next crystal structure and properties: Step $iStep\n" if($PrintStructures);

	my $pHash = {};

	my $line;
	my $iter = 0;
	my ($PrevE, $TotE) = (0);
	my @EConvergence;
if($PrintStructures) {
print " Iteration $iter: Total energy     dE Ry\n";
}
	while(1) {
		$line = $in->ReadLine();
		return undef if(!$line);
		
		if($line =~ /!\s*total energy\s*=\s*([+\-0-9\.eEdD]+)/) {
			$TotE = $1;
			$EConvergence[$iter] = $TotE;
			my $dE = $TotE - $PrevE;
			$PrevE = $TotE;
if($PrintStructures) {
printf "         %04d  %16.8e %10.4e\n", $iter, $TotE, $dE;
}
			$iter++;
		}
		if($line =~ /Forces acting on atoms/) {
			last;
		}
	}
	$pHash->{TotalEnergy} = $TotE;
	$pHash->{pEConvergence} = \@EConvergence;

#	$line = $in->SkipTo('Forces acting on atoms');
#print "line [$line]\n";
	return undef if(!defined $line or $line eq '');

#次のデータがあることを確認してからリストをクリアー
	my @AtomTypeList         = $Crystal->GetCAtomTypeList();
	my $nAtomType            = @AtomTypeList;
	my @AtomSiteList         = $Crystal->GetCAsymmetricAtomSiteList();
	my $nAtomSite            = @AtomSiteList;
	my @ExpandedAtomSiteList = $Crystal->GetCExpandedAtomSiteList();
	my $nTotalExpandedAtomSite = @ExpandedAtomSiteList;

	$Crystal->ClearAtomTypeList();
	$Crystal->ClearAtomSiteList();

	$line = $in->SkipTo('atom ');
	my (@fx, @fy, @fz);
	my $nAtomSite = 0;
	while(1) {
#print "$nAtomSite: line [$line]\n";
		($fx[$nAtomSite], $fy[$nAtomSite], $fz[$nAtomSite])
				= ($line =~ /force\s*=\s*([+-]?[\d\.eEdD]+)\s+([+-]?[\d\.eEdD]+)\s+([+-]?[\d\.eEdD]+)/);
		last if(!defined $fz[$nAtomSite]);

		$fx[$nAtomSite] =~ s/d/e/i;
		$fy[$nAtomSite] =~ s/d/e/i;
		$fz[$nAtomSite] =~ s/d/e/i;
#print "$nAtomSite: F=($fx[$nAtomSite], $fy[$nAtomSite], $fz[$nAtomSite])\n";
		$nAtomSite++;

		$line = $in->ReadLine();
		Utils::DelSpace($line);
		last if($line eq '');
	}
if($PrintStructures) {
print "nAtomSite = $nAtomSite\n";
}

	my $pos = $in->tell();
	my $line = $in->SkipTo('CELL_PARAMETERS');
	if(!$line) {
		$in->seek($pos, 0);
	}
	else {
		my ($alat) = ($line =~ /alat\s*=\s*([+\-0-9\.eEdD]+)/);
#print "line: $line";
#print "alat=$alat\n";
		my ($a11, $a12, $a13) = Utils::Split("\\s+", $in->ReadLine());
		my ($a21, $a22, $a23) = Utils::Split("\\s+", $in->ReadLine());
		my ($a31, $a32, $a33) = Utils::Split("\\s+", $in->ReadLine());
		$a11 *= $a0*1.0e10 * $alat;
		$a12 *= $a0*1.0e10 * $alat;
		$a13 *= $a0*1.0e10 * $alat;
		$a21 *= $a0*1.0e10 * $alat;
		$a22 *= $a0*1.0e10 * $alat;
		$a23 *= $a0*1.0e10 * $alat;
		$a31 *= $a0*1.0e10 * $alat;
		$a32 *= $a0*1.0e10 * $alat;
		$a33 *= $a0*1.0e10 * $alat;
if($PrintStructures) {
print "Base lattice parameter: $alat a.u.\n";
print "  (a11, a12, a13) = ($a11, $a12, $a13)\n";
print "  (a21, a22, a23) = ($a21, $a22, $a23)\n";
print "  (a31, a32, a33) = ($a31, $a32, $a33)\n";
}
		$Crystal->SetLatticeVector($a11, $a12, $a13, $a21, $a22, $a23, $a31, $a32, $a33);
		my ($a, $b, $c, $alpha, $beta, $gamma) = $Crystal->LatticeParameters();
if($PrintStructures) {
print "  Lattice parameters: $a, $b, $c, $alpha, $beta, $gamma\n";
}
	}

	my $line = $in->SkipTo('ATOMIC_POSITIONS');
#print "line: $line";
	my $iatom = 0;
	my %iAtom;
	for(my $i = 0 ; $i < $nAtomSite ; $i++) {
		$line = $in->ReadLine();
		last if(!defined $line);

#print "line [$line]\n";
		my ($atomname, $x, $y, $z) 
			= ($line =~ /(\w+)\s+([+-]?[\d\.eEdD]+)\s+([+-]?[\d\.eEdD]+)\s+([+-]?[\d\.eEdD]+)/);
		$atomname = ucfirst(lc($atomname));
		$iAtom{$atomname}++;
		$Crystal->AddAtomSite($atomname . $iAtom{$atomname}+1, $atomname, 
				$x, $y, $z, 1.0, $fx[$i], $fy[$i], $fz[$i]);
if($PrintStructures) {
print "$i: $atomname ($x, $y, $z) F=($fx[$i], $fy[$i], $fz[$i])\n";
}
		$iatom++;
	}

	$Crystal->CalMetrics();
	$Crystal->ExpandCoordinates();

	$pHash->{pCrystal} = $Crystal;

	return $pHash;
};

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

	my $pCrystalArray = $this->ReadCrystalStructuresFromOutput(
		$App, $InpFile, $WoutFile, $PrintStructures);
	my $n = @$pCrystalArray;
	my $p = $pCrystalArray->[$n-1];
#print "n=$n  p=$p\n";
#my @latt = $p->LatticeParameters();
#print "latt: ", join(', ', @latt), "\n";
	return $p;
}

sub ReadCrystalStructuresFromOutput
{
	my ($this, $App, $InpFile, $WoutFile, $PrintStructures) = @_;
	$PrintStructures = 1 if(!defined $PrintStructures);

	$App->print("Read Crystal Structures From Output:\n");
	$App->print("Input File : $InpFile\n");
	$App->print("Output File: $WoutFile\n");
	$App->print("\n");

	my @CrystalArray;
	my $Crystal = new Crystal;
	$App->print("\nRead initial crystal structure from input [$InpFile].\n");
	$this->ReadCrystalStructureFromInput($App, $InpFile, $Crystal, $PrintStructures);

	$App->print("\nRead relaxed/md crystal structures from output [$WoutFile].\n");
	my $nStep     = $this->GetnFrameInOutput($App, $WoutFile);
	$App->print("\nStep in [$WoutFile]: $nStep\n");

	my @AtomTypeList         = $Crystal->GetCAtomTypeList();
	my $nAtomType            = @AtomTypeList;
	my @AtomSiteList         = $Crystal->GetCAsymmetricAtomSiteList();
	my $nAtomSite            = @AtomSiteList;
	$App->printf("nAtomSites: %d  nAtomType: %d\n", $nAtomSite, $nAtomType);

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

		$CrystalArray[$i] = $NewCrystal;
#		%{$CrystalArray[$i]} = %$NewCrystal;
	}
	$in->Close();

	return \@CrystalArray;
}

sub ReadCrystalStructureFromInput
{
	my ($this, $App, $InpFile, $Crystal, $PrintStructures) = @_;
	$PrintStructures = 1 if(!defined $PrintStructures);
#print "PrintStructures=$PrintStructures\n";
#exit;

print "ReadCrystalStructureFromInput: from [$InpFile]\n";

	my $in  = JFile->new($InpFile, "r") or die "$!: Can not read [$InpFile].\n";

	my $line = $in->SkipTo("title\\s*=");
#print "l:$line\n";
	my ($s) = ($line =~ /=(.*)$/);
	Utils::DelSpace($s);
	Utils::DelQuote($s);
	$Crystal->SetCrystalName($s);
print "title: $s\n";

	my $CoordinateUnit;
	my $UseCELLPARAMETERS = 0;
	my ($ibrav, @latt, @celldm);

print "Search CELL_PARAMETERS:\n";
	$in->rewind();
	my $line = $in->SkipTo("CELL_PARAMETERS");
	if($line) {
		$line =~ /CELL_PARAMETERS\s+(\S+)/i;
		my $unit = $1;
		$CoordinateUnit = $unit if(defined $unit);
		$UseCELLPARAMETERS = 1;
		
		my ($a11, $a12, $a13) = Utils::Split("\\s+", $in->ReadLine());
		my ($a21, $a22, $a23) = Utils::Split("\\s+", $in->ReadLine());
		my ($a31, $a32, $a33) = Utils::Split("\\s+", $in->ReadLine());

		my $k = 1.0;
		if($CoordinateUnit eq 'bohr') {
			$k = $a0 * 1.0e10;
		}
		$a11 *= $k;
		$a12 *= $k;
		$a13 *= $k;
		$a21 *= $k;
		$a22 *= $k;
		$a23 *= $k;
		$a31 *= $k;
		$a32 *= $k;
		$a33 *= $k;
		$Crystal->SetLatticeVector($a11, $a12, $a13, $a21, $a22, $a23, $a31, $a32, $a33);
		($latt[0], $latt[1], $latt[2], $latt[3], $latt[4], $latt[5]) = $Crystal->LatticeParameters();
if($PrintStructures) {
print "  Lattice parameters: $latt[0], $latt[1], $latt[2], $latt[3], $latt[4], $latt[5]\n";
print "  Lattice vectors:\n";
print "    (a11, a12, a13) = ($a11, $a12, $a13)\n";
print "    (a21, a22, a23) = ($a21, $a22, $a23)\n";
print "    (a31, a32, a33) = ($a31, $a32, $a33)\n";
}
	}

if($PrintStructures) {
print "Atom Types:\n";
}
	$in->rewind();
	my $line = $in->SkipTo("ATOMIC_SPECIES");
	my $natomtypes = 0;
	while(1) {
		$line = Utils::DelSpace($in->ReadLine());
#print "l:$line";
		last if(!defined $line or $line eq '');
		last if($line =~ /ATOMIC_POSITIONS/);

#		my ($name, $mass, $ppfile) = ($line =~ /(\w+)\s+([\.0-9eEdD]+)/);
		my ($name, $mass, $ppfile) = ($line =~ /(\w+)\s+([\.0-9eEdD]+)\s+(\S.*)$/);
		last if(!defined $ppfile);
		$name = ucfirst(lc($name));
if($PrintStructures) {
print "$natomtypes: $name M=$mass PP=$ppfile\n";
}
		$Crystal->AddAtomType($name, 1, 1.0);
		$natomtypes++;
	}

	$in->rewind();
	my $line = $in->SkipTo("&SYSTEM");
	while(1) {
		$line = $in->ReadLine();
		last if(!defined $line);
		last if($line =~ /\//);

#print "l:$line";
		if($line =~ /^\s*ibrav\s*=\s*(\S+)\s*,/) {
			$ibrav = $1;
if($PrintStructures) {
print "ibrav=$ibrav\n";
}
		}
		elsif($line =~ /^\s*celldm\((\d+)\)\s*=\s*(\S+)\s*,/) {
			my $idx        = $1;
			my $idx0       = $idx-1;
			$celldm[$idx0] = $2;
if($PrintStructures) {
print "celldm($idx)=$celldm[$idx0]\n";
}
		}
		elsif($line =~ /^\s*A\s*=\s*(\S+)\s*,/) {
			$latt[0] = $1;
#print "A=$1, $latt[0]\n";
		}
		elsif($line =~ /^\s*B\s*=\s*(\S+)\s*,/) {
			$latt[1] = $1;
#print "B=$1, $latt[1]\n";
		}
		elsif($line =~ /^\s*C\s*=\s*(\S+)\s*,/) {
			$latt[2] = $1;
#print "C=$1, $latt[2]\n";
		}
		elsif($line =~ /^\s*cosBC\s*=\s*(\S+)\s*,/) {
			$latt[3] = Utils::Round(Sci::acos($1) * $todeg, 4, 0);
#print "cosBC=$1, $latt[3]\n";
		}
		elsif($line =~ /^\s*cosAC\s*=\s*(\S+)\s*,/) {
			$latt[4] = Utils::Round(Sci::acos($1) * $todeg, 4, 0);
#print "cosAC=$1, $latt[4]\n";
		}
		elsif($line =~ /^\s*cosAB\s*=\s*(\S+)\s*,/) {
			$latt[5] = Utils::Round(Sci::acos($1) * $todeg, 4, 0);
#print "cosAB=$1, $latt[5]\n";
		}
	}
#print      "    latt: ", join(', ', @latt), "\n";
#print      "    celldm: ", join(', ', @celldm), "\n";
#exit;

# celldm()で指定するときは、a.u.。 A,B,Cの時はオングストローム
	if(defined $celldm[0]) {
#printf OUT "    celldm(2) = %10.6f,\n", $b / $a;
#printf OUT "    celldm(3) = %10.6f,\n", $c / $a;
#printf OUT "    celldm(4) = %10.6f,\n", $cosA;
#printf OUT "    celldm(5) = %10.6f,\n", $cosB;
#printf OUT "    celldm(6) = %10.6f,\n", $cosG;
		my $k = $a0 * 1.0e10;
		$latt[0] = $k * $celldm[0];
#print "k=$k  $a0  $latt[0], $celldm[0]\n";
		if(defined $celldm[1]) {
			$latt[1] = $k* $celldm[1] * $celldm[0];
		}
		else {
			$latt[1] = $latt[0];
		}
		if(defined $celldm[2]) {
			$latt[2] = $k* $celldm[2] * $celldm[0];
		}
		else {
			$latt[2] = $latt[0];
		}
		if(defined $celldm[3]) {
			$latt[3] = Utils::Round(Sci::acos($celldm[3]) * $todeg, 4, 0);
		}
		else {
			$latt[3] = 90.0;
		}
		if(defined $celldm[4]) {
			$latt[4] = Utils::Round(Sci::acos($celldm[3]) * $todeg, 4, 0);
		}
		else {
			$latt[4] = 90.0;
		}
		if(defined $celldm[5]) {
			$latt[5] = Utils::Round(Sci::acos($celldm[3]) * $todeg, 4, 0);
		}
		else {
			$latt[5] = 90.0;
		}
		if($ibrav == 4) { #Hexagonal
			$latt[5] = 120.0;
		}
	}

	if(!$UseCELLPARAMETERS) {
		print "Lattice parameters: $latt[0], $latt[1], $latt[2], $latt[3], $latt[4], $latt[5]\n";
		print "\n";
		$Crystal->SetLatticeParameters($latt[0], $latt[1], $latt[2], $latt[3], $latt[4], $latt[5]);
#print "latt: $latt[0], $latt[1], $latt[2], $latt[3], $latt[4], $latt[5]\n";
	}

if($PrintStructures) {
print "Asymmetric Atom Sites:\n";
}
	$in->rewind();
	my $line = $in->SkipTo("ATOMIC_POSITIONS");
	my $natoms = 0;
	my %iSite;
	while(1) {
		$line =~ /ATOMIC_POSITIONS\s+(\S+)/i;
		my $unit = lc($1);
		
		my $k = 1.0;
		if($unit eq 'crystal') {
		}
		elsif($unit eq 'bohr') {
			$k = $a0 * 1.0e10;
		}

		$line = Utils::DelSpace($in->ReadLine());
#print "l:$line";
		last if(!defined $line or $line eq '');
		last if($line =~ /K_POINTS/);

		my ($name, $x, $y, $z) = ($line =~ /(\w+)\s+([+-\.0-9eEdD]+)\s+([+-\.0-9eEdD]+)\s+([+-\.0-9eEdD]+)/);
		last if(!defined $z);
		$name = ucfirst(lc($name));
		$x = Utils::Reduce01($x);
		$y = Utils::Reduce01($y);
		$z = Utils::Reduce01($z);
		if($unit eq 'bohr') {
			$x = $x / ($latt[0] / $k);
			$y = $z / ($latt[0] / $k);
			$z = $y / ($latt[0] / $k);
		}
if($PrintStructures) {
print "$natoms: $name ($x, $y, $z)\n";
}

		$iSite{$name}++;
		$Crystal->AddAtomSite($name . $iSite{$name}+1, $name, $x, $y, $z, 1.0);
		if($ibrav == 2 or $ibrav == 10) { #FCC, FCO
			$Crystal->AddAtomSite($name . ++$iSite{$name}, $name, $x+0.5, $y+0.5, $z, 1.0);
			$Crystal->AddAtomSite($name . ++$iSite{$name}, $name, $x+0.5, $y, $z+0.5, 1.0);
			$Crystal->AddAtomSite($name . ++$iSite{$name}, $name, $x, $y+0.5, $z+0.5, 1.0);
		}
		elsif($ibrav == 3 or $ibrav == 7 or $ibrav == 11) { #BCC, BCT, BCO
			$Crystal->AddAtomSite($name . ++$iSite{$name}, $name, $x+0.5, $y+0.5, $z+0.5, 1.0);
		}
		elsif($ibrav == -9) { #Orth A-centerd
			$Crystal->AddAtomSite($name . ++$iSite{$name}, $name, $x, $y+0.5, $z+0.5, 1.0);
		}
		elsif($ibrav == 9 or $ibrav == 13) { #Orth C-centerd, Mono C-centered
			$Crystal->AddAtomSite($name . ++$iSite{$name}, $name, $x+0.5, $y+0.5, $z, 1.0);
		}
		elsif($ibrav == -13) { #Mono unique axis b
			$Crystal->AddAtomSite($name . ++$iSite{$name}, $name, $x+0.5, $y, $z+0.5, 1.0);
		}
		$natoms++;
	}


	$in->Close();

	$Crystal->CalMetrics();
	$Crystal->ExpandCoordinates();

	return $Crystal;
};

sub ReadOutputToHash
{
	my ($this, $OutFile, $IsPrint) = @_;
	$IsPrint = 0 if(!defined $IsPrint);
#$IsPrint=1;

	my $in = JFile->new($OutFile, "r");
	if(!$in) {
		print "Error in PWSCF::ReadOutputToHash: Can not read [$OutFile].\n";
		return 0;
	}
	
	my %hash;
	while(1) {
		my $line = $in->ReadLine();
#print "l:$line";
		last if(!defined $line);

		if($line =~ /^\s*Title:/) {
			$line = $in->ReadLine();
			$hash{Title} = Utils::DelSpace($line);
		}
		elsif($line =~ /^\s*celldm\(/) {
			$line = $in->ReadLine();
			($hash{celldm1}) = ($line =~ /celldm\(1\)\s*=\s*(\S+)/);
			($hash{celldm2}) = ($line =~ /celldm\(2\)\s*=\s*(\S+)/);
			($hash{celldm3}) = ($line =~ /celldm\(3\)\s*=\s*(\S+)/);
			$line = $in->ReadLine();
			($hash{celldm4}) = ($line =~ /celldm\(4\)\s*=\s*(\S+)/);
			($hash{celldm5}) = ($line =~ /celldm\(5\)\s*=\s*(\S+)/);
			($hash{celldm6}) = ($line =~ /celldm\(6\)\s*=\s*(\S+)/);
		}
		elsif($line =~ /^\s*crystal axes:/) {
			$line = $in->ReadLine();
			($hash{a11}, $hash{a12}, $hash{a13}) = ($line =~ /\(\s*(\S+)\s*(\S+)\s*(\S+)/);
			($hash{a21}, $hash{a22}, $hash{a23}) = ($line =~ /\(\s*(\S+)\s*(\S+)\s*(\S+)/);
			($hash{a31}, $hash{a32}, $hash{a33}) = ($line =~ /\(\s*(\S+)\s*(\S+)\s*(\S+)/);
		}
		elsif($line =~ /^\s*reciprocal axes:/) {
			$line = $in->ReadLine();
			($hash{Ra11}, $hash{Ra12}, $hash{Ra13}) = ($line =~ /\(\s*(\S+)\s*(\S+)\s*(\S+)/);
			($hash{Ra21}, $hash{Ra22}, $hash{Ra23}) = ($line =~ /\(\s*(\S+)\s*(\S+)\s*(\S+)/);
			($hash{Ra31}, $hash{Ra32}, $hash{Ra33}) = ($line =~ /\(\s*(\S+)\s*(\S+)\s*(\S+)/);
		}
		elsif($line =~ /^\s*the Fermi energy is/) {
			($hash{EF}) = ($line =~ /is\s+(\S+)/);
			$hash{EFRy} = $hash{EF} / $RyToeV;
print "EF   = $hash{EF} eV\n" if($IsPrint);
		}
		elsif($line =~ /^!\s*total\s+energy\s*=/) {
			($hash{EtotRy}) = ($line =~ /=\s*(\S+)/);
			$hash{Etot} = $hash{EtotRy} * $RyToeV;
print "Etot = $hash{EtotRy} Ry  $hash{Etot} eV\n" if($IsPrint);
		}
		elsif($line =~ /^\s*total\s+stress\s*\(Ry/) {
			($hash{Pkbar}) = ($line =~ /P\s*=\s*(\S+)/);
			$hash{PGPa} = $hash{Pkbar} * 0.1;
			my ($a1, $a2, $a3);
			($a1, $a2, $a3, $hash{P11kbar}, $hash{P12kbar}, $hash{P13kbar}) 
				= Utils::Split("\\s+", $in->ReadLine());
			($a1, $a2, $a3, $hash{P21kbar}, $hash{P22kbar}, $hash{P23kbar}) 
				= Utils::Split("\\s+", $in->ReadLine());
			($a1, $a2, $a3, $hash{P31kbar}, $hash{P32kbar}, $hash{P33kbar}) 
				= Utils::Split("\\s+", $in->ReadLine());
if($IsPrint) {
print "P = $hash{Pkbar} kbar  $hash{PGPa} GPa\n";
printf "    (%12.4g %12.4g %12.4g)\n", $hash{P11kbar}, $hash{P12kbar}, $hash{P13kbar};
printf "    (%12.4g %12.4g %12.4g)\n", $hash{P21kbar}, $hash{P22kbar}, $hash{P23kbar};
printf "    (%12.4g %12.4g %12.4g)\n", $hash{P31kbar}, $hash{P32kbar}, $hash{P33kbar};
}
		}
 		elsif($line =~ /The maximum number of steps has been reached/) {
			$hash{RelaxationStatus} = "Error: " . Utils::DelSpace($line);
print "Error: Relaxation not converged\n" if($IsPrint);
		}
		elsif($line =~ /convergence has been achieved/) {
			$hash{SCFStatus} = "Completed: converged in "
					  ."$hash{nRelaxationConvergenceStep} steps\n";
			($hash{nRelaxationConvergenceStep}) = ($line =~ /in\s+(\d+)/);
print "Completed: SCF converged in $hash{nRelaxationConvergenceStep} step\n" if($IsPrint);
		}
	}
	$in->Close();

	return \%hash;
}

sub ReadEtotIneV
{
	my ($this, $OutFile) = @_;

	my $in = JFile->new($OutFile, "r") or die "$!: Can not read [$OutFile].\n\n";
	my $Etot;
	while(1) {
		my $line = $in->SkipTo("!\\s*total\\s+energy\\s*=");
		last if(!defined $line);
		($Etot) = ($line =~ /=\s*([+\-\.\d]+)/);
	}
	$in->Close();
#print "Etot=$Etot\n";
	$Etot *= $RyToeV;
	return $Etot;
}

sub UpdateInpFile
{
	my ($this, $App, $InpFile, $NewInpFile, $Crystal) = @_;

	$App->print("Update PWSCF Input File:\n");
	$App->print("Input File : $InpFile\n");
	$App->print("Output File: $NewInpFile\n");
	$App->print("\n");

	my ($a,$b,$c,$alpha,$beta,$gamma) = $Crystal->LatticeParameters();
	my @AtomTypeList           = $Crystal->GetCAtomTypeList();
	my $nAtomType              = @AtomTypeList;
	my @AtomSiteList           = $Crystal->GetCAsymmetricAtomSiteList();
	my $nAsymmetricAtomSite    = @AtomSiteList;
	my @ExpandedAtomSiteList   = $Crystal->GetCExpandedAtomSiteList();
	my $nTotalExpandedAtomSite = $Crystal->nTotalExpandedAtomSite();

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

#ファイル作製開始
	$out->print(" &CONTROL\n");

	while(1) {
		my $line = $in->ReadLine();
		last if(!defined $line);

		if($line =~ /^\s*A\s*=/) {
			$out->printf("    A =  %14.10f,\n", $a);
		}
		elsif($line =~ /^\s*B\s*=/) {
			$out->printf("    B =  %14.10f,\n", $c);
		}
		elsif($line =~ /^\s*C\s*=/) {
			$out->printf("    C =  %14.10f,\n", $c);
		}
		elsif($line =~ /^\s*cosAB\s*=/) {
			$out->printf("    cosAB = %14.10f,\n", cos($torad*$gamma));
		}
		elsif($line =~ /^\s*cosAC\s*=/) {
			$out->printf("    cosAC = %14.10f,\n", cos($torad*$beta));
		}
		elsif($line =~ /^\s*cosBC\s*=/) {
			$out->printf("    cosBC = %14.10f,\n", cos($torad*$alpha));
		}
		elsif($line =~ /^\s*ATOMIC_POSITIONS/) {
			$out->printf("%s", $line);
			for(my $i = 0 ; $i < $nAsymmetricAtomSite ; $i++) {
				$line = $in->ReadLine();

				my $atom        = $AtomSiteList[$i];
				my ($x, $y, $z) = $atom->Position();
				my $atomtype    = $atom->AtomName();
				$out->printf(" %6s  %14.10f  %14.10f  %14.10f\n", $atomtype, $x, $y, $z);
			}
		}
		else {
			$out->printf("%s", $line);
		}
	}
	$out->Close();
	$in->Close();

	return 1;
};


1;