#!/usr/bin/perl
##!d:/Perl/bin/perl

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

use strict;
#use warnings;
#use CGI::Carp qw(fatalsToBrowser);
use Cwd;

use ProgVars;
use Deps;
use Utils;

use MyApplication;
use Sci qw($torad $eVToJ $RyToeV);
use Sci::FrozenPhonon;
use Sci::Optimize;

BEGIN {
#	$| = 1;
	open (STDERR, ">&STDOUT");
}

#===============================================
# Definition
#===============================================


#===============================================
# 大域変数
#===============================================
my $ConfigFile = 'configure.ini';

my $vars = new ProgVars;
my $OS         = Utils::GetOSName();
my $SGROUPPath = ProgVars::SGROUPPath();
my $TEEPath    = ProgVars::TEEPath();

my $ProgramDir     = ProgVars::ProgramDir();
my $TkPerlDir      = ProgVars::TkPerlDir();
my $ElasticDir     = Deps::MakePath($TkPerlDir,  ["ElasticConstant"], 0);
my $TemplateSource = Deps::MakePath($ElasticDir, ["Template"], 0);

my $BaseDir = cwd();
#my $TemplateDir = Deps::MakePath($BaseDir, "Template", 0);

my $TableFile      = "table.txt";
my $SummaryCSV     = "summary.csv";
my $InitScriptPath = 'initialize.bat';
my $RunScriptPath  = 'CalculateFrozenPhonon.bat';
my $DoVASPPath     = 'DoVASP.sh';
if($OS ne 'MSWin32') {
	$InitScriptPath = 'initialize.sh';
	$RunScriptPath  = 'CalculateFrozenPhonon.sh';
}
my $pParams = {};
$pParams = $vars->ReadConfigure($ConfigFile, $pParams, 1, 1);

#===============================================
# Applicationクラス作製
#===========================================
my $App = new MyApplication;
exit if($App->Initialize() < 0);

# 実行プログラム名（デフォルトでは$0から取得）
my $Program = $App->Program();
# アプリケーション名
$App->SetAppName($Program);
# バージョン
my $Version = $App->SetVersion("Ver 0.1");

my $fp = new FrozenPhonon;
$fp->SetApplication($App);

#===============================================
# コマンドラインオプション読み込み
#===============================================
$App->AddArgument("--Action",            "--Action=[MakeInputs|MakeRunScript|cal]: Default ''", '');
$App->AddArgument("--Mode",              "--Mode=[vib|init]: Default vib",              'vib');
$App->AddArgument("--iAtom",             "--iAtom=[index|AtomName|(x,y,z)]: Default 0",    0);
$App->AddArgument("--xdiff",             "--xdiff=val [in angstrom]: Default 0.01",     0.01);
$App->AddArgument("--nStep",             "--nStep=[integer >= 2]: Default 3",              3);
$App->AddArgument("--nOrder",            "--nOrder=[integer <= nStep-1]: Default 2",       2);
$App->AddArgument("--Engine",            "--Engine=[vasp|gulp]: Default ''",          '');
$App->AddArgument("--Library",           "--Library=[lib path]: Default ''",          '');
$App->AddArgument("--OS",                "--OS=[MSWin32|linux]: Default '$OS'",          $OS);


$App->AddArgument("--OriginalDir",       "--OriginalDir=path: Default original",       'original');
$App->AddArgument("--OriginalCIF",       "--OriginalCIF=path: Default original.cif",   'original.cif');
$App->AddArgument("--OriginalInput",     "--OriginalInput=path: Default original.cif", 'original.cif');
$App->AddArgument("--UseBravaisLattice", "--UseBravaisLattice=[0|1]: Default: 1",       1);
$App->AddArgument("--Index",             "--Index=[ijkl]: Default 1111",                  '1111');
$App->AddArgument("--nRelaxationStep",   "--nRelaxationStep=[integer]: Default 100",      100);
$App->AddArgument("--e_tol",             "--e_tol=val     EPS for Etot convergence  : Default 1.0", 0.5e-4);
$App->AddArgument("--F_tol",             "--F_tol=val     EPS for force convergence : Default 1.0", 0.5e-3);
$App->AddArgument("--conv_thr",          "--conv_thr=val  EPS for SCF convergence   : Default 1.0", 0.5e-6);
$App->AddArgument("--pressure_tol",      "--pressure_tol=val in atm. EPS for P convergence: Default 1.0", 1.0);
$App->AddArgument("--help",              "--help    : Show this help");
exit 1 if($App->ReadArgs() != 1);

$App->SetDefault();
$App->AddVars(%$pParams);

#===============================================
# 変数設定
#===============================================
my $iAtom  = $App->var('iAtom');
my $CIFFile = $App->GetGetArg(0);
   $CIFFile = 'ZnO.cif' if($CIFFile eq '');
my $SampleName = $CIFFile;
$SampleName =~ s/\..*?$//;   

#$App->AddVar('OS',                $OS);
$App->AddVar('CIFFile',           $App->GetGetArg(0));
$App->AddVar('SampleName',        $SampleName);
$App->AddVar('TemplateSource',    $TemplateSource);
my @files = $fp->GetFileConfig($TemplateSource, $App->var('Engine'), $App->var('Mode'));
$App->AddVar('RunScriptTemplate', $files[0]);
$App->AddVar('DoVASPTemplate',    $files[1]);
$App->AddVar('RunScriptPath',     $RunScriptPath);
$App->AddVar('DoVASPPath',        $DoVASPPath);
$App->AddVar('InitScriptPath',    $InitScriptPath);
$App->AddVar('RunScriptPath',     $RunScriptPath);
$App->AddVar('TableFile',         $TableFile);
$App->AddVar('SummaryCSV',        $SummaryCSV);

$fp->SetEngine($App->var('Engine'));
$fp->SetpEngine();
$fp->SetMode($App->var('Mode'));

$fp->print("\n");
$App->PrintArgs("%-30s: %s\n", 
	[qw(OS CIFFile Action Engine Mode 
		OriginalCIF OriginalInput OriginalDir
	    UseBravaisLattice Library
	    iAtom
	    xdiff nStep nOrder nRelaxationStep 
	    conv_thr e_tol F_tol pressure_tol), sort $App->ArgKeys()]);

#===============================================
# 実行
#===============================================
print "\n";
if($App->var('Action') =~ /^MakeInputs$/i) {
	$fp->MakeInputs();
}
elsif($App->var('Action') =~ /^MakeRunScript$/i) {
	$fp->MakeRunScript();
}
elsif($App->var('Action') =~ /^cal$/i) {
	$fp->Calculate();
}
else {
	$App->InvalidParameter('Action', 'Main', 1);
	exit 1;
}

exit;

#=========================================
# Subroutines
#=========================================

#export "WIEN2kplCMD=$PerlPATH $TkPerlDir/WIEN2k/WIEN2k.pl"
#MakeSymmetrizedCIFCMD="$WIEN2kplCMD --Action=MakeSymmetrizedStructFile --TOL=0.001 *-final.cif out.struct"

sub MakeSymmetrizedStructFile
{
	my ($App, $CIFile, $StructFile, $TOL, $Prim) = @_;
	$TOL  = 0.0001 if(!defined $TOL or $TOL <= 0.0);
	$Prim = 0 if(!defined $Prim);

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

#ファイル名を（ベース名, ディレクトリ名, 拡張子）に分解
	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();
	unless($CIF->Read($CIFFile)) {
		$App->print("Error: Can not read $CIFFile.\n\n");
		return 0;
	}
	my $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);
	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;
}