#!/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::ElasticConstant;
use Sci::Optimize;

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

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

# eiJ = e_i^j = d a_i / d a_j
# U = 1/2 * Cijkl * eIJ * eKL = 1/2 * C_{ijkl} * e^{ij} * e^{IJ}
#   = 1/2 * CIjKl * eiJ * ekL
# CIjKl = d2U / deiJ / dekL
# U = U0 + a11*dx2 + a22*dy2 + 2a12*dx*dy
# (dx, dy) = ( h,  0) => Up0 =  a11
# (dx, dy) = (-h,  0) => Um0 =  a11
# (dx, dy) = ( 0,  h) => U0p =  a22
# (dx, dy) = ( 0, -h) => U0m =  a22
# (dx, dy) = ( h,  h) => Upp = 2a12
# (dx, dy) = (-h, -h) => Umm = 2a12

#===============================================
# 大域変数
#===============================================
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  = 'CalculateElasticConstant.bat';

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 $ec = new ElasticConstant;
$ec->SetApplication($App);

#===============================================
# コマンドラインオプション読み込み
#===============================================
$App->AddArgument("--Action",            "--Action=[MakeInputs|MakeRunScript|cal]: Default ''", '');
$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("--BurstToP1",         "--BurstToP1=[0|1]: Default: 1",               1);
$App->AddArgument("--UseBravaisLattice", "--UseBravaisLattice=[0|1]: Default: 1",       1);
$App->AddArgument("--Mode",              "--Mode=[Cabc|CXYZ|Vconst|Vconp|init]: Default Vconst", 'Vconst');
$App->AddArgument("--Index",             "--Index=[ijkl]: Default 1111",                  '1111');
$App->AddArgument("--nStep",             "--nStep=[integer >= 2]: Default 3",             3);
$App->AddArgument("--nRelaxationStep",   "--nRelaxationStep=[integer]: Default 100",      100);
$App->AddArgument("--nOrder",            "--nOrder=[integer <= nStep-1]: Default 2",        2);
$App->AddArgument("--pressure0",         "--pressure0=val in atm. Origin of pressure: Default 1.0",           1.0);
$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("--fdiff",             "--fdiff=da/a=dV/V=val: Default 0.01",     0.01);
$App->AddArgument("--fdiff",             "--fdiff=da/a=dV/V=val: Default 0.01",     0.01);
$App->AddArgument("--Engine",            "--Engine=[vasp|gulp]: Default ''",          '');
$App->AddArgument("--Library",           "--Library=[lib path]: Default ''",          '');
$App->AddArgument("--help",              "--help    : Show this help");
exit 1 if($App->ReadArgs() != 1);

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

#===============================================
# 変数設定
#===============================================
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);
$App->AddVar('RunScriptTemplate', $ec->GetFileConfig($TemplateSource, $App->var('Engine'), $App->var('Mode')));
$App->AddVar('RunScriptPath',     $RunScriptPath);
$App->AddVar('InitScriptPath',    $InitScriptPath);
$App->AddVar('RunScriptPath',     $RunScriptPath);
$App->AddVar('TableFile',         $TableFile);
$App->AddVar('SummaryCSV',        $SummaryCSV);

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

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

#===============================================
# 実行
#===============================================
if($App->var('Action') =~ /^MakeInputs$/i) {
	$ec->MakeInputs();
}
elsif($App->var('Action') =~ /^MakeRunScript$/i) {
	$ec->MakeRunScript();
}
elsif($App->var('Action') =~ /^cal$/i) {
	$ec->CalculateElasticTensor();
}
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;
}