#!/usr/bin/perl

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

use strict;
#use warnings;

use ProgVars;
use Math::Matrix;
use Math::MatrixReal;

use Deps;
use Utils;
use JFile;

use MyApplication;
use Sci::Algorism;
use Sci::GeneralFileFormat;

use Crystal::CIF;
use Crystal::Crystal;
use Crystal::SpaceGroup;
use Crystal::AtomType;
use Crystal::WIEN2k;

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

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

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

#===============================================
# 定数
#===============================================
#my $OS = Deps::OS();
my $SGROUPPath = ProgVars::SGROUPPath();
my $TEEPath = ProgVars::TEEPath();

#===============================================
# 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();

#==========================================
# コマンドラインオプション読み込み
#==========================================
$App->AddArgument("--Action",
		"--Action=[MakeStruct|MakeSymmetrizedStruct|MakeSymmetrizedStructFile|MakeInt|MakeInsp|MakeInM|"
			."MakeCIF|ShowScript|MakeScript|"
			."AnalyzeOptStructure|MakeClonDir|BackupStruct]",       '');
$App->AddArgument("--OutFile",          "--OutFile=FileName",       "");
$App->AddArgument("--NewStructFile",    "--NewStructFile=FileName", "");
$App->AddArgument("--FittingFile",      "--Fitting=FileName", "");
$App->AddArgument("--Parameter",        "--Parameter=[Volume|coa] (def: Volume)", "Volume");
$App->AddArgument("--Task",             "--Task=[SCF|DOS|Band] (def: SCF)",   "SCF");
$App->AddArgument("--UseComplex",       "--UseComplex=[0|1] (def: 0)",        0);
$App->AddArgument("--Parallel",         "--Parallel=[0|1] (def: 0)",          0);
$App->AddArgument("--SpinPolarized",    "--SpinPolarized=[0|1] (def: 0)",     0);
$App->AddArgument("--SpinOrbit",        "--SpinOrbit=[0|1] (def: 0)",         0);
$App->AddArgument("--OrbitalPotential", "--OribitalPotential=[0|1] (def: 0)", 0);
$App->AddArgument("--EnergyConvergence", "--EnergyConvergence=val (def: 0.001)", 0.001);
$App->AddArgument("--ForceConvergence",  "--ForceConvergence=val (def: 1)");
$App->AddArgument("--StructOnly",  "--StructOnly=[0|1] (def: 0)");
$App->AddArgument("--TotalPDOSOnly", 
		"--TotalPDOSOnly=[0|1] (def: 0): Specify if decmposed PDOS is stored in *.int file", 0);
$App->AddArgument("--OnceForEachAtomType", 
		"--OnceForEachAtomType=[0|1] (def: 1): Write PDOS data only once for the same atom", 1);
$App->AddArgument("--ForceOverWrite", "--ForceOverWrite=[0|1] (def: 0): Overwrite Directory if it exists", 0);
$App->AddArgument("--DoSGROUP",      "--DoSGROUP=[0|1] (def: 1)");
$App->AddArgument("--PrimitiveCell", "--PrimitiveCell=[0|1] (def: 0)");
$App->AddArgument("--TOL",           "--TOL=val (def: 0.0001)");
$App->AddArgument("--UseSameBasisForSameAtom",  "--UseSameBasisForSameAtom=[0|1] (def: 0)");
$App->AddArgument("--DebugMode", "--DebugMode: Set DebugMode", '');
exit 1 if($App->ReadArgs(0) != 1);
my $Args = $App->Args();
#my $form = new CGI;
#$Args->SetCGIForm($form);
#$Args->parseInput($WebCharCode);

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

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

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

my $TOL = $Args->GetGetArg("TOL");
   $TOL = 0.0001 if(!defined $TOL or $TOL <= 0.0);

my $ret = 0;
if($Action =~ /^MakeSymmetrizedStruct/i) {
	&MakeSymmetrizedStructFile();
}
elsif($Action =~ /MakeStruct/i) {
	&MakeStructFile();
}
elsif($Action =~ /MakeInt/i) {
	&MakeIntFile();
}
elsif($Action =~ /MakeInsp/i) {
	&MakeInspFile();
}
elsif($Action =~ /MakeInM/i) {
	&MakeInMFile();
}
elsif($Action =~ /MakeCIF/i) {
	&MakeCIFFile();
}
elsif($Action =~ /BackupStruct/i) {
	&BackupStructFile();
}
elsif($Action =~ /ShowScript/i) {
	&ShowScript();
}
elsif($Action =~ /MakeScript/i) {
	$ret = 1 if(&MakeScript() <= 0);
}
elsif($Action =~ /AnalyzeOptStructure/i) {
	$ret = 1 if(&AnalyzeOptStructure() <= 0);
}
elsif($Action =~ /MakeClonDir/i) {
	&MakeClonDir();
}
else {
	$App->print("Error: Invald Action: $Action\n");
}

#Utils::EndHTML();

exit $ret;

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

#==========================================
# &Subroutines
#==========================================
sub MakeClonDir
{
	$App->print("\n\n<b>Make Clone Dir: Copy all files with new prefix:</b>\n");

	my $StructOnly = $Args->GetGetArg("StructOnly");
	my $ForceOverWrite = $Args->GetGetArg("ForceOverWrite");
	my $Prefix    = $Args->GetGetArg(0);
	my $CurDir    = $Args->GetGetArg(1);
	my $NewPrefix = $Args->GetGetArg(2);
	my $ParentDir = $CurDir;
	$ParentDir =~ s/[\\\/]$//;
	$ParentDir =~ s/[\\\/][^\\\/]*?$//;
	my $NewDir = Deps::MakePath($ParentDir, $NewPrefix);
	$App->print("  Prefix   : $Prefix\n");
	$App->print("  NewPrefix: $NewPrefix\n");
	$App->print("  CurDir   : $CurDir\n");
	$App->print("  ParentDir: $ParentDir\n");
	$App->print("  NewDir   : $NewDir\n");
	$App->print("  StructOnly    : $StructOnly\n");
	$App->print("  ForceOverWrite: $ForceOverWrite\n");

	if(-f $NewDir) {
		$App->print("  Error: File [$NewDir] exists!\n");
		$App->print("  Program terminated\n");
		return -1;
	}

	if(-d $NewDir) {
		$App->print("  Error: [$NewDir] exists!\n");
		if(!$ForceOverWrite) {
			$App->print("  Program terminated\n");
			return -1;
		}
		else {
			$App->print("  Overwrite it.\n");
		}
	}
	else {
		if(mkdir($NewDir)) {
			$App->print("    Create [$NewDir]\n");
		}
		else {
			$App->print("    [$NewDir] could not be created\n");
			$App->print("  Program terminated\n");
			return -1;
		}
	}

	if($StructOnly) {
		my $fname = Deps::MakePath($CurDir, "$Prefix.struct", 0);
		my $newpath = Deps::MakePath($NewDir, "$NewPrefix.struct", 0);
		$App->print("  cp $fname : $newpath\n");
		system("cp $fname $newpath");
		unless(-f $newpath) {
			$App->print("    Error: $newpath was not created.\n");
		}
		$fname = Deps::MakePath($CurDir, "$Prefix.cif", 0);
		$newpath = Deps::MakePath($NewDir, "$NewPrefix.cif", 0);
		$App->print("  cp $fname : $newpath\n");
		system("cp $fname $newpath");
		unless(-f $newpath) {
			$App->print("    Error: $newpath was not created.\n");
		}
	}
	else {
		my $fname = Deps::MakePath($CurDir, "*", 0);
		my @files = glob($fname);
		$fname = Deps::MakePath($CurDir, ".*", 0);
		@files = (@files, glob($fname));
		my $RegPrefix = Utils::RegExpQuote($Prefix);
		my $RegNewPrefix = Utils::RegExpQuote($NewPrefix);
		foreach my $f (@files) {
			next if(-d $f);
			my ($drive, $dir, $filename, $ext, $lastdir, $filebody)
				= Deps::SplitFilePath($f);
			next if($filename =~ /^(\.|:|tmp|core)/);
			next if($filename =~ /(\.def|\.\d+|_\d+|\.ori?g[^\.]*|\.prev[^\.]*|\.png)$/);
			next if($filename !~ /^$RegPrefix\./ and $filename ne 'xcrysden.klist');
			next if($filename =~ /^${RegPrefix}_/);
	
			$filename =~ s/^$RegPrefix/$RegNewPrefix/;
			my $newpath = Deps::MakePath($NewDir, $filename);
			$App->print("  cp $f : $newpath\n");
			system("cp $f $newpath");
			unless(-f $newpath) {
				$App->print("    Error: $newpath was not created.\n");
			}
		}
	}
	$App->print("\n\n<b>Make Clone Dir: Finished</b>\n\n");
}

sub AnalyzeOptStructure
{
	$App->print("\n\n<b>Analyze Optimized Structure:</b>\n");

	my $OutFile       = $Args->GetGetArg("OutFile");
	my $NewStructFile = $Args->GetGetArg("NewStructFile");
	my $FittingFile   = $Args->GetGetArg("FittingFile");
	my $Parameter = $Args->GetGetArg("Parameter");
	my $Prefix = $Args->GetGetArg(0);
	$App->print("  OutFile      : $OutFile\n");
	$App->print("  NewStructFile: $NewStructFile\n");
	$App->print("  FittingFile  : $FittingFile\n");
	$App->print("  Prefix       : $Prefix\n");
	$App->print("  Parameter    : $Parameter\n");
	my $Debug  = $Args->GetGetArg("DebugMode");
	$Debug = 0 unless(defined $Debug);

	my $fmask = "";
	my $ParamName = "";
	if($Parameter =~ /volume/i) {
		$fmask = "$Prefix*vol*.scf";
		$ParamName = "Volume";
	}
	elsif($Parameter =~ /coa/i) {
		$fmask = "$Prefix*coa*.scf";
		$ParamName = "coa";
	}
	else {
		print "\n***Error: Parameter [$Parameter] is not supported.\n";
		return -1;
	}

	my $out = new JFile($OutFile, "w");
	my $fit = new JFile($FittingFile, "w");
	if($out) {
		$out->Write("Analyze Optimized Structure:\n");
	}

	my $WIEN2k = new WIEN2k;

	$App->print("\n");
	$App->print("  Search for [$fmask]\n");
	my @files = glob($fmask);
	@files = sort { my ($aa) = ($a =~ /([\-\.\d]+)\.[^\.]*$/);
			my ($bb) = ($b =~ /([\-\.\d]+)\.[^\.]*$/);
#			print "$a / $b : aa=$aa / bb=$bb\n";
			return $aa <=> $bb; } @files;
	$App->print("\n");
	for(my $i = 0 ; $i < @files ; $i++) {
		$App->print("    $files[$i]\n");
	}

# 平均値からの差を最少自乗法にかける
	my $toDeg = Sci::todeg();
	my @StructFile;
	my @SCFFile;
	my @Param;
	my @LattStr;
	my @Energy;
	my $AvrEng   = 0.0;
	my $AvrParam = 0.0;
	my $nStruct = 0;
	for(my $i = 0 ; $i < @files ; $i++) {
		my $scf = WIEN2k::GetSCFFile($files[$i]);
		$Energy[$nStruct] = $WIEN2k->ReadSCFFileParameter($scf, "ENE");
#print "$nStruct: E=$Energy[$nStruct]\n";
		next unless(defined $Energy[$nStruct]);
		($Energy[$nStruct]) = ($Energy[$nStruct] =~ /([\-\+\d\.eE]+)\s*$/);
		$Param[$nStruct] = $WIEN2k->ReadSCFFileParameter($scf, "VOL");
		next unless(defined $Param[$nStruct]);
		($Param[$nStruct]) = ($Param[$nStruct] =~ /([\-\+\d\.eE]+)\s*$/);
		$LattStr[$nStruct] = $WIEN2k->ReadSCFFileParameter($scf, "LAT");
		next unless(defined $LattStr[$nStruct]);

		if($Parameter =~ /coa/i) {
			my ($a, $b, $c, $alpha, $beta, $gamma) 
				= ($LattStr[$nStruct] =~ /([\d\.]+)\s+([\d\.]+)\s+([\d\.]+)\s+([\d\.]+)\s+([\d\.]+)\s+([\d\.]+)\s*$/);
			$alpha = Utils::Round($alpha*$toDeg, 4);
			$beta  = Utils::Round($beta*$toDeg,  4);
			$gamma = Utils::Round($gamma*$toDeg, 4);
#print "Latt: $a $b $c  $alpha $beta $gamma\n";
			$Param[$nStruct] = $c / $a;
		}

		$StructFile[$nStruct] = WIEN2k::GetStructFile($files[$i]);
		$SCFFile[$nStruct]    = WIEN2k::GetSCFFile($files[$i]);
		$AvrEng += $Energy[$nStruct];
		$AvrParam += $Param[$nStruct];

		$App->print("  $i: $SCFFile[$i]: V=$Param[$nStruct] a0^3 E=$Energy[$nStruct] Ry\n");
		if($out) {
			$out->Write("  $i: $SCFFile[$i]: V=$Param[$nStruct] a0^3 E=$Energy[$nStruct] Ry\n");
		}
		$nStruct++;
	}

	$AvrEng /= $nStruct;
	$AvrParam /= $nStruct;
	$App->print("  Average $ParamName: $AvrParam\n");
	$App->print("  Average Energy: $AvrEng\n");
	for(my $i = 0 ; $i < $nStruct ; $i++) {
		$Energy[$i] -= $AvrEng;
		$Energy[$i] = Utils::Round($Energy[$i], 6);
		$Param[$i]  -= $AvrParam;
		$Param[$i] = Utils::Round($Param[$i], 6);

#		$App->print("  $i: V=$Param[$i] a0^3 E=$Energy[$i] Ry\n");
	}

# 最少自乗
	$App->print("\n");
	$App->print("  Execute least squre procedure.\n");
	$App->print("    nStruct: $nStruct\n");
	my @Coeff;
	my $nOrder = $nStruct-2;
#	$nOrder = $nStruct-2 if($nStruct <= $nOrder+2);
	$App->print("    nOrder: $nOrder\n");
	if($nOrder < 2) {
		print "\n***Error: # of data [$nStruct] is too fee.\n";
		return -1;
	}
	my $ret = Algorism::lstsq(\@Param, \@Energy, $nStruct, $nOrder, \@Coeff);
	if($ret != 0) {
		print "\n***Error in Algorism::lstsq with ErrorCode=$ret.\n";
		return -1;
	}

# 最少自乗の結果をチェック
	if($fit) {
		$fit->Write("Fitting result\n");
		$fit->Write("$nStruct\n");
	}
	for(my $i = 0 ; $i < $nOrder+1 ; $i++) {
		print "    C[$i]=$Coeff[$i]\n";
		if($out) {
			$out->Write("    C[$i]=$Coeff[$i]\n");
		}
	}
	$App->print("    $ParamName: dE(data) - dE(fit) = dE(diff)\n");
	if($out) {
		$out->Write("    $ParamName: dE(data) - dE(fit) = dE(diff)\n");
	}
	for(my $i = 0 ; $i < $nStruct ; $i++) {
		my $cal = Sci::Polynomial($Param[$i], $nOrder, @Coeff);
		$cal = Utils::Round($cal, 6);
		my $d = $Energy[$i] - $cal;
		$App->printf("    %8.4f: %8.4f - %8.4f = %8.4f\n",
				$Param[$i], $Energy[$i], $cal, $d);
		my $v = $AvrParam + $Param[$i];
		my $e = $AvrEng + $Energy[$i];
		my $c = $AvrEng + $cal;
		if($out) {
			$out->Write("$v\t$e\t$c\t$d\n");
		}
		if($fit) {
#			$fit->Write("$Param[$i]\t$Energy[$i]\t$cal\t$d\n");
			$fit->Write("$v\t$Energy[$i]\t$cal\t$d\n");
		}
	}
#最小点を捜す
	my $DiffFunc = "$Coeff[1]";
	$DiffFunc = "$DiffFunc+$Coeff[2]*2*\$x" if($Coeff[2] >= 0);
	$DiffFunc = "$DiffFunc$Coeff[2]*2*\$x" if($Coeff[2] < 0);
	for(my $i = 3 ; $i < $nOrder+1 ; $i++) {
		my $i1 = $i-1;
		if($Coeff[$i] >= 0) {
			$DiffFunc = "$DiffFunc+$Coeff[$i]*$i*\$x**$i1";
		}
		else {
			$DiffFunc = "$DiffFunc$Coeff[$i]*$i*\$x**$i1";
		}
	}
#print "func: $DiffFunc\n";
	my $minParam;
	my $minEng;
	$minParam = -$Coeff[1] / 2.0 / $Coeff[2];
#print "Init: $minParam\n";
	$minParam = Algorism::SolveByNewton($DiffFunc, $minParam, 0.001, abs($Param[0]*1.0e-6), 1.0e-10, 100);
	$minEng = Sci::Polynomial($minParam, $nOrder, @Coeff);

	my $e = $minEng;
	$minEng = Utils::Round($AvrEng + $minEng, 10);
	my $v = $minParam;
	$minParam = Utils::Round($AvrParam + $minParam, 8);
	$e = Utils::Round($e, 6);
	$v = Utils::Round($v, 6);
	$App->print("     Minimum energy $minEng ($e) at $ParamName=$minParam ($v)\n");
	if($out) {
		$out->Write("     Minimum energy $minEng ($e) at $ParamName=$minParam ($v)\n");
	}

# .structファイルのVolume / c/a / b/aのupdate
	my $Crystal = $WIEN2k->ReadStructFile($StructFile[0], $Debug, tol => $TOL);
	my ($a,$b,$c,$alpha,$beta,$gamma) = $Crystal->LatticeParameters();
	if($Parameter =~ /volume/i) {
		my $V = $Crystal->Volume(1);
		my $k = ($minParam / $V) ** (1.0/3.0);
#print "k: $k = $minVol / $V\n";
		$a = Utils::Round($a*$k, 6);
		$b = Utils::Round($b*$k, 6);
		$c = Utils::Round($c*$k, 6);
		$alpha = Utils::Round($alpha, 6);
		$beta  = Utils::Round($beta, 6);
		$gamma = Utils::Round($gamma, 6);
		$App->print("  New lattice parameters: $a, $b, $c, $alpha, $beta, $gamma\n");
		$Crystal->SetLatticeParameters($a, $b, $c, $alpha, $beta, $gamma);
	}
	elsif($Parameter =~ /coa/i) {
		my $ca = $c * $a;
		my $coa = $minParam;
		$a = sqrt($ca / $coa);
# Tetra.の場合
		$b = $a;
		$c = $ca / $a;
		$a = Utils::Round($a, 6);
		$b = Utils::Round($b, 6);
		$c = Utils::Round($c, 6);
		$alpha = Utils::Round($alpha, 6);
		$beta  = Utils::Round($beta, 6);
		$gamma = Utils::Round($gamma, 6);
		$App->print("  New lattice parameters: $a, $b, $c, $alpha, $beta, $gamma\n");
		$Crystal->SetLatticeParameters($a, $b, $c, $alpha, $beta, $gamma);
	}

#座標の最小点を捜す
	my @wien;
	my @crystal;
	for(my $i = 0 ; $i < $nStruct ; $i++) {
		my $fname = $StructFile[$i];
#		$App->print("  Read [$fname].\n");
		$wien[$i] = new WIEN2k;
		$crystal[$i] = $wien[$i]->ReadStructFile($fname, $Debug, tol => $TOL);
	}
	my @AtomSites = $Crystal->GetCAsymmetricAtomSiteList();
	for(my $i = 0 ; $i < @AtomSites ; $i++) {
		my $IsXChanged = 0;
		my $IsYChanged = 0;
		my $IsZChanged = 0;
		my (@x, @y, @z);
		for(my $is = 0 ; $is < $nStruct ; $is++) {
			my @AsymAtomSites = $crystal[$is]->GetCAsymmetricAtomSiteList();
			my $site = $AsymAtomSites[$i];
			($x[$is], $y[$is], $z[$is]) = $site->Position();
			if($is > 0) {
				$IsXChanged = 1 if(abs($x[$is] - $x[0]) > 0.0001);
				$IsYChanged = 1 if(abs($y[$is] - $y[0]) > 0.0001);
				$IsZChanged = 1 if(abs($z[$is] - $z[0]) > 0.0001);
			}
		}
		my ($xp, $yp, $zp) = $AtomSites[$i]->Position();
		if($IsXChanged) {
			my @Coeff;
			my $ret = Algorism::lstsq(\@Param, \@x, $nStruct, $nOrder, \@Coeff);
			$xp = $Coeff[0];
			for(my $i = 1; $i <= $nOrder ; $i++) {
				$xp += $Coeff[$i] * ($v ** $i);
			}
			$xp = Utils::Round($xp, 6);
			$App->print("    x for atom $i was optimized to x=$xp.\n");
		}
		if($IsYChanged) {
			my @Coeff;
			my $ret = Algorism::lstsq(\@Param, \@y, $nStruct, $nOrder, \@Coeff);
			$yp = $Coeff[0];
			for(my $i = 1; $i <= $nOrder ; $i++) {
				$yp += $Coeff[$i] * ($v ** $i);
			}
			$yp = Utils::Round($yp, 6);
			$App->print("    y for atom $i was optimized to y=$yp.\n");
		}
		if($IsZChanged) {
			my @Coeff;
			my $ret = Algorism::lstsq(\@Param, \@z, $nStruct, $nOrder, \@Coeff);
			$zp = $Coeff[0];
			for(my $i = 1; $i <= $nOrder ; $i++) {
				$zp += $Coeff[$i] * ($v ** $i);
			}
			$zp = Utils::Round($zp, 6);
			$App->print("    z for atom $i was optimized to z=$zp.\n");
		}
		$AtomSites[$i]->SetPosition($xp, $yp, $zp);
	}

#新しい.structファイルを保存する
	$App->print("  Save Optimized Structure to [$NewStructFile].\n");
	unless($WIEN2k->SaveStructFile($Crystal, $NewStructFile, 0, 0)) {
		$App->print("\n***Error: Can not write to [$NewStructFile].\n");
		return -1;
	}

	$out->Close() if($out);
	$fit->Close() if($fit);

	$App->print("\n<b>Analyze Optimized Structure: Finished</b>\n");
	return 1;
}

sub MakeScript
{
	$App->print("\n\n<b>Make WIEN2k Script:</b>\n");

	my $Task       = $Args->GetGetArg("Task");
	my $UseComplex = $Args->GetGetArg("UseComplex");
	my $Parallel   = $Args->GetGetArg("Parallel");
	my $SP         = $Args->GetGetArg("SpinPolarized");
	my $SO         = $Args->GetGetArg("SpinOrbit");
	my $Orb        = $Args->GetGetArg("OrbitalPotential");
#	$UseComplex = 1 if($SP or $SO);
	my $ec         = $Args->GetGetArg("EnergyConvergence");
	$ec = 0.001 unless(defined $ec);
	my $fc         = $Args->GetGetArg("ForceConvergence");
	$fc = 1.0 unless(defined $fc);

	my $FileName   = $Args->GetGetArg(0);
	unless($FileName) {
		$App->print("Script File name should be specified.\n");
		$App->print("  usage: WIEN2k.pl --Action=MakeScript --Task=[SCF|Band|DOS] ScriptFile\n");
		return 0;
	}

	$App->print("  Task: $Task\n");
	$App->print("  UseComplex      : $UseComplex\n");
	$App->print("  Parallel        : $Parallel\n");
	$App->print("  SpinPolarized   : $SP\n");
	$App->print("  SpinOrbit       : $SO\n");
	$App->print("  OrbitalPotential: $Orb\n");
	$App->print("  Script File Name: $FileName\n");
	$App->print("\n");

	my $WIEN2k = new WIEN2k;
	my @lines = $WIEN2k->BuildScript($Task, $Parallel, $UseComplex, $SP, $SO, $Orb);

	my $out = new JFile($FileName, "w");
	unless($out) {
		$App->print("Error: Can not write to [$FileName].\n\n");
		return -1;
	}
	$out->Write("#!/bin/bash\n\n");
	if($Task =~ /SCF/i) {
		my $sp = '';
		$sp = 'sp' if($SP);
		my $opt = '';
		$opt = "$opt -c" if($UseComplex);
		$opt = "$opt -so" if($SO);
		$opt = "$opt -orb" if($Orb);
		$opt = "$opt -i 40";
		$opt = "$opt -fc $fc" if($fc);
		$opt = "$opt -ec $ec";
		$opt = "$opt -p" if($Parallel);
		$out->Write("run${sp}_lapw$opt\n");
		$App->print("run${sp}_lapw$opt\n");
	}
	else {
		for(my $i = 0 ; $i < @lines ; $i++) {
			$out->Write("echo Executing $lines[$i]\n");
			$out->Write("$lines[$i]\n");
			$App->print("$lines[$i]\n");
		}
	}
	$out->Close();

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

sub ShowScript
{
	$App->print("\n\n<b>Show WIEN2k Script:</b>\n");

	my $Task       = $Args->GetGetArg("Task");
	my $UseComplex = $Args->GetGetArg("UseComplex");
	my $Parallel   = $Args->GetGetArg("Parallel");
	my $SP         = $Args->GetGetArg("SpinPolarized");
	my $SO         = $Args->GetGetArg("SpinOrbit");
	my $Orb        = $Args->GetGetArg("OrbitalPotential");
	$UseComplex = 1 if($SP or $SO);
	$App->print("  Task: $Task\n");
	$App->print("  UseComplex      : $UseComplex\n");
	$App->print("  Parallel        : $Parallel\n");
	$App->print("  SpinPolarized   : $SP\n");
	$App->print("  SpinOrbit       : $SO\n");
	$App->print("  OrbitalPotential: $Orb\n");
	$App->print("\n");
	
	my $WIEN2k = new WIEN2k;
	my @lines = $WIEN2k->BuildScript($Task, $Parallel, $UseComplex, $SP, $SO, $Orb);
	for(my $i = 0 ; $i < @lines ; $i++) {
		$App->print("$lines[$i]\n");
	}
	$App->print("\n");
}

sub BackupStructFile
{
	$App->print("\n\n<b>Backup WIEN2k Struct File:</b>\n");

	my $SourceStructFile = $Args->GetGetArg(0);
	my $TargetStructFile = $Args->GetGetArg(1);
	unless($TargetStructFile) {
		$App->print("File names should be specified.\n");
		$App->print("   Usage: WIEN2k.pl --Action=BackupStruct "
			   ."SourceStructFile TargetStructFile\n");
		return 0;
	}

	$App->print("  Read $SourceStructFile.\n");
	$App->print("  Save to $TargetStructFile.\n");

	my $SourceWIEN2k = new WIEN2k;
	my $Crystal = $SourceWIEN2k->ReadStructFile($SourceStructFile, 0, tol => $TOL);
	unless($Crystal) {
		$App->print("Error: Can not read [$SourceStructFile].\n");
		return 0;
	}

	$Crystal->SetOutputMode('asymmetric');
	my $TargetWIEN2k = new WIEN2k;
	$TargetWIEN2k->SaveStructFile($Crystal, $TargetStructFile, 0, 0);

	$App->print("\nBackup WIEN2k Struct File: Finished\n");

	return 1;
}

sub MakeCIFFile
{
	$App->print("\n\n<b>Make CIF File from WIEN2k Struct File:</b>\n");

	my $StructFile = $Args->GetGetArg(0);
	my $CIFFile = $Args->GetGetArg(1);
	unless($CIFFile) {
		$App->print("File names should be specified.\n");
		$App->print("   Usage: WIEN2k.pl --Action=MakeCIF StructFile CIFFile\n");
		return 0;
	}

	$App->print("  Read $StructFile.\n");
	$App->print("  Save to $CIFFile.\n");

	my $WIEN2k = new WIEN2k;
	my $Crystal = $WIEN2k->ReadStructFile($StructFile, 0, tol => $TOL);
	unless($Crystal) {
		$App->print("Error: Can not read [$StructFile].\n");
		return 0;
	}

	my $CIF = new CIF;
	$CIF->CreateCIFFileFromCCrystal($Crystal, $CIFFile, 0, 'unix');

	$App->print("\nMake CIF File from WIEN2k Struct File: Finished\n");

	return 1;
}

sub MakeInMFile
{
	$App->print("\n\n<b>Make WIEN2k InM (Mini input) File:</b>\n");

	my $StructFile = $Args->GetGetArg(0);
	my $InMFile   = $Args->GetGetArg(1);
	unless($InMFile) {
		$App->print("File name should be specified.\n");
		$App->print("   Usage: WIEN2k.pl --Action=MakeInM StructFile InMFile\n");
		return 0;
	}

	my ($drive, $dir, $filename, $ext, $lastdir, $filebody)
		= Deps::SplitFilePath($StructFile);

	$App->print("  Read $StructFile.\n");
	$App->print("  Save to $InMFile.\n");

	my $WIEN2k = new WIEN2k;
	my $Crystal = $WIEN2k->ReadStructFile($StructFile, 0, tol => $TOL);
	unless($Crystal) {
		$App->print("Error: Can not read [$StructFile].\n");
		return 0;
	}

	my $out = new JFile($InMFile, "w");
	unless($out) {
		$App->print("Error: Can not write to [$InMFile].\n");
		return 0;
	}

	$App->print("=== $InMFile ===\n");
	$App->print("NEWT 2.0           # PORT/NEWT;  tolf (geom.opt. stops when F<tolf) (a4,f5.2))\n");
	$out->print("NEWT 2.0           # PORT/NEWT;  tolf (geom.opt. stops when F<tolf) (a4,f5.2))\n");
	my @AsymAtomSites = $Crystal->GetCAsymmetricAtomSiteList();
	for(my $i = 0 ; $i < @AsymAtomSites ; $i++) {
		my $site = $AsymAtomSites[$i];
		my @x = $site->Position();

		my @VariablePos;
		$VariablePos[0] = $Crystal->IsVariablePosition($x[0], $x[1], $x[2], 0);
		$VariablePos[1] = $Crystal->IsVariablePosition($x[0], $x[1], $x[2], 1);
		$VariablePos[2] = $Crystal->IsVariablePosition($x[0], $x[1], $x[2], 2);

		$App->print("$VariablePos[0] $VariablePos[1] $VariablePos[2]  1.0\n");
		$out->print("$VariablePos[0] $VariablePos[1] $VariablePos[2]  1.0\n");
	}
	$out->print("\n\n\n");
	$out->print("#PORT: delta1,2,3=precond., 4:Bondorder\n");
	$out->print("#NEWT: delta1,2,3=timestep, 4:friction (1=MOLD))\n");

	$out->Close();
	
	$App->print("\nMake InM file Finished.\n\n");
	return 1;
}


sub MakeInspFile
{
	$App->print("\n\n<b>Make WIEN2k Insp (spaghetti input) File:</b>\n");

	my $StructFile = $Args->GetGetArg(0);
	my $InspFile   = $Args->GetGetArg(1);
	unless($InspFile) {
		$App->print("File name should be specified.\n");
		return 0;
	}

	my ($drive, $dir, $filename, $ext, $lastdir, $filebody)
		= Deps::SplitFilePath($StructFile);
	my $Scf2File = $filebody . ".scf2";
	my $b = -f $Scf2File;
	$Scf2File = $filebody . ".scf2up" unless(-e $Scf2File);
	$Scf2File = $filebody . ".scf2dn" unless(-e $Scf2File);

	$App->print("  Read $StructFile.\n");
	$App->print("  Read $Scf2File.\n");
	$App->print("  Save to $InspFile.\n");

	my $in = new JFile($Scf2File, "r");
	unless($in) {
		$App->print("Error: Can not read [$Scf2File].\n");
		return 0;
	}
	my $EF;
	while(!$in->eof()) {
		my $line = $in->SkipTo(":FER ");
		last if($in->eof());
$App->DebugPrint("line: $line\n");
		($EF) = ($line =~ /=\s*([\d\+-\.Ee]+)/);
	}
	$in->Close();
$App->print("EF: $EF Ry\n");

	my $out = new JFile($InspFile, "w");
	unless($out) {
		$App->print("Error: Can not write to [$InspFile].\n");
		return 0;
	}

	my $outh = $out->HANDLE();
	print $outh <<EOT;
### Figure configuration
 5.0   3.0                         # paper offset of plot
10.0  15.0                         # xsize,ysize [cm]
 1.0   4                           # major ticks, minor ticks
 1.0   1                           # character height, font switch
 1.1   1    0                      # line width, line switch, color switch
### Data configuration            
-14.0  8.0  2                      # energy range, energy switch (1:Ry, 2:eV)
0      $EF                     # Fermi switch,  Fermi-level (in Ry units)
1   999                            # number of bands for heavier plotting   1,1
0      1    0.2                    # jatom, jtype, size  of heavier plotting   


Fermi switch:
  0...no line
  1...solid line
  2...dashed line
  3...dotted line

Line switch:
  0...dots
  1...lines
  2...lines and open circle
  3...lines and filled circles

Color switch (re-define your colors in defins.f)
  0...black
  1...one-color plot  
  2...three-color plot   
  3...multi-color plot 
  4...multi-color plot, one color for each irr. representations

Font switch:
  0...no text
  1...Times and Symbol
  2...Times, Symbol, and Times-Italic
  3...Helvetica, Symbol, and  Helvetica-Italic
  4...(include your own fonts in defins.f)
EOT
	$out->Close();
	
	$App->print("\nMake Insp file Finished.\n\n");
	return 1;
}

sub MakeIntFile
{
	$App->print("\n\n<b>Make WIEN2k Int File:</b>\n");

	my $StructFile = $Args->GetGetArg(0);
	my $IntFile    = $Args->GetGetArg(1);
	unless($IntFile) {
		$App->print("File name should be specified.\n");
		return 0;
	}

	my $TotalPDOSOnly       = $Args->GetGetArg("TotalPDOSOnly");
	my $OnceForEachAtomType = $Args->GetGetArg("OnceForEachAtomType");

	my ($drive, $dir, $filename, $ext, $lastdir, $filebody)
		= Deps::SplitFilePath($StructFile);
	my $QtlFile = $filebody . ".qtl";
	$QtlFile = $filebody . ".qtlup" if(!-e $QtlFile or -s $QtlFile < 10);
	$QtlFile = $filebody . ".qtldn" if(!-e $QtlFile or -s $QtlFile < 10);

	$App->print("  Read $StructFile.\n");
	$App->print("  Read $QtlFile.\n");
	$App->print("  Save to $IntFile.\n");

	my $WIEN2k = new WIEN2k;
	$WIEN2k->SetApplication($App);

	$App->print("  Reading [$StructFile]...\n");
	my $Crystal = $WIEN2k->ReadStructFile($StructFile, 0, tol => $TOL);
	my $Sample = $Crystal->CrystalName();

	$App->print("\n");
	$App->print("  Reading [$QtlFile]...\n");
	my @OrbitalSymmetries = $WIEN2k->ReadOrbitalSymmetriesFromQtlFile($QtlFile);

	$App->print("  Saving to [$IntFile]...\n");

	my $out = new JFile($IntFile, "w");
	unless($out) {
		$App->print("Error: Can not write to [$IntFile].\n");
		return 0;
	}
	$out->Write("$Sample\n");
	$out->Write(" -0.50 0.002 1.500 0.003   EMIN, DE, EMAX, Gauss-broadening(>de)\n");

	my $nSite = @OrbitalSymmetries;
	my @AsymAtomSites = $Crystal->GetCAsymmetricAtomSiteList();
	my $nOrb = 0;
	my $nPrintOrb = 0;
	my %PrintedAtomType;
	for(my $i = 0 ; $i < $nSite ; $i++) {
		my @a = Utils::Split(",", $OrbitalSymmetries[$i]);
		$nOrb += @a;

		my $site = $AsymAtomSites[$i];
		my $atomname = $site->AtomNameOnly();
		next if($PrintedAtomType{$atomname});
		my @a = Utils::Split(",", $OrbitalSymmetries[$i]);
		for(my $j = 0 ; $j < @a ; $j++) {
			$a[$j] =~ tr/A-Z/a-z/;
			$a[$j] = 's total' if($a[$j] =~ /^\s*0\s*$/);
			$a[$j] = 'p total' if($a[$j] =~ /^\s*1\s*$/);
			$a[$j] = 'd total' if($a[$j] =~ /^\s*2\s*$/);
			$a[$j] = 'f total' if($a[$j] =~ /^\s*3\s*$/);
			next if($TotalPDOSOnly and $a[$j] !~ /tot/i);

			$nPrintOrb++;
			$PrintedAtomType{$atomname}++;
		}
	}
	$nOrb += 2; # TotalとInterstitialを加える
	$nPrintOrb += 2;
	my $line = sprintf(" %4d  (Causion: must be equal to or less than 21: "
			  ."# OF DOS-CASES specified below\n", $nPrintOrb);
	$out->Write($line);
	$out->Write("    0    1   total         atom, case=column in qtl-header, label\n");
	%PrintedAtomType = {};
	for(my $i = 0 ; $i < $nSite ; $i++) {
		my $site = $AsymAtomSites[$i];
		my $atomname = $site->AtomNameOnly();
		next if($PrintedAtomType{$atomname});
		my @a = Utils::Split(",", $OrbitalSymmetries[$i]);
#$App->print("$i: $atomname: $OrbitalSymmetries[$i]\n");
		for(my $j = 0 ; $j < @a ; $j++) {
			$a[$j] =~ tr/A-Z/a-z/;
			$a[$j] = 's total' if($a[$j] =~ /^\s*0\s*$/);
			$a[$j] = 'p total' if($a[$j] =~ /^\s*1\s*$/);
			$a[$j] = 'd total' if($a[$j] =~ /^\s*2\s*$/);
			$a[$j] = 'f total' if($a[$j] =~ /^\s*3\s*$/);
			$line = sprintf(" %4d %4d   %s %s\n", $i+1, $j+1, $atomname, $a[$j]);
			next if($TotalPDOSOnly and $a[$j] !~ /tot/i);

			$out->Write($line);
			$PrintedAtomType{$atomname}++;
		}
	}
	$line = sprintf(" %4d    1   Interstitial\n", $nSite+1);
	$out->Write($line);
	
	$out->Close();

	$App->print("\n$nPrintOrb orbitals out of $nOrb were saved to [$IntFile].\n");
	$App->print("Make Int file Finished.\n\n");

	return 1;
}

sub MakeStructFile
{
	$App->print("\n\n<b>Make WIEN2k Struct File:</b>\n");

	my $DoSGROUP = $Args->GetGetArg("DoSGROUP");
	my $UseSameBasisForSameAtom = $Args->GetGetArg("UseSameBasisForSameAtom");
	my $CIFFile = $Args->GetGetArg(0);
	my $StructFile = $Args->GetGetArg(1);
	unless($StructFile) {
		$App->print("File name should be specified.\n");
		$App->print("  usage: WIEN2k.pl --Action=MakeStruct CIFFile StructFile.\n");
		return 0;
	}

#ファイル名を（ベース名, ディレクトリ名, 拡張子）に分解
	my ($drive, $dir, $filename, $ext, $lastdir, $filebody)
		= Deps::SplitFilePath($CIFFile);
	my $SampleName = $filebody;

	$App->print("  DoSGROUP: $DoSGROUP\n");
	$App->print("  UseSameBasisForSameAtom: $UseSameBasisForSameAtom\n");
	$App->print("  Read $CIFFile.\n");
	$App->print("  Save to $StructFile.\n");

	$App->print("Read [$CIFFile].\n");

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

#CIFクラスの内容から、Crystalクラスを作成
	$Crystal = $CIF->GetCCrystal();
	$Crystal->ExpandCoordinates();

	$App->print("Save to $StructFile.\n");
	my $WIEN = new WIEN2k;
	$WIEN->SetSampleName($SampleName);
	$WIEN->SaveStructFile($Crystal, $StructFile, 0, 0, $UseSameBasisForSameAtom);

	if($DoSGROUP) {
		$App->print("instgen\n");
		system("instgen");
		$App->print("$SGROUPPath -wi $StructFile | $TEEPath $StructFile.out\n");
		system("$SGROUPPath -wi $StructFile | tee $StructFile.out");
		$App->print("mv $StructFile $StructFile.original\n");
		system("mv $StructFile $StructFile.original");
		$App->print("$SGROUPPath -wi -wo $StructFile.original $StructFile\n");
		system("$SGROUPPath -wi -wo $StructFile.original $StructFile");
	}

	$App->print("\nMake Struct file Finished.\n\n");

	return;
}

sub MakeSymmetrizedStructFile
{
	$App->print("\n\n<b>Make Symmetrized WIEN2k Struct File:</b>\n");

	my $CIFFile = $Args->GetGetArg(0);
	my $StructFile = $Args->GetGetArg(1);
	unless($StructFile) {
		$App->print("File name should be specified.\n");
		$App->print("  usage: WIEN2k.pl --Action=MakeSymmetrizedStruct CIFFile StructFile.\n");
		return 0;
	}
	my $TOL = $Args->GetGetArg("TOL");
	$TOL = 0.0001 if(!defined $TOL or $TOL <= 0.0);
	$App->print("Torrelance: $TOL\n");
	my $Prim = $Args->GetGetArg("PrimitiveCell");

#ファイル名を（ベース名, ディレクトリ名, 拡張子）に分解
	my ($drive, $dir, $filename, $ext, $lastdir, $filebody)
		= Deps::SplitFilePath($CIFFile);
	my $SampleName = $filebody;
	my $BasePath = "$dir";
	$BasePath = "$drive$dir" if($drive ne '');

	my $StructPath        = Deps::MakePath($BasePath, "${filebody}.struct");
	my $SymStructPath     = Deps::MakePath($BasePath, "${filebody}-Symmetrized.struct");
	my $SymCIFPath = Deps::MakePath($BasePath, "${filebody}-Symmetrized.cif");
	$App->print("Original(non-symmetric) CIF   : $CIFFile\n");
	$App->print("Original(non-symmetric) struct: $StructPath\n");
	$App->print("Symmetrized struct            : $SymStructPath\n");
	$App->print("Symmetrized CIF               : $SymCIFPath\n");

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

#CIFクラスの内容から、Crystalクラスを作成
	$Crystal = $CIF->GetCCrystal();
	$Crystal->ExpandCoordinates();

	$App->print("Save to $StructPath.\n");
	my $WIEN = new WIEN2k;
	$WIEN->SetSampleName($SampleName);
	$WIEN->SaveStructFile($Crystal, $StructPath, 0, 0, 0);

#	$App->print("instgen\n");
#	system("instgen");
	my $SGROUPCommand = "$SGROUPPath -set-TOL=$TOL -wi $StructPath";
	$SGROUPCommand    = "$SGROUPPath -prim -set-TOL=$TOL -wi $StructPath" if($Prim);
	$App->print("$SGROUPCommand\n");
	system("$SGROUPCommand | $TEEPath $StructPath.out");

	$SGROUPCommand = "$SGROUPPath -set-TOL=$TOL -wi $StructPath -wo $SymStructPath";
	$SGROUPCommand = "$SGROUPPath -prim -set-TOL=$TOL -wi $StructPath -wo $SymStructPath" if($Prim);
	$App->print("$SGROUPCommand\n");
	system("$SGROUPCommand");

	my $WIEN2k2 = new WIEN2k;
	my $SymCrystal = $WIEN2k2->ReadStructFile($SymStructPath, 0, tol => $TOL);
	unless($SymCrystal) {
		$App->print("Error: Can not read [$SymStructPath].\n");
		return 0;
	}
	
	my $SPG = $SymCrystal->GetCSpaceGroup();
	my $SPGName = $SPG->SPGName();
	if($SPGName =~ /^R/i or $SPGName =~ /^P\s*\-?3/i) {
		$SPG->SetSPGName("$SPGName R");
		my ($a,$b,$c,$alpha,$beta,$gamma) = $SymCrystal->LatticeParameters();
		$App->print("  Hexagonal axis: a=$a   c=$c\n");
		my $ar = sqrt(3.0*$a*$a + $c*$c) / 3.0;
		my $alphah = 1.5 / sqrt(3.0 + ($c/$a)*($c/$a));
		$alphah = 2.0 * Sci::asin($alphah) * Sci::todeg();
		$App->print("  Rhombohedral axis: a=$ar   alpha=$alphah\n");
		$SymCrystal->SetLatticeParameters($ar,$ar,$ar,$alphah,$alphah,$alphah);
	}

	my $CIF2 = new CIF;
	if(!$CIF2->CreateCIFFileFromCCrystal($SymCrystal, $SymCIFPath, 0, 'unix')) {
		$App->print("Error: Can not create [$SymCIFPath].\n");
		return 0;
	}

	$App->print("\nMake Symmetrized Struct file Finished.\n\n");

	return;
}