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

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

#use lib "d:/Programs/Perl/lib";

use strict;
use File::Basename;
use Cwd;

use ProgVars;
use Utils;
use GraphData;
use Sci::Science;
use Sci::EnergyBandArray;
#use MyTk::GraphFrameArray;

use Crystal::Crystal;
use Crystal::AtomType;
use Crystal::AtomSite;
use Crystal::XCrySDen;
use Crystal::CIF;
use Crystal::BandStructure;
use Sci qw($e $me $pi $kB $h $hbar $KToeV $JToeV erfc);
use Sci::ChemicalReaction;

#===============================================
# スクリプト大域変数
#===============================================
my $LF = "<br>\n";
my $DirectorySeparator = Deps::DirectorySeparator(); #"\\";

# 擬ポテンシャル、PAWポテンシャルの優先順位
my %PPPriority;
$PPPriority{B} = [qw(B B_h B_s)];
$PPPriority{C} = [qw(C C_h C_s)];
$PPPriority{N} = [qw(N N_h N_s)];
$PPPriority{O} = [qw(O O_h O_s)];
$PPPriority{F} = [qw(F F_h F_s)];
#$PPPriority{Ne} = [qw(Ne)];

$PPPriority{H}  = [qw(H  H_h)];
$PPPriority{Li} = [qw(Li_sv Li)];
$PPPriority{Be} = [qw(Be Be_sv)];
$PPPriority{Na} = [qw(Na_pv Na Na_sv)];
$PPPriority{K}  = [qw(K_sv K_pv)];
$PPPriority{Ca} = [qw(Ca_pv Ca_sv)];
$PPPriority{Rb} = [qw(Rb_sv Rb_pv)];
#$PPPriority{Sr} = [qw(Sr_sv)];
#$PPPriority{Cs} = [qw(Cs_sv)];
#$PPPriority{Ba} = [qw(Ba_sv)];

$PPPriority{Sc} = [qw(Sc_sv Sc)];
$PPPriority{Ti} = [qw(Ti_pv Ti Ti_sv)];
$PPPriority{V}  = [qw(V_pv V V_sv)];
$PPPriority{Cr} = [qw(Cr_pv Cr Cr_sv)];
$PPPriority{Mn} = [qw(Mn_pv Mn Mn_sv)];
#$PPPriority{Y}  = [qw(Y_sv)];
$PPPriority{Zr} = [qw(Zr_sv Zr)];
#$PPPriority{Nb} = [qw(Nb_pv)];
$PPPriority{Mo} = [qw(Mo_pv Mo)];
$PPPriority{Tc} = [qw(Tc_pv Tc)];
$PPPriority{Hf} = [qw(Hf_pv Hf)];
$PPPriority{Ta} = [qw(Ta_pv Ta)];
$PPPriority{W}  = [qw(W_pv W)];
$PPPriority{Re} = [qw(Re Re_pv)];
$PPPriority{Fe} = [qw(Fe Fe_pv Fe_sv)];
$PPPriority{Co} = [qw(Co Co_pv)];
$PPPriority{Ni} = [qw(Ni Ni_pv)];
$PPPriority{Cu} = [qw(Cu Cu_pv)];
#$PPPriority{Zn} = [qw(Zn)];
$PPPriority{Ru} = [qw(Ru Ru_pv)];
$PPPriority{Rh} = [qw(Rh Rh_pv)];
$PPPriority{Pd} = [qw(Pd Pd_pv)];
#$PPPriority{Ag} = [qw(Ag)];
#$PPPriority{Cd} = [qw(Cd)];
$PPPriority{Os} = [qw(Os_pv Os)];
#$PPPriority{Ir} = [qw(Ir)];
$PPPriority{Pt} = [qw(Pt Pt_pv)];
#$PPPriority{Au} = [qw(Au)];
#$PPPriority{Hg} = [qw(Hg)];

#$PPPriority{Al} = [qw(Al)];
#$PPPriority{Si} = [qw(Si)];
$PPPriority{P}  = [qw(P P_h)];
$PPPriority{S}  = [qw(S S_h)];
$PPPriority{Cl} = [qw(Cl Cl_h)];
#$PPPriority{Ar} = [qw(Ar)];
$PPPriority{Ga} = [qw(Ga_d Ga_h Ga)];
$PPPriority{Ge} = [qw(Ge_d Ge_h Ge)];
#$PPPriority{As} = [qw(As)];
#$PPPriority{Se} = [qw(Se)];
#$PPPriority{Br} = [qw(Br)];
#$PPPriority{Kr} = [qw(Kr)];
$PPPriority{In} = [qw(In_d In)];
$PPPriority{Sn} = [qw(Sn_d Sn)];
#$PPPriority{Sb} = [qw(Sb)];
#$PPPriority{Te} = [qw(Te)];
#$PPPriority{I}  = [qw(I)];
#$PPPriority{Xe} = [qw(Xe)];
$PPPriority{Tl} = [qw(Tl_d Tl)];
$PPPriority{Pb} = [qw(Pb_d Pb)];
$PPPriority{Bi} = [qw(Bi_d Bi)];
$PPPriority{Po} = [qw(Po_d Po)];
$PPPriority{At} = [qw(At_d At)];
#$PPPriority{Rn} = [qw(Rn)];

$PPPriority{La} = [qw(La La_s)];
$PPPriority{Ce} = [qw(Ce Ce_3)];
$PPPriority{Pr} = [qw(Pr Pr_3)];
$PPPriority{Nd} = [qw(Nd Nd_3)];
$PPPriority{Pm} = [qw(Pm Pm_3)];
$PPPriority{Sm} = [qw(Sm Sm_3)];
$PPPriority{Eu} = [qw(Eu Eu_2)];
$PPPriority{Gd} = [qw(Gd Gd_3)];
$PPPriority{Tb} = [qw(Tb_3)];
$PPPriority{Dy} = [qw(Dy_3)];
$PPPriority{Ho} = [qw(Ho_3)];
$PPPriority{Tm} = [qw(Tm Tm_3)];
$PPPriority{Yb} = [qw(Yb Yb_2)];
$PPPriority{Lu} = [qw(Lu Lu_3)];
$PPPriority{Ac} = [qw(Ac Ac_s)];
$PPPriority{Th} = [qw(Th Th_s)];
$PPPriority{Pa} = [qw(Pa Pa_s)];
$PPPriority{U}  = [qw(U  U_s)];
$PPPriority{Np} = [qw(Np Np_s)];
$PPPriority{Pu} = [qw(Pu Pu_s)];
$PPPriority{Pa} = [qw(Pa Pa_s)];

#===============================================
# パス（読み込みDB）
# Web関係変数
# CGI の仮想パス
# プログラム名
#===============================================
my $Program        = ProgVars::Program();
my $ProgramDir     = ProgVars::ProgramDir();
#my $ProgramOldDir  = ProgVars::ProgramOldDir();
my $VASPDir        = ProgVars::VASPDir();
#my $VASPProgOldDir = ProgVars::VASPProgOldDir();
my $VASPDBDir      = ProgVars::VASPDBDir();
my $WebRootDir     = ProgVars::WebRootDir();
my $KListDBDir     = ProgVars::KListDBDir();

#============================================================
#　変数等取得、再定義関数
#============================================================
sub ClearAll      { my $this=shift; delete $this->{FileType}; delete $this->{DataArray}; }
sub FileType      { return shift->{'FileType'}; }
sub FileName      { return shift->{'FileName'}; }
sub SetFileName   { my ($this,$f)=@_; return $this->{'FileName'} = $f; }
#sub DataArray     { return shift->{'DataArray'}; }
#sub SetDataArray  { my ($this, $da) = @_; return $this->{'DataArray'} = $da; }
sub SetSampleName { my ($this,$n)=@_; return $this->{'SampleName'} = $n; }
sub SampleName    { return shift->{'SampleName'}; }
sub KListDBDir {
	my $this = shift;
	return $this->{'KListDBDir'} if($this->{'KListDBDir'});
	return $this->{'KListDBDir'} = $KListDBDir;
}
sub SetKListDBDir { 
	my($this, $d) = @_; 
	$this->{'KListDBDir'} = $d;
	return $d;
}
sub GetFileName {
	my($this, $path, $fname) = @_; 
	my $s;
	if(-d $path) {
		$s = $path;
	}
	else {
		my ($drive, $dir, $filename, $ext, $lastdir, $filebody) = Deps::SplitFilePath($path);
#print "[$drive] [$dir] [$filename]\n";
		if($drive) {
			$s = Deps::MakePath($drive, $dir);
		}
		elsif($dir) {
			$s = $dir;
		}
		else {
			$s = '';
		}
	}
	my $path2 = Deps::MakePath($s, $fname);
#print "s=$path2\n";
	return $path2;
}
sub GetINCARFileName    { my($this,$f)=@_; return $this->GetFileName($f, "INCAR"); }
sub GetOUTCARFileName   { my($this,$f)=@_; return $this->GetFileName($f, "OUTCAR"); }
sub GetKPOINTSFileName  { my($this,$f)=@_; return $this->GetFileName($f, "KPOINTS"); }
sub GetPOSCARFileName   { my($this,$f)=@_; return $this->GetFileName($f, "POSCAR"); }
sub GetCONTCARFileName  { my($this,$f)=@_; return $this->GetFileName($f, "CONTCAR"); }
sub GetPOTCARFileName   { my($this,$f)=@_; return $this->GetFileName($f, "POTCAR"); }
sub GetEIGENVALFileName { my($this,$f)=@_; return $this->GetFileName($f, "EIGENVAL"); }
sub GetDOSCARFileName   { my($this,$f)=@_; return $this->GetFileName($f, "DOSCAR"); }
sub GetPROCARFileName   { my($this,$f)=@_; return $this->GetFileName($f, "PROCAR"); }
sub GetCARDir { 
	my($this, $path)=@_;
	return $path if(-d $path);
	my ($drive, $dir, $filename, $ext, $lastdir, $filebody)
		= Deps::SplitFilePath($path);
	if($drive) {
		return Deps::MakePath($drive, $dir);
	}
	return $dir;
}
sub SpinPolarized { my ($this,$s) = @_; return $this->{IsSpinPolarized} = $s; }
sub IsSpinPolarized {
	my ($this) = @_;
	return $this->{IsSpinPolarized} if(defined $this->{IsSpinPolarized});
	return 0;
}
sub IsNonCollinear {
	my ($this) = @_;
	my $this->ReadINCARtoHash();
	return $this->{IsNonCollinear} if(defined $this->{IsNonCollinear});
	return 0;
}
sub DataArray { my ($this)=@_; return $this->{DataArray}; }
sub SetDataArray
{
	my ($this, $da) = @_;
	if($this->{Application} and $this->{Application}->{MainWindow}) {
		return $this->mw()->{DataArray} = $this->{DataArray} = $da;
	}
	else {
		return $this->{DataArray} = $da;
	}
}

#============================================================
#　継承クラスで定義しなおす関数
#============================================================
sub CheckFileType
{
	my ($path) = @_;
	my ($drive, $dir, $filename, $ext, $lastdir, $filebody) 
		= Deps::SplitFilePath($path);
#print "f: $filename\n";
	if($filename =~ /^OUTCAR$/i) {
		return "VASP Output file";
	}
	if($filename =~ /^DOSCAR$/i) {
		my $in = new JFile($path, "r");
		return undef unless($in);
		my $IsBig = 1;
		for(my $i = 0 ; $i < 6 ; $i++) {
			my $line1 = $in->ReadLine();
			if($in->eof()) {
				$IsBig = 0;
				last;
			}
		}
		$in->Close();
		return "VASP DOSCAR file" if($IsBig);
		return undef;
	}
	if($filename =~ /^EIGENVAL(\..*)?$/i) {
		return "VASP EIGENVAL file";
	}
	if($filename =~ /\.csv$/i) {
#print "f: $filename [$path]\n";
		my $in = new JFile($path, "r");
		return undef if(!$in);
		my $line = $in->ReadLine();
#		my @a = split(/\s*,\s*/, $line);
		my @a = Utils::Split("\\s*,\\s*", $line);
		$in->Close();
		return "CSV Band file" if($a[1] =~ /kx/i);
		return undef;
	}
	return undef;
}

#============================================================
#　コンストラクタ、デストラクタ
#============================================================
sub new
{
	my ($module) = @_;
	my $this = {};
	bless $this;
	return $this;
}

sub DESTROY { my $this = shift; }

#============================================================
#　ファイル読み込み関数
#============================================================
sub ReadOUTCAR
{
	my ($this, $path, $target) = @_;

	my $in = new JFile($path, "rb");
	return undef unless($in);

	$this->SetDataArray(new GraphDataArray);
	my $Data = new GraphData;

	my $line = $in->SkipTo("^\\s+POSCAR: ");
	my ($SampleName) = ($line =~ /^\s+POSCAR:\s*(\w*)\s*$/);
	$SampleName = Deps::GetLastDirectory($path) unless($SampleName);
#print "Sample: $SampleName\n";
	$Data->SetTitle($SampleName);

	$in->SkipTo("Dimension of arrays:");
#   k-Points           NKPTS =      1   number of bands    NBANDS=    383
	$line = $in->ReadLine();
	my ($nBands) = ($line =~ /NBANDS\s*=\s*(\d+)/);
#print "nBands=$nBands\n";
#   number of dos      NEDOS =    301   number of ions     NIONS =     72
	$line = $in->ReadLine();
	my ($nIons) = ($line =~ /NIONS\s*=\s*(\d+)/);
#print "nIons=$nIons\n";
# number of electron  639.0000000 magnetization 
	$line = $in->SkipTo("^\\s*number of electron");
	my ($nElectrons) = ($line =~ /number of electron\s+(\S+)/);
	$nElectrons += 0.0;
#print "nElectrons=$nElectrons\n";

	if($target eq 'EnergyConvergence') {
		my @Iteration;
		my @EnergyWithoutEntropy;
		my $iter = 0;
		$in->rewind();
		while(!$in->eof()) {
			$line = $in->SkipTo("--- Iteration ");
#print "line: $line\n";
			last if($in->eof());
			$line = $in->SkipTo(" energy without entropy =");
			my ($eng) = ($line =~ /entropy =\s*(\S+)/);
			$EnergyWithoutEntropy[$iter] = $eng;
			$Iteration[$iter] = $iter+1;
print "i=$iter: E(total)=$eng\n";
			$iter++;
		}
		$Data->{'x0'}      = \@Iteration;
		$Data->{'x0_Name'} = "Iteration";
		$Data->{'y0'}      = \@EnergyWithoutEntropy;
		$Data->{'y0_Name'} = "Energy / Ry";
	}

	if($target =~ /ForceConvergence/i) {
		my @Iteration;
		my @RMSForce;
		my $iter = 0;
		$in->rewind();
		while(!$in->eof()) {
			$line = $in->SkipTo("POSITION\\s+TOTAL-FORCE");
#print "line: $line\n";
			last if($in->eof());
			$in->ReadLine();
			my $n = 0;
			my $F2 = 0.0;
			while(1) {
				$line = $in->ReadLine();
				last if($line =~ /-----/);
				my ($x,$y,$z,$Fx,$Fy,$Fz) = Utils::Split("\\s+", $line);
				$F2 += $Fx*$Fx + $Fy*$Fy + $Fz*$Fz;
				$n++;

				if($target !~ /RMS/i) {
					my ($n0,$n1) = ($target =~ /\(AtomRange:\s*(\d+)\s*-\s*(\d+)\s*\)/);
					$n0 =  0 unless(defined($n0));
					$n1 = -1 unless(defined($n1));
#print "AtomRange: $n0 - $n1\n";
					next if($n1 > 0 and not ($n0 <= $n and $n <= $n1));
					my $ix = 3*($n-1) + 1;
					my $iy = 3*($n-1) + 2;
					my $iz = 3*($n-1) + 3;
					unless(defined $Data->{"x$ix"}) {
						my (@Fx, @Fy, @Fz);
						$Data->{"x$ix"} = \@Iteration;
						$Data->{"x${ix}_Name"} = "Iteration";
						$Data->{"y${ix}"}      = \@Fx;
						$Data->{"y${ix}_Name"} = "Atom ${n} Force x / eV/Angstrom";
						$Data->{"x$iy"} = \@Iteration;
						$Data->{"x${iy}_Name"} = "Iteration";
						$Data->{"y${iy}"}      = \@Fy;
						$Data->{"y${iy}_Name"} = "Atom ${n} Force y / eV/Angstrom";
						$Data->{"x$iz"} = \@Iteration;
						$Data->{"x${iz}_Name"} = "Iteration";
						$Data->{"y${iz}"}      = \@Fz;
						$Data->{"y${iz}_Name"} = "Atom ${n} Force z / eV/Angstrom";
					}
					$Data->{"y${ix}"}->[$iter] = $Fx;
					$Data->{"y${iy}"}->[$iter] = $Fy;
					$Data->{"y${iz}"}->[$iter] = $Fz;
				}
			}
			my $F = 0.0;
			$F = sqrt($F2 / $n) if($n > 0);
			$RMSForce[$iter] = $F;
			$Iteration[$iter] = $iter+1;
print "i=$iter: F(RMS)=$F\n";
			$iter++;
		}
		$Data->{'x0'}      = \@Iteration;
		$Data->{'x0_Name'} = "Iteration";
		$Data->{'y0'}      = \@RMSForce;
		$Data->{'y0_Name'} = "RMS Force / eV/Angstrom";
	}

	if($target eq 'EnergyLevels') {
		$in->rewind();
		my $iter = 0;
		while(!$in->eof()) {
#print "iter: $iter\n";
			$line = $in->SkipTo("  band No.  band");
			last unless(defined $line);

			my @XConstant;
			my @BandEnergy;
			my @Occupation;
			for(my $i = 0 ; $i < $nBands ; $i++) {
				$line = $in->ReadLine();
#print "line: $line";
				my @a = split(/\s+/, $line);
				@a = Utils::RemoveSpaceElement(@a);
#print "i=$i: $a[0], $a[1], $a[2]\n";
				$XConstant[$i]  = 1.0;
				$BandEnergy[$i] = $a[1];
				$Occupation[$i] = $a[2];
			}
			$Data->{'x0'}                 = \@XConstant;
			$Data->{'x0_Name'}            = "Band No";
			$Data->{"y${iter}"}          = \@BandEnergy;
			$Data->{"y${iter}_Name"}     = "Energy (Iter# $iter)";
			$Data->{"Occupation${iter}"} = \@Occupation;
			
			$iter++;
		}
	}
	$in->Close();

	$Data->CalMinMax();
	$this->{'DataArray'}->AddGraphData($Data);
	$this->{'DataArray'}->{'TargetData'} = $target;

	return $path;
}

sub ReadKeyValue
{
	my ($this, $RegExp, $def, $path) = @_;
	my $in = new JFile($path, "rb");
	return 0.0 if(!defined $in);
	my $line = $in->SkipTo($RegExp);
#	my $val = $1;
	my ($val) = ($line =~ /$RegExp/);
	$in->Close();
	$val = $def unless(defined $val);
	return $val;
}

sub ReadFermiEnergy
{
	my ($this, $path) = @_;
	my $in = new JFile($path, "rb");
	return 0.0 if(!defined $in);
	my $line = $in->SkipTo(" E-fermi :");
#print "p[$path]\n";
#print "l:$line\n";
	my $EF;
	($EF) = ($line =~ /E-fermi\s*:\s*(\S+)/) if($line);
	$in->Close();
	$EF = 0.0 unless(defined $EF);
#print "***EF=$EF eV\n";
#exit;
	return $EF;
}

sub ReadOUTCARValues
{
	my ($this, $path) = @_;
	my $in = new JFile($path, "rb");
	return 0.0 if(!defined $in);
	my $energy;
	my $pos;
	while(1) {
		my $l = $in->SkipTo("\\s*free  energy   TOTEN  =");
		last if(!defined $l);
		$pos = $in->tell();
		($energy) = ($l =~ /free  energy   TOTEN  =\s*(\S+)/);
	}
	$in->seek($pos, 0) if(defined $pos);
#print "energy: $energy\n";
	$this->{FreeEnergy} = $energy;

	my $l1 = $in->SkipTo("\\s*total charge");
#print "l1=$l1\n";
	$in->SkipTo("# of ion     s       p       d       tot");
	$in->ReadLine();
	for(my $i = 1 ; ; $i++) {
		my $l = $in->ReadLine();
		Utils::DelSpace($l);
		last if($l =~ /^-----/ or !defined $l or $l eq '');
		my ($iIon, $s, $p, $d, $tot) = Utils::Split("\\s+", $l);
		$this->{"Ion${iIon}_s"} = $s;
		$this->{"Ion${iIon}_p"} = $p;
		$this->{"Ion${iIon}_d"} = $d;
		$this->{"Ion${iIon}_tot"} = $tot;
		$this->{"nIon"} = $iIon;
#print "iIon=$iIon, $s, $p, $d, $tot\n";
	}

	$in->Close();
	return;
}

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

	my $in = new JFile($path, 'r');
	if(!$in) {
		print "Error in VASP::GetOrbitalsFromPROCAR: $!: Can not read [$path].\n";
		return -1;
	}
	if(!$in->SkipTo("band\\s+\\d")) {
		print "Error in VASP::GetOrbitalsFromPROCAR: Can not find 'band' label in [$path].\n";
		return -2;
	}
	$in->ReadLine();
	my $line = $in->ReadLine();
	my ($ion, @a) = Utils::Split("\\s+", $line);
	my $na = @a;
	if($a[$na-1] eq 'tot') {
		delete $a[$na-1];
	}
	
	$in->Close();
	
	return @a;
}

sub CalSele
{
	my ($this, $Emin, $Estep, $nE, $pD, $T, $EF) = @_;

	my $Sele = 0.0;
	for(my $ie = 0 ; $ie < $nE ; $ie++) {
		my $Ei = $Emin + $ie * $Estep;
		my $fe = Sci::FEe($Ei, $EF, $T);
		my $fh = 1.0 - $fe;
		next if($fe < 1.0e-10 or $fh < 1.0e-10);

		$Sele += $pD->[$ie] * ($fe * log($fe) + $fh * log($fh));
	}
	return -$kB * $Sele / $e;
}

sub CalEele
{
	my ($this, $Emin, $Estep, $nE, $pD, $T, $EF) = @_;

	my $Eele = 0.0;
	for(my $ie = 0 ; $ie < $nE ; $ie++) {
		my $Ei = $Emin + $ie * $Estep;
		$Eele += $Ei * $pD->[$ie] * Sci::FEe($Ei, $EF, $T);
	}
	return $Eele;
}

sub CalNe
{
	my ($this, $Dunder, $Emin, $Estep, $nE, $pD, $T, $EF) = @_;

	my $Ne = $Dunder;
	for(my $ie = 0 ; $ie < $nE ; $ie++) {
		my $Ei = $Emin + $ie * $Estep;
		$Ne += $pD->[$ie] * Sci::FEe($Ei, $EF, $T);
	}
	return $Ne;
}

sub CalEF
{
	my ($this, $Dunder, $Emin, $Estep, $nE, $pD, $T, $NELECT, $EFmin, $EFmax, $EFEPS, $NeEPS) = @_;
	$EFEPS = 1.0e-6 if(!defined $EFEPS);
	$NeEPS = 1.0e-4 if(!defined $NeEPS);

	my $dNe0 = $this->CalNe($Dunder, $pD, $T, $EFmin) - $NELECT;
	my $dNe1 = $this->CalNe($Dunder, $pD, $T, $EFmax) - $NELECT;
#print "  dNe($EFmin) = $dNe0  dNe($EFmax) = $dNe1\n";
		if($dNe0 * $dNe1 > 0.0) {
			print "Error in VASP::CalEF: dNe($EFmin)=$dNe0 and dNe($EFmax)=$dNe1 has the same signs.\n";
			print "  Reconsider dEFmin and dEFmax applicable for binary search.\n";
			print "\n";
			exit;
		}

	my ($EFnew, $dNe0, $dNe1, $dNenew);
	my $EF0 = $EFmin;
	my $EF1 = $EFmax;
	my $nMaxIter = 100;
	for(my $i = 0; $i < $nMaxIter ; $i++) {
		$dNe0   = $this->CalNe($Dunder, $Emin, $Estep, $nE, $pD, $T, $EF0) - $NELECT;
		$dNe1   = $this->CalNe($Dunder, $Emin, $Estep, $nE, $pD, $T, $EF1) - $NELECT;
#print "T=$T  EFrange = ($EF0, $EF1)  dNe = ($dNe0, $dNe1)\n";
		if(abs($dNe0) < $NeEPS) {
			return ($EF0, $dNe0);
		}
		if(abs($dNe1) < $NeEPS) {
			return ($EF1, $dNe1);
		}

#		if($EF1 - $EF0 < $EFEPS and abs($dNe0) < $NeEPS and abs($dNe1) < $NeEPS) {
		if($EF1 - $EF0 < $EFEPS) {
			$EFnew = ($EF0 + $EF1) / 2.0;
			$dNenew = $this->CalNe($Dunder, $Emin, $Estep, $nE, $pD, $T, $EFnew) - $NELECT;
#exit;
			return ($EFnew, $dNenew);
		}

		$EFnew  = ($EF0 + $EF1) / 2.0;
		$dNenew = $this->CalNe($Dunder, $Emin, $Estep, $nE, $pD, $T, $EFnew) - $NELECT;
		if($dNe0*$dNenew < 0.0) {
			$EF1 = $EFnew;
		}
		else {
			$EF0 = $EFnew;
		}
	}
	print "Error in VASP::CalEF: Not converged after $nMaxIter iterations\n";

	return (undef, undef);
}

sub CalculateExactDOSFromEIGENVAL
{
	my ($this, $EIGENVALPath, $EF, $Emin, $Emax, $Estep, $nE) = @_;

	my (@D, @Ne);
	my ($Dunder, $Dover) = (0.0, 0.0);

	my $in = JFile->new($EIGENVALPath, 'r') or die "Can not read [$EIGENVALPath].\n";
	my $line;
	for(my $i = 0 ; $i < 5 ; $i++) {
		$in->ReadLine();
	}
#print "line: $line\n";
	my $SampleName = $line;
	$SampleName =~ s/^\s*//;
	$SampleName =~ s/\s*$//;

	$line = $in->ReadLine();
	my ($nn, $nK, $nLevels) = Utils::Split("\\s+", $line);
print "    nn,nK,nLevels=$nn, $nK, $nLevels\n";

	for(my $ie = 0 ; $ie < $nLevels ; $ie++) {
		$D[$ie] = 0.0;
	}

	for(my $i = 0 ; $i < $nK ; $i++) {
		$in->ReadLine();
		$line = $in->ReadLine();
		my ($kx, $ky, $kz, $w) = Utils::Split("\\s+", $line);
		for(my $ie = 0 ; $ie < $nLevels ; $ie++) {
			$line = $in->ReadLine();
			my ($idx, $Ei, $occ) = Utils::Split("\\s+", $line);
			$Ei -= $EF;
			my $iE = int(($Ei - $Emin) / $Estep);
#print "i=$i ie=$ie   Ei = $Ei ($Emin - $Emax)  iE=$iE (< $nE)\n";
			if($iE < 0.0) {
				$Dunder += 2.0 * $w;
			}
			elsif($nE <= $iE) {
				$Dover += 2.0 * $w;
			}
			else {
				$D[$iE] += 2.0 * $w;
#print "   D[$iE] = $D[$iE]\n";
			}
		}
	}
	$in->Close();
	
	$Ne[0] = $Dunder;
	for(my $ie = 0 ; $ie < $nE ; $ie++) {
		$Ne[$ie+1] += $Ne[$ie] + $D[$ie];
	}
	
	return ($Dunder, $Dover, \@D, \@Ne);
}

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

	my $DOSCARPath  = $this->GetDOSCARFileName($path);
print "VASP::ReadDOSCAR: Read [$DOSCARPath].\n";
	my $PROCARPath  = $this->GetPROCARFileName($path);
print "VASP::ReadDOSCAR: Read [$PROCARPath].\n";
	my $POSCAR      = $this->GetPOSCARFileName($path);
print "VASP::ReadPROCAR: Read [$POSCAR].\n";
	my $Crystal = $this->ReadStructureFromCARFiles($POSCAR, 1);

	my @AtomTypeList = $Crystal->GetCAtomTypeList();
	my $nAtomType = @AtomTypeList;
	my %AtomNum;
	for(my $i = 0 ; $i < $nAtomType ; $i++) {
		my $name = $AtomTypeList[$i]->AtomNameOnly();;
		$AtomNum{$name} = 0 if(!defined $AtomNum{$name});
		my $n = ++$AtomNum{$name};
		$AtomTypeList[$i]->SetAtomName("$name$n");
#print "$i: $name$n\n";
	}
	my @ExpandedAtomSiteList = $Crystal->GetCExpandedAtomSiteListByOutputMode();
#	my @M = $Crystal->GetCnMultiplicityExpandedAtomSiteList();
	my $nTotalExpandedAtomSite = @ExpandedAtomSiteList;
	my @AtomNames;
	for(my $i = 0 ; $i < $nTotalExpandedAtomSite ; $i++) {
		my $atomsite = $ExpandedAtomSiteList[$i];
		$AtomNames[$i] = $atomsite->AtomName(); #AtomNameOnly();
#print "Site[$i]: $AtomNames[$i] ($M[$i])\n";
print "Site[$i]: $AtomNames[$i]\n";
	}

	my $OUTCAR = $this->GetOUTCARFileName($path);
	my $EF = $this->ReadFermiEnergy($OUTCAR);
print "EF: $EF eV\n";

print "VASP::ReadPROCAR: Read [$PROCARPath].\n";
	my @Orbitals = $this->GetOrbitalsFromPROCAR($PROCARPath);
	if(@Orbitals == 0) {
		return -1;
	}

	my $in = new JFile($DOSCARPath, "rb");
	return undef unless($in);

	my $line;
	for(my $i = 0 ; $i < 5 ; $i++) {
		$line = $in->ReadLine();
	}
# BaZn2As2-I4mmm [PAW_PBE] (dos)          
	my $SampleName = $line;
	$SampleName =~ s/^\s*//;
	$SampleName =~ s/\s*$//;
print "Sample [$SampleName]\n";

#     15.00000000    -35.00000000 5000      3.50770154      1.00000000
	$line = $in->ReadLine();
	my ($E1, $E0, $nData, $w);
	($E1, $E0, $nData, $EF, $w) = Utils::Split("\\s+", $line);
print "nData=$nData\n";

	my (@TE, %DOS);
	for(my $i = 0 ; $i < $nData ; $i++) {
		$line = $in->ReadLine();
# 指数が3桁の場合、eが削除されるので追加
		$line =~ s/(\d)\-(\d)/$1e\-$2/g;
		my ($e, $tot, $ne) = Utils::Split("\\s+", $line);
		$TE[$i]              = $e;
		$DOS{Total}[$i]        += $tot;
		$DOS{Interstitial}[$i] += $tot;
		$DOS{NE}[$i]           += $ne;
#print "  $i: $e\t$up\t$dn\n";
	}

	my $iAtom = 0;
	while(1) {
		$line = $in->ReadLine();
# 指数が3桁の場合、eが削除されるので追加
		$line =~ s/(\d)\-(\d)/$1e\-$2/g;
		my ($E1, $E0, $nData, $EF, $w) = Utils::Split("\\s+", $line);
		my $atomsite = $ExpandedAtomSiteList[$iAtom];
		my $AtomName = $atomsite->AtomNameOnly();
print "Site[$iAtom]: $AtomName: nData=$nData\n";

		for(my $i = 0 ; $i < $nData ; $i++) {
			$line = $in->ReadLine();
# 指数が3桁の場合、eが削除されるので追加
			$line =~ s/(\d)\-(\d)/$1e\-$2/g;
			my ($e, @a) = Utils::Split("\\s+", $line);
			last if(!$line or @a <= 1);
if($i == 0) {
	Utils::DelSpace($line);
	my $nOrb = scalar @a;
	print "  nO=$nOrb [$line]\n";
}
			for(my $j = 0 ; $j < @a ; $j++) {
				my $orb = $Orbitals[$j];
#print("$j: $orb  ");
				$DOS{Interstitial}[$i]     -= $a[$j];
				$DOS{$AtomName}{Total}[$i] += $a[$j];
				$DOS{$AtomName}{$orb}[$i]  += $a[$j];
			}
		}
		last if($in->eof());

		$iAtom++;
	}

	$in->Close();

	return (\@AtomTypeList, \@ExpandedAtomSiteList, \@Orbitals, $EF, $nData, \@TE, \%DOS);
}

sub ReadDOSCAR
{
	my ($this, $path, $IsSpinPolarized, $IsNonCollinear) = @_;
	$IsSpinPolarized = $this->IsSpinPolarized() if(!defined $IsSpinPolarized);
	if(!defined $IsNonCollinear) {
		my $INCAR  = $this->GetINCARFileName($path);
		my $pINCAR = $this->ReadINCARtoHash($INCAR);
		$IsNonCollinear  = ($pINCAR->{LSORBIT} =~ /[1T]/i or $pINCAR->{LNONCOLLINEAR} =~ /[1T]/i)? 1 : 0;
		$IsSpinPolarized = 0 if($IsNonCollinear);
	}

	my $POSCAR = $this->GetPOSCARFileName($path);
print "VASP::ReadDOSCAR: Read [$POSCAR].\n";
	my $Crystal = $this->ReadStructureFromCARFiles($POSCAR, 1);
	my $pPOS    = $this->ReadPOSCARtoHash($POSCAR);

	my @AtomTypeList = $Crystal->GetCAtomTypeList();
	my $nAtomType = @AtomTypeList;
	my %AtomNum;
	for(my $i = 0 ; $i < $nAtomType ; $i++) {
		my $name = $AtomTypeList[$i]->AtomNameOnly();;
		$AtomNum{$name} = 0 if(!defined $AtomNum{$name});
		my $n = ++$AtomNum{$name};
		$AtomTypeList[$i]->SetAtomName("$name$n");
#print "$i: $name$n\n";
	}
	my @ExpandedAtomSiteList = $Crystal->GetCExpandedAtomSiteListByOutputMode();
#	my @M = $Crystal->GetCnMultiplicityExpandedAtomSiteList();
	my $nTotalExpandedAtomSite = @ExpandedAtomSiteList;
	my @AtomNames;
	for(my $i = 0 ; $i < $nTotalExpandedAtomSite ; $i++) {
		my $atomsite = $ExpandedAtomSiteList[$i];
		$AtomNames[$i] = $atomsite->AtomName(); #AtomNameOnly();
#print "Site[$i]: $AtomNames[$i] ($M[$i])\n";
print "Site[$i]: $AtomNames[$i]\n";
	}

	my $OUTCAR = $this->GetOUTCARFileName($path);
	my $EF = $this->ReadFermiEnergy($OUTCAR);
print "EF: $EF eV\n";

	my $in = new JFile($path, "rb");
	return undef unless($in);

	$this->SetDataArray(new GraphDataArray);
	my $Data0 = new GraphData;
	my $Data1 = new GraphData;
	my $Data2 = new GraphData;
	my $Data3 = new GraphData;

	my $line;
	for(my $i = 0 ; $i < 5 ; $i++) {
		$line = $in->ReadLine();
	}
# NaCl (dos)                              
	my $SampleName = $line;
	$SampleName =~ s/^\s*//;
	$SampleName =~ s/\s*$//;
	$Data0->SetTitle($SampleName);
	$Data1->SetTitle($SampleName);

#     12.29201091    -13.90257885  301      0.14679613      1.00000000
	$line = $in->ReadLine();
	my ($a,$b,$nData) = Utils::Split("\\s+", $line);
print "nData=$nData\n";

	my @Energy;
	my @DOSup;
	my @DOSdn;
	my @Neup;
	my @Nedn;
	for(my $i = 0 ; $i < $nData ; $i++)  {
		$line = $in->ReadLine();
		last if(not defined $line);
#my $f = 0;
#if($line =~ /\d\-\d/) {
#print "i=$i:$line\n";
#$f = 1;
#exit;
#}
# 指数が3桁の場合、eが削除されるので追加
		$line =~ s/(\d)\-(\d)/$1e\-$2/g;
#if($f) {
#print "    $line\n";
#exit;
#}

		my ($e,$du,$dd, $nu, $nd) = Utils::Split("\\s+", $line);
#print "$i ($nData): $e eV\n";
		if(!defined $nu) {
			$this->SpinPolarized(0);
			($du, $nu) = ($du, $dd);
		}
		else {
			$this->SpinPolarized(1);
		}

		last if(not defined $e);
		$Energy[$i] = $e - $EF;
		$DOSup[$i]    = $du + 0.0;
		$Neup[$i]     = $nu + 0.0;
		if($this->IsSpinPolarized()) {
			$DOSdn[$i]    = $dd;
			$Nedn[$i]     = $nd;
		}
		last if($in->eof());
	}
#print "EF: $EF eV\n";

	my @IntDOSup = @DOSup;
	my @IntDOSdn = @DOSdn;

	my $nAtomKinds = 0;
	my %pPDOSupArray;
	my @pPDOStupArray;
	my @pPDOSsupArray;
	my @pPDOSpupArray;
	my @pPDOSdupArray;
	my @pPDOSfupArray;
	my %pPDOSdnArray;
	my @pPDOStdnArray;
	my @pPDOSsdnArray;
	my @pPDOSpdnArray;
	my @pPDOSddnArray;
	my @pPDOSfdnArray;

	my @pPDOSsmxArray;
	my @pPDOSsmyArray;
	my @pPDOSsmzArray;
	my @pPDOSpmxArray;
	my @pPDOSpmyArray;
	my @pPDOSpmzArray;
	my @pPDOSdmxArray;
	my @pPDOSdmyArray;
	my @pPDOSdmzArray;

	my $pAtomType = $pPOS->{pAtomType};
	my $pAtomSite = $pPOS->{pAtomSite};
#print "pD=$pPOS, $pAtomType, $pAtomSite\n";
my $cc = 0;
	while(1) {
		my $site = $pAtomSite->[$cc];
		last if(!$site);
#print "cc=$cc: site=$site\n";
		my $name = $site->AtomName();
		my $occ  = $site->Occupancy();
#print "$i: $name: $occ\n";

$cc++;
		$in->ReadLine();
		my @PDOStup;
		my @PDOSsup;
		my @PDOSpup;
		my @PDOSdup;
		my @PDOSfup;
		my @PDOStdn;
		my @PDOSsdn;
		my @PDOSpdn;
		my @PDOSddn;
		my @PDOSfdn;

		my @PDOSsmx;
		my @PDOSsmy;
		my @PDOSsmz;
		my @PDOSpmx;
		my @PDOSpmy;
		my @PDOSpmz;
		my @PDOSdmx;
		my @PDOSdmy;
		my @PDOSdmz;

		for(my $i = 0 ; $i < $nData ; $i++) {
			$line = $in->ReadLine();
			last if(not defined $line);

# 指数が3桁の場合、eが削除されるので追加
			$line =~ s/(\d)\-(\d)/$1e\-$2/g;

			my ($e, $dt, $dsu, $dsd, $dpu, $dpd, $ddu, $ddd, $dfu, $dfd);
			my ($dsmx, $dsmy, $dsmz, $dpmx, $dpmy, $dpmz, $ddmx, $ddmy, $ddmz);
			my @aa = Utils::Split("\\s+", $line);
#print "n=", scalar @aa, "\n";
			if($IsNonCollinear) {
				($e, $dsu, $dsmx, $dsmy, $dsmz, $dpu, $dpmx, $dpmy, $dpmz, 
				     $ddu, $ddmx, $ddmy, $ddmz) = @aa;
				$dfu = 0;
				$dt = $dsu + $dpu + $ddu + $dfu;
			}
			elsif($IsSpinPolarized) {
				($e, $dsu, $dsd, $dpu, $dpd, $ddu, $ddd, $dfu, $dfd) = @aa;
				$dt = $dsu + $dpu + $ddu + $dfu;
			}
			else {
				if(@aa == 3) {
					($dsu, $dpu, $ddu) = @aa;
				}
				else {
					($e, $dsu, $dpu, $ddu) = @aa;
				}
				$dt = $dsu + $dpu + $ddu + $dfu;
			}
			last if(not defined $e);

			$dfu = 0.0 if(!defined $dfu);
			$dfd = 0.0 if(!defined $dfd);

			$PDOStup[$i] = $dt;
			$PDOSsup[$i] = $dsu + 0.0;
			$PDOSpup[$i] = $dpu + 0.0;
			$PDOSdup[$i] = $ddu + 0.0;
			$PDOSfup[$i] = $dfu + 0.0;
			$IntDOSup[$i] -= $PDOStup[$i] * $occ;
			if($IsSpinPolarized) {
				$PDOStdn[$i] = $dsd + $dpd + $ddd + $dfd;
				$PDOSsdn[$i] = $dsd + 0.0;
				$PDOSpdn[$i] = $dpd + 0.0;
				$PDOSddn[$i] = $ddd + 0.0;
				$PDOSfdn[$i] = $dfd + 0.0;
				$IntDOSdn[$i] -= $PDOStdn[$i] * $occ;
			}
			last if($in->eof());
		}
		$pPDOStupArray[$nAtomKinds] = \@PDOStup;
		$pPDOSsupArray[$nAtomKinds] = \@PDOSsup;
		$pPDOSpupArray[$nAtomKinds] = \@PDOSpup;
		$pPDOSdupArray[$nAtomKinds] = \@PDOSdup;
		$pPDOSfupArray[$nAtomKinds] = \@PDOSfup;
		if($this->IsSpinPolarized()) {
			$pPDOStdnArray[$nAtomKinds] = \@PDOStdn;
			$pPDOSsdnArray[$nAtomKinds] = \@PDOSsdn;
			$pPDOSpdnArray[$nAtomKinds] = \@PDOSpdn;
			$pPDOSddnArray[$nAtomKinds] = \@PDOSddn;
			$pPDOSfdnArray[$nAtomKinds] = \@PDOSfdn;
		}

my $aname = $AtomNames[$nAtomKinds];
$aname =~ s/\d+//;
#print "atom[$cc]:$aname\n";
		my $st = "$aname total";
		my $ss = "$aname s";
		my $sp = "$aname p";
		my $sd = "$aname d";
		my $sf = "$aname f";
		my $ptu = $pPDOSupArray{$st};
		my $psu = $pPDOSupArray{$ss};
		my $ppu = $pPDOSupArray{$sp};
		my $pdu = $pPDOSupArray{$sd};
		my $pfu = $pPDOSupArray{$sf};
		my $ptd = $pPDOSdnArray{$st};
		my $psd = $pPDOSdnArray{$ss};
		my $ppd = $pPDOSdnArray{$sp};
		my $pdd = $pPDOSdnArray{$sd};
		my $pfd = $pPDOSdnArray{$sf};
		if(!defined $psu) {
#print "cc=$cc: initialize for [$st]\n";
			$ptu = $pPDOSupArray{$st} = [];
			$psu = $pPDOSupArray{$ss} = [];
			$ppu = $pPDOSupArray{$sp} = [];
			$pdu = $pPDOSupArray{$sd} = [];
			$pfu = $pPDOSupArray{$sf} = [];
			$ptd = $pPDOSdnArray{$st} = [];
			$psd = $pPDOSdnArray{$ss} = [];
			$ppd = $pPDOSdnArray{$sp} = [];
			$pdd = $pPDOSdnArray{$sd} = [];
			$pfd = $pPDOSdnArray{$sf} = [];
			for(my $i = 0 ; $i < $nData ; $i++) {
				$ptu->[$i] = 0.0;
				$psu->[$i] = 0.0;
				$ppu->[$i] = 0.0;
				$pdu->[$i] = 0.0;
				$pfu->[$i] = 0.0;
				$ptd->[$i] = 0.0;
				$psd->[$i] = 0.0;
				$ppd->[$i] = 0.0;
				$pdd->[$i] = 0.0;
				$pfd->[$i] = 0.0;
			}
		}
else {
#print "cc=$cc: not initialize for [$st]\n";
}
		for(my $i = 0 ; $i < $nData ; $i++) {
			$ptu->[$i] += $PDOStup[$i];
			$psu->[$i] += $PDOSsup[$i];
			$ppu->[$i] += $PDOSpup[$i];
			$pdu->[$i] += $PDOSdup[$i];
			$pfu->[$i] += $PDOSfup[$i];
			if($this->IsSpinPolarized()) {
				$ptd->[$i] += $PDOStdn[$i];
				$psd->[$i] += $PDOSsdn[$i];
				$ppd->[$i] += $PDOSpdn[$i];
				$pdd->[$i] += $PDOSddn[$i];
				$pfd->[$i] += $PDOSfdn[$i];
			}
		}

		$nAtomKinds++;
		last if($in->eof());
	}
	$in->Close();

	$Data1->{"x0"}      = \@Energy;
	$Data1->{"x0_Name"} = "Energy / eV";
	$Data1->{"y0"}      = \@Neup;
	$Data1->{"y0_Name"} = "Ne";
	if($this->IsSpinPolarized()) {
#		$Data1->{'y0_Name'} = "Ne up";
		$Data3->{'x0'}      = \@Energy;
		$Data3->{'x0_Name'} = "Energy / eV";
		$Data3->{'y0'}      = \@Nedn;
		$Data3->{'y0_Name'} = "Ne down";
	}

	my $cup = 0;
	my $cdn = 0;
	$Data0->{"x$cup"}        = \@Energy;
	$Data0->{"x${cup}_Name"} = "Energy / eV";
	$Data0->{"y$cup"}        = \@DOSup;
	$Data0->{"y${cup}_Name"} = "DOS";
	$cup++;
	if($this->IsSpinPolarized()) {
#		$Data0->{"y${cdn}_Name"} = "DOS up";
		$Data2->{"x$cdn"}        = \@Energy;
		$Data2->{"x${cdn}_Name"} = "Energy / eV";
		$Data2->{"y$cdn"}        = \@DOSdn;
		$Data2->{"y${cdn}_Name"} = "DOS down";
		$cdn++;
	}

	$Data0->{"x$cup"}        = \@Energy;
	$Data0->{"x${cup}_Name"} = "Energy / eV";
	$Data0->{"y$cup"}        = \@IntDOSup;
	$Data0->{"y${cup}_Name"} = "Interstitial";
	$cup++;
	if($this->IsSpinPolarized()) {
#		$Data0->{"y${cup}_Name"} = "Interstitial up";
		$Data2->{"x$cdn"}        = \@Energy;
		$Data2->{"x${cdn}_Name"} = "Energy / eV";
		$Data2->{"y$cdn"}        = \@IntDOSdn;
		$Data2->{"y${cdn}_Name"} = "Interstitial down";
		$cdn++;
	}

my $SiteDOS = 0;
	if($SiteDOS) {
		for(my $i = 0 ; $i < $nAtomKinds ; $i++) {
			my $i1 = $i+1;
			my $pPDOStup = $pPDOStupArray[$i];
			$Data0->{"x$cup"} = \@Energy;
			$Data0->{"y$cup"} = $pPDOStup;
			$Data0->{"x${cup}_Name"} = "Energy / eV";
			$Data0->{"y${cup}_Name"} = "Atom$i total";
			$cup++;
			if($this->IsSpinPolarized()) {
				my $pPDOStdn = $pPDOStdnArray[$i];
#				$Data0->{"y${cdn}_Name"} = "Atom$i total up";
				$Data0->{"x$cdn"} = \@Energy;
				$Data0->{"y$cdn"} = $pPDOStdn;
				$Data0->{"x${cdn}_Name"} = "Energy / eV";
				$Data0->{"y${cdn}_Name"} = "Atom$i total down";
				$cdn++;
			}
		}
		for(my $i = 0 ; $i < $nAtomKinds ; $i++) {
			my $i1 = $i+1;
			my $pPDOSsup = $pPDOSsupArray[$i];
			my $pPDOSpup = $pPDOSpupArray[$i];
			my $pPDOSdup = $pPDOSdupArray[$i];
			my $pPDOSfup = $pPDOSfupArray[$i];
			$Data0->{"x$cup"} = \@Energy;
			$Data0->{"y$cup"} = $pPDOSsup;
			$Data0->{"x${cup}_Name"} = "Energy / eV";
			$Data0->{"y${cup}_Name"} = "Atom$i s";
			$cup++;
			$Data0->{"x$cup"} = \@Energy;
			$Data0->{"y$cup"} = $pPDOSpup;
			$Data0->{"x${cup}_Name"} = "Energy / eV";
			$Data0->{"y${cup}_Name"} = "Atom$i p";
			$cup++;
			$Data0->{"x$cup"} = \@Energy;
			$Data0->{"y$cup"} = $pPDOSdup;
			$Data0->{"x${cup}_Name"} = "Energy / eV";
			$Data0->{"y${cup}_Name"} = "Atom$i d";
			$cup++;
			$Data0->{"x$cup"} = \@Energy;
			$Data0->{"y$cup"} = $pPDOSfup;
			$Data0->{"x${cup}_Name"} = "Energy / eV";
			$Data0->{"y${cup}_Name"} = "Atom$i f";
			$cup++;
			if($this->IsSpinPolarized()) {
				my $pPDOSsdn = $pPDOSsdnArray[$i];
				my $pPDOSpdn = $pPDOSpdnArray[$i];
				my $pPDOSddn = $pPDOSddnArray[$i];
				my $pPDOSfdn = $pPDOSfdnArray[$i];
				$Data2->{"x$cdn"} = \@Energy;
				$Data2->{"y$cdn"} = $pPDOSsdn;
				$Data2->{"x${cdn}_Name"} = "Energy / eV";
				$Data2->{"y${cdn}_Name"} = "Atom$i s down";
				$cdn++;
				$Data2->{"x$cdn"} = \@Energy;
				$Data2->{"y$cdn"} = $pPDOSpdn;
				$Data2->{"x${cdn}_Name"} = "Energy / eV";
				$Data2->{"y${cdn}_Name"} = "Atom$i p down";
				$cdn++;
				$Data2->{"x$cdn"} = \@Energy;
				$Data2->{"y$cdn"} = $pPDOSddn;
				$Data2->{"x${cdn}_Name"} = "Energy / eV";
				$Data2->{"y${cdn}_Name"} = "Atom$i d down";
				$cdn++;
				$Data2->{"x$cdn"} = \@Energy;
				$Data2->{"y$cdn"} = $pPDOSfdn;
				$Data2->{"x${cdn}_Name"} = "Energy / eV";
				$Data2->{"y${cdn}_Name"} = "Atom$i f down";
				$cdn++;
			}
		}
	}
#print "c0=$cup $cdn\n";
my $PDOS = 1;
	if($PDOS) {
		for(my $i = 0 ; $i < $nAtomType ; $i++) {
			my $atomtype = $AtomTypeList[$i];
			my $atomname = $atomtype->AtomName(); #AtomNameOnly();
my $aname = $atomname;
$aname =~ s/\d+//;
			my $st = "$aname total";
print "  PDOS for $st\n";
			my $ptup = $pPDOSupArray{$st};
			$Data0->{"x$cup"} = \@Energy;
			$Data0->{"y$cup"} = $pPDOSupArray{$st};
			$Data0->{"x${cup}_Name"} = "Energy / eV";
			$Data0->{"y${cup}_Name"} = "$st";
			$cup++;
			if($this->IsSpinPolarized()) {
				my $ptdn = $pPDOSdnArray{$st};
print "  PDOS for $st down\n";
				$Data2->{"x$cdn"} = \@Energy;
				$Data2->{"y$cdn"} = $pPDOSdnArray{$st};
				$Data2->{"x${cdn}_Name"} = "Energy / eV";
				$Data2->{"y${cdn}_Name"} = "$st down";
				$cdn++;
			}
		}
		for(my $i = 0 ; $i < $nAtomType ; $i++) {
			my $atomtype = $AtomTypeList[$i];
			my $atomname = $atomtype->AtomName(); # AtomNameOnly();
my $aname = $atomname;
$aname =~ s/\d+//;
			my $ss = "$aname s";
			my $sp = "$aname p";
			my $sd = "$aname d";
			my $sf = "$aname f";
print "  PDOS for $ss\n";
print "  PDOS for $sp\n";
print "  PDOS for $sd\n";
print "  PDOS for $sf\n";
			my $psu = $pPDOSupArray{$ss};
			my $ppu = $pPDOSupArray{$sp};
			my $pdu = $pPDOSupArray{$sd};

			$Data0->{"x$cup"} = \@Energy;
			$Data0->{"y$cup"} = $pPDOSupArray{$ss};
			$Data0->{"x${cup}_Name"} = "Energy / eV";
			$Data0->{"y${cup}_Name"} = "$ss";
			$cup++;
			$Data0->{"x$cup"} = \@Energy;
			$Data0->{"y$cup"} = $pPDOSupArray{$sp};
			$Data0->{"x${cup}_Name"} = "Energy / eV";
			$Data0->{"y${cup}_Name"} = "$sp";
			$cup++;
			$Data0->{"x$cup"} = \@Energy;
			$Data0->{"y$cup"} = $pPDOSupArray{$sd};
			$Data0->{"x${cup}_Name"} = "Energy / eV";
			$Data0->{"y${cup}_Name"} = "$sd";
			$cup++;
			$Data0->{"x$cup"} = \@Energy;
			$Data0->{"y$cup"} = $pPDOSupArray{$sf};
			$Data0->{"x${cup}_Name"} = "Energy / eV";
			$Data0->{"y${cup}_Name"} = "$sf";
			$cup++;
			if($this->IsSpinPolarized()) {
				my $psd = $pPDOSdnArray{$ss};
				my $ppd = $pPDOSdnArray{$sp};
				my $pdd = $pPDOSdnArray{$sd};
print "  PDOS for $ss down\n";
print "  PDOS for $sp down\n";
print "  PDOS for $sd down\n";
print "  PDOS for $sf down\n";

				$Data2->{"x$cdn"} = \@Energy;
				$Data2->{"y$cdn"} = $pPDOSdnArray{$ss};
				$Data2->{"x${cdn}_Name"} = "Energy / eV";
				$Data2->{"y${cdn}_Name"} = "$ss down";
				$cdn++;
				$Data2->{"x$cdn"} = \@Energy;
				$Data2->{"y$cdn"} = $pPDOSdnArray{$sp};
				$Data2->{"x${cdn}_Name"} = "Energy / eV";
				$Data2->{"y${cdn}_Name"} = "$sp down";
				$cdn++;
				$Data2->{"x$cdn"} = \@Energy;
				$Data2->{"y$cdn"} = $pPDOSdnArray{$sd};
				$Data2->{"x${cdn}_Name"} = "Energy / eV";
				$Data2->{"y${cdn}_Name"} = "$sd down";
				$cdn++;
				$Data2->{"x$cdn"} = \@Energy;
				$Data2->{"y$cdn"} = $pPDOSdnArray{$sf};
				$Data2->{"x${cdn}_Name"} = "Energy / eV";
				$Data2->{"y${cdn}_Name"} = "$sf down";
				$cdn++;
			}
#print "c=$cup, $cdn\n"
		}
	}
#my $p = $Data0->{x0};
#for(my $i = 0 ; $i < @$p ; $i++) {
#	print "$i: $p->[$i]\n";
#}

	$this->{EFDOS} = $EF;

	$Data0->CalMinMax();
	$this->{'DataArray'}->AddGraphData($Data0);
	$Data1->CalMinMax();
	$this->{'DataArray'}->AddGraphData($Data1);
	if($this->IsSpinPolarized()) {
#print "Add data for down spins\n";
		$Data2->CalMinMax();
		$this->{'DataArray'}->AddGraphData($Data2);
		$Data3->CalMinMax();
		$this->{'DataArray'}->AddGraphData($Data3);
	}

	return $path;
}

sub ReadEIGENVALtoArray
{
	my ($this, $EIGENVALPath, $EF) = @_;
	$EF = 0.0 if(!defined $EF);

	my $pHash = $this->ReadINCARtoHash($EIGENVALPath);
	my $ISPIN = $pHash->{ISPIN};

	my $in = new JFile($EIGENVALPath, "rb");
	return undef unless($in);

	my $line;
	for(my $i = 0 ; $i < 5 ; $i++) {
		$line = $in->ReadLine();
	}
	my $SampleName = $line;

	$line = $in->ReadLine();
	my ($nn, $nK, $nLevels) = Utils::Split("\\s+", $line);

	my @EnergyBands;
	my @KPointInf;
	for(my $iK = 0 ; $iK < $nK ; $iK++) {
#Skip blank line
		$line = $in->ReadLine();
		$line =~ s/^\s+//;
		$line =~ s/\s+$//;
		if($line eq '') {
			$line = $in->ReadLine();
		}
#print "line2: $line\n";
		$line =~ s/^.*?://;
		my ($kx, $ky, $kz, $w) = Utils::Split("\\s+", $line);
		$KPointInf[$iK] = [$kx, $ky, $kz, $w];
#print "k[$iK]: $kx $ky $kz\n";
		for(my $iL = 0 ; $iL < $nLevels ; $iL++) {
			$line = $in->ReadLine();
#print "l: $line";
			my ($ib, $eup, $edn) = Utils::Split("\\s+", $line);
			$EnergyBands[$iL][$iK][0] = $eup - $EF;
			if($ISPIN == 2) {
#print "Spin Polarized\n";
				$this->SpinPolarized(1);
				$EnergyBands[$iL][$iK][1] = $edn - $EF;
#print("VASP.pm: iL=$iL  iK=$iK\n");
			}
			else {
				$this->SpinPolarized(0);
			}
		}
	}
	$in->Close();
	return ($nn, $nK, $nLevels, \@EnergyBands, \@KPointInf);
}

sub ReadCSVBandFile
{
	my ($this, $path) = @_;
print "VASP::ReadCSVBandFile:Read band structure from CSV [$path]\n";

	$this->SpinPolarized(0);
print "  Read EIGENVAL from [$path]\n";
	my $in = new JFile($path, "rb");
	return undef unless($in);

	my $line = $in->ReadLine();
	my ($kname, $kx, $ky, $kz, $ktot, @E) = split(/\s*,\s*/, $line);
	my $nLevels = @E;
print "nLevels = $nLevels\n";
	my $EnergyBandArray = new EnergyBandArray;
#	$EnergyBandupArray->SetTitle($SampleName);
#	$EnergyBandupArray->SetCrystal($Crystal);
	for(my $iL = 0 ; $iL < $nLevels ; $iL++) {
		$EnergyBandArray->AddEnergyBand();
	}

	for(my $iK = 0 ; ; $iK++) {
		last if($in->eof());

		$line = $in->ReadLine();
		($kname, $kx, $ky, $kz, $ktot, @E) = split(/\s*,\s*/, $line);
		next if(!defined $E[0]);
#print "  **k[$iK]: $kx $ky $kz (ktot=$ktot)\n";

		for(my $iL = 0 ; $iL < $nLevels ; $iL++) {
			$EnergyBandArray->AddByKDE($iL, $kx, $ky, $kz, $ktot, $E[$iL]);

print "  **k[$iK]: $kx $ky $kz (ktot=$ktot) E($iL)=$E[$iL]\n";
		}
	}
	$EnergyBandArray->CalMinMax();
	$EnergyBandArray->GetBandBoundary();
	$this->{EnergyBandupArray} = $EnergyBandArray;

	my $DataArray = new GraphDataArray;
	$this->SetDataArray($DataArray);
	$EnergyBandArray->SetGraphDataArray($DataArray);

	return $path;
}

sub ReadEIGENVAL
{
	my ($this, $path, $EF, $SkipNonZerow) = @_;
	$SkipNonZerow = 0 if(!defined $SkipNonZerow);
print "VASP::ReadEIGENVAL:Read EIGENVAL from [$path]\n";

	my $pHash = $this->ReadINCARtoHash($path);
	my $ISPIN = $pHash->{ISPIN};
	my $IsHF  = ($pHash->{LHFCALC} =~ /T/i)? 1 : 0;
print "  IsHF=$IsHF.\n";

	my $POSCAR = $this->GetPOSCARFileName($path);
print "  Read [$POSCAR].\n";
	my $Crystal = $this->ReadStructureFromCARFiles($POSCAR, 1);

	my $OUTCAR = $this->GetOUTCARFileName($path);
print "  Read [$OUTCAR].\n";
	$EF = $this->ReadFermiEnergy($OUTCAR) if(!defined $EF);
print "  EF: $EF eV\n";

	my $KPOINTS = $this->GetKPOINTSFileName($path);
print "  Read [$KPOINTS]\n";
	my ($nKPoints, $pKPOINTSLines, $pKPOINTSList, $iK0) = $this->ReadKListFromKPOINTS($KPOINTS);
	$iK0 = 0 if(!$IsHF and !$SkipNonZerow);
	$this->{iKPoint0} = $iK0;
	$this->{iKPoint1} = $nKPoints;
print "    iK range for band plot = $iK0 - $nKPoints\n";

print "  Read EIGENVAL from [$path]\n";
	my $in = new JFile($path, "rb");
	return undef unless($in);

	my $line;
	for(my $i = 0 ; $i < 5 ; $i++) {
		$line = $in->ReadLine();
	}
#print "line: $line\n";
	my $SampleName = $line;
	$SampleName =~ s/^\s*//;
	$SampleName =~ s/\s*$//;

	$line = $in->ReadLine();
	my ($nn, $nK, $nLevels) = Utils::Split("\\s+", $line);
print "    nn,nK,nLevels=$nn, $nK, $nLevels\n";
	my $EnergyBandupArray = new EnergyBandArray;
	my $EnergyBanddnArray = new EnergyBandArray;
	$EnergyBandupArray->SetTitle($SampleName);
	$EnergyBandupArray->SetCrystal($Crystal);
	$EnergyBanddnArray->SetTitle($SampleName);
	$EnergyBanddnArray->SetCrystal($Crystal);
	for(my $iL = 0 ; $iL < $nLevels ; $iL++) {
		$EnergyBandupArray->AddEnergyBand();
		$EnergyBanddnArray->AddEnergyBand();
	}

	my $Distance = 0.0;
	for(my $iK = 0 ; $iK < $nK ; $iK++) {
#Skip blank line
		$line = $in->ReadLine();
		$line =~ s/^\s+//;
		$line =~ s/\s+$//;
		if($line eq '') {
			$line = $in->ReadLine();
		}
#print "line2: $line\n";
		my ($kx, $ky, $kz, $w) = Utils::Split("\\s+", $line);
print "  **k[$iK/$iK0]: $kx $ky $kz (w=$w)\n";
		for(my $iL = 0 ; $iL < $nLevels ; $iL++) {
			$line = $in->ReadLine();
#print "l: $line";
			my ($ib, $eup, $edn) = Utils::Split("\\s+", $line);
			$EnergyBandupArray->AddByKE($iL, $kx, $ky, $kz, $eup-$EF) if($iK >= $iK0);
			if($ISPIN == 2) {
#print "Spin Polarized\n";
				$this->SpinPolarized(1);
				$EnergyBanddnArray->AddByKE($iL, $kx, $ky, $kz, $edn-$EF) if($iK >= $iK0);
			}
			else {
				$this->SpinPolarized(0);
			}
		}
	}
	$EnergyBandupArray->CalMinMax();
	$EnergyBandupArray->GetBandBoundary();
	$this->{EnergyBandupArray} = $EnergyBandupArray;
	if($this->IsSpinPolarized()) {
		$EnergyBanddnArray->CalMinMax();
		$EnergyBanddnArray->GetBandBoundary();
		$this->{EnergyBanddnArray} = $EnergyBanddnArray;
	}

	my $DataArray = new GraphDataArray;
	$this->SetDataArray($DataArray);

#print "Add up data\n";
	$EnergyBandupArray->SetGraphDataArray($DataArray);
	if($this->IsSpinPolarized()) {
#print "Add down data\n";
		$EnergyBanddnArray->SetGraphDataArray($DataArray);
	}

	return $path;
}

sub ReadFiles
{
	my ($this, $filename, $TargetData) = @_;
	$TargetData = 'EnergyConvergence' unless($TargetData);
	$this->ClearAll();
	my $FileType = $this->{'FileType'} = CheckFileType($filename);
print "Crystal::VASP: ReadFiles: FileType [$FileType]\n";
	$this->SetFileName($filename);
	if($FileType eq "VASP Output file") {
		return $this->ReadOUTCAR($filename, $TargetData);
#		return $this->ReadOUTCAR($filename, "EnergyLevels");
	}
	if($FileType eq "VASP DOSCAR file") {
		return $this->ReadDOSCAR($filename);
	}
	if($FileType eq "VASP EIGENVAL file") {
		return $this->ReadEIGENVAL($filename);
	}
	if($FileType eq "CSV Band file") {
		return $this->ReadCSVBandFile($filename);
	}

	return undef;
}

#============================================================
#　一般関数
#============================================================

sub RWIGSFromINCAR
{
	my ($this, $filename) = @_;
	my $in = new JFile($filename, "r");
	unless($in) {
		print "Error: Can not read [$filename].\n";
		return -1;
	}
	my $line;
	while(1) {
		$line = $in->SkipTo("RWIGS\\s*=");
		last if(!$line);
		next if($line =~ /#\s*RWIGS\s*=/);
		last;
	}

	($line) = ($line =~ /RWIGS\s*=\s*(\S.*\S)\s*$/);
	my @RWIGS = Utils::Split("\\s+", $line);
	$in->Close();
	return @RWIGS;
}

sub ReadAtomTypesFromPOTCAR
{
	my ($this, $filename, $IsPrint, $AllowDuplicate) = @_;
	$IsPrint = 1 unless(defined $IsPrint);
	$AllowDuplicate = 0 if(!defined $AllowDuplicate);
my @chars = qw(a b c d e f g h i j k l m n o p q r s t u v w x y z);
	my $in = new JFile($filename, "r");
	unless($in) {
		print "Error: Can not read [$filename].\n" if($IsPrint);
		return -1;
	}
	my %nAtomTypes;
	my @AtomTypes;
	while(!$in->eof()) {
		my $line = $in->ReadLine();
#print "line: $line\n";
#		my ($atomtype) = ($line =~ /(\S+)\s*?$/);
		my ($atomtype) = ($line =~ /^\s*\S+\s+(\S+)/);
		$atomtype =~ s/_.*//;
##print "at: $atomtype\n";
		if($nAtomTypes{$atomtype}) {
			if($AllowDuplicate) {
				my $key = $chars[$nAtomTypes{$atomtype}];
#print "key: $key\n";
				push(@AtomTypes, "${atomtype}{$key}");
			}
		}
		else {
			push(@AtomTypes, $atomtype);
		}
		$nAtomTypes{$atomtype}++;

		$in->SkipTo("End of Dataset");
	}
	$in->Close();
#for(my $i = 0 ; $i < @AtomTypes ; $i++) {
#	print "$i: $AtomTypes[$i]\n";
#}
	return @AtomTypes;
}

sub ReadStructureFromCARFiles
{
	my ($this, $CARDir, $IsPrint, %hash) = @_;
	$IsPrint = 1 unless(defined $IsPrint);
	my $a0 = Sci::a0();
#print "ReadStructureFromCARFiles\n";

	$CARDir = $this->GetCARDir($CARDir);
	my $INCARPath  = ($hash{INCARPath})?  $hash{INCARPath}  : Deps::MakePath($CARDir, "INCAR", 0);
	my $POSCARPath = ($hash{POSCARPath})? $hash{POSCARPath} : Deps::MakePath($CARDir, "POSCAR", 0);
	my $POTCARPath = Deps::MakePath($CARDir, "POTCAR", 0);

	my $pINCAR = $this->ReadINCARtoHash($INCARPath);
	my @Occ = Utils::Split("\\s+", $pINCAR->{VCA});

	my $crystal = new Crystal();

#print "VASP::ReadStructureFromCARFiles: Read [$POTCARPath]\n" if($IsPrint);
	my @Atomtypes = $this->ReadAtomTypesFromPOTCAR($POTCARPath, 1, 1, \@Occ);
	return -1 if(@Atomtypes == 0 or $Atomtypes[0] eq "-1");

	for(my $i = 0 ; $i < @Atomtypes ; $i++) {
		my $occ = (defined $Occ[$i])? $Occ[$i] : 1.0;
#print "aat occ=$occ\n";
		$crystal->AddAtomType($Atomtypes[$i], 0, $occ);
#print "\n===\nReadStructureFromCARFiles: Atom[$i]=$Atomtypes[$i]\n";
	}

#print "VASP::ReadStructureFromCARFiles: Read [$POSCARPath]\n" if($IsPrint);
	my $in = new JFile($POSCARPath, "r");
	unless($in) {
		print "Error: Can not read [$POSCARPath].\n" if($IsPrint);
		return -1;
	}

	my $line = $in->ReadLine();
	Utils::DelSpace($line);
	$crystal->SetCrystalName($line);
	$crystal->SetSampleName($line);
	$line = $in->ReadLine();
	my $l0 = $line + 0.0;
	$line = $in->ReadLine();
	my ($l11, $l12, $l13) = Utils::Split("\\s+", $line);
	$line = $in->ReadLine();
	my ($l21, $l22, $l23) = Utils::Split("\\s+", $line);
	$line = $in->ReadLine();
	my ($l31, $l32, $l33) = Utils::Split("\\s+", $line);
#print "l0=$l0\n";
#print "l1=$l11, $l12, $l13\n";
#print "l2=$l21, $l22, $l23\n";
#print "l3=$l31, $l32, $l33\n";
#print "l1=", $l11*$l0, $l12*$l0, $l13*$l0, "\n";
#print "l2=", $l21*$l0, $l22*$l0, $l23*$l0, "\n";
#print "l3=", $l31*$l0, $l32*$l0, $l33*$l0, "\n";
	$crystal->SetLatticeVector($l11*$l0, $l12*$l0, $l13*$l0, 
				   $l21*$l0, $l22*$l0, $l23*$l0, 
				   $l31*$l0, $l32*$l0, $l33*$l0);
	my ($a, $b, $c, $alpha, $beta, $gamma) = $crystal->CalculateLatticeConstantFromVector();
	$line = $in->PreReadLine();
	if($line !~ /^[\s\d]+$/) {
		$line = $in->ReadLine();
	}
	$line = $in->ReadLine();
	my @nAtomTypes = Utils::Split("\\s+", $line);
	$in->SkipTo("Direct");
	my %AtomCount;
	for(my $i = 0 ; $i < @nAtomTypes ; $i++) {
		my $atomtype = $Atomtypes[$i];
		my $occ = (defined $Occ[$i])? $Occ[$i] : 1.0;
#		$crystal->AddAtomType($atomtype, 0, $occ);
#print "it=$i:$atomtype occ=$occ\n";
		for(my $j = 0 ; $j < $nAtomTypes[$i] ; $j++) {
#print "  j=$j\n";
			$line = $in->ReadLine();
			my ($x, $y, $z) = Utils::Split("\\s+", $line);
			$AtomCount{$atomtype}++;
			my $Label = $atomtype . $AtomCount{$atomtype};
			$crystal->AddAtomSite($Label, $Label, $x, $y, $z, $occ);
#			$crystal->AddAtomSite($Label, $Label, $x, $y, $z, 1.0);
##			$crystal->AddAtomSite($Label, $atomtype, $x, $y, $z, 1.0);
		}
	}
	$crystal->SetLatticeParameters($a, $b, $c, $alpha, $beta, $gamma);

	$in->Close();
	$crystal->ExpandCoordinates();

	return $crystal;
}

sub GetPOTCARList
{
	my ($this, $name, $Functional, $UseDPP, $UseRecommendedPOTCAR) = @_;
	$UseRecommendedPOTCAR = 0 if(!defined $UseRecommendedPOTCAR);

	$name =~ s/\{[^\}]*\}//;

	my $PotDir;
	my $lcF = lc $Functional;
	if($lcF eq 'gga') {
		$PotDir = Deps::MakePath($VASPDBDir, "pot_GGA", 0);
	}
	elsif($lcF eq 'paw') {
		$PotDir = Deps::MakePath($VASPDBDir, "potpaw", 0);
	}
	elsif($lcF eq 'paw_gga') {
		$PotDir = Deps::MakePath($VASPDBDir, "potpaw_GGA", 0);
	}
	elsif($lcF eq 'paw_pbe') {
		$PotDir = Deps::MakePath($VASPDBDir, "potpaw_PBE", 0);
	}
	elsif($lcF eq 'paw_lda52') {
		$PotDir = Deps::MakePath($VASPDBDir, "potpaw_LDA52", 0);
	}
	elsif($lcF eq 'paw_pbe52') {
		$PotDir = Deps::MakePath($VASPDBDir, "potpaw_PBE52", 0);
	}
	elsif($lcF eq 'paw_lda54') {
		$PotDir = Deps::MakePath($VASPDBDir, "potpaw_LDA54", 0);
	}
	elsif($lcF eq 'paw_pbe54') {
		$PotDir = Deps::MakePath($VASPDBDir, "potpaw_PBE54", 0);
	}
	elsif($lcF eq 'paw_lda64') {
		$PotDir = Deps::MakePath($VASPDBDir, "potpaw_LDA64", 0);
	}
	elsif($lcF eq 'paw_pbe64') {
		$PotDir = Deps::MakePath($VASPDBDir, "potpaw_PBE64", 0);
	}
	else { #us
		$PotDir = Deps::MakePath($VASPDBDir, "pot", 0);
	}
	my $PotDir2 = Deps::MakePath($PotDir, "elements", 0);
	if(-e $PotDir2) {
		$PotDir = $PotDir2;
	}
print("POTCAR Dir [$PotDir] for functional [$Functional]\n");

	my @POTFiles;
	my %POTFileNames;
	if($UseRecommendedPOTCAR =~ /=/) {
		my @a = Utils::Split("\\s*,\\s*", $UseRecommendedPOTCAR);
		for(my $i = 0 ; $i < @a ; $i++) {
			my ($atom, $file) = Utils::Split("=", $a[$i]);
			$POTFileNames{$atom} = $file;
#print "$atom: $file\n";
		}
	}
	elsif($UseRecommendedPOTCAR) {
		if($PPPriority{$name}) {
			my $p = $PPPriority{$name};
			for(my $i = 0 ; $i < @$p ; $i++) {
				my $PotcarFile = Deps::MakePath($PotDir, [$p->[$i], "POTCAR"], 0);
#print "f[$i] [$PotcarFile]\n";
				if(-f $PotcarFile) {
#print " OK\n";
					push(@POTFiles, $PotcarFile);
				}
			}
			return @POTFiles if(@POTFiles > 0);
		}
	}

	my $PotcarDir;
	my $PotcarFile;
	if($POTFileNames{$name}) {
		$PotcarDir  = Deps::MakePath($PotDir,    $POTFileNames{$name}, 0);
		$PotcarFile = Deps::MakePath($PotcarDir, "POTCAR", 0);
#print "F: $PotcarFile\n";
		push(@POTFiles, $PotcarFile);
	}
	else {
		$PotcarDir = Deps::MakePath($PotDir, "${name}_d", 0);
		if($UseDPP eq '') {
			$PotcarDir = Deps::MakePath($PotDir, $name, 0);
		}

		$PotcarFile = Deps::MakePath($PotcarDir, "POTCAR", 0);
		unless(-e $PotcarFile) {
			$PotcarDir  = Deps::MakePath($PotDir,    $name,    0);
			$PotcarFile = Deps::MakePath($PotcarDir, "POTCAR", 0);
		}
		if(-e $PotcarFile) {
			push(@POTFiles, $PotcarFile);
		}
		foreach my $suffix ("_sv_new", "_sv", "_pv_new", "_pv", 
				    "_s_new", "_s", "_h_new", "_h", 
				    "_2_new", "_2", "_3_new", "_3", "1.25_new", "1.25") {
			$PotcarDir  = Deps::MakePath($PotDir, "${name}${suffix}",  0);
			$PotcarFile = Deps::MakePath($PotcarDir, "POTCAR", 0);
#print "f [$PotcarFile]\n";
			if(-e $PotcarFile) {
				push(@POTFiles, $PotcarFile);
			}
		}
	}

	return (@POTFiles);
}

sub MakePOTCARFile
{
	my ($this, $Crystal, $Functional, $UseDPP, $fname, $pParams) = @_;
#print "VASP Old Dir [$VASPProgOldDir]<br>\n";

print "<H1>Make POTCAR</H1>\n";
	my @AtomTypeList = $Crystal->GetCAtomTypeList();
#for(my $i = 0 ; $i < @AtomTypeList ; $i++) {
#	print "MakePOTCARFile: atom[$i] = ", $AtomTypeList[$i]->AtomNameOnly(), "$LF";
#}

	my $nAtomType = @AtomTypeList;

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

	my @POTFiles;
	for(my $i = 0 ; $i < $nAtomType ; $i++) {
		my $atom = $AtomTypeList[$i];
#		my $name = $atom->Label();
		my $name = $atom->AtomNameOnly();

		my @POTRecommendations = $this->GetPOTCARList(
				$name, $Functional, $UseDPP, $pParams->{UseRecommendedPOTCAR});
		if(@POTRecommendations == 0) {
			print "Error in VASP::MakePOTCARFile: Potential for [$name] not found.\n";
			print "     [$Functional][UseDPP=$UseDPP][UseRecommended=$pParams->{UseRecommendedPOTCAR}]\n";
			print("\nPress ENTER to exit>>\n");
			<STDIN>;
			exit();
		}
		my $PotcarFile = $POTRecommendations[0];
		push(@POTFiles, $PotcarFile);
#print "  [$PotcarFile] is used for [$name]\n";
		if(defined $fname) {
			unless(open(IN, "<$PotcarFile")) {
				print "<H2>Error: Can not read [$PotcarFile]!!</H2>\n";
				return -1;
			}
			print "<b>Use $PotcarFile for [$name].</b><br>\n";
			while(<IN>) {
				print OUT $_;
			}
			close(IN);
		}
	}
	if(defined $fname) {
		close(OUT);
	}
	return ($fname, @POTFiles);
}

sub ModifyPOSCARFile
{
	my ($this, $POSCARPath, $A, $a, $b, $c, $alpha, $beta, $gamma, %args) = @_;
	my $in = new JFile();
	my $plines = $in->ReadFileToLines($POSCARPath);
	my @lines = @$plines;
#print join(', ', @lines);

	my $Function = $this->Args()->GetGetArg("Function");
	$Function = 'scf' if($Function eq '');

#	my $out = new JFile("${POSCARPath}.new", "w");
	my $out = new JFile("${POSCARPath}", "w");
	if(!$out) {
		print "Error in VASP::ModifyPOSCARFile: Can not write to [$POSCARPath].\n";
		return undef;
	}
#Si                                      
	$out->print($lines[0]);

	my $Crystal = new Crystal;
	my $A0 = $lines[1] + 0.0;
	my ($a11, $a12, $a13) = Utils::Split("\\s+", $lines[2]);
	my ($a21, $a22, $a23) = Utils::Split("\\s+", $lines[3]);
	my ($a31, $a32, $a33) = Utils::Split("\\s+", $lines[4]);
#print "A: $A0, $a11, $a12, $a13, $a21, $a22, $a23, $a31, $a32, $a33\n";
	$Crystal->SetLatticeVector($A0*$a11, $A0*$a12, $A0*$a13,
				   $A0*$a21, $A0*$a22, $A0*$a23,
				   $A0*$a31, $A0*$a32, $A0*$a33);
	my ($a0, $b0, $c0, $alpha0, $beta0, $gamma0) = $Crystal->LatticeParameters();
#print "($a0, $b0, $c0, $alpha0, $beta0, $gamma0)\n";

print "Positions: $args{Positions}\n";
	my @a = Utils::Split('\s*;\s*', $args{Positions});
	my @NewXYZ;
	for(my $i = 0; $i < @a ; $i++) {
		my ($isite, $xyz, $coord) = Utils::Split('\s*:\s*', $a[$i]);
		if($xyz eq 'x') {
			$xyz = 0;
		}
		elsif($xyz eq 'y') {
			$xyz = 1;
		}
		elsif($xyz eq 'z') {
			$xyz = 2;
		}
		else {
			next;
		}

print "  $i: No.$isite xyz[$xyz] = $coord\n";
		$NewXYZ[$isite][$xyz] = $coord;
	}

	$A0     = ($A)? $A : $A0;
	$a0     = ($a)? $a : $a0;
	$b0     = ($b)? $b : $b0;
	$c0     = ($c)? $c : $c0;
	$alpha0 = ($alpha)? $alpha : $alpha0;
	$beta0  = ($beta)?  $beta  : $beta0;
	$gamma0 = ($gamma)? $gamma : $gamma0;
	$Crystal->SetLatticeParameters($a0, $b0, $c0, $alpha0, $beta0, $gamma0);
	($a11, $a12, $a13, $a21, $a22, $a23, $a31, $a32, $a33)
			= $Crystal->LatticeVectors();
#   5.40640000000000     
#     1.0113565403215667    0.0000000000000000   -0.0000000000000000
#     0.0000000000000000    1.0113565403215667   -0.0000000000000000
#     0.0000000000000000    0.0000000000000000    1.0113565403215667
	$out->print("$A0\n");
	$out->printf("    %21.18f   %21.18f   %21.18f\n", $a11/$A0, $a12/$A0, $a13/$A0);
	$out->printf("    %21.18f   %21.18f   %21.18f\n", $a21/$A0, $a22/$A0, $a23/$A0);
	$out->printf("    %21.18f   %21.18f   %21.18f\n", $a31/$A0, $a32/$A0, $a33/$A0);

	my $iLine = 5;
#   Cl   Na   Sn
	$out->print($lines[$iLine++]) if($lines[$iLine] !~ /^[\s\d]+$/);
#    8    8    1
#Selective dynamics
#Direct
	$out->print($lines[$iLine++]);
	if($Function ne 'phonon' and $lines[$iLine] =~ /dynamics/) {
		$iLine++;
	}
	else {
		$out->print($lines[$iLine++]) 
	}
	$out->print($lines[$iLine++]);
	for(my $is = 0 ; ; $is++) {
		my $idx = $is + $iLine;
		last if(!defined $lines[$idx]);

#print("$is: ", $lines[$idx]);
		my ($x, $y, $z, $idx, $idy, $idz) = Utils::Split('\s+', $lines[$idx]);
		if($idz eq 'T' or $idz eq 'F') {
		}
		else {
			$out->print($lines[$idx]);
			next;
		}

		if(defined $NewXYZ[$is]) {
			$x = $NewXYZ[$is][0] if(defined $NewXYZ[$is][0]);
			$y = $NewXYZ[$is][1] if(defined $NewXYZ[$is][1]);
			$z = $NewXYZ[$is][2] if(defined $NewXYZ[$is][2]);
		}

#		$out->print($lines[$idx]);
#      0.500000000    0.500000000    0.000000000 T T T
		printf("$is: %14.9f %14.9f %14.9f %s %s %s\n", $x, $y, $z, $idx, $idy, $idz);
		$out->printf("   %14.9f %14.9f %14.9f %s %s %s\n", $x, $y, $z, $idx, $idy, $idz);
	}
	$out->Close();
}

sub ModifyPOSCARFile2
{
	my ($this, $Function, $InputFile, $OutputFile, %hash) = @_;
#print "p=$Function, $InputFile, $OutputFile, %hash\n";

	unlink($OutputFile);

	my $plines = JFile->new()->ReadFileToLines($InputFile);
	my @lines = @$plines;
	return undef if(@lines == 0);

	my $out = new JFile;
	$out->Open($OutputFile, "w") or return undef;
	my $iline = 0;
#Si                                      
#   5.40640000000000     
#     1.0113565403215667    0.0000000000000000   -0.0000000000000000
#     0.0000000000000000    1.0113565403215667   -0.0000000000000000
#     0.0000000000000000    0.0000000000000000    1.0113565403215667
	$out->print($lines[$iline++]);
	$out->print($lines[$iline++]);
	$out->print($lines[$iline++]);
	$out->print($lines[$iline++]);
#   Si
#     8
	$out->print($lines[$iline++]);
	$out->print($lines[$iline]);
#print "l**: $lines[$iline]\n";
	if($lines[$iline] =~ /[a-zA-Z]/) {
		$iline++;
		$out->print($lines[$iline]);
	}
#print "l**: $lines[$iline]\n";
	my (@nAtoms) = Utils::Split("\\s+", $lines[$iline++]);
print("VASP::ModifyPOSCARFile: nAtoms: ", join(', ', @nAtoms), "\n");
#Selective dynamics
#Direct
	if($Function =~ /phonon/ and $lines[$iline] =~ /selective/i) {
		$iline++;
	}
	else {
		$out->print($lines[$iline++]) 
	}
	$out->print($lines[$iline++]);
	for(my $i = 0 ; $i < @nAtoms ; $i++) {
		for(my $j = 0 ; $j < $nAtoms[$i] ; $j++) {
			$out->print($lines[$iline++]);
		}
	}
	if(!$hash{InitializeVelocity}) {
		while(1) {
			last if(!defined $lines[$iline]);
			$out->print($lines[$iline++]);
		}
	}
	$out->Close();
}

sub ModifyINCARFile
{
	my ($this, $Function, $InputFile, $OutputFile, $AllSeparated, 
		$MinEnergy, $MaxEnergy, $DOSMeasuredFromEF,
		$FromScratch, $Precision, $KeepSymmetry, $nMesh, 
		$SpinPolarized, $PStress, $pParamHash,
		$NBANDS, $NELECT, $pParams) = @_;
	$Precision     = "Normal" if($Precision eq '');
	$KeepSymmetry  = 0        if($KeepSymmetry eq '');
	if($SpinPolarized ne '') {
		$SpinPolarized = 1 if($SpinPolarized eq 'No');
		$SpinPolarized = 2 if($SpinPolarized eq 'Yes');
	}
	if($pParams->{HybridFunctional} ne '') {
		$pParamHash->{LHFCALC} = ".TRUE.";
	}
#print "0: HybridFunctional=[$pParams->{HybridFunctional}]\n";
#print "0: LHFCALC=[$pParamHash->{LHFCALC}]\n";
	if(!defined $pParams->{METAGGA}) {
		$pParams->{METAGGA} = '';
	}

#	$PStress       = 0        if($PStress eq '');
	Utils::DelSpace($MinEnergy);
	Utils::DelSpace($MaxEnergy);
	$DOSMeasuredFromEF = 1 if(!defined $DOSMeasuredFromEF or $DOSMeasuredFromEF eq '');
	if($pParams->{HybridFunctional}) {
#		$pParams->{IBRION} = 5 if($pParams->{IBRION} == 7);
#		$pParams->{IBRION} = 6 if($pParams->{IBRION} == 8);
	}

#	my $v = $this->ReadKeyValue("NBANDS\\s*=\\s*(\\S+)", undef, $OUTCAR);
	my $POSCAR = $this->GetPOSCARFileName($InputFile);
	if($NBANDS) {
		my $v = $this->EstimateDefaultNBANDS($POSCAR);
		my $v = $this->GetTotalValenceElectrons($POSCAR);
		if($NBANDS =~ /^\-(.*)$/) {
			$NBANDS = $v - $1;
		}
		elsif($NBANDS =~ /^\+(.*)$/) {
			$NBANDS = $v + $1;
		}
		elsif($NBANDS =~ /^x(.*)$/) {
			$NBANDS = $v * $1;
		}
		elsif($NBANDS =~ /^\/(.*)$/) {
			$NBANDS = $v / $1;
		}
		$NBANDS = int($NBANDS);
print "NBANDS changed from $v to $NBANDS\n";
	}

	my $Crystal = $this->ReadStructureFromCARFiles($InputFile, 0);
#print "ak=$pParams->{aKProduct}\n";
	my ($nx, $ny, $nz, $NKREDX, $NKREDY, $NKREDZ) 
		= $this->AnalyzeaKProduct($Crystal, $pParams->{aKProduct});
	if($Function =~ /band/i) {
		$NKREDX = 1;
		$NKREDY = 1;
		$NKREDZ = 1;
	}
print "(nx,ny,nz)=($nx,$ny,$nz)\n";
print "(NKREDX,NKREDY,NKREDZ)=($NKREDX,$NKREDY,$NKREDZ)\n";
	$pParamHash->{NKREDX} = $NKREDX if(defined $NKREDX);
	$pParamHash->{NKREDY} = $NKREDY if(defined $NKREDY);
	$pParamHash->{NKREDZ} = $NKREDZ if(defined $NKREDZ);

	my $OUTCAR  = $this->GetOUTCARFileName($InputFile);
	my $EF = 0;
	if(!$DOSMeasuredFromEF) {
print("  OUTCAR : $OUTCAR\n");
		$EF = $this->ReadFermiEnergy($OUTCAR);
print("    EF = $EF eV\n");
		if($MinEnergy ne '' and $MaxEnergy ne '') {
			$MinEnergy += $EF;
			$MaxEnergy += $EF;
		}
print("    DOS range: $MinEnergy - $MaxEnergy eV\n");
	}

	unless(open(IN, "<$InputFile")) {
		print "Can not read $InputFile.$LF$LF";
		return;
	}
	my @in;
	while(not eof(IN)) {
		my $line = <IN>;
		chomp($line);
		push(@in, $line);
#print "l: $line [", scalar @in, "]\n";
	}
	close(IN);

	unless(open(OUT, ">$OutputFile")) {
		print "Can not write to $OutputFile.$LF$LF";
		return;
	}
	binmode(OUT);

#print "ENCUT=[$pParams->{ENCUT}]\n";
#exit;
	for(my $i = 0 ; $i < @in ; $i++) {
		my $line = $in[$i];
#print "line: $line\n";
		my ($a,$a2,$b,$c);
		my ($head,$key,$equal,$val,$rest);

		if($line =~ /#LDAUTYPE: 1:.*2:the same/) {
			print OUT "  #LDAUTYPE: 1: Rotational invaliant LSDA+U 4:the same as 1, but LDA+U\n";
			print OUT "  #          2: Dudarev's LSDA+U\n";
			$i++;
			print OUT "  #LDAUTYPE = 2\n";
			$i++;
			next;
		}

		if( ($head,$key,$equal,$val,$rest) = ($line =~ /^([\s#]*)([^=\s]*)(\s*=\s*)(\S.*)?$/) ) {
#print "key: $key = $val [$pParamHash->{$key}]\n";
			if(defined $pParamHash->{$key}) {
				if($key eq 'KPAR' and ($pParamHash->{KPAR} eq '' or $Function =~ /eDensity/)) {
					$head =~ s/#//;
					print OUT "$head#$key$equal$pParamHash->{$key}$rest\n";
				}
				elsif($key eq 'NPAR' and $pParamHash->{NPAR} eq '') {
					$head =~ s/#//;
					print OUT "$head#$key$equal$pParamHash->{$key}$rest\n";
				}
				else {
					$head =~ s/#//;
					print OUT "$head$key$equal$pParamHash->{$key}$rest\n";
				}
				next;
			}
		}

		my $AGGAC;
		if( ($a,$b,$c) = ($line =~ /^(\s*SYSTEM\s*=)(.+)$/) ) {
#print "b: $b\n";
			$b =~ s/\([^()]+\)/($Function)/;
			print OUT "$a$b\n";
			next;
		}
		elsif( ($a,$b,$c) = ($line =~ /^([\s#]*PREC\s*=\s*)(\S+)(.*)$/) ) {
			if($Precision !~ /Keep/i) {
				print OUT "${a}$Precision$c\n";
			}
			else {
				print OUT "$a$b$c\n";
			}
			next;
		}
		elsif( ($a,$b,$c) = ($line =~ /^([\s#]*ISPIN\s*=\s*)(\S+)(.*)$/) ) {
			if($SpinPolarized ne '') {
				print OUT "${a}$SpinPolarized$c\n";
			}
			else {
				print OUT "$a$b$c\n";
			}
			next;
		}
		elsif( ($a,$b,$c) = ($line =~ /^([\s#]*PSTRESS\s*=\s*)(\S+)(.*)$/) ) {
			if($PStress ne '') {
				$a =~ s/^(\s*)#/$1/;
				print OUT "${a}$PStress$c\n";
			}
			else {
				print OUT "$a$b$c\n";
			}
			next;
		}
		elsif(($a,$b,$c) = ($line =~ /^([\s#]*ENCUT\s*=\s*)(\S*)(.*)$/) ) {
			if($pParams->{ENCUT} ne '') {
				$a =~ s/^(\s*)#/$1/;
				print OUT "$a$pParams->{ENCUT}$c\n";
			}
			else {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			next;
		}
		elsif(($a,$b,$c) = ($line =~ /^([\s#]*NBANDS\s*=\s*)(\S*)(.*)$/) ) {
print "NBANDS=$NBANDS\n";
			if($NBANDS ne '') {
				$a =~ s/^(\s*)#/$1/;
				print OUT "$a$NBANDS$c\n";
			}
			else {
#				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			next;
		}
		elsif($NELECT ne '' and ($a,$b,$c) = ($line =~ /^([\s#]*NELECT\s*=\s*)(\S*)(.*)$/) ) {
			my $v = $this->ReadKeyValue("NELECT\\s*=\\s*(\\S+)", undef, $OUTCAR);
			if($NELECT =~ /^\-(.*)$/) {
				$NELECT = $v - $1;
			}
			elsif($NELECT =~ /^\+(.*)$/) {
				$NELECT = $v + $1;
			}
			elsif($NELECT =~ /^x(.*)$/) {
				$NELECT = $v * $1;
			}
			elsif($NELECT =~ /^\/(.*)$/) {
				$NELECT = $v / $1;
			}
print "NELECT changed from $v to $NELECT$c\n";
			$a =~ s/^(\s*)#/$1/;
			print OUT "$a$NELECT$c\n";
			next;
		}
		elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LHFCALC\s*=\s*)(\S*)(.*)$/) ) {
#print "LHFCALC=[$pParams->{LHFCALC}]\n";
#exit;
			if(defined $pParams->{LHFCALC}) {
				if($pParams->{LHFCALC} ne '') {
					$a =~ s/^(\s*)#/$1/;
					print OUT "${a}$pParams->{LHFCALC}$c\n";
				}
				else {
					$a =~ s/^(\s*)([^#\s])/$1#$2/;
					print OUT "$a$b$c\n";
#print "1 $a$b$c\n";
#exit;
				}
			}
			elsif(!defined $pParams->{LHFCALC}) {
#				$a =~ s/^(\s*)#/$1/;
#				print OUT "${a}$pParams->{LHFCALC}$c\n";
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
#print "2 $a$b$c\n";
#exit;
			}
			elsif(defined $pParams->{HybridFunctional}) {
				if($pParams->{HybridFunctional} eq '') {
					$a =~ s/^(\s*)([^#\s])/$1#$2/;
				}
				else {
					$a =~ s/^(\s*)#/$1/;
				}
				print OUT "${a}.TRUE.$c\n";
			}
			elsif($pParams->{HybridFunctional} ne '') {
				$a =~ s/^(\s*)#/$1/;
				print OUT "${a}.TRUE.$c\n";
			}
			else {
				$a =~ s/^(\s*)#/$1/;
				print OUT "${a}$pParams->{LHFCALC}$c\n";
			}
			next;
		}
		elsif( ($a,$b,$c) = ($line =~ /^([\s#]*AEXX\s*=\s*)(\S*)(.*)$/) ) {
			if(!defined $pParams->{HybridFunctional} or $pParams->{HybridFunctional} eq '') {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
			}
			else {
				$a =~ s/^(\s*)#/$1/;
			}
			print OUT "$a$b$c\n";
		}
		elsif( ($a,$b,$c) = ($line =~ /^([\s#]*AGGAC\s*=\s*)(\S*)(.*)$/) ) {
			$AGGAC = $b if(!defined $AGGAC);
			if(defined $pParams->{AGGAC}) {
				$AGGAC = $pParams->{AGGAC};
				$a =~ s/^(\s*)#/$1/;
			}
			elsif(!defined $pParams->{HybridFunctional} or $pParams->{HybridFunctional} eq '') {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
			}
			else {
				$a =~ s/^(\s*)#/$1/;
			}
			print OUT "$a$AGGAC$c\n";
		}
		elsif( ($a,$b,$c) = ($line =~ /^([\s#]*ALDAC\s*=\s*)(\S*)(.*)$/) ) {
			if(!defined $pParams->{HybridFunctional} or $pParams->{HybridFunctional} eq '') {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
			}
			else {
				$a =~ s/^(\s*)#/$1/;
			}
			print OUT "$a$b$c\n";
		}
		elsif( ($a,$b,$c) = ($line =~ /^([\s#]*HFSCREEN\s*=\s*)(\S*)(.*)$/) ) {
			if(!defined $pParams->{HybridFunctional} 
				or $pParams->{HybridFunctional} eq ''
				or $pParams->{HybridFunctional} =~ /(PBE0|HF)/) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
			}
			else {
				$a =~ s/^(\s*)#/$1/;
			}
			print OUT "$a$b$c\n";
		}
		elsif( ($a,$b,$c) = ($line =~ /^([\s#]*PRECFOCK\s*=\s*)(\S*)(.*)$/) ) {
			if(!defined $pParams->{HybridFunctional} or $pParams->{HybridFunctional} eq '') {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
			}
			else {
				$a =~ s/^(\s*)#/$1/;
			}
			print OUT "$a$b$c\n";
		}
		elsif( ($a,$b,$c) = ($line =~ /^([\s#]*IMIX\s*=\s*)(\S*)(.*)$/) ) {
print "IMIX=$pParams->{IMIX}\n";
			if(defined $pParams->{IMIX} and $pParams->{IMIX} ne '') {
				$a =~ s/^(\s*)#/$1/;
				print OUT "$a$pParams->{IMIX}$c\n";
			}
			else {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
		}
		elsif( ($a,$b,$c) = ($line =~ /^([\s#]*NELM\s*=\s*)(\S*)(.*)$/) ) {
print "NELM=$pParams->{NELM}\n";
			if(defined $pParams->{NELM}) {
				$a =~ s/^(\s*)#/$1/;
				print OUT "$a$pParams->{NELM}$c\n";
			}
			else {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
		}
		elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LASPH\s*=\s*)(\S*)(.*)$/) ) {
print "LASPH=$pParams->{LASPH}\n";
			if(defined $pParams->{LASPH} and $pParams->{LASPH} eq '.TRUE.') {
				$a =~ s/^(\s*)#/$1/;
				print OUT "$a$pParams->{LASPH}$c\n";
			}
			elsif(defined $pParams->{LASPH}) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$pParams->{LASPH}$c\n";
			}
			else {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
		}
		elsif( ($a,$b,$c) = ($line =~ /^([\s#]*ADDGRID\s*=\s*)(\S*)(.*)$/) ) {
print "ADDGRID=$pParams->{ADDGRID}\n";
			if(defined $pParams->{ADDGRID}) {
				$a =~ s/^(\s*)#/$1/;
				print OUT "$a$pParams->{LASPH}$c\n";
			}
			else {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
		}
		elsif( ($a,$b,$c) = ($line =~ /^([\s#]*METAGGA\s*=\s*)(\S*)(.*)$/) ) {
print "METAGGA=$pParams->{METAGGA}\n";
			if(defined $pParams->{METAGGA} and $pParams->{METAGGA} ne '') {
				$a =~ s/^(\s*)#/$1/;
				print OUT "$a$pParams->{METAGGA}$c\n";
			}
			elsif(defined $pParams->{METAGGA}) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$pParams->{METAGGA}$c\n";
			}
			else {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
		}
		elsif( ($a,$b,$c) = ($line =~ /^([\s#]*CMBJA\s*=\s*)(\S*)(.*)$/) ) {
print "CMBJA=$pParams->{CMBJA}\n";
			if(defined $pParams->{CMBJA} and $pParams->{CMBJA} ne '') {
				$a =~ s/^(\s*)#/$1/;
				print OUT "$a$pParams->{CMBJA}$c\n";
			}
			elsif(defined $pParams->{CMBJA}) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$pParams->{CMBJA}$c\n";
			}
			else {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
		}
		elsif( ($a,$b,$c) = ($line =~ /^([\s#]*CMBJB\s*=\s*)(\S*)(.*)$/) ) {
print "CMBJB=$pParams->{CMBJB}\n";
			if(defined $pParams->{CMBJB} and $pParams->{CMBJB} ne '') {
				$a =~ s/^(\s*)#/$1/;
				print OUT "$a$pParams->{CMBJB}$c\n";
			}
			elsif(defined $pParams->{CMBJA}) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$pParams->{CMBJB}$c\n";
			}
			else {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
		}
		elsif( ($a,$b,$c) = ($line =~ /^([\s#]*CMBJ\s*=\s*)(\S*)(.*)$/) ) {
			$pParams->{CMBJ} =~ s/,/ /g if(defined $pParams->{CMBJ});
print "CMBJ=$pParams->{CMBJ}\n";
			if(defined $pParams->{CMBJ} and $pParams->{CMBJ} ne '') {
				$a =~ s/^(\s*)#/$1/;
				print OUT "$a$pParams->{CMBJ}\n";
			}
			elsif(defined $pParams->{CMBJ}) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$pParams->{CMBJ}\n";
			}
			else {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
		}
		elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LMAXTAU\s*=\s*)(\S*)(.*)$/) ) {
print "LMAXTAU=$pParams->{LMAXTAU}\n";
			if(defined $pParams->{LMAXTAU} and $pParams->{LMAXTAU} ne '') {
				$a =~ s/^(\s*)#/$1/;
				print OUT "$a$pParams->{CMBJ}$c\n";
			}
			elsif(defined $pParams->{LMAXTAU}) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$pParams->{LMAXTAU}$c\n";
			}
			else {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
		}
		elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LMIXTAU\s*=\s*)(\S*)(.*)$/) ) {
print "LMIXTAU=$pParams->{LMIXTAU}\n";
			if(defined $pParams->{LMIXTAU} and $pParams->{LMIXTAU} ne '') {
				$a =~ s/^(\s*)#/$1/;
				print OUT "$a$pParams->{LMAXTAU}$c\n";
			}
			elsif(defined $pParams->{LMIXTAU}) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$pParams->{LMAXTAU}$c\n";
			}
			else {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
		}
		elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LSORBIT\s*=\s*)(\S*)(.*)$/) ) {
print "SO=$pParams->{SpinOrbit}\n";
			if(defined $pParams->{SpinOrbit} and $pParams->{SpinOrbit} =~ /[1T]/i) {
				$a =~ s/^(\s*)#/$1/;
				print OUT "$a$pParams->{SpinOrbit}$c\n";
			}
			elsif(defined $pParams->{SpinOrbit}) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$pParams->{SpinOrbit}$c\n";
			}
			else {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
		}
		elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LNONCOLLINEAR\s*=\s*)(\S*)(.*)$/) ) {
			if(defined $pParams->{NonCollinear} and $pParams->{NonCollinear} =~ /[1T]/i) {
				$a =~ s/^(\s*)#/$1/;
				print OUT "$a$pParams->{NonCollinear}$c\n";
			}
			elsif(defined $pParams->{NonCollinear}) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$pParams->{NonCollinear}$c\n";
			}
			else {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
		}
		elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LMAXMIX\s*=\s*)(\S*)(.*)$/) ) {
print "LMM=$pParams->{LMAXMIX}\n";
			if(!defined $pParams->{LMAXMIX} or $pParams->{LMAXMIX} eq '') {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif(defined $pParams->{LMAXMIX} and $pParams->{LMAXMIX} ne '') {
				$a =~ s/^(\s*)#/$1/;
				print OUT "$a$pParams->{LMAXMIX}$c\n";
			}
			else {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$pParams->{LMAXMIX}$c\n";
			}
		}
		elsif($Function =~ /band/i or $Function =~ /FS/i) {
			if( ($a,$b,$c) = ($line =~ /^([\s#]*ISTART\s*=\s*)(\S+)(.*)$/) ) {
				print OUT "${a}1$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*ICHARG\s*=\s*)(\S+)(.*)$/) ) {
				if($pParams->{HybridFunctional} ne '' or $pParams->{METAGGA} ne '') {
					$a =~ s/^(\s*)#/$1/;
					print OUT "${a}0$c\n";
				}
				else {
					$a =~ s/^(\s*)#/$1/;
#					print OUT "${a}10$c\n";
					print OUT "${a}11$c\n";
				}
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*INIWAV\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*IBRION\s*=\s*)(\S+)(.*)$/)) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LEPSILON\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LRPA\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LOPTICS\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*CSHIFT\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*NSW\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*TIME\s*=\s*)(\S+)(.*)$/) ) {
				if($pParams->{HybridFunctional} ne '') {
					$a =~ s/^(\s*)([^#\s])/$1$2/;
					print OUT "$a$b$c\n";
				}
				else {
					$a =~ s/^(\s*)([^#\s])/$1#$2/;
					print OUT "$a$b$c\n";
				}
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*POTIM\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*ISIF\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*MDALGO\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LANGEVIN_GAMMA\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LANGEVIN_GAMMA_L\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*PMASS\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*ISYM\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/#//;
				print OUT "${a}1$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*SMASS\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*ISMEAR\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/#//;
# tetrahedron method can not be used for strings of k-points
				print OUT "${a}-1$c\n";
#				print OUT "${a}$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*SIGMA\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/#//;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$a2,$b,$c) 
				= ($line =~ /^([\s#]*(LPARD|NBMOD|EINT|IBNAD|KPUSE|LSEPB|LSEPK)\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			else {
				print OUT "$line\n";
			}
		}
		elsif($Function =~ /scf/i or $Function =~ /phonon/i) {
			if( ($a,$b,$c) = ($line =~ /^([\s#]*ISTART\s*=\s*)(\S+)(.*)$/) ) {
				if($FromScratch) {
					print OUT "${a}0$c\n";
				}
				else {
					print OUT "${a}1$c\n";
				}
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*ICHARG\s*=\s*)(\S+)(.*)$/) ) {
				if($FromScratch) {
					$a =~ s/^(\s*)([^#\s])/$1#$2/;
					print OUT "$a$b$c\n";
				}
				else {
					$a =~ s/#//;
					print OUT "${a}0$c\n";
				}
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*INIWAV\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*IBRION\s*=\s*)(\S+)(.*)$/)) {
				if($Function =~ /phonon/ or $pParams->{CalVibration}) {
					$a =~ s/#//;
					$b = (defined $pParams->{IBRION})? $pParams->{IBRION} : 8;
				}
				else {
					$a =~ s/^(\s*)([^#\s])/$1#$2/;
				}
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*NPAR\s*=\s*)(\S+)(.*)$/)) {
				if($Function =~ /phonon/ or $pParams->{CalVibration}) {
					$a =~ s/^(\s*)([^#\s])/$1#$2/;
				}
				elsif(defined $pParams->{NPAR}) {
					$a =~ s/#//;
					$b = $pParams->{NPAR};
				}
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*MDALGO\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LANGEVIN_GAMMA\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LANGEVIN_GAMMA_L\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*PMASS\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*NWRITE\s*=\s*)(\S+)(.*)$/) ) {
				if($Function =~ /phonon/ or $pParams->{CalVibration}) {
					$a =~ s/#//;
					$b = 3;
				}
				else {
					$a =~ s/^(\s*)([^#\s])/$1#$2/;
				}
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LEPSILON\s*=\s*)(\S+)(.*)$/) ) {
				if($Function =~ /phonon/ or $pParams->{CalVibration}) {
					$a =~ s/#//;
					$b = '.TRUE.';
				}
				else {
					$a =~ s/^(\s*)([^#\s])/$1#$2/;
				}
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LRPA\s*=\s*)(\S+)(.*)$/) ) {
				if($pParams->{CalVibration}) {
					$a =~ s/#//;
					$b = (defined $pParams->{LRPA})? $pParams->{LRPA} : $b;
				}
				else {
					$a =~ s/^(\s*)([^#\s])/$1#$2/;
				}
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LOPTICS\s*=\s*)(\S+)(.*)$/) ) {
				if($pParams->{CalOptics}) {
					$a =~ s/#//;
					$b = '.TRUE.';
				}
				else {
					$a =~ s/^(\s*)([^#\s])/$1#$2/;
				}
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*CSHIFT\s*=\s*)(\S+)(.*)$/) ) {
				if($pParams->{CalOptics}) {
					$a =~ s/#//;
					$b = (defined $pParams->{CSHIFT})? $pParams->{CSHIFT} : $b;
				}
				else {
					$a =~ s/^(\s*)([^#\s])/$1#$2/;
				}
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*NSW\s*=\s*)(\S+)(.*)$/) ) {
				if($Function =~ /phonon/) {
					$a =~ s/#//;
					$b = 1;
				}
				else {
					$a =~ s/^(\s*)([^#\s])/$1#$2/;
				}
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*POTIM\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*ISIF\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*IBRION\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*ISYM\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/#//;
				print OUT "${a}1$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*SMASS\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*ISMEAR\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/#//;
				print OUT "${a}$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*SIGMA\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/#//;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$a2,$b,$c) 
				= ($line =~ /^([\s#]*(LPARD|NBMOD|EINT|IBNAD|KPUSE|LSEPB|LSEPK)\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			else {
				print OUT "$line\n";
			}
		}
		elsif($Function =~ /dos/i) {
			if( ($a,$b,$c) = ($line =~ /^([\s#]*ISTART\s*=\s*)(\S+)(.*)$/) ) {
				print OUT "${a}1$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*ICHARG\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/#//;
				if( $pParams->{METAGGA} eq '' ) {
#					print OUT "${a}10$c\n";
					print OUT "${a}11$c\n";
				}
				else {
					print OUT "${a}0$c\n";
				}
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*INIWAV\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*SIGMA\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/#//;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*IBRION\s*=\s*)(\S+)(.*)$/)) {
				if($pParams->{CalVibration}) {
					$a =~ s/#//;
					$b = (defined $pParams->{IBRION})? $pParams->{IBRION} : 8;
				}
				else {
					$a =~ s/^(\s*)([^#\s])/$1#$2/;
				}
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*MDALGO\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LANGEVIN_GAMMA\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LANGEVIN_GAMMA_L\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*PMASS\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LEPSILON\s*=\s*)(\S+)(.*)$/) ) {
				if($pParams->{CalVibration}) {
					$a =~ s/#//;
					$b = '.TRUE.';
				}
				else {
					$a =~ s/^(\s*)([^#\s])/$1#$2/;
				}
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LRPA\s*=\s*)(\S+)(.*)$/) ) {
				if($pParams->{CalVibration}) {
					$a =~ s/#//;
					$b = (defined $pParams->{LRPA})? $pParams->{LRPA} : $b;
				}
				else {
					$a =~ s/^(\s*)([^#\s])/$1#$2/;
				}
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LOPTICS\s*=\s*)(\S+)(.*)$/) ) {
				if($pParams->{CalOptics}) {
					$a =~ s/#//;
					$b = '.TRUE.';
				}
				else {
					$a =~ s/^(\s*)([^#\s])/$1#$2/;
				}
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*CSHIFT\s*=\s*)(\S+)(.*)$/) ) {
				if($pParams->{CalOptics}) {
					$a =~ s/#//;
					$b = (defined $pParams->{CSHIFT})? $pParams->{CSHIFT} : $b;
				}
				else {
					$a =~ s/^(\s*)([^#\s])/$1#$2/;
				}
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*NSW\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*POTIM\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*ISIF\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*IBRION\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*ISYM\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/#//;
				print OUT "${a}1$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*SMASS\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*ISMEAR\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/#//;
				print OUT "${a}$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*SIGMA\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/#//;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*EMIN\s*=\s*)(\S+)(.*)$/) ) {
#				if($MinEnergy ne '' and $MaxEnergy ne '' and $nMesh ne '') {
				if($MinEnergy ne '' and $MaxEnergy ne '') {
					$a =~ s/#//;
					print OUT "$a$MinEnergy$c\n";
				}
				else {
					print OUT "$a$b$c\n";
				}
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*EMAX\s*=\s*)(\S+)(.*)$/) ) {
#				if($MinEnergy ne '' and $MaxEnergy ne '' and $nMesh ne '') {
				if($MinEnergy ne '' and $MaxEnergy ne '') {
					$a =~ s/#//;
					print OUT "$a$MaxEnergy$c\n";
				}
				else {
					print OUT "$a$b$c\n";
				}
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*NEDOS\s*=\s*)(\S+)(.*)$/) ) {
#				if($MinEnergy ne '' and $MaxEnergy ne '' and $nMesh ne '') {
				if($nMesh ne '') {
					$a =~ s/#//;
					print OUT "$a$nMesh$c\n";
				}
				else {
					print OUT "$a$b$c\n";
				}
			}
			elsif( ($a,$a2,$b,$c) 
				= ($line =~ /^([\s#]*(LPARD|NBMOD|EINT|IBNAD|KPUSE|LSEPB|LSEPK)\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			else {
				print OUT "$line\n";
			}
		}
		elsif($Function =~ /eDensity/i) {
			if( ($a,$b,$c) = ($line =~ /^([\s#]*ISTART\s*=\s*)(\S+)(.*)$/) ) {
				print OUT "${a}1$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*KPAR\s*=\s*)(\S+)(.*)$/) ) {
# Never entered to this line if --Param:KPAR= is active
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "#$a$b$c\n";
			}
#			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*PREC\s*=\s*)(\S+)(.*)$/) ) {
#				print OUT "${a}Normal$c\n";
#			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*ICHARG\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/#//;
				if( $pParams->{METAGGA} eq '') {
#					print OUT "${a}10$c\n";
					print OUT "${a}11$c\n";
				}
				else {
					print OUT "${a}0$c\n";
				}
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*INIWAV\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*IBRION\s*=\s*)(\S+)(.*)$/)) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*MDALGO\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LANGEVIN_GAMMA\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LANGEVIN_GAMMA_L\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*PMASS\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LEPSILON\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LRPA\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LOPTICS\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*CSHIFT\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*SIGMA\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/#//;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*NSW\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*POTIM\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*ISIF\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*ISYM\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/#//;
				print OUT "${a}1$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*SMASS\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*ISMEAR\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/#//;
				print OUT "${a}$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*SIGMA\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/#//;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LPARD\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/#//;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*NBMOD\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/#//;
				if($AllSeparated and ($MinEnergy eq '' or $MaxEnergy eq '')) {
					print OUT "${a}0$c\n";
				}
				else {
					if($DOSMeasuredFromEF) {
						print OUT "${a}-3$c\n";
					}
					else {
						print OUT "${a}-2$c\n";
					}
				}
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*EINT\s*=\s*)(\S+\s+\S+)(.*)$/) ) {
				if($AllSeparated and ($MinEnergy eq '' or $MaxEnergy eq '')) {
					$a =~ s/^(\s*)([^#\s])/$1#$2/;
					print OUT "$a$b$c\n";
				}
				else {
					$a =~ s/#//;
					print OUT "${a} $MinEnergy $MaxEnergy$c\n";
				}
			}
			elsif( ($a,$a2,$b,$c) = ($line =~ /^([\s#]*(LSEPB|LSEPK)\s*=\s*)(\S+)(.*)$/) ) {
				if($AllSeparated) {
					$a =~ s/#//;
				}
				else {
					$a =~ s/^(\s*)([^#\s])/$1#$2/;
				}
				print OUT "$a$b$c\n";
			}
			else {
				print OUT "$line\n";
			}
		}
		elsif($Function =~ /md/i) {
			if( ($a,$b,$c) = ($line =~ /^([\s#]*ISTART\s*=\s*)(\S+)(.*)$/) ) {
				if($FromScratch) {
					print OUT "${a}0$c\n";
				}
				else {
					print OUT "${a}1$c\n";
				}
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*ICHARG\s*=\s*)(\S+)(.*)$/) ) {
				if($FromScratch) {
					$a =~ s/^(\s*)([^#\s])/$1#$2/;
					print OUT "$a$b$c\n";
				}
				else {
					$a =~ s/#//;
					print OUT "${a}0$c\n";
				}
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*INIWAV\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*ISIF\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/#//;
				if($Function =~ /vc/i) {
					print OUT "${a}3$c\n";
				}
				else {
					print OUT "${a}0$c\n";
				}
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*IBRION\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/#//;
				print OUT "${a}0$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*NSW\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/#//;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*MDALGO\s*=\s*)(\S+)(.*)$/) ) {
				if($Function =~ /vc/i) {
					$a =~ s/#//;
					print OUT "${a}3\n";
				}
				else {
					$a =~ s/^(\s*)([^#\s])/$1#$2/;
					print OUT "$a$b$c\n";
				}
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LANGEVIN_GAMMA\s*=\s*)(\S+)(.*)$/) ) {
				if($Function =~ /vc/i) {
					$a =~ s/#//;
					print OUT "$a$b$c\n";
				}
				else {
					$a =~ s/^(\s*)([^#\s])/$1#$2/;
					print OUT "$a$b$c\n";
				}
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LANGEVIN_GAMMA_L\s*=\s*)(\S+)(.*)$/) ) {
				if($Function =~ /vc/i) {
					$a =~ s/#//;
					print OUT "$a$b$c\n";
				}
				else {
					$a =~ s/^(\s*)([^#\s])/$1#$2/;
					print OUT "$a$b$c\n";
				}
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*PMASS\s*=\s*)(\S+)(.*)$/) ) {
				if($Function =~ /vc/i) {
					$a =~ s/#//;
					print OUT "$a$b$c\n";
				}
				else {
					$a =~ s/^(\s*)([^#\s])/$1#$2/;
					print OUT "$a$b$c\n";
				}
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*ISYM\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/#//;
				print OUT "${a}$KeepSymmetry$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*POTIM\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/#//;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*SMASS\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/#//;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*TEBEG\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/#//;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*TEEND\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/#//;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*PSTRESS\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/#//;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*POMASS\s*=\s*)(\S+)(.*)$/) ) {
#				$a =~ s/#//;
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*ZVAL\s*=\s*)(\S+)(.*)$/) ) {
#				$a =~ s/#//;
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			else {
				print OUT "$line\n";
			}
		}
		elsif($Function =~ /relax/i) {
			if( ($a,$b,$c) = ($line =~ /^([\s#]*ISTART\s*=\s*)(\S+)(.*)$/) ) {
				if($FromScratch) {
					print OUT "${a}0$c\n";
				}
				else {
					print OUT "${a}1$c\n";
				}
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*ICHARG\s*=\s*)(\S+)(.*)$/) ) {
				if($FromScratch) {
					$a =~ s/^(\s*)([^#\s])/$1#$2/;
					print OUT "$a$b$c\n";
				}
				else {
					$a =~ s/#//;
					print OUT "${a}0$c\n";
				}
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*INIWAV\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*NSW\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/#//;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*ISIF\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/#//;
				if($Function =~ /vc-relax/i) {
					print OUT "${a}3$c\n";
				}
				else {
					print OUT "${a}2$c\n";
				}
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*MDALGO\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
#				$a =~ s/#//;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LANGEVIN_GAMMA\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
#				$a =~ s/#//;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LANGEVIN_GAMMA_L\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
#				$a =~ s/#//;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*PMASS\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
#				$a =~ s/#//;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*POMASS\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
#				$a =~ s/#//;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*ZVAL\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
#				$a =~ s/#//;
				print OUT "$a$b$c\n";
			}
#			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*IBRION\s*=\s*)(\S+)(.*)$/)) {
#				$a =~ s/#//;
##				$a =~ s/^(\s*)([^#\s])/$1#$2/;
#				print OUT "$a$b$c\n";
#			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*IBRION\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/#//;
				print OUT "${a}2$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*NBLOCK\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
#				$a =~ s/#//;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*TEBEG\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
#				$a =~ s/#//;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*TEEND\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
#				$a =~ s/#//;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*POTIM\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/#//;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*ISYM\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/#//;
				print OUT "${a}$KeepSymmetry$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*SMASS\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/#//;
				print OUT "${a}-3$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LEPSILON\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LRPA\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*LOPTICS\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*CSHIFT\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*ISMEAR\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/#//;
#				if($pParams->{HybridFunctional} ne '') {
#					print OUT "${a}$pParams->{ISMEARHybridFunctional}$c\n";
#				}
#				else {
					print OUT "${a}$b$c\n";
#				}
#				$a =~ s/^(\s*)([^#\s])/$1#$2/;
#				print OUT "$a$b$c\n";
			}
			elsif( ($a,$b,$c) = ($line =~ /^([\s#]*SIGMA\s*=\s*)(\S+)(.*)$/) ) {
#				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				$a =~ s/#//;
#				if($pParams->{HybridFunctional} ne '') {
#					print OUT "${a}$pParams->{SIGMAHybridFunctionl}$c\n";
#				}
#				else {
					print OUT "$a$b$c\n";
#				}
			}
			elsif( ($a,$a2,$b,$c) 
				= ($line =~ /^([\s#]*(LPARD|NBMOD|EINT|IBAND|KPUSE|LSEPB|LSEPK)\s*=\s*)(\S+)(.*)$/) ) {
				$a =~ s/^(\s*)([^#\s])/$1#$2/;
				print OUT "$a$b$c\n";
			}
			else {
				print OUT "$line\n";
			}
		}
	}

	close(OUT);
}



sub ReadElementFromPOT
{
	my ($this, $potfile) = @_;
	unless(open(IN,"<$potfile")) {
		print "ReadElementFromPOT: Can not read [$potfile].$LF$LF";
		return -1;
	}
	binmode(IN);
	my $line = <IN>;
	my ($val) = ($line =~ /US\s*(\S+)/i);
	($val) = ($line =~ /PAW_PBE\s+(\S+)/i) if($val eq '');
	($val) = ($line =~ /PAW\s+(\S+)/i) if($val eq '');
	$val =~ s/\_.*$//;

	close(IN);
	return $val;
}

sub ReadRWIGSFromPOT
{
	my ($this, $potfile) = @_;
	unless(open(IN,"<$potfile")) {
		print "ReadRWIGSFromPOT: Can not read [$potfile].$LF$LF";
		return -1;
	}
	binmode(IN);
	my $val;
	while(!eof(IN)) {
		my $line = <IN>;
		($val) = ($line =~ /RWIGS\s*=\s*(\S+)\s+?wigner/i);
		last if(defined $val);
	}
	close(IN);
	return $val;
}

sub ReadSystemName
{
	my ($this, $infile) = @_;
	$infile = $this->GetINCARFileName($infile);
#print "in: $infile\n";
	my $pHash = $this->ReadINCARtoHash($infile);
	return undef if(!defined $pHash->{SYSTEM});
	$pHash->{SYSTEM} =~ s/(\s+\[.*$)//;
	Utils::DelSpace($pHash->{SYSTEM});
	return $pHash->{SYSTEM};
}

sub GetFileNameFromSystemName
{
	my ($this, $dir, $ext) = @_;
	my $fname = $this->ReadSystemName($dir);
	return undef if(!defined $fname);
#print "f: $fname\n";
	my ($drive, $directory, $filename, $ext1, $lastdir, $filebody) = Deps::SplitFilePath($dir);
	$ext = $ext1 if(!defined $ext);
	return Utils::ConvertToValidFileName(Deps::MakePath("$drive$directory", "$fname$ext", 0));
}

sub LatticeConversionForDensityFile
{
	my ($this, $InputVASPFile, $ConvertedVASPFile, $NewCrystal, $T, $tRT) = @_;

	my $Tinv = $tRT->inverse();
#&PrintMatrix("FCC to Primitive coordinate matrix:\n", $Tinv);

	my $pDensity = $this->ReadDensityFileToHash($InputVASPFile);
print "Title: [$pDensity->{Title}]\n";

	my ($a11,$a12,$a13,$a21,$a22,$a23,$a31,$a32,$a33) = $NewCrystal->LatticeVectors();
	$pDensity->{Base} = $a11;
	($pDensity->{a11}, $pDensity->{a12}, $pDensity->{a13}) = ($a11/$a11, $a12/$a11, $a13/$a11);
	($pDensity->{a21}, $pDensity->{a22}, $pDensity->{a23}) = ($a21/$a11, $a22/$a11, $a23/$a11);
	($pDensity->{a31}, $pDensity->{a32}, $pDensity->{a33}) = ($a31/$a11, $a32/$a11, $a33/$a11);

	my @AtomType = $NewCrystal->GetCAtomTypeList();
	my @site = $NewCrystal->GetCAsymmetricAtomSiteList();
	my %nAtoms;
	for(my $i = 0 ; $i < @site ; $i++) {
		my $type = lc $site[$i]->AtomNameOnly();
		$nAtoms{$type}++;
#print "i=$i: $type\n";
	}
	my @nAtoms;
	my @pos;
	my $ic = 0;
	for(my $i = 0 ; $i < @AtomType ; $i++) {
		my $name = lc $AtomType[$i]->AtomNameOnly();
		$nAtoms[$i] = $nAtoms{$name};
#print "i=$i: $name\n";
		for(my $is = 0 ; $is < @site ; $is++) {
			my $name1 = lc $site[$is]->AtomNameOnly();
			next if($name ne $name1);

			my ($x, $y, $z) = $site[$is]->Position();
			$pos[$ic] = [$x, $y, $z];
			$ic++;
		}
	}
	$pDensity->{pnAtoms}  = \@nAtoms;
	$pDensity->{pAtomPos} = \@pos;

	my $pPrimValues = $pDensity->{pValues};
	my @ConvValues;
	my $v = new Math::MatrixReal(3, 1);
	for(my $ix = 0 ; $ix < $pDensity->{nx} ; $ix++) {
		for(my $iy = 0 ; $iy < $pDensity->{ny} ; $iy++) {
			for(my $iz = 0 ; $iz < $pDensity->{nz} ; $iz++) {
				my $x = $ix / $pDensity->{nx};
				my $y = $iy / $pDensity->{ny};
				my $z = $iz / $pDensity->{nz};
				$v->assign(1, 1, $x);
				$v->assign(2, 1, $y);
				$v->assign(3, 1, $z);
				my $v2 = $T * $v;
				my ($x2, $y2, $z2) = ($v2->element(1,1), $v2->element(2,1), $v2->element(3,1));
				my $v = Algorism::InterpolateValuesInCube(
						$pDensity->{nx}, $pDensity->{ny},$pDensity->{nz},
						$pPrimValues, $x2, $y2, $z2);
				$ConvValues[$ix][$iy][$iz] = $v;
#printf "$ix,$iy,$iz: (%g,%g,%g) => (%g,%g,%g)\n", $x, $y, $z, $x2, $y2, $z2;
#print "($ix,$iy,$iz) => ($ix2,$iy2,$iz2)\n";
			}
		}
	}
	$pDensity->{pValues} = \@ConvValues;

	$this->SaveDensityFileFromHash($ConvertedVASPFile, $pDensity);
}

# Charge density * Vcell is returned in $hash{pValues}
sub ReadDensityFileToHash
{
	my ($this, $path, $CoreChargePath) = @_;

	my $pCore;
	if($CoreChargePath) {
print "Read Core Charge from [$CoreChargePath]\n";
		$pCore = $this->ReadDensityFileToHash($CoreChargePath);
		return undef if(!$pCore);
	}

print "Read Charge from [$path]\n";

	my %hash;

	my $in = new JFile;
	return undef if(!$in->Open($path, "r"));

	my $line = $in->ReadLine();
	$hash{Title} = Utils::DelSpace($line);
	$line = $in->ReadLine();
	$hash{Base} = $line + 0.0;
	($hash{a11}, $hash{a12}, $hash{a13}) = Utils::Split("\\s+", $in->ReadLine());
	($hash{a21}, $hash{a22}, $hash{a23}) = Utils::Split("\\s+", $in->ReadLine());
	($hash{a31}, $hash{a32}, $hash{a33}) = Utils::Split("\\s+", $in->ReadLine());
	$line = $in->ReadLine();
	if($line !~ /^(\d\s)$/) {
		$line =~ s/[\r\n]+$//;
		$hash{Atoms} = $line;
		$line = $in->ReadLine() 
	}
	my @nAtoms = Utils::Split("\\s+", $line);
	$hash{pnAtoms} = \@nAtoms;
	$line = $in->ReadLine();
	$hash{LatticeMode} = Utils::DelSpace($line);
	
	my @AtomPos;
	my $ic = 0;
	for(my $i = 0 ; $i < @nAtoms ; $i++) {
		for(my $j = 0 ; $j < $nAtoms[$i] ; $j++) {
			my ($x, $y, $z) = Utils::Split("\\s+", $in->ReadLine());
			$AtomPos[$ic] = [$x, $y, $z];
			$ic++;
		}
	}
	$hash{pAtomPos} = \@AtomPos;

	$in->ReadLine();
	($hash{nx}, $hash{ny}, $hash{nz}) = Utils::Split("\\s+", $in->ReadLine());
	my ($ix, $iy, $iz) = (0, 0, 0);
	my @t;
	my @v;
	while(1) {
		$line = $in->ReadLine();
		last if(!defined $line);
		my @a = Utils::Split("\\s+", $line);
		last if(@a == 0);
		
		push(@t, @a);
	}

	my $c = 0;
	my $sum = 0.0;
	for(my $iz = 0 ; $iz < $hash{nz} ; $iz++) {
#	for(my $ix = 0 ; $ix < $hash{nx} ; $ix++) {
	for(my $iy = 0 ; $iy < $hash{ny} ; $iy++) {
	for(my $ix = 0 ; $ix < $hash{nx} ; $ix++) {
#	for(my $iz = 0 ; $iz < $hash{nz} ; $iz++) {
		$v[$ix][$iy][$iz] = $t[$c] + 0.0;
		$sum += $v[$ix][$iy][$iz];
		$c++;
	}
	}
	}
	$hash{pValues} = \@v;

	$in->Close();

	if($pCore) {
print "Merge Valence and Core Charges\n";
		$sum = 0.0;
		for(my $iz = 0 ; $iz < $hash{nz} ; $iz++) {
		for(my $iy = 0 ; $iy < $hash{ny} ; $iy++) {
		for(my $ix = 0 ; $ix < $hash{nx} ; $ix++) {
#			$sum += $v[$ix][$iy][$iz];
			$v[$ix][$iy][$iz] += $pCore->{pValues}[$ix][$iy][$iz];
			$sum += $v[$ix][$iy][$iz];
		}
		}
		}
	}
	$hash{Sum} = $sum;
#	$hash{Sum} = Algorism::Integrate3D($hash{nx}, $hash{ny}, $hash{nz}, $hash{pValues});
#print "sum=$sum / $hash{Sum}\n";
#exit;

	return \%hash;
}

sub SaveDensityFileFromHash
{
	my ($this, $path, $phash) = @_;

	my $out = new JFile;
	return undef if(!$out->Open($path, "w"));

	$out->print("$phash->{Title}\n");
	$out->printf("%6.4f\n", $phash->{Base});
	$out->printf("   %10.6f  %10.6f  %10.6f\n", $phash->{a11}, $phash->{a12}, $phash->{a13});
	$out->printf("   %10.6f  %10.6f  %10.6f\n", $phash->{a12}, $phash->{a22}, $phash->{a23});
	$out->printf("   %10.6f  %10.6f  %10.6f\n", $phash->{a13}, $phash->{a32}, $phash->{a33});

	my $pnAtoms = $phash->{pnAtoms};
	for(my $i = 0 ; $i < @$pnAtoms ; $i++) {
		$out->printf(" %4d", $pnAtoms->[$i]);
	}
	$out->print("\n");
	$out->print("$phash->{LatticeMode}\n");

	my $pAtomPos = $phash->{pAtomPos};
	for(my $i = 0 ; $i < @$pAtomPos ; $i++) {
		$out->printf(" %9.6f %9.6f %9.6f\n", $pAtomPos->[$i][0],  $pAtomPos->[$i][1],  $pAtomPos->[$i][2]);
	}
	$out->print("\n");

	$out->printf(" %d %d %d\n", $phash->{nx}, $phash->{ny}, $phash->{nz});
	my $pv = $phash->{pValues};
	my $ic = 0;
	for(my $ix = 0 ; $ix < $phash->{nx} ; $ix++) {
		for(my $iy = 0 ; $iy < $phash->{ny} ; $iy++) {
			for(my $iz = 0 ; $iz < $phash->{nz} ; $iz++) {
				$out->printf(" %11.4g", $pv->[$ix][$iy][$iz]);
				$ic++;
				if($ic % 10 == 0) {
					$out->print("\n");
				}
			}
		}
	}

	$out->Close();
}

sub ReadINCARtoHash
{
	my ($this, $infile) = @_;

	$infile = $this->GetINCARFileName($infile);

	my $in = new JFile;
	if(!$in->Open($infile, "r")) {
		print("Error in ReadINCARtoHash: Can not read [$infile].\n");
		return undef;
	}
	my %hash;
	while(1) {
		my $line = $in->ReadLine(1);
		last if(!defined $line);
		next if($line =~ /^#/);
		$line =~ s/#.*$//;
		if($line =~ /^Source:\s*(\S.*)\s*/) {
			$hash{Source} = $1;
			next;
		}
		my ($key, $val) = ($line =~ /^\s*(\S+)\s*=\s*(.*)$/);
#		my ($key, $val) = ($line =~ /^\s*(\S+)\s*=\s*(\S*)/);
		$val =~ s/#.*$// if(defined $val);
		Utils::DelSpace($val);
		next if(!defined $key);
		$hash{$key} = $val;
#print("$key: $val\n");
	}
	$in->Close();

	if(defined $hash{ISPIN} and $hash{ISPIN} == 2) {
		$this->{IsSpinPolarized} = 1;
	}
	else {
		$this->{IsSpinPolarized} = 0;
	}
		return \%hash;
}

sub ReadKPOINTStoHash
{
	my ($this, $infile) = @_;
print "VASP::ReadKPOINTStoHash from [$infile]\n";

	my $pHash = $this->ReadINCARtoHash($infile);
	my $ISPIN = $pHash->{ISPIN};
	my $IsHF  = ($pHash->{LHFCALC} =~ /T/i)? 1 : 0;
print "VASP::ReadKPOINTStoHash: IsHF=$IsHF.\n";

	my $in = new JFile;
	if(!$in->Open($infile, "r")) {
		print("Error in ReadKPOINTStoHash: Can not read [$infile].\n");
		return undef;
	}

	my %hash;
	my $s = '';

	my $line = $in->ReadLine();
	$s .= $line;
	Utils::DelSpace($line);
	$hash{Title}    = $line;

	$line = $in->ReadLine();
	$s .= $line;
	Utils::DelSpace($line);
	$hash{nKPoints} = $line;

	$line = $in->ReadLine();
	$s .= $line;
	Utils::DelSpace($line);
	$hash{KPointMode} = $line;
	
	if($hash{KPointMode} =~ /^rec/i or $hash{KPointMode} =~ /^dir/i) {
		$hash{LatticeMode} = $line;
	}
	else {
		$line = $in->ReadLine();
		$s .= $line;
		Utils::DelSpace($line);
		$hash{LatticeMode} = $line;
	}

	while(1) {
		my $line = $in->ReadLine();
		last if(!defined $line);
		$s .= $line;
	}
	$in->Close();
	
	$hash{KPOINTS} = $s;
	return \%hash;
}

sub ReadPOTCARtoHash
{
	my ($this, $infile) = @_;

	$infile = $this->GetPOTCARFileName($infile);

	my $in = new JFile;
	if(!$in->Open($infile, "r")) {
		print("Error in ReadPOTCARtoHash: Can not read [$infile].\n");
		return undef;
	}
	my %hash = (FileName => $infile);
	my $count = 0;
	while(1) {
		my $PPKey = "PseudoPotential$count";

		my $line = $in->ReadLine();
		Utils::DelSpace($line);

		$hash{$PPKey}          = $line;
		$hash{PseudoPotential} = $line . "\n";

		my ($AtomName) = ($line =~ / ([A-Z][a-z]?)/);

		my $PPHashKey = "PseudoPotential${count}Hash";
		$line = $in->SkipTo("VRHFIN\\s*=\\s*(\\w.*)\$");
		my ($VRHFIN) = ($line =~ /VRHFIN\s*=\s*(\w.*)$/);
		$hash{$PPHashKey} = {AtomName => $AtomName, VRHFIN => $VRHFIN};

		while(1) {
			$line = $in->ReadLine();
			Utils::DelSpace($line);
			if($line eq 'Description') {
				last;
			}
#print "l: $line\n";
			if($line =~ /(\w+)\s*=\s*([\+\-\.\w]+)\s*;\s*(\w+)\s*=\s*([\+\-\.\w]+)/) {
				my ($key1, $val1, $key2, $val2) = ($1, $2, $3, $4);
				if($key1 eq $key2) {
					$key1 .= '1';
					$key2 .= '2';
				}
#print "$key1=$val1; $key2=$val2\n";
				$hash{$PPHashKey}{$key1} = $val1;
				$hash{$PPHashKey}{$key2} = $val2;
			}
			elsif($line =~ /(\w+)\s*=\s*([\+\-\.\w]+)/) {
				my ($key1, $val1) = ($1, $2);
#print "$key1=$val1\n";
				$hash{$PPHashKey}{$key1} = $val1;
			}
		}
		$line = $in->SkipTo("\\s*End of Dataset");
		last if(!defined $line);

		$count++;
	}
	$in->Close();
	return \%hash;
}

sub EstimateDefaultNBANDS
{
	my ($this, $infile) = @_;
	
	my $NIONS = $this->GetNIONS($infile, 1);
	my $NELE  = $this->GetTotalValenceElectrons($infile);
	my $NBANDS = int( ($NIONS + $NELE + 0.5) / 2.0 );
	$NBANDS = 8 if($NBANDS < 8);
	return $NBANDS;
}

sub GetNIONS
{
	my ($this, $infile, $ConsiderOccupancy) = @_;
	$ConsiderOccupancy = 0 if(!defined $ConsiderOccupancy);

	my $Crystal      = $this->ReadStructureFromCARFiles($infile);
#	my @AtomTypeList = $Crystal->GetCAtomTypeList();
	my @AtomSiteList = $Crystal->GetCExpandedAtomSiteList();
	return @AtomSiteList if(!$ConsiderOccupancy);
	
	my $n = 0;
	for(my $i = 0 ; $i < @AtomSiteList ; $i++) {
		my $site = $AtomSiteList[$i];
		$n += $site->Occupancy();
	}
	return $n;
}

sub GetTotalValenceElectrons
{
	my ($this, $infile) = @_;
	return undef if(!-f $infile);

#print "in:$infile GetTotalValenceElectrons\n";
	my $POTCAR = $this->GetPOTCARFileName($infile);
	my $POSCAR = $this->GetPOSCARFileName($infile);
#print "POSCAR=[$POSCAR]\n";
#print "POTCAR=[$POTCAR]\n";
	my $pPOT = $this->ReadPOTCARtoHash($POTCAR);
	my $pPOS = $this->ReadPOSCARtoHash($POSCAR);
	my $pn   = $pPOS->{pnAtomsForEachSite};
	my $pAtomType = $pPOS->{pAtomType};
#	my $pAtomSite = $pPOS->{pAtomSite};

	my $nVEL = 0;
	for(my $i = 0 ; $i < @$pn ; $i++) {
		my $atom = $pAtomType->[$i];
		my $occ  = $atom->Occupancy();
		my $key = "PseudoPotential$i";
		my $p = $pPOT->{"${key}Hash"};
		$nVEL += $pn->[$i] * $p->{ZVAL} * $occ;
#print "VASP:pm: GetTotalValenceElectrons: Atom#$i: n(site)=$pn->[$i]: "
#	."$p->{AtomName}: ZVAL = $p->{ZVAL}: nVEL = $nVEL\n";
	}
	return $nVEL;
}

sub GetFinalIonCharges
{
	my ($this, $infile) = @_;

	my $OUTCARPath     = $this->GetOUTCARFileName($infile);
	my $pOUT           = $this->ReadOUTCARtoHash($OUTCARPath);
	my $pCharge        = $pOUT->{pCharge};
	my $pMagnetization = $pOUT->{pMagnetization};
	my $pPOT           = $this->ReadPOTCARtoHash($infile);
	my $pPOS           = $this->ReadPOSCARtoHash($infile);
	my $pn             = $pPOS->{pnAtomsForEachSite};
	my $isite = 0;
	my (@AtomName, @SiteIndex, @Charge, @Magnetization, @ZVAL);
	for(my $i = 0 ; $i < @$pn ; $i++) {
		my $key = "PseudoPotential$i";
		my $p = $pPOT->{"${key}Hash"};
		for(my $j = 0 ; $j < $pn->[$i] ; $j++) {
print "$i: $pn->[$i]: $p->{AtomName}: nVEL = $p->{ZVAL} ($pCharge->[$isite])\n";
			$AtomName[$isite]      = $p->{AtomName};
			$ZVAL[$isite]          = $p->{ZVAL};
			$Charge[$isite]        = $p->{ZVAL} - $pCharge->[$isite] if(defined $pCharge->[$isite]);
			$Magnetization[$isite] = $pMagnetization->[$isite]       if(defined $pMagnetization->[$isite]);
			$SiteIndex[$isite]     = $i;
			$isite++;
		}
	}
	return (\@AtomName, \@SiteIndex, \@Charge, \@Magnetization, \@ZVAL);
}

sub ReadBaderOuttoHash
{

	my ($this, $dir) = @_;

	$dir = Deps::ExtractDirectory($dir) if(-f $dir);
	my $InFile = Deps::MakePath($dir, "Bader.out", 0);
print "Bader.out [$InFile]\n";
	return undef if(!-f $InFile);

	my $in = JFile->new($InFile, "r") or return undef;
	my $line = $in->SkipTo("AtomType:");
	my ($nAtomType) = ($line =~ /:\s*(\d+)/);
print "nAtomType=$nAtomType\n";
	my @AtomType;
	for(my $i = 0 ; $i < $nAtomType ; $i++) {
		$line = $in->ReadLine();
		my @a = Utils::Split(
		"\\s+", $line);
		$AtomType[$i]->{Name} = $a[1];
#print "  $i: $AtomType[$i]->{Name}\n";
	}

	$line = $in->SkipTo("Site:");
	my ($nAtomSite) = ($line =~ /:\s*(\d+)/);
print "nAtomSite=$nAtomSite\n";
	my @AtomSite;
	for(my $i = 0 ; $i < $nAtomSite ; $i++) {
		$line = $in->ReadLine();
		my @a = Utils::Split("[\\s\\(\\),]+", $line);
		$AtomSite[$i] = {
			Label => $a[1],
			x     => $a[2]+0,
			y     => $a[3]+0,
			z     => $a[4]+0,
			};
#print "  $i: $AtomSite[$i]->{Label} ($AtomSite[$i]->{x}, $AtomSite[$i]->{y}, $AtomSite[$i]->{z})\n";
	}

	$in->SkipTo("Ion charges");
	for(my $i = 0 ; $i < $nAtomSite ; $i++) {
		$line = $in->ReadLine();
		my ($Z) = ($line =~ /Charge\s*=\s*(\S+)/);
		$AtomSite[$i]->{BaderCharge} = $Z+0;
print "  $i: $AtomSite[$i]->{Label} ($AtomSite[$i]->{x}, $AtomSite[$i]->{y}, $AtomSite[$i]->{z}) Z=$AtomSite[$i]->{BaderCharge}\n";
	}
	$in->SkipTo("Ion charges re-adjusted");
	for(my $i = 0 ; $i < $nAtomType ; $i++) {
		$line = $in->ReadLine();
		my ($Z) = ($line =~ /Charge\s*=\s*(\S+)/);
		$AtomType[$i]->{BaderCharge} = $Z+0;
print "  $i: $AtomType[$i]->{Name} Z=$AtomType[$i]->{BaderCharge}\n";
	}

	my %hash = (
		DirName => $dir,
		nAtomType => $nAtomType,
		nAtomSite => $nAtomSite,
		pAtomType => \@AtomType,
		pAtomSite => \@AtomSite,
		);
	$hash{BaderCharges} = '';
	for(my $i = 0 ; $i < $nAtomSite ; $i++) {
		$hash{BaderCharges} .= "$AtomSite[$i]->{Label} "
			."($AtomSite[$i]->{x}, $AtomSite[$i]->{y}, $AtomSite[$i]->{z}): "
			. Utils::Round($AtomSite[$i]->{BaderCharge}, 4) . "\n";
	}
	for(my $i = 0 ; $i < $nAtomType ; $i++) {
		$hash{BaderCharges} .= "$AtomType[$i]->{Name}: " . Utils::Round($AtomType[$i]->{BaderCharge}, 4) . "\n";
	}

	return \%hash;
}

sub ReadMadelungPotentialtoHash
{
	my ($this, $dir) = @_;
	$dir = Deps::ExtractDirectory($dir) if(-f $dir);
	my $InFile1 = Deps::MakePath($dir, "LDEnergy.out", 0);
	my $InFile2 = Deps::MakePath($dir, "LDEnergy-modified.out", 0);
	my $p1 = $this->ReadMadelungPotentialtoHash1($InFile1);
	my $p2 = $this->ReadMadelungPotentialtoHash1($InFile2);

	my %hash = (
		DirName => $dir,
		pFormalCharge => $p1,
		pBaderCharge  => $p2,
		);
	if($p1) {
		my $key = 'MP(Formal Z)';
		$hash{$key} = '';
		my $nAtom = $p1->{nAtom};
		my $p     = $p1->{pAtomSite};
		for(my $i = 0 ; $i < $nAtom ; $i++) {
			my $p2 = $p->[$i];
			$hash{$key} .= "$p2->{Name}: " . Utils::Round($p2->{MP}, 4) . "\n";
		}
	}
	if($p2) {
		my $key = 'MP(Bader Z)';
		$hash{$key} = '';
		my $nAtom = $p2->{nAtom};
		my $p     = $p2->{pAtomSite};
		for(my $i = 0 ; $i < $nAtom ; $i++) {
			my $p2 = $p->[$i];
			$hash{$key} .= "$p2->{Name}: " . Utils::Round($p2->{MP}, 4) . "\n";
#print "$p2->{Name}: $p2->{MP} eV\n";
		}
	}

	return \%hash;
}

sub ReadMadelungPotentialtoHash1
{
	my ($this, $InFile) = @_;

print "Energy.out [$InFile]\n";
	return undef if(!-f $InFile);

	my $in = JFile->new($InFile, "r") or return undef;
	my $line = $in->SkipTo("<<ENERGY CALCULATION>>");
	my $nAtom = 0;
	my @AtomSite;
	while(1) {
		$line = $in->SkipTo("No\\.");
		last if(!defined $line);
		my ($idx, $Name) = ($line =~ /No\.\s*(\d+)\s+(\S+)/);
		$line = $in->SkipTo("Madelung potential:");
		my ($MP) = ($line =~ /:\s+(\S+)/);
		last if(!defined $MP);

		$AtomSite[$nAtom] = {
			idx  => $idx+0,
			Name => $Name,
			MP   => $MP+0,
			};
print "idx=$idx: Name=$Name  MP=$MP eV\n";
		$nAtom++;
	}

	my %hash = (
		nAtom     => $nAtom,
		pAtomSite => \@AtomSite,
		);

	return \%hash;
}

sub ReadPOSCARtoHash
{
	my ($this, $dir) = @_;
	$dir = Deps::ExtractDirectory($dir) if(-f $dir);
#print "dir:[$dir]\n";

	my $Crystal = $this->ReadStructureFromCARFiles($dir, 0);
#print "C[$Crystal]\n";
	if(!$Crystal or $Crystal <= 0) {
		print("Error in ReadPOSCARtoHash: Can not read [$dir].\n");
		return undef;
	}
	my %hash = (DirName => $dir);
	my ($a,$b,$c,$alpha,$beta,$gamma) = $Crystal->LatticeParametersByOutputMode(0);
	$hash{pCrystal} = $Crystal;
	$hash{a} = $a;
	$hash{b} = $b;
	$hash{c} = $c;
	$hash{alpha} = $alpha;
	$hash{beta}  = $beta;
	$hash{gamma} = $gamma;
	$hash{Latt}  = "$a\n$b\n$c\n$alpha\n$beta\n$gamma";
	$hash{Volume}    = $Crystal->Volume();
	$hash{Density}   = $Crystal->Density();
	$hash{FormulaUnit}            = $Crystal->FormulaUnit();
	$hash{SumChemicalComposition} = $Crystal->SumChemicalComposition();
	$hash{ChemicalComposition}    = $Crystal->ChemicalComposition();
	my @AtomTypeList = $Crystal->GetCAtomTypeList();
	my @AtomSiteList = $Crystal->GetCExpandedAtomSiteList();
	$hash{pAtomType} = \@AtomTypeList;
	$hash{AtomTypes} = join("\n", @AtomTypeList);
	$hash{nAtomType} = scalar @AtomTypeList;
	$hash{pAtomSite} = \@AtomSiteList;
	$hash{nAtomSite} = scalar @AtomSiteList;
	$hash{AtomSites} = '';
	for(my $i = 0 ; $i < @AtomSiteList ; $i++) {
		my $site = $AtomSiteList[$i];
		my $name = $site->AtomNameOnly();
		my ($x, $y, $z) = $site->Position();
		my $i1 = $i;
		$hash{AtomSites} .= ($hash{AtomSites} eq '')?
				sprintf("%03d: %2s:(%4.3g,%4.3g,%4.3g)",   $i1, $name, $x, $y, $z) :
				sprintf("\n%03d: %2s:(%4.3g,%4.3g,%4.3g)", $i1, $name, $x, $y, $z);
#printf("%03d: %2s:(%4.3g,%4.3g,%4.3g)",   $i1, $name, $x, $y, $z);
	}
	$hash{nAtomSite} = scalar @AtomSiteList;

	my $POSCARPath = $this->GetPOSCARFileName("$dir/POSCAR");
	my $in = new JFile($POSCARPath, "r");
	if($in) {
		$in->ReadLine();
		$in->ReadLine();
		$in->ReadLine();
		$in->ReadLine();
		$in->ReadLine();
		my $line = $in->ReadLine();
		if($line !~ /^[\s\d]+$/) {
			$line = $in->ReadLine();
		}
		my @n = Utils::Split("\\s+", $line);
		$hash{pnAtomsForEachSite} = \@n;
		$in->Close();
	}

#print "ReadPOSCARtoHash out\n";
	return \%hash;
}

sub ReadStructuresFromOUTCAR
{
	my ($this, $infile, $IsPrint) = @_;
	$IsPrint = 1 if(!defined $IsPrint);

	my $Crystal = $this->ReadStructureFromCARFiles($infile);
	my @AtomSiteList = $Crystal->GetCExpandedAtomSiteList();

	my $in = new JFile;
	$in->Open($infile, "r") or return ();

	my $line;
	my @Crystals;
	my $iCrystal = 0;
	while(1) {
		my $pos = $in->tell();

#  external pressure =      -80.85 kB  Pullay stress =        0.00 kB
		$line = $in->SkipTo("  external pressure =");
		my ($Pext, $PPullay) = ($line =~ /external pressure =\s*(\S+)\s*kB\s+Pullay stress =\s*(\S+)\s+kB/);
#print "*Pext,pullay: $Pext, $PPullay\n";
#  kinetic energy (ideal gas correction) =     39.58 kB
		$line = $in->SkipTo("  kinetic energy");
		my ($Pkin) = ($line =~ /=\s*(\S+)/);
#print "Pkin: $Pkin\n";
#  Total+kin.   -31.487     -46.886     -45.454      -6.709     -51.500     -11.152
		$line = $in->SkipTo("  Total\\+kin\\.");
		my ($P11, $P22, $P33, $P23, $P13, $P12) 
				= ($line =~ /kin\.\s*(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)/);
		my ($P32, $P31, $P21) = ($P23, $P13, $P12);
#print "Pij: $P11,$P22,$P33\n";

		$in->seek($pos);

# VOLUME and BASIS-vectors are now :
		$line = $in->SkipTo(" VOLUME and");
#		$line = $in->SkipTo(" VOLUME and BASIS-vectors are");
#print "**i=$iCrystal: l[$line]\n";
		if(!defined $line) {
print "Can not find \" VOLUME and BASIS-vectors are now \"\n";
			last;
		}
		$line = $in->SkipTo(" direct lattice vectors ");
		if(!defined $line) {
print "Can not find \" direct lattice vectors \"\n";
			last;
		}
		$Crystals[$iCrystal] = new Crystal;
		$Crystals[$iCrystal]->{Pexternal} = $Pext;
		$Crystals[$iCrystal]->{PPullay}   = $Pext;
		$Crystals[$iCrystal]->{Pkin}      = $Pkin;
		$Crystals[$iCrystal]->{P11}       = $P11;
		$Crystals[$iCrystal]->{P12}       = $P12;
		$Crystals[$iCrystal]->{P13}       = $P13;
		$Crystals[$iCrystal]->{P21}       = $P21;
		$Crystals[$iCrystal]->{P22}       = $P22;
		$Crystals[$iCrystal]->{P23}       = $P23;
		$Crystals[$iCrystal]->{P31}       = $P31;
		$Crystals[$iCrystal]->{P32}       = $P32;
		$Crystals[$iCrystal]->{P33}       = $P33;
#		$Crystals[$iCrystal]->{TMD}       = $TMD;
		
		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());
print "$iCrystal: Vector($a11, $a12, $a13, $a11, $a22, $a23, $a31, $a32, $a33)\n" if($IsPrint);
		$Crystals[$iCrystal]->SetLatticeVector($a11, $a12, $a13, $a21, $a22, $a23, $a31, $a32, $a33);

		$line = $in->SkipTo(" POSITION ");
		$in->ReadLine();
		for(my $i = 0 ; $i < @AtomSiteList ; $i++) {
			my $site = $AtomSiteList[$i];
			$line = $in->ReadLine(1);
			last if($line =~ /-----/);
			my ($xc, $yc, $zc, $fx, $fy, $fz) = Utils::Split("\\s+", $line);
			my ($x, $y, $z) = $Crystals[$iCrystal]->CartesianToFractional($xc, $yc, $zc);
print "    Pos=($x, $y, $z) Force=($fx, $fy, $fz)\n" if($IsPrint);
			$Crystals[$iCrystal]->AddAtomSite($site->Label(), $site->AtomNameOnly(), 
					$x, $y, $z, $site->Occupancy, $fx, $fy, $fz);
		}
		$Crystals[$iCrystal]->ExpandCoordinates();

		$pos = $in->tell();
#  FREE ENERGIE OF THE ION-ELECTRON SYSTEM (eV)
		$line = $in->SkipTo("\\s*FREE ENERGIE OF THE ION\\-ELECTRON SYSTEM");
		if($line) {
#  free  energy   TOTEN  =      -526.50050340 eV
			$line = $in->SkipTo("TOTEN\\s*=");
			($Crystals[$iCrystal]->{TOTEN}) = ($line =~ /TOTEN\s*=\s*(\S+)/);
#  energy  without entropy=     -526.52119211  energy(sigma->0) =     -526.50739963
			$line = $in->SkipTo("energy  without entropy=");
			($Crystals[$iCrystal]->{TOTENwithEntropy}) = ($line =~ /energy\(signma-\>0\)\s*=\s*(\S+)/);
		}
		else {
			$in->seek($pos);
		}

		$pos = $in->tell();
		$line = $in->SkipTo("\\s*ENERGY OF THE ELECTRON-ION-THERMOSTAT");
#print "l1[$line]\n";
		$line = $in->SkipTo(".*temperature\\s+(.*)\\s*K");
#print "l2[$line]\n";
		if($line) {
			my ($T) = ($line =~ /temperature\s+(.*?)\s*?K/);
			$Crystals[$iCrystal]->SetTemperature($T);
		}
		else {
			$in->seek($pos);
		}
#print "T=$T\n";
#  ENERGY OF THE ELECTRON-ION-THERMOSTAT SYSTEM (eV)
#  kin. lattice  EKIN_LAT=         0.777034  (temperature 2000.38 K)

		$iCrystal++;
	}
	$in->Close();
	return (@Crystals);
}

sub ReadOUTCARtoHash
{
	my ($this, $infile) = @_;

	my $in = new JFile;
	if(!$in->Open($infile, "r")) {
		print("Error in ReadOUTCARtoHash: Can not read [$infile].\n");
		return undef;
	}
	my (@charge, @magnetization, @Charge, @Magnetization, @PhononMode);
	my %hash = (FileName => $infile, pCharge => \@Charge, pMagnetization => \@Magnetization);
	while(1) {
		my $line = $in->ReadLine();
		last if(!defined $line);
#print "l: $line\n";

		if($line =~ /!!!/ and $line !~ /g\s+!!!/i) {
			for(my $i = 0 ; $i < 5 ; $i++) {
				my $line = $in->ReadLine();
			}
			my $s = '';
			for(my $i = 0 ; $i < 3 ; $i++) {
				my $s2 = $in->ReadLine();
				$s2 =~ s/^[\|\s]*//;
				$s2 =~ s/[\|\s]*$/ /;
				$s .= $s2;
			}
			Utils::DelSpace($s);
			$hash{Warning} = $s;
		}

		if($line =~ /static.*point symmetry\s+(\S+)/) {
			$hash{StaticSymmetry} = $1;
		}
		elsif($line =~ /Fermi energy:\s*(\S+);/) {
			$hash{EF} = $1;
		}
# E-fermi :  -2.4956     XC(G=0):  -3.7271     alpha+bet : -4.0683
		elsif($line =~ /E\-fermi\s*:\s*(\S+)/) {
			$hash{EF} = $1;
		}
		elsif($line =~ /dynamic.*point symmetry\s+(\S+)/) {
			$hash{DynamicSymmetry} = $1;
		}
		elsif($line =~ /constrained.*point symmetry\s+(\S+)/) {
			$hash{ConstrainedSymmetry} = $1;
		}
		elsif($line =~ /ENCUT\s*=\s*([\d\.Ee+\-]+)/) {
			$hash{ENCUT}  = $1;
		}
		elsif($line =~ /NBANDS\s*=\s*(\d+)/) {
			$hash{NBANDS}  = $1;
		}
		elsif($line =~ /NKPTS\s*=\s*(\d+).*NBANDS\s*=\s*(\d+)/) {
			$hash{NKPTS}  = $1;
			$hash{NBANDS} = $2;
		}
		elsif($line =~ /NEDOS\s*=\s*(\d+).*NIONS\s*=\s*(\d+)/) {
			$hash{NEDOS} = $1;
			$hash{NIONS} = $2;
		}
		elsif($line =~ /LDIM\s*=\s*(\d+).*LMDIM\s*=\s*(\d+)/) {
			$hash{LDIM}  = $1;
			$hash{LMDIM} = $2;
		}
		elsif($line =~ /NPLWV\s*=\s*(\d+)/) {
			$hash{NPLWV} = $1;
		}
		elsif($line =~ /dimension x,y,z NGX\s*=\s*(\d+).*NGY\s*=\s*(\d+).*NGZ\s*=\s*(\d+)/) {
			$hash{CurrNGXYZ} = "$1,$2,$3";
		}
		elsif($line =~ /I would recommend the setting:/) {
			$line = $in->ReadLine;
			if($line =~ /dimension x,y,z NGX\s*=\s*(\d+).*NGY\s*=\s*(\d+).*NGZ\s*=\s*(\d+)/) {
				$hash{RecommendNGXYZ} = "$1,$2,$3";
			}
		}
		elsif($line =~ /TEIN\s*=\s*(\d+)/) {
			$hash{TEIN} = $1;
		}
		elsif($line =~ /TEBEG\s*=\s*(\S+);.*TEEND\s*=\s*(\S+)/) {
			$hash{TEBEG} = $1;
			$hash{TEEND} = $2;
		}
		elsif($line =~ /POMASS\s*=\s*(.*?)\s*$/) {
			$hash{POMASS} = $1;
		}
		elsif($line =~ /ZVAL\s*=\s*(.*?)\s*$/) {
			$hash{ZVAL} = $1;
		}
		elsif($line =~ /NELECT\s*=\s*([+\-\d+\.eE]+)/) {
			$hash{NELECT} = $1;
		}
		elsif($line =~ /NUPDOWN\s*=\s*([\-\d]+)/) {
			$hash{NUPDOWN} = $1;
		}
		elsif($line =~ /IALGO\s*=\s*(\d+)/) {
			$hash{IALGO} = $1;
		}
		elsif($line =~ /LDIAG\s*=\s*(\S+)/) {
			$hash{LDIAG} = $1;
		}
		elsif($line =~ /IMIX\s*=\s*(\d+)/) {
			$hash{IMIX} = $1;
		}
		elsif($line =~ /AMIX\s*=\s*([\d\.]+)/) {
			$hash{AMIX} = $1;
			$line =~ /BMIX\s*=\s*([\d\.]+)/;
			$hash{BMIX} = $1;
		}
		elsif($line =~ /AMIX\_MAG\s*=\s*([\d\.]+)/) {
			$hash{AMIX_MAG} = $1;
			$line =~ /BMIX_MAG\s*=\s*([\d\.]+)/;
			$hash{BMIX_MAG} = $1;
		}
		elsif($line =~ /AMIN\s*=\s*([\d\.]+)/) {
			$hash{AMIN} = $1;
		}
		elsif($line =~ /TIME\s*=\s*(\d+)/) {
			$hash{TIME} = $1;
		}
		elsif($line =~ /LSECVAR\s*=\s*(\S+)/) {
			$hash{LSECVAR} = $1;
		}
		elsif($line =~ /IDIPOL\s*=\s*(\S+)/) {
			$hash{IDIPOL} = $1;
		}
		elsif($line =~ /LWAVE\s*=\s*(\S+)/) {
			$hash{LWAVE} = $1;
		}
		elsif($line =~ /LCHARG\s*=\s*(\S+)/) {
			$hash{LCHARG} = $1;
		}
		elsif($line =~ /LVTOT\s*=\s*(\S+)/) {
			$hash{LVTOT} = $1;
		}
		elsif($line =~ /LVHAR\s*=\s*(\S+)/) {
			$hash{LVHAR} = $1;
		}
		elsif($line =~ /LELF\s*=\s*(\S+)/) {
			$hash{LELF} = $1;
		}
		elsif($line =~ /LORBIT\s*=\s*(\S+)/) {
			$hash{LORBIT} = $1;
		}
		elsif($line =~ /LMONO\s*=\s*(\S+)/) {
			$hash{LMONO} = $1;
		}
		elsif($line =~ /LDIPOL\s*=\s*(\S+)/) {
			$hash{LDIPOL} = $1;
		}
		elsif($line =~ /IDIPOL\s*=\s*(\d+)/) {
			$hash{IDIPOL} = $1;
		}
		elsif($line =~ /EPSILON\s*=\s*([\d\.Ee+\-]+)/) {
			$hash{EPSILON} = $1;
		}
		elsif($line =~ /LEPSILON\s*=\s*(\S+)/) {
			$hash{LEPSILON} = $1;
		}
		elsif($line =~ /LRPA\s*=\s*(\S+)/) {
			$hash{LRPA} = $1;
		}
		elsif($line =~ /CSHIFT\s*=\s*([\d\.+\-Ee]+)/) {
			$hash{CSHIFT} = $1;
		}
		elsif($line =~ /GGA\s*=\s*(\S+)/) {
			$hash{EPSILON} = $1;
		}
		elsif($line =~ /LEXCH\s*=\s*(\d+)/) {
			$hash{LEXCH} = $1;
		}
		elsif($line =~ /LHFCALC\s*=\s*(\S+)/) {
			$hash{LHFCALC} = $1;
		}
		elsif($line =~ /PRECFOCK\s*=\s*(\S+)/) {
			$hash{PRECFOCK} = $1;
		}
		elsif($line =~ /LTHOMAS\s*=\s*(\S+)/) {
			$hash{LTHOMAS} = $1;
		}
		elsif($line =~ /LMAXFOCK\s*=\s*(\S+)/) {
			$hash{LMAXFOCK} = $1;
		}
		elsif($line =~ /AEXX\s*=\s*([\d+\-\.Ee]+)/) {
			$hash{AEXX} = $1;
		}
		elsif($line =~ /HFSCREEN\s*=\s*([\d+\-\.Ee]+)/) {
			$hash{HFSCREEN} = $1;
		}
		elsif($line =~ /ALDAX\s*=\s*([\d+\-\.Ee]+)/) {
			$hash{ALDAX} = $1;
		}
		elsif($line =~ /AGGAX\s*=\s*([\d+\-\.Ee]+)/) {
			$hash{AGGAX} = $1;
		}
		elsif($line =~ /ALDAC\s*=\s*([\d+\-\.Ee]+)/) {
			$hash{ALDAC} = $1;
		}
		elsif($line =~ /AGGAC\s*=\s*([\d+\-\.Ee]+)/) {
			$hash{AGGAC} = $1;
		}
		elsif($line =~ /ENCUT\s*=\s*(\S+)/) {
#   ENCUT  =  400.0 eV  29.40 Ry    5.42 a.u.   5.29  5.29  8.47*2*pi/ulx,y,z
			$hash{ENCUT} = $1;
		}
		elsif($line =~ /TOTEN\s*=\s*(\S+)/) {
			$hash{TOTEN} = $1;
		}
		elsif($line =~ /energy\s+without\s+entropy\s*=\s*(\S+).*sigma-\>0\)\s*=\s*(\S+)/) {
			$hash{TOTENwoS}    = $1;
			$hash{TOTENwoSlim} = $2;
		}
		elsif($line =~ /Total CPU time.*:\s*(\S+)/) {
			$hash{TotCPUTime} = $1;
		}
		elsif($line =~ /User time.*:\s*(\S+)/) {
			$hash{UserTime} = $1;
		}
		elsif($line =~ /System time.*:\s*(\S+)/) {
			$hash{SystemTime} = $1;
		}
		elsif($line =~ /Elapsed time.*:\s*(\S+)/) {
			$hash{ElapsedTime} = $1;
		}
		elsif($line =~ /(abort.*?)[\-\s]*$/) {
			$hash{AbortMessage} = $1;
		}
		elsif($line =~ /^(.*reached required accuracy.*)$/) {
			$hash{RequiredAccuracyMessage} = $1;
		}
		elsif($line =~ /^ *total charge +$/) {
#print "l2: [$line]\n";
			$hash{pChargeLines} = \@charge;
			$in->SkipTo("-----");
			my $nCharge = 0;
			while(1) {
				$line = $in->ReadLine();
				Utils::DelSpace($line);
				last if(!defined $line or $line eq '');
				last if($line =~ /-----/);
#print "l3: [$line]\n";
				my @a = Utils::Split("\\s+", $line);
#print "l3: [$line]: ", scalar @a, "\n";
				$charge[$nCharge]   = $line;
				$Charge[$nCharge++] = $a[@a-1];

			}
			$line = $in->ReadLine();
			$hash{TotalChargeLine} = $line;
		}
		elsif($line =~ /^ *magnetization \(x\) *$/) {
			$hash{pMagnetizationLines} = \@magnetization;
			$in->SkipTo("-----");
			my $nCharge = 0;
			while(1) {
				$line = $in->ReadLine();
				Utils::DelSpace($line);
				last if(!defined $line or $line eq '');
				last if($line =~ /-----/);

				my @a = Utils::Split("\\s+", $line);
				$magnetization[$nCharge]   = $line;
				$Magnetization[$nCharge++] = $a[@a-1];
			}
			$line = $in->ReadLine();
			$hash{TotalMagnetizationLine} = $line;
		}
		elsif($line =~ /Eigenvectors/) {
#		elsif($line =~ /^ *Eigenvectors and eigenvalues of the dynamical matrix *$/) {
			$hash{pPhononModeLines} = \@PhononMode;
			$line = $in->ReadLine();
			my $nPhononMode = 0;
			my $iline = 0;
			while(1) {
				$line = $in->ReadLine(1);
				last if(!defined $line);
				next if($line eq '');
				last if($line =~ /-----/);

				$PhononMode[$iline++] = $line;
#print "$nPhononMode: $line\n";
				if($line =~ /\d+\s+f(\/i)?\s*=\s+\d/) {
					$nPhononMode++;
				}
			}
			$hash{nPhononMode} = $nPhononMode;
		}

	}
	$in->Close();
	return \%hash;
}

#===================================================================
# DOS/OUTCARからエネルギー準位を読み取り、EVB, ECB, Egを読み取る
#===================================================================
sub ReadEgFromDOSOUTCAR
{
	my ($this, $DOSOUTCAR) = @_;

#print "\n=== From $DOSOUTCAR\n";
	my $in = new JFile;
	if(!$in->Open($DOSOUTCAR, "r")) {
		print "Error: Can not read [$DOSOUTCAR].\n";
		return ();
	}

	my $line = $in->SkipTo("E-fermi :");
	my ($EF) = ($line =~ /:\s*([\d+-\.eE]+)/);
#print "line: [$line] [EF = $EF]\n";
#print "   EF = $EF eV\n";

	my ($VBM, $CBM);
	while(1) {
		$line = $in->SkipTo("k-point\\s+\\d+\\s:");
		last if(!$line);

#print "line: $line";
		$line = $in->ReadLine();
		for(my $i = 0 ; ; $i++) {
			$line = $in->ReadLine();
			my ($iband, $E, $occ) = Utils::Split("\\s+", $line); #, "\\s+");
			last if(!defined $occ);
#print "  $iband: $E, $occ\n";

			if($occ > 0.1 and $E <= $EF) {
				if(!defined $VBM or $VBM < $E) {
					$VBM = $E;
#print "  VB: $iband: $E, $occ\n";
				}
			}
			elsif($occ < 0.1 and $E >= $EF) {
				if(!defined $CBM or $CBM > $E) {
					$CBM = $E;
#print "  CB: $iband: $E, $occ\n";
				}
			}
		}
	}
	my $Eg = $CBM - $VBM;

#print "   VBM=$VBM   CBM=$CBM   Eg=$Eg eV\n";

	$in->Close();
	return ($EF, $VBM, $CBM, $Eg);
}

#===================================================================
# DOS/*.csvからDOSを読み取り、EVB, ECB, Eg, Ne, Nspinを読み取る
#===================================================================
sub ReadFromDOSCSV
{
	my ($this, $Dir) = @_;

	my @DOSupFiles = glob(Deps::MakePath2($Dir, ["DOS", "*-1.csv"], 0));
	if(@DOSupFiles == 0) {
		@DOSupFiles = glob(Deps::MakePath2($Dir, ["DOS", "DOS-up.csv"], 0));
	}
	my @DOSdnFiles = glob(Deps::MakePath2($Dir, ["DOS", "*-2.csv"], 0));
	if(@DOSdnFiles == 0) {
		@DOSdnFiles = glob(Deps::MakePath2($Dir, ["DOS", "DOS-dn.csv"], 0));
	}

	my $DOSup = $DOSupFiles[0];
	my $DOSdn = $DOSdnFiles[0];
	my ($DEFup, $DEFdn, $NVBup, $NVBdn, $Ntotup, $Ntotdn);
	my ($EVBLowestup, $EVBLowestdn);
	my ($pXup, $pYup, $pXdn, $pYdn, $iXVBLowestup, $iXVBLowestdn);
	foreach my $file ($DOSup, $DOSdn) {
		next if(!-e $file);
#print "\n=== From $file\n";

		my $csv = new CSV;
		$csv->Read($file, 1);
		my $pLabelArray = $csv->LabelArray();
		my $pDataArray  = $csv->DataArray();
		my $nData = $csv->nData();
#print "  nData: $nData\n";
		my $pX = $pDataArray->[0];
		my $pY = $pDataArray->[1];
		last if(!defined $pY);
#print @$pLabelArray, "\n";
		my ($YMin, $YMax) = Utils::CalMinMax($pY);
		my $Yth = $YMax * 1.0e-3;
#print "  Y: $YMin - $YMax  Yth=$Yth\n";

		my $iXZero;
		for(my $i = 0 ; $i < $nData ; $i++) {
#print "$i: $pX->[$i]\n";
			if($pX->[$i] >= -1.0e-4) {
				$iXZero = $i;
#print "iXZero=$iXZero [$nData]\n";
				last;
			}
		}
		my $nzero = 0;
		my $iXVBLowest;
		for(my $i = $iXZero ; $i >= 0 ; $i--) {
			if($pY->[$i] < $Yth) {
				$nzero++;
#print "i=$i: nzero=$nzero\n";
				if($nzero > 10) {
					$iXVBLowest = $i;
					last;
				};
			}
			else {
				$nzero = 0;
			}
		}
#print "iXVBLowest=$iXVBLowest\n";
		if($file eq $DOSup) {
			$EVBLowestup = $pX->[$iXVBLowest];
			$iXVBLowestup = $iXVBLowest;
#print "iXVBLowestup=$iXVBLowestup\n";
			$pXup = $pX;
			$pYup = $pY;
		}
		else {
			$EVBLowestdn = $pX->[$iXVBLowest];
			$iXVBLowestdn = $iXVBLowest;
#print "iXVBLowestdn=$iXVBLowestdn\n";
			$pXdn = $pX;
			$pYdn = $pY;
		}
	}
	my $iXVBLowest = (!defined $iXVBLowestdn or $iXVBLowestup < $iXVBLowestdn)? $iXVBLowestup : $iXVBLowestdn;
	my $EVBLowest  = ($EVBLowestup < $EVBLowestdn)? $EVBLowestup : $EVBLowestdn;
#print "EVBLowest: $EVBLowest [$iXVBLowest] (up=$EVBLowestup [$iXVBLowestup]  dn=$EVBLowestdn [$iXVBLowestdn])\n";

	my ($pSup, $pSVBLowestup, $pSEFup);
	if($pYup) {
		my $csv = new CSV;
		$pSup         = Algorism::IntegrateByConstantStepSimpson(scalar @$pXup, $pXup->[1] - $pXup->[0], $pYup);
		$pSVBLowestup = $pSup->[$iXVBLowest];
		$pSEFup       = $csv->YVal($pXup, $pSup, 0.0);
#print "diff: $pSEFup - $pSVBLowestup\n";
		$NVBup        = $pSEFup - $pSVBLowestup;
		$DEFup        = $csv->YVal($pXup, $pYup, 0.0);
#print "  Total electrons in VB: up=$NVBup e\n";
#print "  D(EF):                 up=$DEFup e/eV\n";
	}

	my ($pSdn, $pSVBLowestdn, $pSEFdn);
	if($pYdn) {
		my $csv = new CSV;
		$pSdn         = Algorism::IntegrateByConstantStepSimpson(scalar @$pXdn, $pXdn->[1] - $pXdn->[0], $pYdn);
		$pSVBLowestdn = $pSdn->[$iXVBLowest];
		$pSEFdn       = $csv->YVal($pXdn, $pSdn, 0.0);
		$NVBdn        = $pSEFdn - $pSVBLowestdn;
		$DEFdn        = $csv->YVal($pXdn, $pYdn, 0.0);
#print "  Total electrons in VB: dn=$NVBdn e\n";
#print "  D(EF):                 dn=$DEFdn e/eV\n";
	}

	if($pYup and $pYdn) {
		my $Ntottotal  = $Ntotup + $Ntotdn;
		my $Ntotspin   = $Ntotup - $Ntotdn;
		my $NVBtotal   = $NVBup + $NVBdn;
		my $NVBspin    = $NVBup - $NVBdn;
		my $DEFtotal   = $DEFup + $DEFdn;
		my $DEFspin    = $DEFup - $DEFdn;
#printf "  Ntot(total)=%9.6f  NVB(total)=%9.6f  DEF(total)=%9.6f\n", $Ntottotal, $NVBtotal, $DEFtotal;
#printf "  Ntot(spin) =%9.6f  NVB(spin) =%9.6f  DEF(spin) =%9.6f\n", $Ntotspin, $NVBspin, $DEFspin;
	}
}

sub MakeINCARFile
{
	my ($this, $Crystal, $Function, $Functional, $UseDPP, $SpinPolarized, $fname, $pParams) = @_;

print "<H1>Make INCAR</H1>\n";
	if($pParams->{HybridFunctional}) {
#		$pParams->{IBRION} = 5 if($pParams->{IBRION} == 7);
#		$pParams->{IBRION} = 6 if($pParams->{IBRION} == 8);
	}
#print "f=$pParams->{SpinOrbit},$pParams->{NonCollinear},$pParams->{LMAXMIX}\n";

	my ($f, @POTFiles) = $this->MakePOTCARFile($Crystal, $Functional, $UseDPP, undef, $pParams);
	my @Element;
	for(my $i = 0 ; $i < @POTFiles ; $i++) {
		my $s = $this->ReadElementFromPOT($POTFiles[$i]);
#print "  <b>Element</b> from @POTFiles[$i]: $s$LF";
		push(@Element, $s);
	}
	my @RWIGS;
	for(my $i = 0 ; $i < @POTFiles ; $i++) {
		my $s = $this->ReadRWIGSFromPOT($POTFiles[$i]);
#print "  <b>RWIGS</b> for $Element[$i] from $POTFiles[$i]: $s$LF";
		push(@RWIGS, $s);
	}

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

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

	my $Sample = $this->SampleName();
	my $Volume = $Crystal->Volume();
	if(defined $Crystal->{'CIFFileName'} and $Crystal->{'CIFFileName'} ne '') {
		print OUT "Source: $Crystal->{'CIFFileName'}\n";
	}
	print OUT "SYSTEM = $Sample [$Functional] ($Function)\n";
	for(my $i = 0 ; $i < @POTFiles ; $i++) {
		print OUT "   #POTFiles: $POTFiles[$i]\n";
	}
	print OUT "\n";
	print OUT "Start parameters for this Run:\n";
	print OUT "  #ISTART: 0:from scratch 1:restart with the same Ecut\n";
	print OUT "  #   2:restart with the same basis set 3:full restart\n";
	print OUT "    ISTART = 0\n";
	print OUT "  #ICHARG: 0:Charge from initial wavef 1: CHGCAR 2:Atomic charge\n";
	print OUT "  #       +10: Non-selfconsistent\n";


	if($Function =~ /band|dos/i) {
		print OUT "    ICHARG = 11\n";
#		print OUT "    ICHARG = 10\n";
	}
	else {
		print OUT "    #ICHARG = 1\n";
	}
	print OUT "  #INIWAV: Only for ISTART=0\n";
	print OUT "    INIWAV = 1\n";
	print OUT "  NWRITE = 2 # Larger, more verbose. Must be 3 for phonon.\n";
	print OUT "  #PREC: High|Accurate|Normal|Medium|Low\n";
	if($Function =~ /vc-relax/i) {
		print OUT "    PREC = Accurate\n";
	}
	else {
		print OUT "    PREC = Normal\n";
	}
	print OUT "  #NBANDS: Set NBANDS to double for collinear calculation\n";
	print OUT "  #NBANDS = \n";
	if(defined $pParams->{ADDGRID}) {
		print OUT "  ADDGRID = $pParams->{ADDGRID}\n";
	}
	else {
		print OUT "  #ADDGRID = .TRUE.\n";
	}
	print OUT "\n";
	print OUT "Electronic Relaxation\n";
	print OUT "  #NELECT = \n";
	print OUT "  #Energy correction\n";
	print OUT "  EDIFF = 1e-04  stopping-criterion for ELM\n";
	print OUT "  # Critera for relaxation. In energy for positive value, force for negative value\n";
	print OUT "  EDIFFG = 1.0e-3\n";
	print OUT "  #LREAL: .TRUE.|.FALSE.|A|O\n";
	if($Volume > 7.0*7.0*7.0) {
		print OUT "  LREAL = .TRUE.\n";
	}
	else {
		print OUT "  LREAL = .FALSE.\n";
	}
	print OUT "  #ENCUT = 200.00 eV\n";
	print OUT "  #Algorithm 8:CG, 38: Davidson, 48:RMM=DIIS, 8,48: for NPAR=1, 48: for NPAR > 1\n";
	print OUT "  #Algorithm Normal  Damped  All  Diag\n";
	print OUT "  #  For hybrid: All, Diag, or Damped. ISMEAR=-5 is not recommended for All\n";
	print OUT "  #IALGO = 48\n";

	print OUT "  # DFT functional\n";
	print OUT "  # PE: PBE RE: revPBE RP: RPBE PS: PBEsol AM: AM05\n";
	print OUT "  # B3: B3LYP with VWN3  B5: B3LYP with VWN5\n";
	if(defined $pParams->{GGA} and $pParams->{GGA} ne '') {
		print OUT "    GGA=$pParams->{GGA}\n";
	}
	else {
		print OUT "    #GGA=PE\n";
	}
	print OUT "\n";

	my $HFALGO = 'All';
#	my $HFALGO = ($Function =~ /(dos|band|eDensity)/i)? '53' : 'Auto';
	print OUT "  #  For wannnier90: Use ALGO = None\n";
	print OUT "  #LWANNIER90_RUN = .FALSE.\n";
	if($pParams->{HybridFunctional} eq 'HF') {
		print OUT "  ALGO      = $HFALGO\n";
		print OUT "  #LDIAG     = .TRUE.\n";
		print OUT "  LHFCALC   = .TRUE.\n";
		print OUT "  AEXX      = 1.0\n";
		print OUT "  AGGAC     = 0.0\n";
		print OUT "  ALDAC     = 0.0\n";
		print OUT "  #HFSCREEN  = 0.2\n";
		print OUT "  #ENCUTFOCK = 0 #No lognger used in vasp.5.2.x\n";
		print OUT "  PRECFOCK  = N #PrecFock=L|M|F|N|A\n";
		print OUT "  #LMAXFOCK  = 4\n";
		print OUT "  TIME      = 0.4\n";
	}
	elsif($pParams->{HybridFunctional} =~ /^HSE/) {
		print OUT "  ALGO      = $HFALGO\n";
		print OUT "  #LDIAG     = .TRUE.\n";
		print OUT "  LHFCALC   = .TRUE.\n";
		print OUT "  AEXX      = 0.25\n";
		print OUT "  AGGAC     = 1.0\n";
		print OUT "  ALDAC     = 1.0\n";
		print OUT "  HFSCREEN  = 0.2\n";
		print OUT "  #ENCUTFOCK = 0\n";
		print OUT "  PRECFOCK  = N #PrecFock=L|M|F|N|A\n";
		print OUT "  #LMAXFOCK  = 4\n";
		print OUT "  TIME      = 0.4\n";
	}
	elsif($pParams->{HybridFunctional} eq 'PBE0') {
		print OUT "  ALGO      = $HFALGO\n";
		print OUT "  #LDIAG     = .TRUE.\n";
		print OUT "  LHFCALC   = .TRUE.\n";
		print OUT "  AEXX      = 0.25\n";
		print OUT "  AGGAC     = 1.0\n";
		print OUT "  ALDAC     = 1.0\n";
		print OUT "  #HFSCREEN  = 0.2\n";
		print OUT "  #ENCUTFOCK = 0\n";
		print OUT "  PRECFOCK  = N #PrecFock=L|M|F|N|A\n";
		print OUT "  #LMAXFOCK  = 4\n";
		print OUT "  TIME      = 0.4\n";
	}
	else {
		print OUT "  ALGO      = $HFALGO\n";
		print OUT "  #LDIAG     = .TRUE.\n";
		print OUT "  #LHFCALC   = .TRUE.\n";
		print OUT "  #AEXX      = 0.25\n";
		print OUT "  #AGGAC     = 1.0\n";
		print OUT "  #ALDAC     = 1.0\n";
		print OUT "  #HFSCREEN  = 0.2\n";
		print OUT "  #ENCUTFOCK = 0 #No lognger used in vasp.5.2.x\n";
		print OUT "  #PRECFOCK  = N #PrecFock=L|M|F|N|A\n";
		print OUT "  #LMAXFOCK  = 4\n";
		print OUT "  #TIME      = 0.4\n";
	}
	print OUT "  #NKRED     = 2\n";
	print OUT "  #NKREDX    =\n";
	print OUT "  #NKREDY    =\n";
	print OUT "  #NKREDZ    =\n";
	if(defined $pParams->{IMIX}) {
		print OUT "  #IMIX      = $pParams->{IMIX}\n";
	}
	else {
		print OUT "  #IMIX      = 4\n";
	}
	print OUT "  #AMIX      = 0.4\n";
	print OUT "  #BMIX      = 1.0\n";

	print OUT "  #NSIM = 4\n";
	print OUT "  #NELM: Number of electronic steps\n";
	if($Function =~ /band/i) {
		print OUT "    NELMIN = 0 # 5 at least for hybrid band calc\n";
		print OUT "    NELM   = 0\n";
		print OUT "    NELMDL = 3 # of ELM steps";
	}
	elsif($Function =~ /wannier/i) {
		print OUT "    #NELMIN = 0 # 5 at least for hybrid band calc\n";
		print OUT "    NELM   = 1\n";
		print OUT "    #NELMDL = 3 # of ELM steps";
	}
	else {
		print OUT "    NELMIN = 0 # 5 at least for hybrid band calc\n";
		if(defined $pParams->{NELM}) {
			print OUT "    #NELM   = $pParams->{NELM}\n";
		}
		else {
			print OUT "    #NELM   = 100\n";
		}
		print OUT "    #NELMDL = 3 # of ELM steps\n";
	}
	print OUT "\n";
	print OUT "  #IDIPOL = 4\n";
	print OUT "  #DIPOL  = 0.00000 0.00000 0.00000\n";
	print OUT "  #LDIPOL = .TRUE.\n";

	print OUT "  #Spin-orbit\n";
	if(defined $pParams->{SpinOrbit} and $pParams->{SpinOrbit} =~ /[1T]/i) {
		print OUT "    LSORBIT = $pParams->{SpinOrbit}\n";
		print OUT "    SAXIS = 0 0 1\n";
	}
	else {
		print OUT "    #LSORBIT = .TRUE.\n";
		print OUT "    #SAXIS = 0 0 1\n";
	}
	print OUT "  #Non-collinear setup\n";
	if(defined $pParams->{NonCollinear} and $pParams->{NonCollinear} ne '') {
		print OUT "    LNONCOLLINEAR = $pParams->{NonCollinear}\n";
	}
	else {
		print OUT "    #LNONCOLLINEAR = .TRUE.\n";
	}
	print OUT "    #MAGMOM = ";
	for(my $ia = 0 ; $ia < $nAtomType ; $ia++) {
		my $atomtype = $AtomTypeList[$ia];
		my $typename = $atomtype->AtomNameOnly();
		for(my $i = 0 ; $i < @AtomSiteList ; $i++) {
			my $atom      = $AtomSiteList[$i];
			my $atomname  = $atom->AtomNameOnly();
			next if(uc $typename ne uc $atomname);

			print OUT "0 0 1   ";
		}
	}
	print OUT "\n";
	print OUT "    #NBANDS: Set NBANDS to double of collinear calculation\n";
	print OUT "    #GGA_COMPAT = .FALSE.\n";
	print OUT "  #ISPIN: 1:Non-polarized 2:Spin-polarized\n";
	if($SpinPolarized =~ /Yes/i or $SpinPolarized >= 2) {
		print OUT "  ISPIN = 2\n";
	}
	else {
		print OUT "  ISPIN = 1\n";
	}
	print OUT "  #mmmmmm = ";
	for(my $ia = 0 ; $ia < $nAtomType ; $ia++) {
		my $atomtype = $AtomTypeList[$ia];
		my $typename = $atomtype->AtomNameOnly();
		for(my $i = 0 ; $i < @AtomSiteList ; $i++) {
			my $atom      = $AtomSiteList[$i];
			my $atomname  = $atom->AtomNameOnly();
			next if(uc $typename ne uc $atomname);

			printf OUT "%2s", $atomname;
		}
	}
	print OUT "\n";
	print OUT "  #MAGMOM = ";
	for(my $ia = 0 ; $ia < $nAtomType ; $ia++) {
		my $atomtype = $AtomTypeList[$ia];
		my $typename = $atomtype->AtomNameOnly();
		for(my $i = 0 ; $i < @AtomSiteList ; $i++) {
			my $atom      = $AtomSiteList[$i];
			my $atomname  = $atom->AtomNameOnly();
			next if(uc $typename ne uc $atomname);

			printf OUT " 0";
		}
	}
	print OUT "\n";
	print OUT "  #NUPDOWN = 0\n";

	print OUT "\n";
	print OUT "  #PAW Control:\n";
	print OUT "  #LMAXMIX: Default: 2\n";
	print OUT "  #      Use 4 and 6 for d and f orbitals in NSCF/SOC/NC/L(S)DA+U calculations\n";
	if(defined $pParams->{LMAXMIX} and $pParams->{LMAXMIX} ne '') {
		print OUT "  LMAXMIX = $pParams->{LMAXMIX}\n";
	}
	else {
		print OUT "  #LMAXMIX = 4\n";
	}
	print OUT "  #LMAXPAW: Default: 2*lmax\n";
	print OUT "  #         -1: No on-site correction  0: Only spherical terms\n";
	print OUT "  #LMAXPAW = 0\n";

	print OUT "\n";
	print OUT "Functional\n";
	if(defined $pParams->{LASPH}) {
		print OUT "  LASPH = $pParams->{LASPH}\n";
	}
	else {
		print OUT "  #LASPH = .TRUE.\n";
	}
	print OUT "\n";
	print OUT "  # van der Waals correction by DFT-D method of Grimme\n";
	print OUT "  #LVDW = .TRUE.\n";
	print OUT "  # Other vdW-DF\n";
	print OUT "  #LUSE_VDW = .TRUE.\n";
	print OUT "  #   vdWDFT (Dion): AGGAC=0, GGA=RE, no Zab_vdW\n";
	print OUT "  #   vdW2   : AGGAC=0, GGA=ML, no Zab_vdW=-1.8867)\n";
	print OUT "  #   optPBE : AGGAC=0, GGA=OR, no Zab_vdW)\n";
	print OUT "  #   optB88 : AGGAC=0, GGA=BO, PARAM1=0.1833333333, PARAM2=0.22, no Zab_vdW)\n";
	print OUT "  #   optB86b: AGGAC=0, GGA=MK, PARAM1=0.1234, PARAM2=1.0, no Zab_vdW)\n";
#	print OUT "  #AGGAC = 0.000\n";
	print OUT "  #PARAM1 = 0.1234\n";
	print OUT "  #PARAM2 = 1.0000\n";
	print OUT "  #Zab_vdW = -1.8867\n";

	print OUT "\n";
	print OUT "  #METAGGA: TPSS|RTPSS|M06L|MBJ\n";
	if(defined $pParams->{METAGGA} and $pParams->{METAGGA} ne '') {
		print OUT "  METAGGA=$pParams->{METAGGA}\n";
	}
	else {
		print OUT "  #METAGGA=MBJ\n";
	}
	print OUT "  #CMBJA=-0.012\n";
	print OUT "  #CMBJB=1.023 # bohr\n";
	print OUT "  #CMBJ: c parameters: One ently per atomic type: c_1 c_2 .. c_n\n";
	print OUT "  #CMBJ=";
	for(my $i = 0 ; $i < $nAtomType ; $i++) {
		print OUT " c_$i";
	}
	print OUT "\n";
	print OUT "  #LMAXTAU =  0\n";
	print OUT "  #LMIXTAU =  F\n";
	print OUT "\n";
	print OUT "  #LDA+U setup\n";
	print OUT "  #LDAU = .TRUE.\n";
	print OUT "  #LDAUTYPE: 1: Rotational invaliant LSDA+U 4:the same as 1, but LDA+U\n";
	print OUT "  #          2: Dudarev's LSDA+U\n";
	print OUT "  #LDAUTYPE = 2\n";
#	print OUT "  #LDAUL = 2 -1 -1\n";
	print OUT "  #LDAUL = 2 ";
	for(my $i = 1 ; $i < $nAtomType ; $i++) {
		print OUT "-1 ";
	}
	print OUT "\n";
#	print OUT "  #LDAUU = 4.0 0.0 0.0\n";
	print OUT "  #LDAUU = 4.0 ";
	for(my $i = 1 ; $i < $nAtomType ; $i++) {
		print OUT "0.0 ";
	}
	print OUT "\n";
#	print OUT "  #LDAUJ = 0.0 0.0 0.0\n";
	print OUT "  #LDAUJ = ";
	for(my $i = 0 ; $i < $nAtomType ; $i++) {
		print OUT "0.0 ";
	}
	print OUT "\n";
	print OUT "  #LDAUPRINT: 0:silent  1:medium  2:detail\n";
	print OUT "  #LDAUPRINT = 2\n";

	print OUT "\n";
	print OUT "Ionic Relaxation\n";
	if($Function =~ /relax/i) {
		print OUT "  #ISIF: 0:Ion relax only 2:Cell fixed/Calc.Pressure 3:VC-Relax/VC-MD\n";
		print OUT "  #      4: Variable shape fixed V (not implemented)  7: fixed shape variable volume (not implemented)\n";
		if($Function =~ /vc/i) {
			print OUT "  ISIF = 3\n";
		}
		else {
			print OUT "  ISIF = 0\n";
		}
		print OUT "  #NSW: Number of ion steps.Should be 0 for phonon.\n";
		print OUT "    NSW = 100\n";
		print OUT "  # Invoke Langevin dynamics (MDALGO=3)\n";
		print OUT "  #MDALGO=3\n";
		print OUT "  # Friction coefficients for atomic degrees of freedom, one for each species defiend in POSCAR\n";
		print OUT "  #    (in ps^-1, typically 0(no thermostating) - 100 ps^-1)\n";
		print OUT "  #    See e.g. book of Allen and Tildesley: Computer simulations of liquids\n";
		print OUT "  #LANGEVIN_GAMMA=";
		for(my $i = 0 ; $i < $nAtomType ; $i++) {
			print OUT "100 ";
		}
		print OUT "\n";
		print OUT "  # Friction coefficients for lattice degrees of freedom\n";
		print OUT "  #LANGEVIN_GAMMA_L=100\n";
		print OUT "  # Fictitious mass for lattice degrees of freedom (in amu) used\n";
		print OUT "  #   in Parrinello-Rahman dynamics (Default 1000, typically 1-10?)\n";
		print OUT "  #PMASS=10\n";
		print OUT "  #IBRION: -1:Ion position fiexed 0:MD\n";
		print OUT "            1: Quasi-Newton 2:Conjugate Gradient\n";
		print OUT "            3: Damped MD, use dumping factor(SMASS/POTIM)\n";
		print OUT "            5,6: Hessian matrix calculated by finite-difference. 6 considers symmetry\n";
		print OUT "            7,8: Hessian matrix calculated analytically. 8 considers symmetry\n";
		print OUT "  IBRION = 2\n";
		print OUT "  # Nose mass definition: -3: a micro canonical ensemble (constant energy molecular dynamics)\n";
		print OUT "  #     -2: the initial velocities are kept constant  -1: the velocities are scaled each NBLOCK  step\n";
		print OUT "  #    >=0: a canonical ensemble by Nose algorism\n";
		print OUT "  SMASS = -3\n";
		print OUT "  #NBLOCK = 50\n";
		print OUT "  POTIM = 0.5 fs  time-step for ion-motion\n";
		print OUT "  TEBEG = 0000.0\n";
		print OUT "  TEEND = 0000.0\n";
		print OUT "  PSTRESS = 0.0 kBar\n";
		print OUT "  #POMASS and ZVAL are read from POTCAR for default\n";
		print OUT "  #POMASS = 102.91\n";
		print OUT "  #ZVAL = 11.0\n";
		print OUT "  #ISYM: 0:Break symmetry 1,2:Keep symmetry\n";
		print OUT "         2:More effecient memory use\n";
		print OUT "  ISYM  = 0\n";
		print OUT "  #SYMPREC: deault 1e-5\n";
		print OUT "  #SYMPREC = 1e-5\n";
	}
	elsif($Function =~ /md/i) {
		print OUT "  #ISIF: 0:Ion relax only 2:Cell fixed/Calc.Pressure 3:VC-Relax/VC-MD\n";
		print OUT "  #      4: Variable shape fixed V (not implemented)  7: fixed shape variable volume (not implemented)\n";
		if($Function =~ /vc/i) {
			print OUT "  ISIF = 3\n";
		}
		else {
			print OUT "  ISIF = 0\n";
		}
		print OUT "  #NSW: Number of ion steps\n";
		print OUT "    NSW = 100\n";
		print OUT "  # Invoke Langevin dynamics (MDALGO=3)\n";
		if($Function =~ /vc/i) {
			print OUT "  MDALGO=3\n";
		}
		else {
			print OUT "  #MDALGO=3\n";
		}
		print OUT "  # Friction coefficients for atomic degrees of freedom, one for each species defiend in POSCAR\n";
		print OUT "  #    (in ps^-1, typically 0(no thermostating) - 100 ps^-1)\n";
		print OUT "  #    See e.g. book of Allen and Tildesley: Computer simulations of liquids\n";
		if($Function =~ /vc/i) {
			print OUT "  #LANGEVIN_GAMMA=";
			for(my $i = 0 ; $i < $nAtomType ; $i++) {
				print OUT "100 ";
			}
			print OUT "\n";
		print OUT "  # Friction coefficients for lattice degrees of freedom\n";
			print OUT "  LANGEVIN_GAMMA_L=100\n";
		}
		else {
			print OUT "  #LANGEVIN_GAMMA=";
			for(my $i = 0 ; $i < $nAtomType ; $i++) {
				print OUT "100 ";
			}
			print OUT "\n";
		print OUT "  # Friction coefficients for lattice degrees of freedom\n";
			print OUT "  #LANGEVIN_GAMMA_L=100\n";
		}
		print OUT "  # Fictitious mass for lattice degrees of freedom (in amu) used\n";
		print OUT "  #   in Parrinello-Rahman dynamics (Default 1000, typically 1-10?)\n";
		if($Function =~ /vc/i) {
			print OUT "  PMASS=10\n";
		}
		else {
			print OUT "  #PMASS=10\n";
		}
		print OUT "  #IBRION: -1:Ion position fiexed 0:MD\n";
		print OUT "            1: Quasi-Newton 2:Conjugate Gradient\n";
		print OUT "            3: Damped MD, use dumping factor(SMASS/POTIM)\n";
		print OUT "            5,6: Hessian matrix calculated by finite-difference. 6 considers symmetry\n";
		print OUT "            7,8: Hessian matrix calculated analytically. 8 considers symmetry\n";
		print OUT "  IBRION = 0\n";
		print OUT "  # Nose mass definition: -3: a micro canonical ensemble (constant energy molecular dynamics)\n";
		print OUT "  #     -2: the initial velocities are kept constant  -1: the velocities are scaled each NBLOCK  step\n";
		print OUT "  #    >=0: a canonical ensemble by Nose algorism\n";
		print OUT "  SMASS = -3\n";
		print OUT "  #NBLOCK = 50\n";
		print OUT "  POTIM = 2.0 fs time-step for ion-motion\n";
		print OUT "  TEBEG = 300.0\n";
		print OUT "  TEEND = 300.0\n";
		print OUT "  PSTRESS = 0.0 kBar\n";
		print OUT "  #POMASS and ZVAL are read from POTCAR for default\n";
		print OUT "  #POMASS = 102.91\n";
		print OUT "  #ZVAL = 11.0\n";
		print OUT "  #ISYM: 0:Break symmetry 1,2:Keep symmetry\n";
		print OUT "         2:More effecient memory use\n";
		print OUT "  ISYM = 0\n";
		print OUT "  #SYMPREC: deault 1e-5\n";
		print OUT "  #SYMPREC = 1e-5\n";
	}
	else {
		print OUT "  #ISIF: 0:Ion relax only 2:Cell fixed/Calc.Pressure 3:VC-Relax/VC-MD\n";
		print OUT "  #      4: Variable shape fixed V (not implemented)  7: fixed shape variable volume (not implemented)\n";
		print OUT "  #ISIF = 3\n";
		print OUT "  #NSW: Number of ion steps\n";
		print OUT "  #NSW = 100\n";
		print OUT "  # Invoke Langevin dynamics (MDALGO=3)\n";
		print OUT "  #MDALGO=3\n";
		print OUT "  # Friction coefficients for atomic degrees of freedom, one for each species defiend in POSCAR\n";
		print OUT "  #    (in ps^-1, typically 0(no thermostating) - 100 ps^-1)\n";
		print OUT "  #    See e.g. book of Allen and Tildesley: Computer simulations of liquids\n";
		print OUT "  #LANGEVIN_GAMMA=";
		for(my $i = 0 ; $i < $nAtomType ; $i++) {
			print OUT "100 ";
		}
		print OUT "\n";
		print OUT "  # Friction coefficients for lattice degrees of freedom\n";
		print OUT "  #LANGEVIN_GAMMA_L=100\n";
		print OUT "  # Fictitious mass for lattice degrees of freedom (in amu) used\n";
		print OUT "  #   in Parrinello-Rahman dynamics (Default 1000, typically 1-10?)\n";
		print OUT "  #PMASS=10\n";
		print OUT "  #IBRION: -1:Ion position fiexed 0:MD\n";
		print OUT "            1: Quasi-Newton 2:Conjugate Gradient\n";
		print OUT "            3: Damped MD, use dumping factor(SMASS/POTIM)\n";
		print OUT "            5,6: Hessian matrix calculated by finite-difference. 6 considers symmetry\n";
		print OUT "            7,8: Hessian matrix calculated analytically. 8 considers symmetry\n";
		print OUT "  IBRION = -1\n";
		print OUT "  # Nose mass definition: -3: a micro canonical ensemble (constant energy molecular dynamics)\n";
		print OUT "  #     -2: the initial velocities are kept constant  -1: the velocities are scaled each NBLOCK  step\n";
		print OUT "  #    >=0: a canonical ensemble by Nose algorism\n";
		print OUT "  #SMASS = -3\n";
		print OUT "  #NBLOCK = 50\n";
		print OUT "  #POTIM = 2.0 fs time-step for ion-motion\n";
		print OUT "  #TEBEG = 300.0\n";
		print OUT "  #TEEND = 300.0\n";
		print OUT "  #PSTRESS = 0.0 kBar\n";
		print OUT "  #POMASS and ZVAL are read from POTCAR for default\n";
		print OUT "  #POMASS = 102.91\n";
		print OUT "  #ZVAL = 11.0\n";
		print OUT "  #ISYM: 0:Break symmetry 1,2:Keep symmetry\n";
		print OUT "         2:More effecient memory use\n";
		print OUT "    ISYM = 1\n";
		print OUT "  #SYMPREC: deault 1e-5\n";
		print OUT "  #SYMPREC = 1e-5\n";
	}

	print OUT "\n";
	print OUT "# EFERMI: LEGACY, MIDGAP, real: Position of Efermi in OUTCAR for semiconductor\n";
	print OUT "  EFERMI=MIDGAP\n";

	print OUT "\n";
	print OUT "  #LEPSILON=.TRUE.\n";
	print OUT "  #LRPA=.FALSE.\n";
	print OUT "  #LOPTICS=.TRUE.\n";
	print OUT "  #CSHIFT=0.1\n";

	print OUT "\n";
	print OUT "File write controls:\n";
	print OUT "  LWAVE  = .TRUE.\n";
	print OUT "  LCHARG = .TRUE.\n";
	print OUT "  # LVHAR, LVTOT: Write Hartree (LVHAR=.T.) or Total (LVTOT=.T.)local potential to LOCPOT in vasp.5.x\n";
	print OUT "  #LVHAR  = .TRUE.\n";
	print OUT "  LVTOT  = .TRUE.\n";
	if(defined $pParams->{LKPOINTS_OPT} and $pParams->{LKPOINTS_OPT} ne '') {
		print OUT "  LKPOINTS_OPT = $pParams->{LKPOINTS_OPT}\n";
	}
	else {
		print OUT "  #LKPOINTS_OPT = .TRUE.\n";
	}
	print OUT "  #LORBIT  = 2\n";
	print OUT "  #LELF    = .TRUE.\n";
	print OUT "  # LAECHG: Write core charge to AECCAR0, valance charge to AECCAR2 e.g. for Bader analysis\n";
	print OUT "  LAECHG  = .TRUE.\n";
	print OUT "\n";

	print OUT "Parallelization configuration:\n";
	print OUT "  #Parallelizaion over bands and plane wave coeff.(the latter is default)\n";
	print OUT "  # For dense cluster: NPAR ~ sqrt(Number of cores) or Number of cores per computer node\n";
	print OUT "  NPAR = 1\n";
	print OUT "  #NCORE = 1 # NCORE = Total number of cores / NPAR\n";
	print OUT "  #LPLANE = .TRUE.\n"; 
	print OUT "  #KPAR = 1\n";
	print OUT "\n";

	print OUT "DOS related values:\n";
	my $IsPartial = 0;
	for(my $i = 0 ; $i < @AtomSiteList ; $i++) {
		my $atom      = $AtomSiteList[$i];
		my $occupancy = $atom->Occupancy();
		if(abs($occupancy-1.0) > 1.0e-4) {
			$IsPartial = 1;
			last;
		}
	}
	if($IsPartial) {
		print OUT "  #Elements: ", join(" ", @Element), "\n";
		my %Done;
		print OUT "  VCA = ";
		my %Done;
		for(my $ia = 0 ; $ia < $nAtomType ; $ia++) {
			my $atomtype    = $AtomTypeList[$ia];
			my $typenameraw = $atomtype->Label();
			my $typename    = $atomtype->AtomNameOnly();
			my $occupancy   = $atomtype->Occupancy();
print "$ia: $typenameraw\n";
			next if($Done{$typenameraw});
			$Done{$typenameraw}++;
			printf OUT "%6.4f ", $occupancy;
		}
		print OUT "\n";
	}

	print OUT "  #Elements: ", join(" ", @Element), "\n";
	print OUT "  RWIGS = ", join(" ", @RWIGS), "\n";
	print OUT "  #ISMEAR: -1: Fermi-smearing  0: Gaussian smearing  N: Methfessel-Paxton order N\n";
	print OUT "  #        -4 tetrahedron method without Blochl corrections\n";
	print OUT "  #        -5 tetrahedron method with Blochl corrections\n";
	print OUT "  #    use -5 for semiconductors/insulators, 1/2 for metals\n";
	print OUT "  #    use -1/0/1/2 for Band\n";
	print OUT "  ISMEAR = 1\n";
	print OUT "  SIGMA = 0.2 broad. in eV\n";
	print OUT "  #DOS range and # of mesh\n";
	print OUT "  #EMIN  = -35\n";
	print OUT "  #EMAX  = 15\n";
	print OUT "  #NEDOS = 5000\n";
	print OUT "\n";

	print OUT "Band decomposed charge density:\n";
	print OUT "  #LPARD = .TRUE.\n";
	print OUT "  #Mode: 0: Total charge density incl. unoccupied bands.\n";
	print OUT "  #     -1: Usual total charge density (default).\n";
	print OUT "  #     -2: Partial charge density within energies specified by EINT.\n";
	print OUT "  #     -3: The same as -2, but the energy is measured from the Fermi energy.\n";
	print OUT "  #NBMOD = 0\n";
	print OUT "  #Energy range to be incorporated\n";
	print OUT "  #EINT = -5.0 0.0\n";
	print OUT "  #Indexes of bands to be used (see EIGENVAL or OUTCAR)\n";
	print OUT "  #IBAND = 20 21 22 23\n";
	print OUT "  #Indexes of K points to be used (see IBZKPT)\n";
	print OUT "  #KPUSE = 1 2 3 4\n";
	print OUT "  #Write band decomposed density to PARCHG.nb.* separately if .TRUE.\n";
	print OUT "  #LSEPB = .TRUE.\n";
	print OUT "  #Write K-point decomposed density to PARCHG.*.nk separately if .TRUE.\n";
	print OUT "  #LSEPK = .TRUE.\n";
	print OUT "\n";

	close(OUT);
	return 1;
}

sub MakePOSCARFile
{
	my ($this, $Crystal, $Function, $fname, $UseConventionalCell, $pParams) = @_;

print "<H1>Make POSCAR</H1>\n";
#ファイル作製開始
	unless(open(OUT,">$fname")) {
		print "Can not write to $fname.$LF$LF";
		return 0;
	}
	binmode(OUT);

	my $Sample  = $this->SampleName();
	my $SPG     = $Crystal->GetCSpaceGroup();
	my $SPGName = $SPG->SPGName();
#print "SPG in VASP.pm:[$this->{ConventionalSPGName}][$SPGName]\n";
#exit;

	print OUT "$Sample\n";
	my ($a, $b, $c, $alpha, $beta, $gamma);
#print "C[$Crystal->{ConventionalSPGName}]n";
	if(defined $Crystal->{ConventionalSPGName} and $Crystal->{ConventionalSPGName} =~ /^\s*[FIABC]/i) {
		my $pl = $Crystal->{ConventionalLatticeParameters};
#print "pl[$pl]\n";
		($a, $b, $c, $alpha, $beta, $gamma) = @$pl;
	}
	else {
		($a,$b,$c,$alpha,$beta,$gamma) = $Crystal->LatticeParametersByOutputMode(0);
	}
print "***L [Conv cell SPG: $Crystal->{ConventionalSPGName}]: $a, $b, $c, $alpha, $beta, $gamma\n";

	printf OUT "   %12.6f\n", $a;
	my ($ax,$ay,$az,$bx,$by,$bz,$cx,$cy,$cz) = $Crystal->LatticeVectorsByOutputMode(0);
	printf OUT "   %10.6f %10.6f %10.6f\n", $ax/$a, $ay/$a, $az/$a;
	printf OUT "   %10.6f %10.6f %10.6f\n", $bx/$a, $by/$a, $bz/$a;
	printf OUT "   %10.6f %10.6f %10.6f\n", $cx/$a, $cy/$a, $cz/$a;

	my @AtomTypeList = $Crystal->GetCAtomTypeList();
	my $nAtomType = @AtomTypeList;
	my @ExpandedAtomSiteList = $Crystal->GetCExpandedAtomSiteListByOutputMode();
	my $nTotalExpandedAtomSite = @ExpandedAtomSiteList;
#for(my $i = 0 ; $i < @AtomTypeList ; $i++) {
#	print "MakePOSCARFile: atom[$i] = ", $AtomTypeList[$i]->AtomNameOnly(), "$LF";
#}

	for(my $i = 0 ; $i < $nAtomType ; $i++) {
		my $atom = $AtomTypeList[$i];
		my $name = $atom->AtomNameOnly();
		printf OUT " %4s", $name;
	}
	print OUT "\n";

	my %nAtomTypes;
	for(my $i = 0 ; $i < @ExpandedAtomSiteList ; $i++) {
		my $atom = $ExpandedAtomSiteList[$i];
		my $name = $atom->AtomName();
#		my $name = $atom->AtomNameOnly();
		$nAtomTypes{$name}++;
	}
	my %Done;
	for(my $i = 0 ; $i < $nAtomType ; $i++) {
		my $atom = $AtomTypeList[$i];
		my $name = $atom->Label();
#		my $name = $atom->AtomNameOnly();
		next if($Done{$name});
		$Done{$name}++;
		printf OUT " %4d", $nAtomTypes{$name};
	}
	print OUT "\n";
	
#	print OUT "Selective dynamics # Necessary to specify T or F for coordinates, but removed for phonon\n";
	print OUT "Direct\n";
	my $iCount = 0;
	%Done = {};
	for(my $ia = 0 ; $ia < $nAtomType ; $ia++) {
		my $atomtype = $AtomTypeList[$ia];
#		my $typename = $atomtype->AtomName();
		my $typename = $atomtype->Label();
#		my $typename = $atomtype->AtomNameOnly();
		next if($Done{$typename});
		$Done{$typename}++;
#print "ia=$ia atom=$typename$LF";
		for(my $i = 0 ; $i < @ExpandedAtomSiteList ; $i++) {
			my $atom      = $ExpandedAtomSiteList[$i];
			my $atomname  = $atom->AtomName();
#print "name: $atomname [type $typename]\n";
#			my $atomname  = $atom->AtomNameOnly();
#print "  i=$i atom=$atomname$LF";
			next if(uc $typename ne uc $atomname);
#print "passed\n";
			my $charge    = $atom->Charge();
			my ($x,$y,$z) = $atom->Position(1);
			my $occupancy = $atom->Occupancy();

#Occupancyが1のサイト
#			my ($xc,$yc,$zc) = $Crystal->FractionalToCartesian($x,$y,$z);
#			printf "%-2s %14.9f %14.9f %14.9f T T T\n", $atomname, $x, $y, $z;
			if($occupancy >= 0.9999) {
				if($iCount == 0) {
					printf OUT "   %14.9f %14.9f %14.9f T T T\n", $x, $y, $z;
				}
				else {
					printf OUT "   %14.9f %14.9f %14.9f T T T\n", $x, $y, $z;
				}
			}
			else {
				my $i1 = $i+1;
				if($iCount == 0) {
					printf OUT "   %14.9f %14.9f %14.9f T T T\n", $x, $y, $z;
#					printf OUT "PAtom%d   %14.9f %14.9f %14.9f T T T\n", $i1, $x, $y, $z;
				}
				else {
					printf OUT "   %14.9f %14.9f %14.9f T T T\n", $x, $y, $z;
#					printf OUT "PAtom%d   %14.9f %14.9f %14.9f T T T\n", $i1, $x, $y, $z;
				}
			}
			$iCount++;
		}
	}
	close(OUT);

	return 1;
}

sub SaveVASPFiles
{
	my ($this, $Crystal, $CellType, $UseConventionalCell, $Function, $Functional, 
		$UseDPP, $SpinPolarized,
		$Dir, $IsChooseRandomly, $pParams) = @_;
	my $LF = "<br>\n";

	my $INCAR   = Deps::MakePath($Dir, "INCAR", 0);
	if(unlink($INCAR)) {
		print "$INCAR was deleted.$LF";
	}
	my $POSCAR  = Deps::MakePath($Dir, "POSCAR", 0);
	if(unlink($POSCAR)) {
		print "$POSCAR was deleted.$LF";
	}
	my $POTCAR  = Deps::MakePath($Dir, "POTCAR", 0);
	if(unlink($POTCAR)) {
		print "$POTCAR was deleted.$LF";
	}
	my $KPOINTS = Deps::MakePath($Dir, "KPOINTS", 0);
	if(unlink($KPOINTS)) {
		print "$KPOINTS was deleted.$LF";
	}

	my $ret = $this->MakeINCARFile($Crystal, $Function, $Functional, $UseDPP, $SpinPolarized, $INCAR, $pParams);
	return () if($ret <= 0);
	$ret = $this->MakePOSCARFile($Crystal, $Function, $POSCAR, $UseConventionalCell, $pParams);
	return ($INCAR) if($ret <= 0);
	$ret = $this->MakePOTCARFile($Crystal, $Functional, $UseDPP, $POTCAR, $pParams);
	return ($INCAR, $POSCAR) if($ret <= 0);
	my $KListFile = $this->MakeKPOINTSFile($Crystal, $CellType, $Function, $KPOINTS, $pParams);
#print "KListFile: $KListFile<br>\n";
	return ($INCAR, $POSCAR, $POTCAR) if(!defined $KListFile);

	return ($KListFile, $INCAR, $POSCAR, $POTCAR, $KPOINTS);
}

sub ReadVASPResults
{
	my ($this, $RootDir, $IsPrint) = @_;

	my ($drive, $directory, $filename, $ext1, $lastdir, $filebody) = Deps::SplitFilePath($RootDir);

	my $VCRelaxDir = Deps::MakePath($RootDir, "VCRelax", 0);
	my $SCFDir     = Deps::MakePath($RootDir, "SCF", 0);
	my $DOSDir     = Deps::MakePath($RootDir, "DOS", 0);

	my $INCARHash          = $this->ReadINCARtoHash(Utils::MakePath($SCFDir, 'INCAR', '/', 0));
	my $KPOINTSHash        = $this->ReadKPOINTStoHash(Utils::MakePath($SCFDir, 'KPOINTS', '/', 0));
	my $VCRelaxKPOINTSHash = $this->ReadKPOINTStoHash(Utils::MakePath($VCRelaxDir, 'KPOINTS', '/', 0));
	my $POTCARHash         = $this->ReadPOTCARtoHash(Utils::MakePath($SCFDir, 'POTCAR', '/', 0));
	my $POSCARHash         = $this->ReadPOSCARtoHash(Utils::MakePath($SCFDir, 'POSCAR', '/', 0));
	my $SCFOUTCARHash      = $this->ReadOUTCARtoHash(Utils::MakePath($SCFDir, 'OUTCAR', '/', 0));
	my $VCROUTCARHash      = $this->ReadOUTCARtoHash(Utils::MakePath($VCRelaxDir, 'OUTCAR', '/', 0));
	my $DOSCARPath         = $this->ReadDOSCAR(Utils::MakePath($DOSDir, 'DOSCAR', '/', 0));
	my $pDataArray = $this->{'DataArray'};
	my $pDOSupArray = $pDataArray->GetGraphData(0);
	my $pNeupArray  = $pDataArray->GetGraphData(1);
	my $pDOSdnArray = $pDataArray->GetGraphData(2);
	my $pNednArray  = $pDataArray->GetGraphData(3);
#print("this->{'DataArray'}: [$pDataArray] [$pDOSupArray] [$pNeupArray] [$pDOSdnArray] [$pNednArray]\n");
#foreach my $key (%$pDOSupArray) {
#	print("key: $key\n");
#}
	my $nData = $pDOSupArray->nData();
	my $title = $pDOSupArray->Title();
	my $pEDOS = $pDOSupArray->GetDataArray("x0");
#	last if(not defined $pX0);
	$nData = @$pEDOS;
	my $xlabel = $pDOSupArray->{"x0_Name"};
#print("Title: $title\n");
#print("    X: $xlabel\n");
#print("nData in DOS(up): $nData\n");
	for(my $id = 0 ; ; $id++) {
		my $pY0   = $pDOSupArray->GetDataArray("y$id");
		last if(not defined $pY0);
		my $label0 = $pDOSupArray->{"y${id}_Name"};
#print("    Y: $label0\n");
	}
#	for(my $i = 0 ; $i < $nData ; $i++) {
#		my $x = $pX0->[$i];
#		$out1->print("$x");
#		for(my $id = 0 ; ; $id++) {
#			my $pY0 = $Data0->GetDataArray("y$id");
#			last if(!defined $pY0);
#			my $y = $pY0->[$i];
#			$out1->print(",$y");
#		}
#		$out1->print("\n");
#	}

	my $Status = '';
	if(-d $VCRelaxDir and $VCROUTCARHash->{RequiredAccuracyMessage} =~ /^\s*$/) {
		$Status = 'Relaxation not converged(1)';
	}
	elsif(-d $VCRelaxDir and $VCROUTCARHash->{AbortMessage} !~ /is\s+reached/) {
		$Status = 'Relaxation not converged(2)';
	}
	elsif(-d $VCRelaxDir and $VCROUTCARHash->{TotCPUTime} eq '') {
		$Status = 'VCR aborted(1)';
	}
	elsif($SCFOUTCARHash->{AbortMessage} =~ /^\s*$/) {
		$Status = 'SCF not converged(1)';
	}
	elsif($SCFOUTCARHash->{TotCPUTime} eq '') {
		$Status = 'SCF aborted(1)';
	}
	if($Status ne '') {
		print("Error in ReadDOSCSV: Abort due to [$Status].\n");
		return ('Error', $Status);
	}
	
	my ($EF, $VBM, $CBM, $Eg, $pEDOS, $pDOS) = $this->ReadEgFromDOSOUTCAR(Deps::MakePath($DOSDir, 'OUTCAR', 0));

	print("ISPIN=$INCARHash->{ISPIN}\n");
	print("PREC=$INCARHash->{PREC}\n");
	print("EDIFF=$INCARHash->{EDIFF}\n");
	print("EDIFFG=$INCARHash->{EDIFFG}\n");
	print("LDAU=$INCARHash->{LDAU}\n");
	print("LDAUTYPE=$INCARHash->{LDAUTYPE}\n");
	print("LDAUU=$INCARHash->{LDAUU}\n");
	print("LDAUJ=$INCARHash->{LDAUJ}\n");
	print("ISMEAR=$INCARHash->{ISMEAR}\n");
	print("SIGMA=$INCARHash->{SIGMA}\n");
	$POTCARHash->{PseudoPotential} =~ s/ .*$//s;
	print("PseudoPotential=$POTCARHash->{PseudoPotential}\n");

	$VCRelaxKPOINTSHash->{KPOINTS} =~ s/^[^\n]*\n[^\n]*\n[^\n]*\n//m;
	$VCRelaxKPOINTSHash->{KPOINTS} =~ s/\n.*$//m;
	Utils::DelSpace($VCRelaxKPOINTSHash->{KPOINTS});
	my @nKVCRelax = Utils::Split("\\s+", $VCRelaxKPOINTSHash->{KPOINTS});

	$KPOINTSHash->{KPOINTS} =~ s/^[^\n]*\n[^\n]*\n[^\n]*\n//m;
	$KPOINTSHash->{KPOINTS} =~ s/\n.*$//m;
	Utils::DelSpace($KPOINTSHash->{KPOINTS});
	my @nK = Utils::Split("\\s+", $KPOINTSHash->{KPOINTS});

	print("KPOINTS=$KPOINTSHash->{KPOINTS}\n");
	print("ISPIN=$INCARHash->{ISPIN}\n");
	print("Etot=$SCFOUTCARHash->{TOTEN}\n");
	print("EF(SCF)=$SCFOUTCARHash->{EF}\n");
	print("EF(DOS)=$EF\n");
	print("EVBM(DOS)=$VBM\n");
	print("ECBM(DOS)=$CBM\n");
	print("Eg(DOS)=$Eg\n");
	print("SumChemicalComposition=$POSCARHash->{SumChemicalComposition}\n");
	print("ChemicalComposition=$POSCARHash->{ChemicalComposition}\n");
	print("FormulaUnit=$POSCARHash->{FormulaUnit}\n");

	$POSCARHash->{Latt} =~ s/\n/,/g;
	my @latt = Utils::Split(",+", $POSCARHash->{Latt});

	my @aKVCRelax = ($latt[0]*$nKVCRelax[0], $latt[1]*$nKVCRelax[1], $latt[2]*$nKVCRelax[2]);
	my ($aKminVCRelax, $aKmaxVCRelax) = Utils::CalMinMax(\@aKVCRelax);
	$aKminVCRelax = sprintf("%4.1f", $aKminVCRelax);

	my @aK = ($latt[0]*$nK[0], $latt[1]*$nK[1], $latt[2]*$nK[2]);
	my ($aKmin, $aKmax) = Utils::CalMinMax(\@aK);
	$aKmin = sprintf("%4.1f", $aKmin);

	print("Latt=$POSCARHash->{Latt}\n");
	print("aKMax=$aKmin-$aKmax\n");
	print("Volume=$POSCARHash->{Volume}\n");
	print("Density=$POSCARHash->{Density}\n");

	my $R = new ChemicalReaction;

	my $Ref = Deps::MakePath(cwd(), $RootDir, 0);
	$Ref =~ s/$ENV{HOME}//;
	my $Eavrg = $SCFOUTCARHash->{TOTEN} / $POSCARHash->{FormulaUnit};
	my $Compound = $POSCARHash->{ChemicalComposition};
#	my $Compuond = $POSCARHash->{SumChemicalComposition};
	my $SortedComposition = $R->SortChemicalFormula($POSCARHash->{ChemicalComposition});
	$Compound =~ s/\s//g;

	$Status .= " / SCFWarning:$SCFOUTCARHash->{Warning}" if($SCFOUTCARHash->{Warning} ne '');
	$Status .= " / VCRWarning:$VCROUTCARHash->{Warning}" if($VCROUTCARHash->{Warning} ne '');

	if($Status ne '') {
		print("\n\nWarning in ReadDOSCSV: [$Status].\n\n");
	}

#	$out->print("$filebody,$Compound,$SortedComposition,$POSCARHash->{SumChemicalComposition},,"
#			   ."$POSCARHash->{FormulaUnit},$SCFOUTCARHash->{TOTEN},$Eavrg,$SCFOUTCARHash->{EF},"
#			   ."$EF,$VBM,$CBM,$Eg,"
#			   ."$POSCARHash->{Density},$latt[0],$latt[1],$latt[2],$POSCARHash->{Volume},"
#			   ."$POTCARHash->{PseudoPotential},"
#			   ."$INCARHash->{PREC},$aKmin,$aKminVCRelax,$INCARHash->{EDIFF},$INCARHash->{EDIFFG},"
#			   ."$INCARHash->{NBANDS},"
#			   ."$INCARHash->{ISPIN},$INCARHash->{NUPDOWN},$INCARHash->{LDAU},$INCARHash->{LDAUU},$INCARHash->{LDAUJ},"
#			   ."$INCARHash->{ISMEAR},$INCARHash->{SIGMA},"
#			   ."$Status,$Ref\n");

	return  ($EF, $VBM, $CBM, $Eg, $pDOSupArray, $pDOSdnArray,
			$INCARHash, $KPOINTSHash, $VCRelaxKPOINTSHash, $POTCARHash, $POSCARHash, $SCFOUTCARHash, $VCROUTCARHash);
}

sub AnalyzeaKProduct
{
	my ($this, $Crystal, $aKProduct) = @_;
	if(!$Crystal) {
		print "Error in VASP::AnalyzeaKProduct: Can not get Crystal object\n";
		exit;
	}

	my ($a,$b,$c,$alpha,$beta,$gamma) 
		= $Crystal->LatticeParametersByOutputMode(0);

	my ($nx, $ny, $nz);
	my ($NKREDX, $NKREDY, $NKREDZ) = ();
	if($aKProduct =~ /^fix:(.*)$/i) {
		my ($a1, $a2, $a3) = Utils::Split(",", $1);
		($nx, $NKREDX) = Utils::Split("[/_]", $a1);
		($ny, $NKREDY) = Utils::Split("[/_]", $a2);
		($nz, $NKREDZ) = Utils::Split("[/_]", $a3);
	}
	else {
		my $an = $aKProduct * 10.0;
		$nx = int($an / $a) + 1;
		$ny = int($an / $b) + 1;
		$nz = int($an / $c) + 1;
	}
	return ($nx, $ny, $nz, $NKREDX, $NKREDY, $NKREDZ);
}

sub ReadIBZKPTToArray
{
	my ($this, $IBZKPTFile) = @_;

	my $in = JFile->new($IBZKPTFile, 'r') or die "$!: Can not read [$IBZKPTFile].\n";
	$in->ReadLine();
	my $nK = $in->ReadLine() + 0;
	$in->ReadLine();
	my @k;
	for(my $i = 0 ; $i < $nK ; $i++) {
		my $line = $in->ReadLine();
		last if(!defined $line);

		my ($kx, $ky, $kz, $w) = Utils::Split("\\s+", $line);
		$k[$i] = {
			kx => $kx,
			ky => $ky,
			kz => $kz,
			w  => $w,
			};
	}
	$in->Close();

	return ($nK, \@k);
}

sub ReadIBZKPT
{
	my ($this, $IBZKPTFile) = @_;

	open(IN, $IBZKPTFile) or die "$!: Can not read [$IBZKPTFile].\n";
	my (@lines, @SCFKList);
	my $c = 0;
	$lines[$c++] = <IN>; #Automatically generated mesh
	my $nK = $lines[$c++] = <IN>; #       6
	$lines[$c++] = <IN>; #Reciprocal lattice
	for(my $i = 0 ; $i < $nK ; $i++) {
		my $line = <IN>;
		
		$lines[$c++] = $line;
		chomp($line);
		$SCFKList[$i] = $line;
	}
	close(IN);
	return ($nK, \@lines, \@SCFKList);
}

sub ReadKListFromklist
{
	my ($this, $klistFile, $nKPointDiv) = @_;

	my $IsLineMode = 0;
	my @klists;
	my @lines;
	my $c = 0;
	my $in = JFile->new($klistFile, "r") or die "$!: Can not read [$klistFile].\n";
#k-points along high symmetry linesA
	my $line = $in->ReadLine();
#15
	$line = $in->ReadLine();
#Line-mode
	$line = $in->ReadLine();
	if($line =~ /Line-mode/i) {
		$IsLineMode = 1;
#Reciprocal
		$line = $in->ReadLine();
		my ($kx0, $ky0, $kz0);
		while(1) {
			my $line = $in->ReadLine();
			last if(!defined $line);
			Utils::DelSpace($line);
			next if($line eq '');

#print "l:[$line";
			my ($kx, $ky, $kz, $comment, $name) = Utils::Split("\\s+", $line);
			next if(!defined $kz);
			next if(abs($kx-$kx0) < 1.0e-4 and abs($ky-$ky0) < 1.0e-4 and abs($kz-$kz0) < 1.0e-4);

#print "GLH: ", $GLabelHash{a}, "\n";
			$name = '' if(!defined $name);
#print "l:[$name, $kx, $ky, $kz]\n";

			push(@klists, sprintf("%10.8f  %10.8f  %10.8f    0.0", $kx, $ky, $kz));
			$c++;
			$kx0 = $kx;
			$ky0 = $ky;
			$kz0 = $kz;
		}
	}
	else {
		$in->rewind();
		while(1) {
			my $line = $in->ReadLine();
			last if(!defined $line);

			$lines[$c++] = $line;
			Utils::DelSpace($line);
			last if($line =~ /^END/i);
			$line =~ s/^\s*[a-zA-Z]+\s+//;
			my ($kx, $ky, $kz, $ndiv) = Utils::Split("\\s+", $line);
			if($ndiv ne '') {
				push(@klists, sprintf("%10.8f  %10.8f  %10.8f    0.0", $kx/$ndiv, $ky/$ndiv, $kz/$ndiv));
			}
		}
	}
	close(IN);

	if($IsLineMode) {
		my @a = @klists;
		@klists = ();
		my $nd = int($nKPointDiv / @a);
		$nd = 1 if($nd == 0);
		for(my $i = 0 ; $i < @a - 1 ; $i++) {
			my ($x0, $y0, $z0) = Utils::Split("\\s+", $a[$i]);
			my ($x1, $y1, $z1) = Utils::Split("\\s+", $a[$i+1]);
			for(my $j = 0 ; $j < $nd ; $j++) {
				my $kx = $x0 + ($x1 - $x0) / $nd * $j;
				my $ky = $y0 + ($y1 - $y0) / $nd * $j;
				my $kz = $z0 + ($z1 - $z0) / $nd * $j;
				push(@klists, sprintf("%10.8f  %10.8f  %10.8f    0.0", $kx, $ky, $kz));
			}
		}
		push(@klists, $a[@a-1]);
	}
	my $nKPoints = @klists;

	return ($nKPoints, \@lines, \@klists);
}

sub ReadKListFromKPOINTS
{
	my ($this, $klistFile, $SkipNonZerow) = @_;
	$SkipNonZerow = 0 if(!defined $SkipNonZerow);

	my @klists;
	my @lines;
	my $c = 0;

	open(IN, $klistFile) or die "$!: Can not read [$klistFile].\n";
	my $line0 = $lines[$c++] = <IN>; #kpoints along high symmetry linesA
	my $nK    = $lines[$c++] = <IN>; #       11
	$nK += 0;
	my $c = 0;
	my $mode;
	if($line0 =~ /^Explicit/i) {
		$mode = $line0;
	}
	else {
		$mode = $lines[$c++] = <IN>; #       Line-mode or Reciprocal
	}
	my $Space;
	my $iK0 = 0;

print "mode: $mode\n";
	if($mode =~ /^Line-mode/i) {
		$Space = $lines[$c++] = <IN>; # Reciprocal
		while(1) {
			my $line1 = <IN>;
#print "line1(1): $line1\n";
			while(1) {
				last if(eof(IN));
				Utils::DelSpace($line1);
				last if($line1 ne '');
				$line1 = <IN>;
			}
			last if(eof(IN));
			my $line2 = <IN>;
			last if(eof(IN));

			my ($kx0, $ky0, $kz0) = Utils::Split("\\s+", $line1);
			my ($kx1, $ky1, $kz1) = Utils::Split("\\s+", $line2);

			for(my $i = 0 ; $i < $nK ; $i++) {
				my $f = 1.0 * $i / ($nK-1);
				my $kx = $kx0 + ($kx1 - $kx0) * $f;
				my $ky = $ky0 + ($ky1 - $ky0) * $f;
				my $kz = $kz0 + ($kz1 - $kz0) * $f;
				$lines[$c]  = sprintf("%10.8f  %10.8f  %10.8f    0.0", $kx, $ky, $kz);
#				$klists[$c] = [$kx, $ky, $kz, 0];
				push(@klists, sprintf("%10.8f  %10.8f  %10.8f    0.0", $kx, $ky, $kz));
#print "k[$c]=($kx,$ky,$kz)\n";
				$c++;
			}
		}
	}
	elsif($mode =~ /^Explicit/i) {
		$Space = $lines[$c++] = <IN>; # Reciprocal
		for(my $i = 0 ; $i < $nK ; $i++) {
			my $line1 = <IN>;
			last if(!defined $line1);
#			last if(eof(IN));

			my ($kx, $ky, $kz, $w) = Utils::Split("\\s+", $line1);
			$lines[$c]  = sprintf("%10.8f  %10.8f  %10.8f    %10.8f", $kx, $ky, $kz, $w);
#			$klists[$c] = [$kx, $ky, $kz, 0];
			push(@klists, sprintf("%10.8f  %10.8f  %10.8f    0.0", $kx, $ky, $kz));
			$c++;
			$iK0 = $i if($iK0 == 0 and $w == 0.0);
		}
	}
	else {
		$Space = $mode;
		$mode = '';
		while(1) {
			my $line1 = <IN>;
			last if(eof(IN));
			my ($kx, $ky, $kz, $w) = Utils::Split("\\s+", $line1);
			next if($kz eq '');
			$lines[$c]  = sprintf("%10.8f  %10.8f  %10.8f    %10.8f", $kx, $ky, $kz, $w);
#			$klists[$c] = [$kx, $ky, $kz, 0];
			push(@klists, sprintf("%10.8f  %10.8f  %10.8f    0.0", $kx, $ky, $kz));
			$c++;
		}
	}

	close(IN);
#	my $nKPoints = @lines;
	my $nKPoints = @klists;
#print "nKPoints=$nKPoints/$iK0\n";
#exit;

	return ($nKPoints, \@lines, \@klists, $iK0);
}

sub MakeKPOINTSForFS
{
	my ($this, $KListPath, $CheckPath, $Crystal, $Title, $UseSymmetry, $nx, $ny, $nz, $K, $IsPrint) = @_;
	$UseSymmetry = 0   if(!defined $UseSymmetry);
	$nx          = 5   if(!defined $nx);
	$ny          = 5   if(!defined $ny);
	$nz          = 5   if(!defined $nz);
	$K           = 1.0 if(!defined $K);
	$IsPrint     = 1   if(!defined $IsPrint);

	$Crystal->SetLatticeSystem(undef, 1, 1.0e-3, 1.0e-2);

	my $check;
if($CheckPath) {
	$check = new JFile();
	if(!$check->Open($CheckPath, "w")) {
		print("$!: Can not write to [$CheckPath].\n") if($IsPrint);
	}
}

	my @KIndex;
	my $count = 0;
	my $BandName = "";
	my @lines;
print("kx range: (0 - $nx)\n") if($IsPrint);
	for(my $ix = 0 ; $ix <= $nx ; $ix++) {
		for(my $iy = 0 ; $iy <= $ny ; $iy++) {
			my $BandNamePrinted = 0;
			for(my $iz = 0 ; $iz <= $nz ; $iz++) {
				my ($IsConverted, $cix, $ciy, $ciz) = $Crystal->ConvertLatticeIndexByLatticeSystem($ix, $iy, $iz);
#print("VASP.pm: $ix-$iy-$iz: $IsConverted, $cix, $ciy, $ciz\n");

				if($UseSymmetry and $IsConverted) {
					next;
				}

				my $x = 1.0 * $K * $ix / $nx;
				my $y = 1.0 * $K * $iy / $ny;
				my $z = 1.0 * $K * $iz / $nz;
				my $weight = 1.0;
				if(!$BandNamePrinted) {
					$BandNamePrinted = 1;
					$lines[$count] = sprintf("%8.4f %8.4f %8.4f %8.4f ! B%02d%02d\n", $x, $y, $z, $weight, $ix+1, $iy+1);
				}
				else {
					$lines[$count] = sprintf("%8.4f %8.4f %8.4f %8.4f\n", $x, $y, $z, $weight);
				}
if($check) {
	$KIndex[$count] = sprintf("%03d%03d%03d", $ix, $iy, $iz);
	$check->printf("%05d: %2d %2d %2d [%s]\n", $count, $ix, $iy, $iz, $KIndex[$count]);
}
				$count++;
			}
		}
	}
if($check) {
	$check->Close();
}

print("Save to [$KListPath]\n") if($IsPrint);
	my $out = new JFile;
	if(!$out->Open($KListPath, "w")) {
		print("Error: $!: Can not write to [$KListPath]\n") if($IsPrint);
		return (0);
	}
	$out->print("$Title\n");
	$out->print("$count\n"); #($nx+1) * ($ny+1) * ($nz+1), "\n");
#	$out->print("Cartesian\n");
	$out->print("Reciprocal\n");
	for(my $i = 0 ; $i < $count ; $i++) {
		$out->print($lines[$i]);
	}
	$out->printf("END\n\n");
	$out->Close();
	
	return (1, @KIndex);
}

sub MakeBandKPOINTSFileFromKLISTForHF
{
	my ($this, $Function, $OutputName, $klistFile, $UseGammaOnly, 
			$nKPointDiv, $KPoints, $aKProduct, $pParams) = @_;
	$UseGammaOnly = 0 if(!defined $UseGammaOnly);
	$nKPointDiv = 11 if(!defined $nKPointDiv or $nKPointDiv <= 0);
	$KPoints = 'File' if(!defined $KPoints);
	$aKProduct = 2.0 if(!defined $aKProduct);

	my $IBZKPTFile = "SCF/IBZKPT";
	$IBZKPTFile = "../SCF/IBZKPT" if(!-f $IBZKPTFile);
print "*Read [$IBZKPTFile] for HF$LF";
	my ($nK, $pHeader, $pSCFKList) = $this->ReadIBZKPT($IBZKPTFile);;
	my @SCFKList = @$pSCFKList;

#klistファイル読み込み
	my @klists;
	if($KPoints =~ /File/i and $klistFile =~ /\.KPOINTS$/) {
#		print "KPOINTS: Copy [$klistFile] to [KPOINTS].\n";
#		Utils::CopyFile($klistFile, $OutputName);
#		return;
print "*Read [$klistFile] for HF$LF";
		my ($nKPoints, $plines, $pklists) = $this->ReadKListFromKPOINTS($klistFile);
		@klists = @$pklists;
	}
	elsif($KPoints =~ /File/i and $klistFile ne '') {
print "*Read [$klistFile] for HF$LF";
		my ($nKPoints, $plines, $pklists) = $this->ReadKListFromklist($klistFile, $nKPointDiv);
		@klists = @$pklists;
	}

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

	my $nKTotal = @SCFKList + @klists;
	print OUT "Explicit k-points\n";
	print OUT "$nKTotal\n";
	print OUT "Reciprocal\n";
	for(my $i = 0 ; $i < @SCFKList ; $i++) {
		print OUT "$SCFKList[$i]\n";
	}

	for(my $i = 0 ; $i < @klists ; $i++) {
print "klist[$i]=$klists[$i]\n";
my $i1 = $i+1;
print "klist[$i1]=$klists[$i1]\n" if($i == @klists-2);
		printf OUT "$klists[$i]\n";
	}
	close(OUT);
}

sub MakeBandKPOINTSFileForHF
{
	my ($this, $Crystal, $CellType, $Function, $fname, $UseGammaOnly, 
			$nKPointDiv, $KPoints, $aKProduct, $pParams) = @_;

	my $IBZKPTFile = "SCF/IBZKPT";
print "*Read [$IBZKPTFile]$LF";
	my ($nK, $pHeader, $pSCFKList) = $this->ReadIBZKPT($IBZKPTFile);;
	my @SCFKList = @$pSCFKList;

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

	my ($a,$b,$c,$alpha,$beta,$gamma) = $Crystal->LatticeParametersByOutputMode(0);
	my ($nx, $ny, $nz, $NKREDX, $NKREDY, $NKREDZ) = $this->AnalyzeaKProduct($Crystal, $aKProduct);

	my $XC = new XCrySDen;
	my $KListDBDir = $XC->SetKListDBDir($this->KListDBDir);
	my $SPGName = $Crystal->SPGName();
	my $iSPG    = $Crystal->iSPG();
	my ($KListFilePath, @KList) = $XC->ReadKList($iSPG, $SPGName, $CellType);
	my $nKList = @KList;

	my $nKTotal = @SCFKList + @KList;
	print OUT "Explicit k-points\n";
	print OUT "$nKTotal\n";
	print OUT "Reciprocal\n";
	for(my $i = 0 ; $i < @SCFKList ; $i++) {
		print OUT "$SCFKList[$i]\n";
	}

	for(my $i = 1 ; $i < $nKList+1 ; $i++) {
		my $pk1    = $KList[$i-1];
		my $pk2    = $KList[$i];
		my $label1 = $pk1->{'Label'};
		my $kx1    = $pk1->{'kx'};
		my $ky1    = $pk1->{'ky'};
		my $kz1    = $pk1->{'kz'};
		my $label2 = $pk2->{'Label'};
		my $kx2    = $pk2->{'kx'};
		my $ky2    = $pk2->{'ky'};
		my $kz2    = $pk2->{'kz'};
		printf OUT "%10.4f %10.4f %10.4f  0.0\n", $kx1, $ky1, $kz1;
		printf OUT "%10.4f %10.4f %10.4f  0.0\n", $kx2, $ky2, $kz2;
		print OUT "\n";
	}
}

sub MakeKPOINTSFileFromKLIST
{
	my ($this, $Function, $OutputName, $LKPOINTS_OPT, $klistFile, $UseGammaOnly, 
			$nKPointDiv, $KPoints, $aKProduct, $pParams) = @_;
	$UseGammaOnly = 0 if(!defined $UseGammaOnly);
	$nKPointDiv = 11 if(!defined $nKPointDiv or $nKPointDiv <= 0);
	$KPoints = 'File' if(!defined $KPoints);
	$aKProduct = 2.0 if(!defined $aKProduct);
	if(!defined $pParams->{METAGGA}) {
		$pParams->{METAGGA} = '';
	}

print "<H1>6113: Make KPOINTS from KLIST [$klistFile]</H1>\n";
	if($Function =~ /band/) {
#klistファイル読み込み
print "*Read [$klistFile]$LF";
		if($pParams->{HybridFunctional} ne '' or $pParams->{METAGGA} ne '') {
			print("Hybrid : [$pParams->{HybridFunctional}]\n");
			print("METAGGA: [$pParams->{METAGGA}]\n");
			print("LKPOINTS_OPT: [$LKPOINTS_OPT]\n");
			if($LKPOINTS_OPT ne '.TRUE.') {
				print("Run MakeBandKPOINTSFileFromKLISTForHF\n");
				return $this->MakeBandKPOINTSFileFromKLISTForHF(
					$Function, $OutputName, $klistFile, $UseGammaOnly, 
					$nKPointDiv, $KPoints, $aKProduct, $pParams);
			}
		}

		my @klists;
		if($KPoints =~ /File/i and $klistFile =~ /\.KPOINTS$/) {
			print "KPOINTS: Copy [$klistFile] to [KPOINTS].\n";
			Utils::CopyFile($klistFile, $OutputName);
			return;
		}
		elsif($KPoints =~ /File/i and $klistFile ne '') {
			my $in = JFile->new($klistFile, "r") or die "$!: Can not read [$klistFile].\n";
#k-points along high symmetry linesA
			my $line = $in->ReadLine();
#15
			$line = $in->ReadLine();
#Line-mode
			$line = $in->ReadLine();
			if($line =~ /Line-mode/i) {
#Reciprocal
				$line = $in->ReadLine();
				my ($kx0, $ky0, $kz0);
				while(1) {
					$line = $in->ReadLine();
					last if(!defined $line);

					my ($kx, $ky, $kz, $comm, $name) = Utils::Split("\\s+", $line);
					next if(!defined $kz);
					next if(abs($kx-$kx0) < 1.0e-3 and abs($ky-$ky0) < 1.0e-3 and abs($kz-$kz0) < 1.0e-3);

					my $s = sprintf("%-7s %f %f %f 1", $name, $kx, $ky, $kz);
					push(@klists, $s);
					$kx0 = $kx;
					$ky0 = $ky;
					$kz0 = $kz;
				}
			}
			else {
				$in->rewind();
				while(1) {
					$line = $in->ReadLine();;
					last if(!defined $line);
					chomp($line);
					my $s = substr($line, 0, 10);
					$s =~ s/^\s*//;
					$s =~ s/\s*$//;
					last if($s =~ /^END/i);
					next if($s =~ /^DELTA/i);
					next if($s =~ /^LAMBDA/i);
					next if($s =~ /^SIGMA/i);
					if($s ne '') {
						my ($xyz) = ($line =~ /^(\S+\s+\S+\s+\S+\s+\S+\s+\S+)/);
						push(@klists, $xyz);
print "   k points: $xyz$LF";
					}
				}
			}
			$in->Close();
		}
		elsif($KPoints =~ /File/i) {
			print "Can not read $klistFile.$LF";
			print "Use default klists.$LF";
			push(@klists, "R     1 1 1 2");
			push(@klists, "GAMMA 0 0 0 2");
			push(@klists, "X     1 0 0 2");
			push(@klists, "Z     0 0 1 2");
			push(@klists, "M     1 1 0 2");
			push(@klists, "GAMMA 0 0 0 2");
		}
		else {
			if(length($KPoints) == 1) {
				$KPoints = "$KPoints$KPoints";
			}
			for(my $i = 0 ; $i < length($KPoints) ; $i++) {
				my $c0 = substr($KPoints, $i, 1);
#print "$i: $c0\n";
				if($c0 =~ /G/i) {
					push(@klists, "GAMMA     0 0 0 2");
				}
				elsif($c0 =~ /X/i) {
					push(@klists, "X         1 0 0 2");
				}
				elsif($c0 =~ /Y/i) {
					push(@klists, "Y         0 1 0 2");
				}
				elsif($c0 =~ /Z/i) {
					push(@klists, "Z         0 0 1 2");
				}
				elsif($c0 =~ /R/i) {
					push(@klists, "R         1 1 1 2");
				}
				elsif($c0 =~ /M/i) {
					push(@klists, "M         1 1 0 2");
				}
			}
		}
#ファイル作製開始
print "*Create [$OutputName]$LF";
		unless(open(OUT,">$OutputName")) {
			print "Can not write to $OutputName.$LF$LF";
			return undef;
		}
		binmode(OUT);
	
		print OUT "k-points along high symmetry linesA\n";
		print OUT "$nKPointDiv\n";
		print OUT "Line-mode\n";
		print OUT "Reciprocal\n";
		for(my $i = 0 ; $i < @klists-1 ; $i++) {
print "klist[$i]=$klists[$i]\n";
my $i1 = $i+1;
print "klist[$i1]=$klists[$i1]\n" if($i == @klists-2);
			my ($name0,$x0,$y0,$z0,$m0) 
				= ($klists[$i] =~ /(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)/);
			my ($name1,$x1,$y1,$z1,$m1) 
				= ($klists[$i+1] =~ /(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)/);
			printf OUT "  %8.4f   %8.4f   %8.4f   ! %s\n", 
				$x0/$m0, $y0/$m0, $z0/$m0, $name0;
			printf OUT "  %8.4f   %8.4f   %8.4f   ! %s\n", 
				$x1/$m1, $y1/$m1, $z1/$m1, $name1;
			print OUT "\n";
		}
		close(OUT);
	}
	else {
		my ($drive, $dir, $filename, $ext, $lastdir, $filebody) 
				= Deps::SplitFilePath($OutputName);
		my $CARDir;
		if($drive) {
			$CARDir = "$drive$dir";
		}
		else {
			$CARDir = $dir;
		}
		my $Crystal = $this->ReadStructureFromCARFiles($CARDir, 0);
		my ($a,$b,$c,$alpha,$beta,$gamma) 
				= $Crystal->LatticeParametersByOutputMode(0);

		my ($nx, $ny, $nz, $NKREDX, $NKREDY, $NKREDZ) 
			= $this->AnalyzeaKProduct($Crystal, $aKProduct);
print "(nx,ny,nz)=($nx,$ny,$nz)\n";

#		if($pParams->{HybridFunctional} =~ /^HSE/ and $pParams->{NKRED} > 1) {
#if(0) {
#			$nx = int( ($nx+1) / $pParams->{NKRED}) * $pParams->{NKRED};
#			$ny = int( ($ny+1) / $pParams->{NKRED}) * $pParams->{NKRED};
#			$nz = int( ($nz+1) / $pParams->{NKRED}) * $pParams->{NKRED};
#			$nx = 4 if($nx > 4);
#			$ny = 4 if($ny > 4);
#			$nz = 4 if($nz > 4);
#}
#		}
#		if($pParams->{CalVibration}) {
#if(0) {
#			$nx = 3 if($nx > 3);
#			$ny = 3 if($ny > 3);
#			$nz = 3 if($nz > 3);
#}
#		}
		unless(open(OUT,">$OutputName")) {
			print "Can not write to $OutputName.$LF$LF";
			return undef;
		}
		binmode(OUT);

		if(!$UseGammaOnly) {
			print OUT "Automatic mesh\n";
			print OUT "0\n";
			print OUT "Gamma\n";
			print OUT "$nx  $ny  $nz\n";
			print OUT "0. 0. 0.\n";
		}
		else {
			print OUT "Automatic mesh\n";
			print OUT "0\n";
			print OUT "Gamma\n";
			print OUT "1  1  1\n";
			print OUT "0. 0. 0.\n";
#			print OUT "Automatic mesh\n";
#			print OUT "0\n";
#			print OUT "Automatic\n";
#			print OUT $nx*$ny*$nz, "\n";
		}
		
		close(OUT);
	}
}

sub MakeKPOINTSForMass
{
	my ($this, $OutputName, $klistFile, $pParams) = @_;

#klistファイル読み込み
print "*Read k points from [$klistFile]\n";
	my @klist = BandStructure->ReadKListFile($klistFile, 0);
	my $nk    = @klist;

#ファイル作製開始
print "*Create [$OutputName]$LF";
	my $out = JFile->new($OutputName, 'w');
	if(!$out) {
		print "Can not write to $OutputName.$LF$LF";
		return undef;
	}

#	if($pParams->{HybridFunctional} eq '') {
#		$out->print("k-points along high symmetry linesA\n");
#		$out->printf("%d\n", $nk);
#		$out->print("Line-mode\n");
#		$out->print("Reciprocal\n");
#		for(my $i = 0 ; $i < $nk ; $i++) {
#			my $p0 = $klist[$i];
#			$out->printf("%10.7f  %10.7f  %10.7f   0.0\n", $p0->{kx}, $p0->{ky}, $p0->{kz});
#		}
#	}
#	else {
	{
		my $IBZKPTFile = "../DFTSCF/IBZKPT";
		$IBZKPTFile = "./DFTSCF/IBZKPT" if(!-f $IBZKPTFile);
		$IBZKPTFile = "../SCF/IBZKPT"   if(!-f $IBZKPTFile);
		$IBZKPTFile = "./SCF/IBZKPT"    if(!-f $IBZKPTFile);
print "*Read IBZKPT from [$IBZKPTFile].\n";
		my ($nKIBZKPT, $pIBZKPTLines, $pIBZKPTSCFKList) = $this->ReadIBZKPT($IBZKPTFile);

		$out->print("Explicit k-points\n");
		$out->printf("%d\n", $nKIBZKPT + $nk);
		$out->print("Reciprocal\n");
		for(my $i = 0 ; $i < $nKIBZKPT ; $i++) {
			$out->print("$pIBZKPTSCFKList->[$i]\n");
		}
		for(my $i = 0 ; $i < $nk ; $i++) {
			my $p0 = $klist[$i];
			$out->printf("%10.7f  %10.7f  %10.7f   0.0\n", $p0->{kx}, $p0->{ky}, $p0->{kz});
		}
	}

	$out->Close();
}

sub MakeKPOINTSForBand
{
	my ($this, $POSCAR, $Function, $OutputName, 
		$klistFile, $nKPointDiv, $KPoints, $pParams) = @_;
#	$UseGammaOnly = 0 if(!defined $UseGammaOnly);
	$nKPointDiv   = 11 if(!defined $nKPointDiv or $nKPointDiv <= 0);
	$KPoints      = 'File' if(!defined $KPoints);
#	$aKProduct    = 2.0 if(!defined $aKProduct);

	if($Function ne 'band') {
		print("Error in VASP.pm::MakeKPOINTSForBand: Function [$Function] is not 'band'.\n");
		return;
	}

#klistファイル読み込み
	my @klist;
	if($KPoints =~ /File/i) {
print "*Read k points from [$klistFile]\n";
		@klist = BandStructure->ReadKListFile($klistFile);
	}

#ファイル作製開始
print "*Create [$OutputName]$LF";
	my $out = JFile->new($OutputName, 'w');
	if(!$out) {
		print "Can not write to $OutputName.$LF$LF";
		return undef;
	}

	if($pParams->{HybridFunctional} eq '') {
		if(@klist == 0 and $klistFile ne 'File' and !-f $klistFile) {
			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});
				}
			}
		}
		elsif(@klist == 0) {
			print "Can not read $klistFile.$LF";
			print "Use default klist.$LF";
			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});
		}

		$out->print("k-points along high symmetry linesA\n");
		$out->print("$nKPointDiv\n");
		$out->print("Line-mode\n");
		$out->print("Reciprocal\n");
		for(my $i = 0 ; $i < @klist-1 ; $i++) {
			my $p0 = $klist[$i];
			my $p1 = $klist[$i+1];
			$out->printf("  %8.4f   %8.4f   %8.4f   ! %s\n", $p0->{kx}, $p0->{ky}, $p0->{kz}, $p0->{name});
			$out->printf("  %8.4f   %8.4f   %8.4f   ! %s\n", $p1->{kx}, $p1->{ky}, $p1->{kz}, $p1->{name});
			$out->printf("\n");
		}
	}
	else {
print "*Read Structure from [$POSCAR] etc.\n";
		my $Crystal = $this->ReadStructureFromCARFiles($POSCAR, 1);
#print "Crystal=$Crystal\n";
		my ($nk, $pk) = BandStructure->DivideKList($Crystal, \@klist, $nKPointDiv);

		my $IBZKPTFile = "../DFTSCF/IBZKPT";
		$IBZKPTFile = "./DFTSCF/IBZKPT" if(!-f $IBZKPTFile);
		$IBZKPTFile = "../SCF/IBZKPT"   if(!-f $IBZKPTFile);
		$IBZKPTFile = "./SCF/IBZKPT"    if(!-f $IBZKPTFile);
print "*Read IBZKPT from [$IBZKPTFile].\n";
		my ($nKIBZKPT, $pIBZKPTLines, $pIBZKPTSCFKList) = $this->ReadIBZKPT($IBZKPTFile);

		$out->print("Explicit k-points\n");
		$out->printf("%d\n", $nKIBZKPT + $nk);
		$out->print("Reciprocal\n");
		for(my $i = 0 ; $i < $nKIBZKPT ; $i++) {
			$out->print("$pIBZKPTSCFKList->[$i]\n");
		}
		for(my $i = 0 ; $i < $nk ; $i++) {
			my $p0 = $pk->[$i];
			$out->printf("%10.7f  %10.7f  %10.7f   0.0\n", $p0->{kx}, $p0->{ky}, $p0->{kz});
		}
	}

	$out->Close();
}

sub MakeKPOINTSFile
{
	my ($this, $Crystal, $CellType, $Function, $fname, $UseGammaOnly, 
			$nKPointDiv, $KPoints, $aKProduct, $pParams) = @_;

	$aKProduct    = 2.0 if($aKProduct eq '');
	$UseGammaOnly = 0 if(!defined $UseGammaOnly);
	$nKPointDiv   = 9 if(!defined $nKPointDiv or $nKPointDiv <= 0);
	if( !defined $pParams->{METAGGA} ) {
		$pParams->{METAGGA} = '';
	}

print "<H1>Make KPOINTS [$fname]</H1>\n";
#ファイル作製開始
	unless(open(OUT,">$fname")) {
		print "Can not write to $fname.$LF$LF";
		return undef;
	}
	binmode(OUT);

	my ($a,$b,$c,$alpha,$beta,$gamma) 
		= $Crystal->LatticeParametersByOutputMode(0);

	my ($nx, $ny, $nz, $NKREDX, $NKREDY, $NKREDZ) 
			= $this->AnalyzeaKProduct($Crystal, $aKProduct);
#print "(nx,ny,nz)=($nx,$ny,$nz)\n";

	my $XC = new XCrySDen;
	my $KListDBDir = $XC->SetKListDBDir($this->KListDBDir());
print "klist DB dir: $KListDBDir$LF";
	my $SPGName = $Crystal->SPGName();
	my $iSPG    = $Crystal->iSPG();
	my ($KListFilePath, @KList) = $XC->ReadKList($iSPG, $SPGName, $CellType);
	my $nKList = @KList;
print "KList DB file: $KListFilePath$LF\n";

	if($Function =~ /(band|nscf)/i) {
		if($pParams->{HybridFunctional} ne '' or $pParams->{METAGGA} ne '') {
			return $this->MakeBandKPOINTSFileForHF($Crystal, $CellType, 
					$Function, $fname, $UseGammaOnly, $nKPointDiv,
					$KPoints, $aKProduct, $pParams);
		}

		print OUT "k-points along high symmetry linesA\n";
		print OUT "$nKPointDiv\n";
#		print OUT "10\n";
		print OUT "Line-mode\n";
		print OUT "Reciprocal\n";
#		print OUT "0.0 0.0 1.0\n";
#		print OUT "0.0 0.0 0.0\n";
#		print OUT "\n";
#		print OUT "0.0 0.0 0.0\n";
#		print OUT "0.0 1.0 1.0\n";
#		print OUT "\n";
#		print OUT "0.0 1.0 1.0\n";
#		print OUT "1.0 1.0 1.0\n";
#		print OUT "\n";
#		print OUT "1.0 1.0 1.0\n";
#		print OUT "0.0 0.0 0.0\n";
#		print OUT "\n";
		for(my $i = 1 ; $i < $nKList+1 ; $i++) {
			my $pk1 = $KList[$i-1];
			my $pk2 = $KList[$i];
			my $label1 = $pk1->{'Label'};
			my $kx1    = $pk1->{'kx'};
			my $ky1    = $pk1->{'ky'};
			my $kz1    = $pk1->{'kz'};
			my $label2 = $pk2->{'Label'};
			my $kx2    = $pk2->{'kx'};
			my $ky2    = $pk2->{'ky'};
			my $kz2    = $pk2->{'kz'};
#			my $nPoint = $pk2->{'nPoint'};
			printf OUT "%10.4f %10.4f %10.4f\n", $kx1, $ky1, $kz1;
			printf OUT "%10.4f %10.4f %10.4f\n", $kx2, $ky2, $kz2;
			print OUT "\n";
		}
	}
	elsif($Function =~ /relax/i) {
		if($nx == 1 and $ny == 1 and $nz == 1) {
			print OUT "Gamma only\n";
			print OUT "1\n";
			print OUT "Cartesian\n";
			print OUT "0.0  0.0  0.0   1.\n";
		}
		else {
			if($UseGammaOnly) {
				print OUT "Automatic mesh\n";
				print OUT "0\n";
				print OUT "Gamma\n";
				print OUT "$nx  $ny  $nz\n";
				print OUT "0. 0. 0.\n";
			}
			else {
				print OUT "Automatic mesh\n";
				print OUT "0\n";
				print OUT "Automatic\n";
				print OUT $nx*$ny*$nz, "\n";
			}
		}
	}
	else {
		if($UseGammaOnly) {
			print OUT "Automatic mesh\n";
			print OUT "0\n";
			print OUT "Gamma\n";
			print OUT "$nx  $ny  $nz\n";
			print OUT "0. 0. 0.\n";
		}
		else {
			print OUT "Automatic mesh\n";
			print OUT "0\n";
			print OUT "Automatic\n";
			print OUT $nx*$ny*$nz, "\n";
		}
	}
	close(IN);
	
	return $KListFilePath;
}

1;