#!/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 Deps;
use Utils;

use MyApplication;
use Crystal::CIF;
use Crystal::Crystal;
use Sci qw($torad);

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

#===============================================
# 大域変数
#===============================================
# eiJを計算するための格子定数の歪率
my $fda = 0.01;

my $OriginalCIF         = 'original.cif';
my $StrainedCIFTemplate = 'e{i}{j}-{sign}.cif';
my $TableFile           = "table.txt";
my $SummaryCSV          = "summary.csv";

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

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

#===============================================
# コマンドラインオプション読み込み
#===============================================
$App->AddArgument("--Action", "--Action=[MakeInputs|run]: Default ''");
$App->AddArgument("--Engine", "--Engine=[vasp|gulp\]: Default ''");
$App->AddArgument("--help",   "--help    : Show this help");
exit 1 if($App->ReadArgs() != 1);

#===============================================
# 変数設定
#===============================================
my $Action  = $App->GetGetArg('Action');
   $Action  = '' if(!$Action);
my $Engine  = $App->GetGetArg('Engine');
   $Engine  = '' if(!$Engine);
my $CIFFile = $App->GetGetArg(0);
   $CIFFile = 'ZnO.cif' if($CIFFile eq '');

unless($CIFFile) {
	$App->Usage();
	exit 1;
}

$App->print("Action: $Action\n");
$App->print("Engine: $Engine\n");
my (@CopyList, $RunScriptPath);
if($Engine =~ /^gulp$/i) {
	@CopyList = ("GoGULP.bat", "runForDielectricConstant.bat");
	$RunScriptPath = "runForDielectricConstant.bat"
}
elsif($Engine =~ /^vasp$/i) {
	@CopyList = ("DoVASP.sh");
	$RunScriptPath = "bash DoVASP.sh"
}
$App->print("Copy files [", join(', ', @CopyList), "] from [$TemplateDir].\\n");
$App->print("Run script: $RunScriptPath\n");

#===============================================
# 実行
#===============================================
if($Action =~ /^MakeInputs$/i) {
	&MakeInputs();
}
elsif($Action =~ /^run$/i) {
	&run();
}
else {
	$App->print("\nError: Invalid Action [$Action].\n\n");
	$App->Usage();
	$App->print("\n");
	exit 1;
}

exit;

#=========================================
# Subroutines
#=========================================
sub ReadEtot
{
	my ($Engine, $file) = @_;

	my ($Etot);
	my $dir = $file;
	$dir =~ s/\..*?$//;
	if($Engine =~ /^gulp$/i) {
		my $OutFile = "$dir.out";
		$App->print("  Read Etot from [$OutFile].\n");
		my $in = JFile->new($OutFile, "r") or die "$!: Can not read [$OutFile].\n\n";
		my $line = $in->SkipTo("Total lattice enthalpy\\s*=");
		($Etot) = ($line =~ /=\s*([+\-\.\d]+)/);
	}
	elsif($Engine =~ /^vasp$/i) {
		my $OutFile = "SCF/OUTCAR";
		$App->print("  Read Etot from [$OutFile].\n");
		my $in = JFile->new($OutFile, "r") or die "$!: Can not read [$OutFile].\n\n";
		while(1) {
			my $line = $in->SkipTo("free\\s+energy\\s+TOTEN\\s*=");
			last if(!defined $line);
			($Etot) = ($line =~ /=\s*([+\-\.\d]+)/);
		}
	}
	else {
		$App->print("\nErrorin ReadEtot: Invalid Engine [$Engine].\n\n");
		exit;
	}
	return $Etot;
}

sub run
{
	my $table = JFile->new($TableFile, "r") or die "$!: Can not read [$TableFile].\n";
	my $out   = JFile->new($SummaryCSV, "w") or die "$!: Can not write to [$SummaryCSV].\n";
	$out->print("i,j,eiJ,Etot(eV),file\n");
	
	$App->print("\n");
	$App->print("Calculate total energies\n");
	$App->print("Base directory: $BaseDir\n");
	$App->print("\n");

	while(1) {
		my $line = $table->ReadLine();
		last if(!defined $line);
		
		my ($li, $lj, $eiJ, $file) = Utils::Split("\\s+", $line);
		$App->print(" for e$li$lj=$eiJ from [$file].\n");
		
		my $dir = $file;
		$dir =~ s/\..*?$//;
		$App->print("  Create dir [$dir].\n");
		Utils::CreateDirectory($dir);
#my $WorkDir = cwd();
#$App->print("  Working dir: [$WorkDir]\n");
#system("ls");
		$App->print("  Copy [$file, ", join(', ', @CopyList), "] to [$dir].\n");
		if(!Utils::CopyFile($file, $dir)) {
			$App->print("\nError: Can not copy [$file] to [$dir]\n\n");
			exit;
		}
		foreach my $f (@CopyList) {
			my $path = "$TemplateDir/$f";
			if(!Utils::CopyFile($path, $dir)) {
				$App->print("\nError: Can not copy [$path] to [$dir]\n\n");
				exit;
			}
		}

		$App->print("  cd [$dir]\n");
		if(!chdir($dir)) {
			$App->print("\nError: Can not change dir to [$dir]\n\n");
			exit;
		}
		my $WorkDir = cwd();
		$App->print("  Working dir: [$WorkDir]\n");
		my $cmd;
		if($Engine =~ /^gulp$/i) {
			$cmd = "$RunScriptPath $dir";
		}
		elsif($Engine =~ /^vasp$/i) {
			$cmd = "$RunScriptPath";
		}
		$App->print("  run [$cmd]\n");
#		system("chmod +x $RunScriptPath");
		system($cmd);

		my $Etot = &ReadEtot($Engine, $file);
		$App->print("  Etot = $Etot eV\n");
		$out->print("$li,$lj,$eiJ,$Etot,$file\n");

		$App->print("  cd [$BaseDir]\n");
		if(!chdir($BaseDir)) {
			$App->print("\nError: Can not change dir to [$dir]\n\n");
			exit;
		}

		$App->print("\n");
	}
	
	$out->Close();
	$table->Close();
}

sub MakeInputs
{
	$App->print("\n");
	$App->print("Make input files for calculating dielectric constant by finite method\n");
	$App->print("\n");

	$App->print("CIF File: $CIFFile\n");
	$App->print("\n");

#CIFファイル名を引数に与えて読み込む
	my $CIF = $App->{'CIF'} = new CIF;
	unless($CIF->Read($CIFFile, 0)) {
		$App->print("Error: Can not read [$CIFFile]\n\n");
		exit 2;
	}
	$App->print("\n");

	my $Crystal = $CIF->GetCCrystal();
#$Crystal->SplitPartialSites();
	$Crystal->ExpandCoordinates();
#pCrystal->m_ASFDCDBFile  = ASFDCPath;
#pCrystal->m_NonrelDBFile = NONRELPath;
#pCrystal->m_SPGDBFile    = SPGDBPath;
#pCrystal->SetFileName(InputFile);
#(pCrystal->SpaceGroup).SetNoSymmetry();
	my $CrystalName = $Crystal->CrystalName();
	$App->print("CrystalName: $CrystalName\n");
	my ($a,$b,$c,$alpha,$beta,$gamma) = $Crystal->LatticeParameters();
	$App->printf("cell: %9.6f %9.6f %9.6f   %9.6f %9.6f %9.6f \n", $a, $b, $c, $alpha, $beta, $gamma);
#my ($a11, $a12, $a13, $a21, $a22, $a23, $a31, $a32, $a33) = $Crystal->LatticeVectors();
#my ($g11, $g12, $g13, $g21, $g22, $g23, $g31, $g32, $g33) = $Crystal->Metrics();
#$my $vol = $Crystal->Volume();
#my $Density = $Crystal->Density();
#my ($Ra,$Rb,$Rc,$Ralpha,$Rbeta,$Rgamma) = $Crystal->ReciprocalLatticeParameters();
#my ($Ra11, $Ra12, $Ra13, $Ra21, $Ra22, $Ra23, $Ra31, $Ra32, $Ra33) = $Crystal->ReciprocalLatticeVectors();
#my ($Rg11, $Rg12, $Rg13, $Rg21, $Rg22, $Rg23, $Rg31, $Rg32, $Rg33) = $Crystal->RMetrics();
#my $Rvol = $Crystal->CalculateReciprocalVolume();
#my $SPG = $Crystal->GetCSpaceGroup();
#my $SPGName = $SPG->SPGName();
#my $iSPG    = $SPG->iSPG();
#my $LatticeSystem = $SPG->LatticeSystem();

	my @AtomTypeList     = $Crystal->GetCAtomTypeList();
	my $nAtomType        = @AtomTypeList;
	my @AtomAsymSiteList = $Crystal->GetCAsymmetricAtomSiteList();
	my $nAtomAsymSites   = @AtomAsymSiteList;
	my @AtomSiteList     = $Crystal->GetCExpandedAtomSiteList();
	my $nAtomSites       = $Crystal->nTotalExpandedAtomSite();

#===============================================
# eiJ計算
#===============================================

	my $table = JFile->new($TableFile, "w") or die "$!: Can not write to [$TableFile].\n";

	print "\n";
	print "Save original CIF to [$OriginalCIF]\n";
	my $NewCrystal = $Crystal->Duplicate();
	$NewCrystal->BurstToP1();
	&SaveCIF($NewCrystal, $OriginalCIF, 0, 'unix');
	$table->printf("%2d %2d %10.6f %s\n", -1, -1, 0.0, $OriginalCIF);

	for(my $i = 0 ; $i < 3 ; $i++) {
		for(my $j = $i ; $j < 3 ; $j++) {

			my $OutCIF = $StrainedCIFTemplate;
			$OutCIF =~ s/{i}/$i/g;
			$OutCIF =~ s/{j}/$j/g;
			$OutCIF =~ s/{sign}/p/g;
			print "Save strained CIF to [$OutCIF]\n";
			my ($eiJ) = &SaveStrainedStructure($Crystal, $i, $j, +$fda, $OutCIF);
			$table->printf("%2d %2d %10.6f %s\n", $i, $j, $eiJ, $OutCIF);
		
			$OutCIF = $StrainedCIFTemplate;
			$OutCIF =~ s/{i}/$i/g;
			$OutCIF =~ s/{j}/$j/g;
			$OutCIF =~ s/{sign}/m/g;
			print "Save strained CIF to [$OutCIF]\n";
			my ($eiJ) = &SaveStrainedStructure($Crystal, $i, $j, -$fda, $OutCIF);
			$table->printf("%2d %2d %10.6f %s\n", $i, $j, $eiJ, $OutCIF);
		}
	}

	$table->Close();
}

sub SaveCIF
{
	my ($Crystal, $OutCIF, $IsChooseRandomly, $strCRLF) = @_;

	my $IsPrint = 1;
	my $cif = new CIF;
	$cif->CreateCIFFileFromCCrystal($Crystal, $OutCIF, $IsChooseRandomly, $strCRLF, $IsPrint);
}

sub SaveStrainedStructure
{
	my ($Crystal, $li, $lj, $fda, $OutCIF) = @_;
	
	my $eiJ = $fda;

	my @v;
	my @latt = $Crystal->LatticeParameters();
	($v[0][0], $v[0][1], $v[0][2], $v[1][0], $v[1][1], $v[1][2], $v[2][0], $v[2][1], $v[2][2])
		= $Crystal->LatticeVectors();
	my $NewCrystal = $Crystal->Duplicate();
	$NewCrystal->BurstToP1();
#	$NewCrystal->SetSpaceGroup("P 1", 1, 1);
#	my @ExpandedAtomSiteList = $Crystal->GetCExpandedAtomSiteList();
#	$NewCrystal->SetCAsymmetricAtomSiteList(@ExpandedAtomSiteList);
#	$NewCrystal->ExpandCoordinates();

	$v[$li][0] += 0.5 * $fda * $v[$lj][0] * $latt[$li] / $latt[$lj];
	$v[$li][1] += 0.5 * $fda * $v[$lj][1] * $latt[$li] / $latt[$lj];
	$v[$li][2] += 0.5 * $fda * $v[$lj][2] * $latt[$li] / $latt[$lj];
	$v[$lj][0] += 0.5 * $fda * $v[$li][0] * $latt[$lj] / $latt[$li];
	$v[$lj][1] += 0.5 * $fda * $v[$li][1] * $latt[$lj] / $latt[$li];
	$v[$lj][2] += 0.5 * $fda * $v[$li][2] * $latt[$lj] / $latt[$li];

	$NewCrystal->SetLatticeVectors(
		$v[0][0], $v[0][1], $v[0][2], $v[1][0], $v[1][1], $v[1][2], $v[2][0], $v[2][1], $v[2][2]
		);

	&SaveCIF($NewCrystal, $OutCIF, 0, 'unix');

	return ($eiJ, $OutCIF);
}