#!/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 Utils;
use MyApplication;
use Crystal::WIEN2k;
use Crystal::XCrySDen;

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

#===============================================
# スクリプト大域変数
#===============================================
# for klist & bxsf
my $div = 360;
my ($nx, $ny, $nz)  = (20, 20, 40);
my $Target = "RHO"; # RHO/DIFF/OVER

#==========================================
# コマンドラインオプション読み込み
#==========================================
$App->AddArgument("--Action",    "--Action=[all]", '');
$App->AddArgument("--nx",        "--nx=num", '');
$App->AddArgument("--ny",        "--ny=num", '');
$App->AddArgument("--nz",        "--nz=num", '');
$App->AddArgument("--Target",    "--Target=[RHO]", '');
$App->AddArgument("--Action",    "--Action=[all]", '');
$App->AddArgument("--DebugMode", "--DebugMode: Set DebugMode", '');
exit 1 if($App->ReadArgs(0) != 1);
my %Args = $App->Args()->GetArgHash();

$App->SetDebug($Args{DebugMode}) if(defined $Args{DebugMode});

$nx     = $Args{nx} if($Args{nx} > 0);
$ny     = $Args{ny} if($Args{ny} > 0);
$nz     = $Args{nz} if($Args{nz} > 0);
$Target = $Args{Target} if($Args{Target} > 0);

#==========================================
# メイン関数スタート
#==========================================
my $Action  = $Args{Action};
my $DataDir = $App->Args()->GetGetArg(0);
$DataDir = "." if($DataDir eq '');
$DataDir = Utils::ReduceDirectory($DataDir);
my $StructFile = &GetStructPath($DataDir);
if(!-e $StructFile) {
	$App->print("Error: StructFile [$StructFile] could not obtained.\n");
	exit 0;
}
my $Title    = &GetTitle($StructFile);
my $In5Path  = Deps::ReplaceExtension($StructFile, ".in5");
my $RhoPath  = Deps::ReplaceExtension($StructFile, ".rho");
my $VASPPath = Deps::ReplaceExtension($StructFile, ".${Target}{spin}.vasp");
my $PrefPath = Deps::MakePath($DataDir, ".prefs-scf");

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

my $config = new IniFile();
$config->ReadAll($PrefPath);

$App->print("Struct: $StructFile\n");
$App->print("Lattice System: ", $Crystal->LatticeSystem(), "\n");
$App->print("Space Group   : ", $Crystal->SPGName(), "\n");
if($config->{spin} eq '') {
	$App->print("Spin polarized: no\n");
}
else {
	$App->print("Spin polarized: yes\n");
}

my $ret = 0;
if($Action eq 'all') {
	if($config->{spin} eq '') {
		&Execute($Crystal, $config, '',   $Title, $Target, $VASPPath, $In5Path, $RhoPath);
	}
	else {
		&Execute($Crystal, $config, 'up', $Title, $Target, $VASPPath, $In5Path, $RhoPath);
		&Execute($Crystal, $config, 'dn', $Title, $Target, $VASPPath, $In5Path, $RhoPath);
	}
}
else {
	$App->print("Error: Invald Action: [$Action]\n");
	$App->Usage();
}

exit $ret;

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

#==========================================
# &Subroutines
#==========================================
sub GetTitle
{
	my ($path) = @_;

	my ($drive, $directory, $filename, $ext, $lastdir, $filebody) = Deps::SplitFilePath($path);
	my $Title   = $filebody;
	$Title =~ s/\s/_/g;
	return $Title;
}

sub GetStructPath
{
	my ($DataDir) = @_;

	my ($drive, $directory, $filename, $ext, $lastdir, $filebody) = Deps::SplitFilePath($DataDir);
	my $StructFile = Deps::MakePath($DataDir, "$filename.struct", 0);
	if(!-e $StructFile) {
		my $fmask = Deps::MakePath($DataDir, "*.struct");
		my @files = glob($fmask);
		if($files[0] =~ /setrmt\.struct$/) {
			$StructFile = $files[1];
		}
		else {
			$StructFile = $files[0];
		}
		unless($StructFile) {
			print("Struct file is not found from Data Directory [$DataDir].\n");
			return undef;
		}
	}
	return $StructFile;
}

sub Execute
{
	my ($Crystal, $config, $spin, $Title, $Target, $VASPPath, $In5Path, $RhoPath) = @_;

	$App->print("\nMake Electron density .vasp file for spin=[$spin]\n");
	$App->print("VASP Path: $VASPPath\n");

	unlink($VASPPath);
	my $nx1 = $nx+1;
	my $ny1 = $ny+1;
	my $nz1 = $nz+1;
	my @Rho;
	for(my $iz = 0 ; $iz <= $nz ; $iz++) {
$App->print("iz=$iz / $nz\n");
		unlink($RhoPath);
		open(IN5, ">$In5Path") or die "Error: $!: Can not write to [$In5Path].\n";

		my $iz1 = int($iz * $div / $nz + 0.01);
		print IN5 <<EOT;
   0    0  $iz1  $div   # origin of plot: x, y,z, denominator
   0  $div  $iz1  $div   # y-end of plot
 $div    0  $iz1  $div   # x-end of plot
   3    3    3        # x, y, z nshells (of unit cells)
 $nx1  $ny1             # nx, ny
$Target                   # RHO/DIFF/OVER; ADD/SUB or blank
ANG VAL NODEBUG       # ANG/ATU, VAL/TOT, DEBUG/NODEBUG
NONORTHO              # optional line: ORTHO|NONORTHO
EOT
		close(IN5);
		
		my $cmd = ($spin eq '')? "x lapw5" : "x lapw5 -$spin";
		if(Utils::Execute($cmd)) {
			$App->print("Error: Could not execute [$cmd].\n");
			exit 0;
		}
		
		my $in = new JFile();
		if(!$in->Open($RhoPath, "r")) {
			$App->print("Error: Could not read [$RhoPath].\n");
			exit 0;
		}
		my $line = $in->ReadLine();
		while(1) {
			$line = $in->ReadLine();
			last if(!defined $line);

			push(@Rho, Utils::Split("\\s+", $line))
		}
		$in->Close();
	}
	unlink($RhoPath);

	my $OutFile = $VASPPath;
	$OutFile =~ s/\{spin\}/$spin/;
	my $out = new JFile();
	if(!$out->Open($OutFile, "w")) {
		$App->print("Error: Could not write to [$OutFile].\n");
		exit 0;
	}

	my ($a11, $a12, $a13, $a21, $a22, $a23, $a31, $a32, $a33) = $Crystal->LatticeVectors();
	$out->print("$Title\n");
	$out->printf("%6.4f\n", $a11);
	$out->printf("   %10.6f  %10.6f  %10.6f\n", 1.0,       $a12/$a11, $a13/$a11);
	$out->printf("   %10.6f  %10.6f  %10.6f\n", $a21/$a11, $a22/$a11, $a23/$a11);
	$out->printf("   %10.6f  %10.6f  %10.6f\n", $a31/$a11, $a32/$a11, $a33/$a11);
	$out->print("    1\n");
	$out->print("Direct\n");
	$out->print("  0.000000  0.000000  0.000000\n");
	$out->print("\n");
	$out->printf(" %d %d %d\n", $nx1, $ny1, $nz1);
	my $count = 0;
	for(my $iz = 0 ; $iz <= $nz ; $iz++) {
		for(my $iy = 0 ; $iy <= $ny ; $iy++) {
			for(my $ix = 0 ; $ix <= $nx ; $ix++) {
				$out->printf("%8.5f    ", $Rho[$count]);
				$count++;
				if($count % 10 == 0) {
					$out->print("\n");
				}
			}
		}
	}
	$out->Close();

}