#!/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 MyHTMLApplication;
use Crystal::VASP;
#use Sci::ChemicalReaction;
use Sci::Optimize;

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

#===============================================
# 大域変数
#===============================================
my $ScriptCharCode     = 'utf8';
my $FileSystemCharCode = 'sjis';

my $nLevel = 0;
my $KeyDir = 'SCF';
my $KeyDir = 'VCRelax';
   $KeyDir = 'VCRelax1' if(!-d $KeyDir);
my $h    = 0.1;
my $maxd = 0.3;
my $EEPS = 0.03;
my $dEPS = 0.001;

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

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

$App->SetOutputMode("console");

#==========================================
# コマンドラインオプション読み込み
#==========================================
$App->AddArgument("--Action", "--Action=[minimize]", '');
exit 1 if($App->ReadArgs(1, "sjis", 0) != 1);

my $Args = $App->Args();
#my $form = new CGI;
#$Args->SetCGIForm($form);
#$Args->parseInput($WebCharCode);

#$App->SetDebug($Debug);
$App->SetAppName('Relax.pl');
$App->SetUsage("perl Relax.pl --Action=[init|minimize|inf|summarize] InFile OutFile");

#my $Debug   = $Args->GetGetArg("DebugMode");
my $Action  = $Args->GetGetArg('Action');

print "Action : $Action\n";

#==========================================
# メイン関数スタート
#==========================================

#Utils::InitHTML("Research", $WebCharSet, "_self");
my $Action = $Args->GetGetArg("Action");
my $InFile  = $Args->GetGetArg(0);
my $OutFile = $Args->GetGetArg(1);
print "Action : $Action\n";
print "InFile : $InFile\n";
print "OutFile: $OutFile\n";

my $vasp = new VASP;

my $ret = 0;
if($Action =~ /init/i) {
#	$ret = &Initialize($InFile, $OutFile);
	$ret = &minimize($vasp);
}
elsif($Action =~ /minimize/i) {
	$ret = &minimize($vasp);
}
elsif($Action =~ /inf/i) {
	my @Dirs = &SearchDirs();
	my ($imin, $cmin, $Emin, $pr) = &GetInf($vasp, \@Dirs);
	$ret = 1;
}
elsif($Action =~ /summarize/i) {
	$ret = 1;
	my @Dirs = &SearchDirs();
	my ($imin, $cmin, $Emin, $pr) = &GetInf($vasp, \@Dirs);
	$ret = &SaveSummary('Relax.csv', $pr, 0);
}
else {
	$App->Usage();
	$ret = -1;
}
#Utils::EndHTML();

exit $ret;

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

#==========================================
# &Subroutines
#==========================================
sub Initialize
{
	my ($InFile, $OutFile) = @_;

	my $params = $Args->GetGetArg("Parameters");
	print "Parameters : $params\n";
	my ($a, $b, $c) = Utils::Split("\\s*,\\s*", $params);
	print "  a=$a  b=$b  c=$c\n";
	&UpdatePOSCAR($InFile, $OutFile, $a, $b, $c);
}

sub UpdatePOSCAR
{
	my ($InFile, $OutFile, $a, $b, $c) = @_;
	
	if(!defined $OutFile) {
		$App->Usage();
		return;
	}

	my $params = $Args->GetGetArg("Parameters");
	print "Parameters : $params\n";
	my ($a1, $b1, $c1) = Utils::Split("\\s*,\\s*", $params);
	$a = $a1 if(!defined $a);
	$b = $b1 if(!defined $b);
	$c = $c1 if(!defined $c);
	
	print "\n";
	print "Update POSCAR from [$InFile] to [$OutFile] for a=$a b=$b c=$c\n";
	
	my $in  = JFile->new($InFile, 'r')  or return -1;
	my $out = JFile->new($OutFile, 'w') or return -1;
#CsCl                                    
	$out->print($in->ReadLine());
#   1.00000000000000     
	my $line = $in->ReadLine();
	my $a0 = $line+0.0;
	$out->print("   1.00000000000000\n");
	my ($a11, $a12, $a13) = Utils::Split("\\s+", $in->ReadLine());
	my ($a21, $a22, $a23) = Utils::Split("\\s+", $in->ReadLine());
	my ($a31, $a32, $a33) = Utils::Split("\\s+", $in->ReadLine());
	$a11 *= $a0;
	$a12 *= $a0;
	$a13 *= $a0;
	$a21 *= $a0;
	$a22 *= $a0;
	$a23 *= $a0;
	$a31 *= $a0;
	$a32 *= $a0;
	$a33 *= $a0;
	$a11 = $a if(defined $a);
	$a22 = $b if(defined $b);
	$a33 = $c if(defined $c);
	$out->printf("   %20.16f   %20.16f   %20.16f\n", $a11, $a12, $a13);
	$out->printf("   %20.16f   %20.16f   %20.16f\n", $a21, $a22, $a23);
	$out->printf("   %20.16f   %20.16f   %20.16f\n", $a31, $a32, $a33);
	while(1) {
		my $line = $in->ReadLine();
		last if(!defined $line);
		$out->print($line);
	}
	$out->Close();
	$in->Close();
}

sub SearchDirs
{
	my @Dirs;
	@Dirs = Utils::SearchFilesRecursive2(".", "*", $nLevel, \@Dirs, 1,
			sub { 
				my ($path) = @_;
				my $KeyDir  = Deps::MakePath($path, $KeyDir, 0);
#				print "dir: $path [$KeyDir] Level=$nLevel\n"; 
				return 0 if(!-d $path);
				if(-d $KeyDir) {
					print "dir: $path [$KeyDir]\n"; 
					return 1;
				}
#				return 1;# if(-d $SCFDir);
				return 0;
			},
			SearchDir => 1,
#			Debug     => 1,
			);

	return @Dirs;
}

sub GetInf
{
	my ($vasp, $pDirs) = @_;

	my @r;
	print("\nInformations\n");
	my $c = 0;
	for(my $i = 0 ; $i < @$pDirs ; $i++) {
		my $RootDir = $pDirs->[$i];
		my $KeyDir  = Deps::MakePath($RootDir, $KeyDir, 0);

		next if(!-d $KeyDir);
		next if($RootDir =~ /junk/i);
#next if($i <= 1);
#next if($i > 0);

#		my $ret = &MakeVASPSummary($App, $out, $pDirs->[$i]);
		my ($drive, $directory, $filename, $ext1, $lastdir, $filebody) = Deps::SplitFilePath($RootDir);
		my $VCRelaxDir    = Deps::MakePath($RootDir, "VCRelax",  0);
		   $VCRelaxDir    = Deps::MakePath($RootDir, "VCRelax1",  0) if(!-d $VCRelaxDir);
		my $SCFDir        = Deps::MakePath($RootDir, "SCF",      0);
		$SCFDir = $VCRelaxDir if(!-d $SCFDir);
#		$App->print("\n$i: $RootDir\n");
		$App->print("\n$i SCF: $SCFDir\n");
		$App->print("\n$i VCRelax: $VCRelaxDir\n");

		my $POSCARHash    = $vasp->ReadPOSCARtoHash( Utils::MakePath($SCFDir,      'POSCAR',  '/', 0));
		my $VCROUTCARHash = $vasp->ReadOUTCARtoHash( Utils::MakePath($VCRelaxDir,  'OUTCAR',  '/', 0));
		my $SCFOUTCARHash = $vasp->ReadOUTCARtoHash( Utils::MakePath($SCFDir,      'OUTCAR',  '/', 0));
#		my $Crystal = $VASP->ReadStructureFromCARFiles($VCROUTCARHash);
#		unless($Crystal) {
#			$App->print("Error: Can not read from [$CARDir].\n");
#			return 0;
#		}

		$POSCARHash->{Latt} =~ s/\n/,/g;
		my @latt = Utils::Split(",+", $POSCARHash->{Latt});

		$App->print("Latt=$POSCARHash->{Latt}\n");
		$App->print("Volume=$POSCARHash->{Volume}\n");
		$App->print("Density=$POSCARHash->{Density}\n");
		$App->print("Etot=$VCROUTCARHash->{TOTEN}\n");
		$r[$c] = {
			Dir => $SCFDir,
			a => $latt[0],
			b => $latt[1],
			c => $latt[2],
			alpha  => $latt[3],
			beta  => $latt[4],
			gamma => $latt[0],
			E => $SCFOUTCARHash->{TOTEN},
			};
		$c++;
	}

	print "\n";
	print "Sorted\n";
	my $Emin = 1.0e10;
	my ($cmin, $imin, $Dirmin);
	@r = sort { $a->{c} <=> $b->{c} } @r;
	for(my $i = 0 ; $i < @r ; $i++) {
		print "  $i: a=$r[$i]->{a}  b=$r[$i]->{b}  c=$r[$i]->{c}  E=$r[$i]->{E}\n";
		if($Emin > $r[$i]->{E}) {
			$imin = $i;
			$cmin = $r[$i]->{c};
			$Emin = $r[$i]->{E};
			$Dirmin = $r[$i]->{Dir};
		}
	}
	print "Minimum: i=$imin  c=$cmin  E=$Emin  in [$Dirmin]\n";
	print "\n";

	return ($imin, $cmin, $Emin, \@r);
}

sub SaveSummary
{
	my ($fname, $pr, $IsPrint) = @_;
	$IsPrint = 1 if(!defined $IsPrint);

	print "\nSave summary to [$fname].\n";
	@$pr = sort { $a->{E} <=> $b->{E} } @$pr;

	my $params = $Args->GetGetArg("Parameters");
	print "Parameters : $params\n" if($IsPrint);
	my ($al, $bl, $cl) = Utils::Split("\\s*,\\s*", $params);
	print "  a=$al  b=$bl  c=$cl\n"; # if($IsPrint);
	
	my $out = JFile->new($fname, 'w') or return 0;
	$out->print("a,b,c,E,dir\n");
	print("a,b,c,E,dir\n");
	for(my $i = 0 ; $i < @$pr ; $i++) {
		$out->print("$al,$bl,$pr->[$i]{c},$pr->[$i]{E},$pr->[$i]{Dir}\n");
		print("$al,$bl,$pr->[$i]{c},$pr->[$i]{E},$pr->[$i]{Dir}\n");
	}
	$out->Close();

	return 1;
}

sub minimize
{
	my ($vasp) = @_;

	my $opt  = new Optimize;

	my @Dirs = &SearchDirs();
	my ($imin, $cmin, $Emin, $pr) = &GetInf($vasp, \@Dirs);
	&SaveSummary('Relax.csv', $pr);

	my @r = @$pr;

	if(@r == 0) {
		print "\n";
		print "Error: No valid E is available\n";
		$ret = &Initialize($InFile, $OutFile);
		return 4;
	}
	if(@r == 1) {
		print "\n";
		print "Warning: Only one E is available\n";
		my $x1 = $r[0]->{c} + $h;
		print "  Next c will be given by h=$h: x1 = $x1\n";
		&UpdatePOSCAR($InFile, $OutFile, $a, $b, $x1);
		return 2;
	}
	if(@r == 2) {
		print "\n";
		print "Warning: Only two E are available\n";
		my $dx = $r[1]->{c} - $r[0]->{c};
		$dx = $h if($dx < 1.0e-6);
		my $x1 = $r[1]->{c} + $dx;
		print "  Next c will be given by c0=$r[0]->{c} and c1=$r[1]->{c}: x1 = $x1\n";
		&UpdatePOSCAR($InFile, $OutFile, $a, $b, $x1);
		return 3;
	}

	my (@x, @y, $i0);
	if($imin < 1 or $imin >= @r - 1) {
		print "  Warning: Minimum not found (imin = $imin)\n";
		$i0 = ($imin == 0)? 0 : @r - 3;
	}
	else {
		print "Minimum at i=$imin, c=$cmin, E=$Emin\n";
		$i0 = $imin - 1;
	}

	print "\n";
	print "fit to E(c) = a0 + a1 * c + a2 * c^2\n";
	@x = ($r[$i0]->{c}, $r[$i0+1]->{c}, $r[$i0+2]->{c});
	@y = ($r[$i0]->{E}, $r[$i0+1]->{E}, $r[$i0+2]->{E});
	print "Use:\n";
	print "  ", $i0,   ": ($r[$i0]->{c}, $r[$i0]->{E})\n";
	print "  ", $i0+1, ": ($r[$i0+1]->{c}, $r[$i0+1]->{E})\n";
	print "  ", $i0+2, ": ($r[$i0+2]->{c}, $r[$i0+2]->{E})\n";

	print "\n";
	my $pc = $opt->mlsq(\@x, \@y, 2, 1);
	print "  E(c) = $pc->[0] + $pc->[1] * c + $pc->[2] * c^2\n";
#print "  c = $x[0]  y = ", $opt->mlsqYCal($x[0]), "\n";
	my $x1 = -$pc->[1] / 2.0 / $pc->[2];
	my $y1 = $opt->mlsqYCal($x1);
	print "\n";
	print "*** Minimum would be at c=$x1 (E=$y1)\n";
		print "\n";

	my $dd = $x1 - $cmin;
	my $dE = $y1 - $Emin;
	print "  dd = $dd (dEPS=$dEPS)  dE = $dE (EEPS=$EEPS)\n";
	if(abs($dd) < $dEPS and abs($dE) < $EEPS) {
		print "Convergency (dd=$dd < dEPS=$dEPS, dE=$dE < EEPS=$EEPS) reached\n";
		&UpdatePOSCAR($InFile, $OutFile, undef, undef, ($x1 + $cmin) / 2.0);
		return 0;
	}

	if(abs($dd) > $maxd) {
		print "  Displacement for x1=$x1 from cmin=$cmin ($dd) is larger than maxd=$maxd\n";
		if($dd < 0) {
			$x1 = $cmin - $maxd;
		}
		else {
			$x1 = $cmin + $maxd;
		}
		print "    Next c will be $x1\n";
		&UpdatePOSCAR($InFile, $OutFile, undef, undef, $x1);
		return 1;
	}
	print "\n";

	print "not converged (dd=$dd (dEPS=$dEPS)  dE=$dE (EEPS=$EEPS)).\n";
	print "    Next c will be $x1\n";
	&UpdatePOSCAR($InFile, $OutFile, undef, undef, $x1);

	return 1;
}