#!/usr/bin/perl

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

use strict;
use Utils;
use JFile;
use CSV;

use MyHTMLApplication;
use Crystal::CIF;
use Crystal::Crystal;
use Crystal::VASP;
use Sci::Algorism;

my $Dir = $ARGV[0];
print "Dir: $Dir\n";

my @Components = ($Dir =~ /^([A-Z][a-z]?)([A-Z][a-z]?)([A-Z][a-z]?)([A-Z][a-z]?)/);
print "Components: ", join(', ', @Components), "\n";

my $App = new MyHTMLApplication;


my @Dirs = $App->CollectDirs($App->pParams(), $Dir, ["*DOS*", "*dos*"]);
for(my $i = 0 ; $i < @Dirs ; $i++) {
my $d = $Dirs[$i];
print("\n");
print("Search from [$d]\n");
my @DOSupFiles = glob(Deps::MakePath2($d, ["*-1.csv"], 0));
my @DOSdnFiles = glob(Deps::MakePath2($d, ["*-2.csv"], 0));
my $DOSup = $DOSupFiles[0];
my $DOSdn = $DOSdnFiles[0];
my ($DEFup, $DEFdn, $NVBup, $NVBdn, $Ntotup, $Ntotdn);
my ($EVBLowestup, $EVBLowestdn);
my ($pXup, $pYup, $pXdn, $pYdn, $iXVBLowestup, $iXVBLowestdn);
foreach my $file ($DOSup, $DOSdn) {
	print "  From [$file]\n";
	my $csv = new CSV;
	$csv->Read($file, 1);
	my $pLabelArray = $csv->LabelArray();
	my $pDataArray  = $csv->DataArray();
	my $nData = $csv->nData();
#print "  nData: $nData\n";
	my $pX = $pDataArray->[0];
	my $pY = $pDataArray->[1];
	last if(!defined $pY);
#print @$pLabelArray, "\n";
	my ($YMin, $YMax) = Utils::CalMinMax($pY);
	my $Yth = $YMax * 1.0e-3;
	print "  Y: $YMin - $YMax  Yth=$Yth\n";
	my $iXZero;
	for(my $i = 0 ; $i < $nData ; $i++) {
#print "$i: $pX->[$i]\n";
		if($pX->[$i] >= -1.0e-4) {
			$iXZero = $i;
			last;
		}
	}
	my $nzero = 0;
	my $iXVBLowest;
	for(my $i = $iXZero ; $i >= 0 ; $i--) {
		if($pY->[$i] < $Yth) {
			$nzero++;
			if($nzero > 10) {
				$iXVBLowest = $i;
				last;
			};
		}
		else {
			$nzero = 0;
		}
	}
	if($file eq $DOSup) {
		$EVBLowestup = $pX->[$iXVBLowest];
		$iXVBLowestup = $iXVBLowest;
		$pXup = $pX;
		$pYup = $pY;
	}
	else {
		$EVBLowestdn = $pX->[$iXVBLowest];
		$iXVBLowestdn = $iXVBLowest;
		$pXdn = $pX;
		$pYdn = $pY;
	}
}
my $iXVBLowest = ($iXVBLowestup < $iXVBLowestdn)? $iXVBLowestup : $iXVBLowestdn;
my $EVBLowest  = ($EVBLowestup < $EVBLowestdn)? $EVBLowestup : $EVBLowestdn;
print "EVBLowest: $EVBLowest [$iXVBLowest] (up=$EVBLowestup [$iXVBLowestup]  dn=$EVBLowestdn [$iXVBLowestdn])\n";

if($pYup and $pYdn) {
my $csv = new CSV;
my $pSup         = Algorism::IntegrateByConstantStepSimpson(scalar @$pXup, $pXup->[1] - $pXup->[0], $pYup);
my $pSVBLowestup = $pSup->[$iXVBLowest];
my $pSEFup       = $csv->YVal($pXup, $pSup, 0.0);
print "diff: $pSEFup - $pSVBLowestup\n";
my $NVBup        = $pSEFup - $pSVBLowestup;
my $DEFup        = $csv->YVal($pXup, $pYup, 0.0);

my $pSdn         = Algorism::IntegrateByConstantStepSimpson(scalar @$pXdn, $pXdn->[1] - $pXdn->[0], $pYdn);
my $pSVBLowestdn = $pSdn->[$iXVBLowest];
my $pSEFdn       = $csv->YVal($pXdn, $pSdn, 0.0);
my $NVBdn        = $pSEFdn - $pSVBLowestdn;
my $DEFdn        = $csv->YVal($pXdn, $pYdn, 0.0);

print "  Total electrons in VB: up=$NVBup e      dn=$NVBdn e\n";
print "  D(EF):                 up=$DEFup e/eV   dn=$DEFdn e/eV\n";

my $Ntottotal  = $Ntotup + $Ntotdn;
my $Ntotspin   = $Ntotup - $Ntotdn;
my $NVBtotal   = $NVBup + $NVBdn;
my $NVBspin    = $NVBup - $NVBdn;
my $DEFtotal   = $DEFup + $DEFdn;
my $DEFspin    = $DEFup - $DEFdn;
printf "  Ntot(total)=%9.6f  NVB(total)=%9.6f  DEF(total)=%9.6f\n", $Ntottotal, $NVBtotal, $DEFtotal;
printf "  Ntot(spin) =%9.6f  NVB(spin) =%9.6f  DEF(spin) =%9.6f\n", $Ntotspin, $NVBspin, $DEFspin;
}
}


my @Dirs = $App->CollectDirs($App->pParams(), $Dir, ["*SCF*", "*scf*"]);
for(my $i = 0 ; $i < @Dirs ; $i++) {
my $d = $Dirs[$i];
print("");
print("Search from [$d]\n");
my $INCAR = Deps::MakePath2($d, ["INCAR"], 0);
print "*From INCAR [$INCAR]\n";
my $in = new JFile($INCAR, "r");
my $SYSTEM = $in->SkipTo("SYSTEM\\s*=");
#print "SYSTEM: $SYSTEM\n";
my @AtomNames;
while(1) {
	my $line = $in->ReadLine();
	Utils::DelSpace($line);
	last if($line eq '');
	my ($AtomName) = ($line =~ /([^\\\/]+)[\\\/][^\\\/]*$/);
	last if(!defined $AtomName);
#print "a:$AtomName\n";
	push(@AtomNames, $AtomName);
}
print "Atoms: ", join(", ", @AtomNames), "\n";

my ($LDAUL, $LDAUU, $LDAUJ);
my @list = qw(PREC EDIFF EDIFFG ISPIN LDAU LDAUTYPE LDAUL LDAUU LDAUJ LMAXMIX ISIF PSTRESS ISMEAR SIGMA);
foreach my $key (@list) {
	my $val = $in->SkipTo("#?$key\\s*=");
	Utils::DelSpace($val);
	$LDAUL = $val if($key eq 'LDAUL');
	$LDAUU = $val if($key eq 'LDAUU');
	$LDAUJ = $val if($key eq 'LDAUJ');
	print "$val\n" if(defined $val and $val ne '');
}
if($LDAUL !~ /^#/) {
	my ($Lline) = ($LDAUL =~ /=\s*(.*)\s*$/);
	my @L = Utils::Split("\\s+", $Lline);
	my ($Uline) = ($LDAUU =~ /=\s*(.*)\s*$/);
	my @U = Utils::Split("\\s+", $Uline);
	my ($Jline) = ($LDAUL =~ /=\s*(.*)\s*$/);
	my @J = Utils::Split("\\s+", $Jline);
	for(my $i = 0 ; $i < @AtomNames ; $i++) {
		print "LDA: $AtomNames[$i]: l=$L[$i] U=$U[$i] J=$U[$i]\n";
	}
}
$in->Close();

my $OUTCAR = Deps::MakePath($d, "OUTCAR", 0);
print "\nFrom [$OUTCAR]\n";
my $in = new JFile($OUTCAR, "r");
my ($Magnetization, $FreeEnergy, $ExternalPressure);
while(1) {
	my $line = $in->SkipTo(".*number of electron\\s.*\\smagnetization\\s");
	last if(!defined $line);
	$Magnetization = $line;
	Utils::DelSpace($Magnetization);

	$line = $in->SkipTo("free energy\\s+TOTEN\\s*=");
	last if(!defined $line);
	$FreeEnergy = $line;
	Utils::DelSpace($FreeEnergy);

#	$line = $in->SkipTo("\\s*external pressure\\s*=");
#	last if(!defined $line);
#	$ExternalPressure = $line;
#	Utils::DelSpace($ExternalPressure);
}
print "$Magnetization\n";
print "$FreeEnergy\n";
#print "$ExternalPressure\n";

$in->rewind();
my $line = $in->SkipTo(".*TOTAL-FORCE");
print $line;
$in->ReadLine();
while(1) {
	$line = $in->ReadLine();
	last if(!defined $line or $line =~ /^\s*---/);
	$line =~ s/   / /g;
	print $line;
}

$line = $in->SkipTo(".*total charge");
print $line;
$in->ReadLine();
$line = $in->ReadLine();
print "ion   s     p     d      f     tot\n"; #$line;
$in->ReadLine();
while(1) {
	$line = $in->ReadLine();
	last if(!defined $line or $line =~ /^---/);
	my ($i) = ($line =~ /(\d+)/);
	$i = int(($i-0.99) / 2);
	$line =~ s/^\s*\d+\s*//g;
	$line =~ s/   / /g;
	print "$AtomNames[$i]: $line";
}

$line = $in->SkipTo(".*magnetization");
print $line;
$in->ReadLine();
$line = $in->ReadLine();
print "ion  s     p     d      f     tot\n"; #$line;
$in->ReadLine();
while(1) {
	$line = $in->ReadLine();
	last if(!defined $line or $line =~ /^---/);
	$line =~ s/   / /g;
	my ($i) = ($line =~ /(\d+)/);
	$i = int(($i-0.99) / 2);
	$line =~ s/^\s*\d+\s*//g;
	$line =~ s/   / /g;
	print "$AtomNames[$i]: $line";
}

$in->Close();
}


my @Dirs = $App->CollectDirs($App->pParams(), $Dir, ["*Relax*", "*relax*"]);
for(my $i = 0 ; $i < @Dirs ; $i++) {
my $d = $Dirs[$i];
my @ExpandedAtomSiteList;
my $nTotalExpandedAtomSite;
my @nMultiplicityExpandedAtomSiteList;
my $RelaxDir = Deps::MakePath($d, "*-final.cif", 0);
my @RelaxFiles = glob($RelaxDir);
#print @RelaxFiles, "\n";
my $RelaxFile = $RelaxFiles[0];
print "\nFrom [$RelaxFile]\n";
if($RelaxFile ne '') {
	my $cif = new CIF();
	$cif->Read($RelaxFile);
	my $Crystal = $cif->GetCCrystal();
	my $CrystalName = $Crystal->CrystalName();
print "Crystal: $CrystalName\n";
	my $Formula = $cif->Formula();
print "Formula: $Formula\n";
	my ($a,$b,$c,$alpha,$beta,$gamma) = $Crystal->LatticeParameters();
print "Lattice: $a $b $c A  $alpha $beta $gamma\n";
	$a = $b = ($a + $b) / 2.0;
	my $vol = $Crystal->Volume();
print "Volume: $vol A3\n";
	$Crystal->ExpandCoordinates();
	my $Density = $Crystal->Density();
print "Density: $Density g/cm3\n";
	@ExpandedAtomSiteList = $Crystal->GetCExpandedAtomSiteList();
	$nTotalExpandedAtomSite = $Crystal->nTotalExpandedAtomSite();
	@nMultiplicityExpandedAtomSiteList = $Crystal->GetCnMultiplicityExpandedAtomSiteList();
#print "nTotalExpandedAtomSite: $nTotalExpandedAtomSite\n";
	my ($zLn, $zPn) = (0.0, 0.0);
	for(my $i = 0 ; $i < $nTotalExpandedAtomSite ; $i++) {
		my $label     = $ExpandedAtomSiteList[$i]->Label();
		my $type      = $ExpandedAtomSiteList[$i]->AtomName();
		my ($x,$y,$z) = $ExpandedAtomSiteList[$i]->Position(1);
		my $occupancy = $ExpandedAtomSiteList[$i]->Occupancy();
		my $id = $ExpandedAtomSiteList[$i]->IdAsymmetricAtomSite();
		my $mult = $nMultiplicityExpandedAtomSiteList[$id-1];
		my $i1 = $i+1;
#print "   $i1: $label ($type): ($x, $y, $z) [$occupancy][$mult]\n";
		if($type eq $Components[0]) {
			if($z < 0.5) {
				$zLn += $z;
			}
			else {
				$zLn += 1.0 - $z;
			}
		}
		if($type eq $Components[3]) {
			if($z > 0.5) {
				$zPn += $z;
			}
			else {
				$zPn += 1.0 - $z;
			}
		}
	}
	$zLn /= 2.0;
	$zPn /= 2.0;
print "z($Components[0])=$zLn   z($Components[3])=$zPn\n";
	my $dLnO = $Crystal->GetInterAtomicDistance(0.75, 0.75, 1.0-$zLn, 0.75, 0.25, 1.0);
printf "d(Ln-O)=%8.5f A\n", $dLnO;
	my $dMPn = $Crystal->GetInterAtomicDistance(0.25, 0.25, $zPn, 0.25, 0.75, 0.5);
printf "d(M-Pn)=%8.5f A\n", $dMPn;
	my $dLnPn = $Crystal->GetInterAtomicDistance(1.75, 0.75, 1.0-$zLn, 1.25, 0.25, $zPn);
printf "d(Ln-Pn)=%8.5f A\n", $dLnPn;
	my $dLnM = $Crystal->GetInterAtomicDistance(1.25, 0.25, $zLn, 1.75, 0.25, 0.5);
printf "d(Ln-M)=%8.5f A\n", $dLnM;
	my $aPnMPnNormal = $Crystal->GetInterAtomicAngle(0.75, 0.25, 0.5, 0.25, 0.25, $zPn, 0.75, -0.25, 1.0-$zPn);
printf "a(Pn-M-Pn,normal to a-b plane)=%8.5f deg\n", $aPnMPnNormal;
	my $aPnMPnIn = $Crystal->GetInterAtomicAngle(0.75, 0.25, 0.5, 1.25, 0.25, $zPn, 0.75, -0.25, 1.0-$zPn);
printf "a(Pn-M-Pn,in a-b plane)=%8.5f deg\n", $aPnMPnIn;
	my $zPnPn = $c * (2.0 * $zPn - 1.0);
	my $xPnPn = $Crystal->GetInterAtomicDistance(0.25, 0.25, $zPn, 1.25, 0.25, $zPn);
printf "MPn4 tetrahedron: height=%8.5g A   width=%8.5g  oblateness=%8.5g\n",
			$zPnPn, $xPnPn, $zPnPn / $xPnPn / 0.707106781;
}
}