#!/usr/bin/perl

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

use strict;
#use warnings;

use MyApplication;
use Utils;

use Crystal::VASP;

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

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

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

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

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

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

#==========================================
# コマンドラインオプション読み込み
#==========================================
$App->AddArgument("--Action",
		 "--Action=[MakePhonopyINCAR|MakeGOnlyKPOINTS|MakeKPOINTSforSPOSCAR|MakeMeshConf|MakePDOSConf|MakeBandConf"
		."|MakeBORN]", '');
$App->AddArgument("--PhonopyTarget", "--PhonopyTarget=[VASP|VASPDFPT] (Def:VASP)", "VASP");
$App->AddArgument("--ENCUT",         "--ENCUT=val (Def: 500)",          "500");
$App->AddArgument("--EDIFF",         "--EDIFF=val (Def: 1.0e-8)",       "1.0e-8");
$App->AddArgument("--IBRION",        "--IBRION=[-1|8] (Def:-1)",        "-1");
$App->AddArgument("--ISMEAR",        "--ISMEAR=[0|-5] (Def:-1)",        "0");
$App->AddArgument("--SIGMA",         "--SIGMA=val (Def: 0.01)",         "0.01");
$App->AddArgument("--aKProduct",     "--aKProduct=val (Def: 2.0)",      "2.0");
$App->AddArgument("--Precision",     "--Precision=val (Def: Accurate)", "Accurate");
$App->AddArgument("--nSupercell",    "--nSupercell=val (Def: 2 2 2)",   "2 2 2");
$App->AddArgument("--MeshPoints",    "--MeshPoints=val (Def: 8 8 8)",   "8 8 8");
$App->AddArgument("--BandPoints",    "--BandPoints=val (Def: )",        "");
$App->AddArgument("--fConversion",   "--fConversion=val (Def: )",       "");
$App->AddArgument("--CARDir",        "--CARDir=val (Def: .)",           ".");
exit 1 if($App->ReadArgs(1, "sjis", 0) != 1);
my $Args = $App->Args();
#my $form = new CGI;
#$Args->SetCGIForm($form);
#$Args->parseInput($WebCharCode);

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

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

#Utils::InitHTML("Research", $WebCharSet, "_self");

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

my $ret = 0;
if($Action =~ /MakePhonopyINCAR/i) {
	&MakePhonopyINCAR();
}
elsif($Action =~ /MakeGOnlyKPOINTS/i) {
	&MakeGOnlyKPOINTS();
}
elsif($Action =~ /MakeKPOINTSforSPOSCAR/i) {
	&MakeKPOINTSforSPOSCAR();
}
elsif($Action =~ /MakeMeshConf/i) {
	&MakeMeshConf();
}
elsif($Action =~ /MakePDOSConf/i) {
	&MakePDOSConf();
}
elsif($Action =~ /MakeBandConf/i) {
	&MakeBandConf();
}
elsif($Action =~ /MakeBORN/i) {
	&MakeBORN();
}
else {
	$App->print("Error: Invalid Action: [$Action]\n");
	$App->Usage();
}

#Utils::EndHTML();

exit $ret;

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

#==========================================
# &Subroutines
#==========================================
sub ReadNBANDS
{
	$App->print(VASP->new()->ReadKeyValue("NBANDS\\s*=\\s*(\\S+)", 0, $Args->GetGetArg(0)) . "\n");
	$App->print(VASP->new()->ReadKeyValue("NELECT\\s*=\\s*(\\S+)", 0, $Args->GetGetArg(0)) . "\n");
	$App->print(VASP->new()->ReadKeyValue("NUPDOWN\\s*=\\s*(\\S+)", 0, $Args->GetGetArg(0)) . "\n");
}

sub ReadAtomTypeList
{
	my ($CARDir) = @_;

	my $IsPrint = 1;
	my $vasp    = new VASP;
	my $Crystal = $vasp->ReadStructureFromCARFiles($CARDir, $IsPrint);
	my @AtomTypeList = $Crystal->GetCAtomTypeList();
#for(my $i = 0 ; $i < @AtomTypeList ; $i++) {
#	print "ReadAtomTypeList: atom[$i]: $AtomTypeList[$i]\n"
#}
	
	return @AtomTypeList;
}

sub BuildAtomTypeListLine
{
	my ($CARDir) = @_;

	my @AtomTypeList = &ReadAtomTypeList($CARDir);
	my $s = '';
	for(my $i = 0 ; $i < @AtomTypeList ; $i++) {
		my $atom = $AtomTypeList[$i];
		$s .= $atom->AtomNameOnly() . " ";
	}
	return $s;
}

sub GetStructureInfo
{
	my ($CARDir) = @_;

	my $IsPrint = 1;
	my $vasp    = new VASP;
	my $Crystal = $vasp->ReadStructureFromCARFiles($CARDir, $IsPrint);

	my @AtomTypeList = $Crystal->GetCAtomTypeList();
#for(my $i = 0 ; $i < @AtomTypeList ; $i++) {
#	print "GetStructureInfo: atom[$i]: $AtomTypeList[$i]\n"
#}
	my @ExpandedAtomSiteList = $Crystal->GetCExpandedAtomSiteListByOutputMode();

	my $nTotalExpandedAtomSite = @ExpandedAtomSiteList;
	my @AtomNames;
	for(my $i = 0 ; $i < @AtomTypeList ; $i++) {
		my $atomtype = $AtomTypeList[$i];
		$AtomNames[$i] = $atomtype->AtomNameOnly();
	}

	my %nAtoms;
	for(my $i = 0 ; $i < $nTotalExpandedAtomSite ; $i++) {
		my $atomsite = $ExpandedAtomSiteList[$i];
#		$AtomNames[$i] = $atomsite->AtomNameOnly();
		$nAtoms{$AtomNames[$i]}++;
	}
	
	my $AtomTypeListLine = join(' ', @AtomNames);
	my $SiteIndexList = "1";
	for(my $i = 2 ; $i <= @ExpandedAtomSiteList ; $i++) {
		$SiteIndexList .= ", $i";
	}
	
	return ($AtomTypeListLine, $SiteIndexList, $Crystal, \@AtomTypeList, \@ExpandedAtomSiteList);
}

sub GetMatrixConventionalToPrimitiveCell
{
	my ($BL) = @_;

#print "BL = $BL\n";
	if($BL eq 'F') {
		return [ [ 0.5, 0.5,   0], 
			 [ 0,   0.5, 0.5],
			 [ 0.5,   0, 0.5] ]
	}
	elsif($BL eq 'I') {
		return [ [-0.5, 0.5, 0.5], 
			 [ 0.5,-0.5, 0.5],
			 [ 0.5, 0.5,-0.5] ];
	}
	elsif($BL eq 'A') {
		return [ [ 1.0, 0.0, 0.0], 
			 [ 0.0, 0.5, 0.5],
			 [ 0.0,-0.5, 0.5] ];
	}
	elsif($BL eq 'B') {
		return [ [ 0.5, 0.0, 0.5], 
			 [ 0.0, 1.0, 0.0],
			 [-0.5, 0.0, 0.5] ];
	}
	elsif($BL eq 'C') {
		return [ [ 0.5, 0.5, 0.0], 
			 [-0.5, 0.5, 0.0],
			 [ 0.0, 0.0, 1.0] ];
	}
	return [ [ 1.0, 0.0, 0.0], 
		 [ 0.0, 1.0, 0.0],
		 [ 0.0, 0.0, 1.0] ];
}

sub MakeBORN
{
	$App->print("\n<b>Make BORN File from OUTCAR with NWRITE=3 for phonopy Non-Analytical Term Correction:</b>\n");

	my $VerboseFile   = $Args->GetGetArg(0);
	$VerboseFile      = "phonopy-verbose.out.txt" if($VerboseFile eq '');
	my $OUTCAR        = $Args->GetGetArg(1);
	$OUTCAR = 'Phonon/OUTCAR' if($OUTCAR eq '');
	my $OutFile       = $Args->GetGetArg(2);
	$OutFile          = 'BORN' if($OutFile eq '');
	my $ReadSuperCell = $Args->GetGetArg('ReadSuperCell');
$App->print("Verbose File : [$VerboseFile]\n");
$App->print("OUTCAR       : [$OUTCAR]\n");
$App->print("Out File     : [$OutFile]\n");
$App->print("ReadSuperCell: [$ReadSuperCell]\n");

	$App->print("\nRead non-equivalent sites from [$VerboseFile].\n");
	my $in = JFile->new($VerboseFile, "r") or die "$!: Cannot read [$VerboseFile].\n";
#-------------------------------- unit cell ---------------------------------
#Lattice vectors:
#Atomic positions (fractional):
#   *1 Na  0.00000000000000  0.00000000000000  0.00000000000000  22.990
#-------------------------------- super cell --------------------------------
#Lattice vectors:
#Atomic positions (fractional):
#   *1 Na  0.00000000000000  0.00000000000000  0.00000000000000  22.990 > 1
	if($ReadSuperCell) {
		$in->SkipTo(" super cell ");
	}
	$in->SkipTo("Atomic positions");
#   *1 Na  0.00000000000000  0.00000000000000  0.00000000000000  22.990 > 1
	my @atoms;
	my $c = 0;
	while(1) {
		my $line = $in->ReadLine();
		last if(!defined $line or $line =~ /^\s*---/);

		my ($i, $name, $x, $y, $z, $mass, $lt, $ieqsite) = Utils::Split("\\s+", $line);
		last if(!defined $z);

		$x += 0.0; $y += 0.0; $z += 0.0; $mass += 0.0;
		my $IsIndependent = 0;
		if($i =~ /\*(\d+)/) {
			$i = $1;
			$IsIndependent = 1;
		}

		$atoms[$c++] = {
			IsIndependent   => $IsIndependent,
			iSite           => $i,
			AtomName        => $name,
			x               => $x,
			y               => $y,
			z               => $z,
			Mass            => $mass,
			iEquivalentSite => $ieqsite,
			};
if($IsIndependent) {
	$App->print("Indep: $name($i) ($x, $y, $z) m=$mass\n");
}
else {
	$App->print("   equivalent to $ieqsite: $name($i) ($x, $y, $z) m=$mass\n");
}
	}
	$in->Close();

	$App->print("\nRead dielectric and BEC tensors from [$OUTCAR].\n");
	my $in = JFile->new($OUTCAR, "r") or die "$!: Cannot read [$OUTCAR].\n";
# MACROSCOPIC STATIC DIELECTRIC TENSOR (including local field effects in DFT)
	$in->SkipTo("MACROSCOPIC STATIC DIELECTRIC TENSOR");
# ------------------------------------------------------
	$in->ReadLine();
#           2.462755    -0.000000    -0.000000
	my @a1 = Utils::Split("\\s+", $in->ReadLine());
	my @a2 = Utils::Split("\\s+", $in->ReadLine());
	my @a3 = Utils::Split("\\s+", $in->ReadLine());
	my @eps = (@a1, @a2, @a3);
	$App->print("Dielectric tensor: [", join(', ', @eps), "]\n");

# BORN EFFECTIVE CHARGES (in e, cummulative output)
	$in->SkipTo("BORN EFFECTIVE CHARGES");
# -------------------------------------------------
	$in->ReadLine();
	my @atoms2;
	my $c = 0;
	while(1) {
# ion    1
		$in->ReadLine();
#    1     1.07590     0.00000     0.00000
#    2     0.00000     1.07590     0.00000
#    3     0.00000     0.00000     1.07590
		my ($i1, $i2, $i3);
		my @Z;
		($i1, $Z[0], $Z[1], $Z[2]) = Utils::Split("\\s+", $in->ReadLine());
		last if(!defined $Z[2]);
		($i2, $Z[3], $Z[4], $Z[5]) = Utils::Split("\\s+", $in->ReadLine());
		last if(!defined $Z[5]);
		($i3, $Z[6], $Z[7], $Z[8]) = Utils::Split("\\s+", $in->ReadLine());
		last if(!defined $Z[8]);

		for(my $i = 0 ; $i < 9 ; $i++) {
			$Z[$i] += 0.0;
		}
		$atoms2[$c] = {
			iSite => $c,
			pZij  => \@Z,
			};
		$App->print("BEC for atom $c: [", join(', ', @Z), "]\n");
		$c++;
	}
	$in->Close();

	$App->print("\nSave BORN file to [$OutFile].\n");
	my $out = JFile->new($OutFile, "w") or die "$!: Cannot write to [$OutFile].\n";
	$out->print("default value\n");
	print("default value\n");
	$out->printf("%8.5f %8.5f %8.5f  %8.5f %8.5f %8.5f  %8.5f %8.5f %8.5f\n", @eps);
	      printf("%8.5f %8.5f %8.5f  %8.5f %8.5f %8.5f  %8.5f %8.5f %8.5f\n", @eps);
	for(my $i = 0 ; $i < @atoms ; $i++) {
		my $pZij = $atoms2[$i]{pZij};

		if($atoms[$i]{IsIndependent}) {
			$out->printf("%8.5f %8.5f %8.5f   %8.5f %8.5f %8.5f   %8.5f %8.5f %8.5f\n", @$pZij);
			      printf("%8.5f %8.5f %8.5f   %8.5f %8.5f %8.5f   %8.5f %8.5f %8.5f\n", @$pZij);
		}
	}
	$out->Close();

	$App->print("\n<b>Finished: Make BORN File from OUTCAR with NWRITE=3 for phonopy Non-Analytical Term Correction</b>\n");
}

sub MakeGOnlyKPOINTS
{
	$App->print("\n<b>Make VASP KPOINTS File for G-point only calc:</b>\n");

	my $KPOINTSFile     = $Args->GetGetArg(0);
$App->print("KPOINTS : [$KPOINTSFile]\n");

	my $out = JFile->new($KPOINTSFile, "w") or die "$!: Cannot write to [$KPOINTSFile].\n";
	$out->print("Automatic mesh\n");
	$out->print("0\n");
	$out->print("Gamma\n");
	$out->print("1  1  1\n");
	$out->print("0. 0. 0.\n");
	$out->Close();

	$App->print("\n<b>Finished: Make VASP KPOINTS File for G-point only calc</b>\n");
}

sub ReadKList
{
	my ($KListDir) = @_;

my %GLabelHash = (
	gamma  => '\Gamma',
	lambda => '\Lambda',
	delta  => '\Delta',
	sigma  => '\Sigma',
	);

	my @files = glob(Deps::MakePath($KListDir, '*.klist', 0));
	return undef if(@files == 0);

	my $f = ''; 
	for(my $i = 0 ; $i < @files ; $i++) {
#print "$i: $files[$i]\n";
		if(-f $files[$i]) {
			$f = $files[$i];
			last;
		}
	}

	my $in = JFile->new($f, "r") or die "$!: Can not read [$f].\n";
#k-points along high symmetry linesA
	my $line = $in->ReadLine();
#15
	$line = $in->ReadLine();
#Line-mode
	$line = $in->ReadLine();
	my @kl;
	my $c = 0;
	if($line =~ /Line-mode/i) {
#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";
			my $sname = lc $name;
			$name = $GLabelHash{$sname} if($GLabelHash{$sname});
			$kl[$c] = {
				name => $name,
				kx   => $kx,
				ky   => $ky,
				kz   => $kz,
				};
#print "$c: $kl[$c]{name}: $kl[$c]{kx}, $kl[$c]{ky}, $kl[$c]{kz}\n";
			$c++;
			$kx0 = $kx;
			$ky0 = $ky;
			$kz0 = $kz;
		}
	}
	else {
		$in->rewind();
		while(1) {
			$line = $in->ReadLine();
			last if(!defined $line);
			last if($line =~ /^\s*END/i);

#print "l:[$line";
			my ($name, $kx, $ky, $kz, $n, $w) = Utils::Split("\\s+", $line);
			next if(!defined $w);
#print "l:[$name, $kx, $ky, $kz, $n, $w]\n";

#print "GLH: ", $GLabelHash{a}, "\n";
			my $sname = lc $name;
			$name = $GLabelHash{$sname} if($GLabelHash{$sname});
			$kl[$c] = {
				name => $name,
				kx   => ($kx+0.0) / $n,
				ky   => ($ky+0.0) / $n,
				kz   => ($kz+0.0) / $n,
				w    => $w
				};
#print "$c: $kl[$c]{name}: $kl[$c]{kx}, $kl[$c]{ky}, $kl[$c]{kz}\n";
			$c++;
		}
	}
	$in->Close();

	return \@kl;
}

sub MakeBandConf
{
	$App->print("\n<b>Make band.conf for phonopy:</b>\n");

	my $CARDir        = $Args->GetGetArg('CARDir');
	my $KListDir      = $Args->GetGetArg('KListDir');
	my $OutFile       = $Args->GetGetArg(0);
	my $PhonopyTarget = $Args->GetGetArg("PhonopyTarget");
	$PhonopyTarget = 'VASP' if(!defined $PhonopyTarget);
	my $nSupercell    = $Args->GetGetArg("nSupercell");
	$nSupercell = "2 2 2" if(!defined $nSupercell);
	my $BandPoints    = $Args->GetGetArg("BandPoints");
	my $fConversion   = $Args->GetGetArg("fConversion");
$App->print("Out File    : [$OutFile]\n");
$App->print("CAR Dir     : [$CARDir]\n");
$App->print("KList Dir   : [$KListDir]\n");
$App->print("Target      : $PhonopyTarget\n");
$App->print("nSupercell  : [$nSupercell]\n");
$App->print("BandPoints  : [$BandPoints]\n");
$App->print("fConversion : [$fConversion]\n");
#	my $AtomTypeList = &BuildAtomTypeListLine($CARDir);
	my ($AtomTypeList, $SiteIndex, $Crystal) = &GetStructureInfo($CARDir);
$App->print("Atom Types  : [$AtomTypeList]\n");
$App->print("Site Indexes: [$SiteIndex]\n");

	my ($CrystalSystem, $BravaisLattice) = $Crystal->GuessBravaisLattice(0.001, 0.01, 0.02);
$App->print("Bravais Lattice: [$BravaisLattice - $CrystalSystem]\n");

	my $pm = &GetMatrixConventionalToPrimitiveCell($BravaisLattice);
#print "pm: $pm->[0][0] $pm->[0][1] $pm->[0][2]\n";
#print "pm: $pm->[1][0] $pm->[1][1] $pm->[1][2]\n";
#print "pm: $pm->[2][0] $pm->[2][1] $pm->[2][2]\n";

	my $pKList = &ReadKList($KListDir);
	my $kpts   = '';
	my $labels = '';
	for(my $i = 0 ; $i < @$pKList ; $i++) {
		my $s = $pKList->[$i];
		$kpts   .= "$s->{kx} $s->{ky} $s->{kz}   ";
		$labels .= "$s->{name} ";
	}
	
	my $out = JFile->new($OutFile, "w") or die "$!: Cannot write to [$OutFile].\n";
	$out->print("ATOM_NAME = $AtomTypeList\n");
	$out->print("DIM = $nSupercell\n");
	$out->print("PRIMITIVE_AXIS = ");
	for(my $i = 0 ; $i < 3 ; $i++) {
		for(my $j = 0 ; $j < 3 ; $j++) {
			$out->printf("%3.1f ", $pm->[$i][$j]);
		}
		$out->print("  ");
	}
	$out->print("\n");
	$out->print("BAND = $kpts\n");
#	$out->print("BAND = 0.5 0.5 0.5  0.0 0.0 0.0  0.5 0.5 0.0  0.0 0.5 0.0\n");
	if(defined $BandPoints) {
		$out->print("BAND_POINTS = $BandPoints\n");
	}
	else {
		$out->print("#BAND_POINTS = 51\n");
	}
	$out->print("BAND_LABELS = $labels\n");
#	$out->print("#BAND_LABELS = X \\Gamma L\n");
	$out->print("#BAND_CONNECTION = .TRUE.\n");
	if($fConversion) {
		$out->print("FREQUENCY_CONVERSION_FACTOR = $fConversion\n");
	}
	else {
		$out->print("#FREQUENCY_CONVERSION_FACTOR = 521.471\n");
	}
	my $fConversion   = $Args->GetGetArg("fConversion");
	if($PhonopyTarget eq 'VASP') {
		$out->print("#FORCE_CONSTANTS = READ\n");
	}
	else {
		$out->print("FORCE_CONSTANTS = READ\n");
	}
	$out->Close();

	$App->print("\n<b>Finished: Make band.conf for phonopy:</b>\n");
}

sub MakePDOSConf
{
	$App->print("\n<b>Make pdos.conf for phonopy:</b>\n");

	my $CARDir        = $Args->GetGetArg('CARDir');
	my $OutFile       = $Args->GetGetArg(0);
	my $PhonopyTarget = $Args->GetGetArg("PhonopyTarget");
	$PhonopyTarget = 'VASP' if(!defined $PhonopyTarget);
	my $nSupercell    = $Args->GetGetArg("nSupercell");
	$nSupercell = "2 2 2" if(!defined $nSupercell);
	my $MeshPoints    = $Args->GetGetArg("MeshPoints");
	$MeshPoints = "8 8 8" if(!defined $MeshPoints);
	my $fConversion   = $Args->GetGetArg("fConversion");
$App->print("Out File  : [$OutFile]\n");
$App->print("CAR Dir   : [$CARDir]\n");
$App->print("Target    : $PhonopyTarget\n");
$App->print("nSupercell: [$nSupercell]\n");
$App->print("MeshPoints: [$MeshPoints]\n");
$App->print("fConversion : [$fConversion]\n");
#	my $AtomTypeList = &BuildAtomTypeListLine($CARDir);
	my ($AtomTypeList, $SiteIndex) = &GetStructureInfo($CARDir);
$App->print("Atom Types  : [$AtomTypeList]\n");
$App->print("Site Indexes: [$SiteIndex]\n");

	my $out = JFile->new($OutFile, "w") or die "$!: Cannot write to [$OutFile].\n";
	$out->print("ATOM_NAME = $AtomTypeList\n");
	$out->print("DIM = $nSupercell\n");
	$out->print("MP = $MeshPoints\n");
	$out->print("PDOS = $SiteIndex\n");
	$out->print("#DOS_RANGE = 0 40 0.1\n");
	$out->print("#SIGMA = 0.1\n");
	$out->print("#TETRAHEDRON\n");
	$out->print("EIGENVECTORS = .TRUE.\n");
	$out->print("#PROJECTION_DIRECTION = -1 1 1\n");
	$out->print("#DEBYE_MODEL = .TRUE.\n");
	if($fConversion) {
		$out->print("FREQUENCY_CONVERSION_FACTOR = $fConversion\n");
	}
	else {
		$out->print("#FREQUENCY_CONVERSION_FACTOR = 521.471\n");
	}
	if($PhonopyTarget eq 'VASP') {
		$out->print("#FORCE_CONSTANTS = READ\n");
	}
	else {
		$out->print("FORCE_CONSTANTS = READ\n");
	}
	$out->Close();

	$App->print("\n<b>Finished: Make pdos.conf for phonopy:</b>\n");
}

sub MakeMeshConf
{
	$App->print("\n<b>Make mesh.conf for phonopy:</b>\n");

	my $CARDir        = $Args->GetGetArg('CARDir');
	my $OutFile       = $Args->GetGetArg(0);
	my $PhonopyTarget = $Args->GetGetArg("PhonopyTarget");
	$PhonopyTarget = 'VASP' if(!defined $PhonopyTarget);
	my $nSupercell    = $Args->GetGetArg("nSupercell");
	$nSupercell = "2 2 2" if(!defined $nSupercell);
	my $MeshPoints    = $Args->GetGetArg("MeshPoints");
	$MeshPoints = "8 8 8" if(!defined $MeshPoints);
	my $fConversion   = $Args->GetGetArg("fConversion");
$App->print("Out File  : [$OutFile]\n");
$App->print("CAR Dir   : [$CARDir]\n");
$App->print("Target    : $PhonopyTarget\n");
$App->print("nSupercell: [$nSupercell]\n");
$App->print("MeshPoints: [$MeshPoints]\n");
$App->print("fConversion : [$fConversion]\n");
#	my $AtomTypeList = &BuildAtomTypeListLine($CARDir);
	my ($AtomTypeList, $SiteIndex) = &GetStructureInfo($CARDir);
$App->print("Atom Types  : [$AtomTypeList]\n");
$App->print("Site Indexes: [$SiteIndex]\n");

	my $out = JFile->new($OutFile, "w") or die "$!: Cannot write to [$OutFile].\n";
	$out->print("ATOM_NAME = $AtomTypeList\n");
	$out->print("DIM = $nSupercell\n");
	$out->print("MP = $MeshPoints\n");
	$out->print("#WRITE_MESH = .FALSE.\n");
	$out->print("#GAMMA_CENTER = .TRUE.\n");
	$out->print("#MP_SHIFT = 1/2 1/2 1/2\n");
	if($fConversion) {
		$out->print("FREQUENCY_CONVERSION_FACTOR = $fConversion\n");
	}
	else {
		$out->print("#FREQUENCY_CONVERSION_FACTOR = 521.471\n");
	}
	if($PhonopyTarget eq 'VASP') {
		$out->print("#FORCE_CONSTANTS = READ\n");
	}
	else {
		$out->print("FORCE_CONSTANTS = READ\n");
	}
	$out->print("\n");
	$out->print("# for thermal properties\n");
	$out->print("#TPROP = .TRUE.\n");
	$out->print("#TMIN = 0\n");
	$out->print("#TMAX = 2000\n");
	$out->print("#TSTEP = 100\n");
	$out->print("\n");
	$out->print("# for thermal displacement\n");
	$out->print("#TDISP = .TRUE.\n");
	$out->print("#PROJECTION_DIRECTION = 1 1 0\n");
	$out->print("#CUTOFF_FREQUENCY\n");
	$out->Close();

	$App->print("\n<b>Finished: Make mesh.conf for phonopy:</b>\n");
}

sub MakeKPOINTSforSPOSCAR
{
	$App->print("\n<b>Make VASP KPOINTS File for SPOSCAR made by phonopy:</b>\n");

	my $CARDir      = $Args->GetGetArg('CARDir');
	my $IsPrint     = 1;
	my $KPOINTSFile = $Args->GetGetArg(0);
	my $aKProduct   = $Args->GetGetArg('aKProduct');
$App->print("KPOINTS  : [$KPOINTSFile]\n");
$App->print("CAR Dir  : [$CARDir]\n");
$App->print("aKProduct: $aKProduct\n");

	my $vasp = new VASP;
	my $Crystal = $vasp->ReadStructureFromCARFiles($CARDir, $IsPrint);
	my ($nx, $ny, $nz, $NKREDX, $NKREDY, $NKREDZ) = $vasp->AnalyzeaKProduct($Crystal, $aKProduct);
$App->print("[nx ny nz]: [$nx $ny $nz]\n");

	my $out = JFile->new($KPOINTSFile, "w") or die "$!: Cannot write to [$KPOINTSFile].\n";
	$out->print("Automatic mesh\n");
	$out->print("0\n");
	$out->print("Gamma\n");
	$out->print("$nx $ny $nz\n");
	$out->print("0. 0. 0.\n");
	$out->Close();

	$App->print("\n<b>Finished: Make VASP KPOINTS File for SPOSCAR made by phonopy</b>\n");
}

sub MakePhonopyINCAR
{
	$App->print("\n<b>Make VASP INCAR File for phonopy:</b>\n");

	my $INCARFile     = $Args->GetGetArg(0);
	my $PhonopyTarget = $Args->GetGetArg("PhonopyTarget");
	$PhonopyTarget = 'VASP' if(!defined $PhonopyTarget);
	my $Precision     = $Args->GetGetArg("Precision");
	$Precision = 'Accurate' if(!defined $Precision);
	my $ENCUT         = $Args->GetGetArg("ENCUT");
	$ENCUT = 500 if(!defined $ENCUT);
	my $EDIFF         = $Args->GetGetArg("EDIFF");
	$EDIFF = 1.0e-8 if(!defined $EDIFF);
	my $IBRION        = ($PhonopyTarget eq 'VASP')? -1 : 8;
	my $ISMEAR        = $Args->GetGetArg("ISMEAR");
	$ISMEAR = 0 if(!defined $ISMEAR);
	my $SIGMA         = $Args->GetGetArg("SIGMA");
	$SIGMA = 0.01 if(!defined $SIGMA);

$App->print("INCAR    : [$INCARFile]\n");
$App->print("Target   : $PhonopyTarget\n");
$App->print("IBRION   : $IBRION\n");
$App->print("Precision: $Precision\n");
$App->print("ENCUT    : $ENCUT\n");
$App->print("EDIFF    : $EDIFF\n");
$App->print("ISMEAR   : $ISMEAR\n");
$App->print("SIGMA    : $SIGMA\n");

	my $out = JFile->new($INCARFile, "w") or die "$!: Cannot write to [$INCARFile].\n";

	$out->print(" PREC    = $Precision\n");
	$out->print(" ENCUT   = $ENCUT\n");
	$out->print(" IBRION  = $IBRION\n");
	$out->print(" EDIFF   = $EDIFF\n");
	$out->print(" IALGO   = 38\n");
	$out->print(" ISMEAR  = $ISMEAR\n");
	$out->print(" SIGMA   = $SIGMA\n");
	$out->print(" LREAL   = .FALSE.\n");
	$out->print(" ADDGRID = .TRUE.\n");
	$out->print(" LWAVE   = .FALSE.\n");
	$out->print(" LCHARG  = .FALSE.\n");
	$out->Close();
 
	$App->print("\n<b>Finished: Make VASP INCAR File for phonopy</b>\n");
}