#!/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 warnings;
use File::Path;
use File::Basename;
use File::Find;

use Deps;
use Utils;
use JFile;

use MyApplication;

use Crystal::CIF;
use Crystal::Crystal;
use Crystal::SpaceGroup;
use Crystal::AtomType;
use Crystal::VASP;

#===============================================
# デバッグ関係変数
#===============================================
#$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();

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

#$App->SetLF("<br>\n");
#$App->SetPrintCharCode("sjis");
#$App->SetDebug($Debug);
$App->SetDeleteHTMLFlag(1);

#===============================================
# スクリプト大域変数
#===============================================
my $RootDir     = ($ARGV[0])? $ARGV[0] : 'SCF'; #'.';

my $CIFFile     = 'structure.cif';
my $MaxDistance = 8.0;#3.4; # Angstrom
my $nMesh       = 100.0;

my $VASP = new VASP();
my $CHGCARPath     = "$RootDir/CHGCAR";
my $CoreCHGCARPath = '';
if(-e "$RootDir/AECCAR0" and -e "$RootDir/AECCAR2") {
	$CHGCARPath     = "$RootDir/AECCAR2";
	$CoreCHGCARPath = "$RootDir/AECCAR0";
}
#my $CSVPath     = $VASP->GetFileName($RootDir, 'DOSlm.csv');

my $pCHGHash = $VASP->ReadDensityFileToHash($CHGCARPath, $CoreCHGCARPath);
if(!$pCHGHash) {
	print "Error: Can not read [$CHGCARPath][$CoreCHGCARPath]\n";
	exit;
}
#	($hash{a11}, $hash{a12}, $hash{a13}) = Utils::Split("\\s+", $in->ReadLine());
#	($hash{a21}, $hash{a22}, $hash{a23}) = Utils::Split("\\s+", $in->ReadLine());
#	($hash{a31}, $hash{a32}, $hash{a33}) = Utils::Split("\\s+", $in->ReadLine());
#			$AtomPos[$ic] = [$x, $y, $z];
#	$hash{pAtomPos} = \@AtomPos;
#	$hash{pValues} = \@v;
print "Lattice Mode: $pCHGHash->{LatticeMode}\n";
my $pAtomPos = $pCHGHash->{pAtomPos};
print "nAtomPos: ", scalar @$pAtomPos, "\n";
for(my $i = 0 ; $i < @$pAtomPos ; $i++) {
	print "  i=$i: ($pAtomPos->[$i][0], $pAtomPos->[$i][1], $pAtomPos->[$i][2])\n";
}
print "n = ($pCHGHash->{nx}, $pCHGHash->{ny}, $pCHGHash->{nz})\n";
#for(my $ix = 0 ; $ix < $pCHGHash->{nx} ; $ix++) {
#	for(my $iy = 0 ; $iy < $pCHGHash->{ny} ; $iy++) {
#		for(my $iz = 0 ; $iz < $pCHGHash->{nz} ; $iz++) {
#			print "  i=($ix,$iy,$iz): $pCHGHash->{pValues}[$ix][$iy][$iz]\n";
#		}
#	}
#}

my $crystal = $VASP->ReadStructureFromCARFiles($RootDir);
unless($crystal) {
	print("Error: Can not read the crystal structure from [$RootDir].\n");
	exit;
}
my $Volume = $crystal->Volume();
print "Volume: $Volume A^3\n";
my $Density = $crystal->Density();
print "Density: $Density g/cm3\n";

my $TotalCharge = $pCHGHash->{Sum} / ($pCHGHash->{nx} * $pCHGHash->{ny} * $pCHGHash->{nz});
#my $TotalCharge = $pCHGHash->{Sum} * $Volume / $pCHGHash->{nx} / $pCHGHash->{ny} / $pCHGHash->{nz};
print "Sum: $pCHGHash->{Sum}\n";
print "Total Charge: $TotalCharge\n";
exit;

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

my @AtomTypeList           = $crystal->GetCAtomTypeList();
my $nAtomType              = @AtomTypeList;
my @ExpandedAtomSiteList   = $crystal->GetCExpandedAtomSiteList();
my $nTotalExpandedAtomSite = $crystal->nTotalExpandedAtomSite();
print "\n";
print "nTotalExpandedAtomSite: $nTotalExpandedAtomSite\n";
for(my $i = 0 ; $i < $nTotalExpandedAtomSite ; $i++) {
	my $atom      = $ExpandedAtomSiteList[$i];
	my $name      = $atom->AtomName();;
	my ($x,$y,$z) = $atom->Position(1);
	my $occ       = $atom->Occupancy();
	my $id        = $atom->Id();
	my $mult      = $atom->Multiplicity();
	
	my $dens = &GetDensity($x, $y, $z, $pCHGHash->{nx}, $pCHGHash->{ny}, $pCHGHash->{nz}, $pCHGHash->{pValues});
	print "  $name: ($x, $y, $z) dens = $dens\n";
}
my %nAtoms;
for(my $i = 0 ; $i < $nTotalExpandedAtomSite ; $i++) {
	my $atomsite = $ExpandedAtomSiteList[$i];
	my $atomname = $atomsite->AtomNameOnly();
	$nAtoms{$atomname}++;
}

my @Pair;
my $ip = 0;

print "\n";
print "Distances\n";
for(my $i = 0 ; $i < 1 ; $i++) {
#for(my $i = 0 ; $i < $nTotalExpandedAtomSite ; $i++) {
	my $atom      = $ExpandedAtomSiteList[$i];
	my $name      = $atom->AtomName();;
	my ($x,$y,$z) = $atom->Position(1);
	my $occ       = $atom->Occupancy();
	my $id        = $atom->Id();
	my $mult      = $atom->Multiplicity();

#	my $dens = &GetDensity($x, $y, $z, $pCHGHash->{nx}, $pCHGHash->{ny}, $pCHGHash->{nz}, $pCHGHash->{pValues});
#	print "$name: ($x, $y, $z) dens = $dens\n";

	for(my $j = $i+1 ; $j < $nTotalExpandedAtomSite ; $j++) {
#	for(my $j = 0 ; $j < $nTotalExpandedAtomSite ; $j++) {
		my $atom2        = $ExpandedAtomSiteList[$j];
		my $name2        = $atom2->AtomName();;
		my ($x2,$y2,$z2) = $atom2->Position(1);
		my $occ2         = $atom2->Occupancy();
		my $id2          = $atom2->Id();
		my $mult2        = $atom2->Multiplicity();
#	my $dens2 = &GetDensity($x2, $y2, $z2, $pCHGHash->{nx}, $pCHGHash->{ny}, $pCHGHash->{nz}, $pCHGHash->{pValues});
#	print "  $name: ($x2, $y2, $z2) dens = $dens2\n";

		for(my $ix = -1 ; $ix <= 1 ; $ix++) {
		for(my $iy = -1 ; $iy <= 1 ; $iy++) {
		for(my $iz = -1 ; $iz <= 1 ; $iz++) {
#			my $dis = $crystal->GetNearestInterAtomicDistance($x, $y, $z, $x2+$ix, $y2+$iy, $z2+$iz);
			my $dis = $crystal->GetInterAtomicDistance($x, $y, $z, $x2+$ix, $y2+$iy, $z2+$iz);
			next if($dis > $MaxDistance);

			print "  $name ($i) - $name2 ($j) ($ix,$iy,$iz): $dis\n";

			my $MinR     = 0.0;
			my $PrevDens = 1.0e100;
			print "distance from $name ($i)\tdensity\tx\ty\tz\n";
			for(my $im = 0 ; $im <= $nMesh ; $im++) {
				my $k = 1.0 * $im / $nMesh;
				my $xm = $x + $k * ($x2+$ix - $x);
				my $ym = $y + $k * ($y2+$iy - $y);
				my $zm = $z + $k * ($z2+$iz - $z);
				my $dens = &GetDensity($xm, $ym, $zm, $pCHGHash->{nx}, $pCHGHash->{ny}, $pCHGHash->{nz}, $pCHGHash->{pValues});
				my $dism = $dis * $k;
				print "$dism\t$dens\t$xm\t$ym\t$zm\n";
				if($PrevDens > $dens) {
					$MinR     = $dism;
					$PrevDens = $dens;
				}
			}
			print "\n";

			$Pair[$ip]{Name1}    = $name;
			$Pair[$ip]{Name2}    = $name2;
			$Pair[$ip]{Name}     = "$name ($i) - $name2 ($j)";
			$Pair[$ip]{Distance} = $dis;
			$Pair[$ip]{R1} = $MinR;
			$Pair[$ip]{R2} = $dis - $MinR;
			$ip++;
		}
		}
		}
	}
}

print "\n";
print "Suggested atomic radius\n";
my %R;
for(my $i = 0 ; $i < @Pair ; $i++) {
	print "$Pair[$i]{Name}:\n";
	print "  Distance = $Pair[$i]{Distance} = $Pair[$i]{R1} + $Pair[$i]{R2}\n";
	my $name1 = $Pair[$i]{Name1};
	$name1 =~ s/\s*\d*$//;
	my $name2 = $Pair[$i]{Name2};
	$name2 =~ s/\s*\d*$//;
	$R{$name1} = $Pair[$i]{R1} if(!defined $R{$name1} or $R{$name1} > $Pair[$i]{R1});
	$R{$name2} = $Pair[$i]{R2} if(!defined $R{$name2} or $R{$name2} > $Pair[$i]{R2});
}

print "\n";
print "Minimum radius\n";
foreach my $name (sort keys %R) {
	print "  *** R($name) = $R{$name} A\n";
}

print "\n";
print "Filling factor\n";
print("  Unit cell volume: $Volume A^3\n");
my $VolumeSum = 0.0;
for(my $i = 0 ; $i < $nAtomType ; $i++) {
	my $atomtype = $AtomTypeList[$i];
	my $atomname = $atomtype->AtomNameOnly();
	my $RWIG0    = $R{$atomname};
	my $n        = $nAtoms{$atomname};
	my $v        = 4.0/3.0 * Sci::pi() * $RWIG0*$RWIG0*$RWIG0;
	$VolumeSum += $v * $n;
$App->printf("  $atomname: r = %8.4f A: %10.4f x %d = %10.4f A^3\n", $RWIG0, $v, $n, $v*$n);
	}
my $f = $VolumeSum / $Volume * 100.0;
$App->printf("  Volume sum: %10.4f A^3\n", $VolumeSum);
$App->printf("  %6.2f %% of the Unit cell\n", $f);

print "\n";
print "For 80% Filling factor:\n";
my $kR = (80.0 / $f)**(1.0/3.0);
print "  kR = $kR\n";
for(my $i = 0 ; $i < $nAtomType ; $i++) {
	my $atomtype = $AtomTypeList[$i];
	my $atomname = $atomtype->AtomNameOnly();
	my $RWIG0    = $R{$atomname};
print "  $atomname: ", $kR * $RWIG0, " A\n";
	}


exit;


my ($nData, $pTE, $pDOS) = &MakeDOSCSVfromPRODOSCARs($RootDir);

exit;

sub GetDensity
{
	my ($x, $y, $z, $nx, $ny, $nz, $pV) = @_;

	return Algorism::InterpolateValuesInCube($nx, $ny, $nz, $pV, $x, $y, $z);




	$x = Utils::Reduce01($x);
	$y = Utils::Reduce01($y);
	$z = Utils::Reduce01($z);
	my $ix = int($x * $nx);
	my $iy = int($y * $ny);
	my $iz = int($z * $nz);
# d = a000 + a100 * x + a010 * y + a001 * z + a110 * xy + a101 * xz + a011 * yz + a111 * xyz
# d[0][0][0] = a000
# d[1][0][0] = a000 + a100
# d[0][1][0] = a000 + a010
# d[0][0][1] = a000 + a010
# d[1][1][0] = a000 + a100 + a010 + a110
# d[1][0][1] = a000 + a100 + a001 + a101
# d[0][1][1] = a000 + a010 + a001 + a011
# d[1][1][1] = a000 + a100 + a010 + a001 + a110 + a101 + a011 + a111
	my $a000 = $pV->[$ix][$iy][$iz];
	my $a100 = $pV->[$ix+1][$iy][$iz]     - $a000;
	my $a010 = $pV->[$ix][$iy+1][$iz]     - $a000;
	my $a001 = $pV->[$ix][$iy][$iz+1]     - $a000;
	my $a110 = $pV->[$ix+1][$iy+1][$iz]   - $a000 - $a100 - $a010;
	my $a101 = $pV->[$ix+1][$iy][$iz+1]   - $a000 - $a100 - $a001;
	my $a011 = $pV->[$ix][$iy+1][$iz+1]   - $a000 - $a010 - $a001;
	my $a111 = $pV->[$ix+1][$iy+1][$iz+1] - $a000 - $a100 - $a010 - $a001 - $a110 - $a101 - $a011;

#	my $dens = $pV->[$ix][$iy][$iz];
	my $dx = $x * $nx - $ix;
	my $dy = $y * $ny - $iy;
	my $dz = $z * $nz - $iz;
	my $dens = $a000 + $a100 * $dx + $a010 * $dy + $a001 * $dz
	         + $a110 * $x*$y + $a101 * $x*$z + $a011 * $y*$z + $a111 * $x*$y*$z;

	return $dens;
}

sub MakeDOSCSVfromPRODOSCARs
{
	my ($RootDir) = @_;

	$App->print("\n<b>Make lm-decomposed DOS CSV file from PROCAR/DOSCAR:</b>\n");

	my $VASP = new VASP();
	my $CSVPath     = $VASP->GetFileName($RootDir, 'DOSlm.csv');
	$App->print("CSV           : $CSVPath\n");

#	my $Width = $Args->GetGetArg("Width");
#	$Width = 0 if(!defined $Width);
#	$App->print("Width: $Width eV\n");
#	my $IgnoreZero = $Args->GetGetArg("IgnoreZero");
#	$IgnoreZero = 0 if(!defined $IgnoreZero);
#	$App->print("$IgnoreZero: $IgnoreZero\n");

	&ReadDOSCAR($RootDir);

	$App->print("\n<b>Make lm-decomposed DOS CSV file from PROCAR/DOSCAR: finished</b>\n");
}

sub ReadDOSCAR
{
	my ($RootDir) = @_;

	my $VASP = new VASP();

	my $PROCARPath  = $VASP->GetFileName($RootDir, 'PROCAR');
print "ReadDOSCAR: Read [$PROCARPath].\n";
	my @Orbitals = &GetOrbitalsFromPROCAR($PROCARPath);
	if(@Orbitals == 0) {
		return -1;
	}

	my $POSCAR = $VASP->GetPOSCARFileName($RootDir);
print "ReadDOSCAR: Read [$POSCAR].\n";
	my $Crystal = $VASP->ReadStructureFromCARFiles($POSCAR, 1);

	my @AtomTypeList = $Crystal->GetCAtomTypeList();
	my $nAtomType = @AtomTypeList;
	my %AtomNum;
	for(my $i = 0 ; $i < $nAtomType ; $i++) {
		my $name = $AtomTypeList[$i]->AtomNameOnly();;
		$AtomNum{$name} = 0 if(!defined $AtomNum{$name});
		my $n = ++$AtomNum{$name};
		$AtomTypeList[$i]->SetAtomName("$name$n");
#print "$i: $name$n\n";
	}
	my @ExpandedAtomSiteList = $Crystal->GetCExpandedAtomSiteListByOutputMode();
#	my @M = $Crystal->GetCnMultiplicityExpandedAtomSiteList();
	my $nTotalExpandedAtomSite = @ExpandedAtomSiteList;
	my @AtomNames;
	for(my $i = 0 ; $i < $nTotalExpandedAtomSite ; $i++) {
		my $atomsite = $ExpandedAtomSiteList[$i];
		$AtomNames[$i] = $atomsite->AtomName(); #AtomNameOnly();
#print "Site[$i]: $AtomNames[$i] ($M[$i])\n";
print "Site[$i]: $AtomNames[$i]\n";
	}

	my $OUTCAR = $VASP->GetOUTCARFileName($RootDir);
	my $EF     = $VASP->ReadFermiEnergy($OUTCAR);
print "EF: $EF eV\n";

	my $DOSCARPath  = $VASP->GetFileName($RootDir, 'DOSCAR');
print "ReadDOSCAR: Read [$DOSCARPath].\n";

	my $in = new JFile($DOSCARPath, 'r');
	if(!$in) {
		print "Error in ReadDOSCAR: $!: Can not read [$DOSCARPath].\n";
		return -1;
	}

	my $line;
	for(my $i = 0 ; $i < 5 ; $i++) {
		$line = $in->ReadLine();
	}
# BaZn2As2-I4mmm [PAW_PBE] (dos)          
	my $SampleName = $line;
	$SampleName =~ s/^\s*//;
	$SampleName =~ s/\s*$//;
print "Sample [$SampleName]\n";

#     15.00000000    -35.00000000 5000      3.50770154      1.00000000
	$line = $in->ReadLine();
	my ($E1, $E0, $nData, $EF, $w) = Utils::Split("\\s+", $line);
print "nData=$nData\n";

	my (@TE, %DOS);
	for(my $i = 0 ; $i < $nData ; $i++) {
		$line = $in->ReadLine();
		my ($e, $tot, $ne) = Utils::Split("\\s+", $line);
		$TE[$i]              = $e;
		$DOS{Total}[$i]        += $tot;
		$DOS{Interstitial}[$i] += $tot;
		$DOS{NE}[$i]           += $ne;
#print "  $i: $e\t$up\t$dn\n";
	}

	my $iAtom = 0;
	for(my $i = 0 ; $i < @ExpandedAtomSiteList ; $i++) {
#	while(1) {
		$line = $in->ReadLine();
		my ($E1, $E0, $nData, $EF, $w) = Utils::Split("\\s+", $line);
		my $atomsite = $ExpandedAtomSiteList[$iAtom];
print "iAtom=$iAtom\n";
		my $AtomName = $atomsite->AtomNameOnly();
print "Site[$iAtom]: $AtomName: nData=$nData\n";

		for(my $i = 0 ; $i < $nData ; $i++) {
			$line = $in->ReadLine();
			my ($e, @a) = Utils::Split("\\s+", $line);
			last if(!$line or @a <= 1);
if($i == 0) {
	Utils::DelSpace($line);
	my $nOrb = scalar @a;
	print "  nO=$nOrb [$line]\n";
}
			for(my $j = 0 ; $j < @a ; $j++) {
				my $orb = $Orbitals[$j];
#print("$j: $orb  ");
				$DOS{Interstitial}[$i]     -= $a[$j];
				$DOS{$AtomName}{Total}[$i] += $a[$j];
				$DOS{$AtomName}{$orb}[$i]  += $a[$j];
			}
		}
		last if($in->eof());

		$iAtom++;
	}

	$in->Close();

	my $CSVPath  = $VASP->GetFileName($RootDir, 'DOSlm.csv');
print "ReadDOSCAR: Write to [$CSVPath].\n";
	my $csv = new JFile($CSVPath, 'w');
	if(!$csv) {
		print "Error in MakeDOSCAR:: $!: Can not write to [$CSVPath].\n";
		return -2;
	}

	$csv->print("E(eV),Total,Interstitial,NE");
	my $nAtomType = @AtomTypeList;
	for(my $i = 0 ; $i < $nAtomType ; $i++) {
		my $atom     = $AtomTypeList[$i];
		my $atomname = $atom->AtomNameOnly();
		$csv->print(",$atomname(Total)");
		for(my $j = 0 ; $j < @Orbitals ; $j++) {
			my $orb = $Orbitals[$j];
			$csv->print(",$atomname($orb)");
		}
	}
	$csv->print("\n");
	
	for(my $i = 0 ; $i < $nData ; $i++) {
		$csv->print("$TE[$i],$DOS{Total}[$i],$DOS{Interstitial}[$i],$DOS{TotalNE}[$i]");
		for(my $it = 0 ; $it < $nAtomType ; $it++) {
			my $atom     = $AtomTypeList[$it];
			my $atomname = $atom->AtomNameOnly();
			$csv->print(",$DOS{$atomname}{Total}[$i]");
			for(my $j = 0 ; $j < @Orbitals ; $j++) {
				my $orb = $Orbitals[$j];
				$csv->print(",$DOS{$atomname}{$orb}[$i]");
			}
		}
		$csv->print("\n");
	}

	$csv->Close();
	
	return ($nData, \@TE, \%DOS);
}


sub GetOrbitalsFromPROCAR
{
	my ($path) = @_;

	my $in = new JFile($path, 'r');
	if(!$in) {
		print "Error in GetOrbitalsFromPROCAR: $!: Can not read [$path].\n";
		return -1;
	}
	if(!$in->SkipTo("band\\s+\\d")) {
		print "Error in GetOrbitalsFromPROCAR: Can not find 'band' label in [$path].\n";
		return -2;
	}
	$in->ReadLine();
	my $line = $in->ReadLine();
	my ($ion, @a) = Utils::Split("\\s+", $line);
	my $na = @a;
	if($a[$na-1] eq 'tot') {
		delete $a[$na-1];
	}
	
	$in->Close();
	
	return @a;
}