#!/usr/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 File::Basename;

use Sci qw($pi $e $e0 erfc erf);
use Sci::Algorism;

use Utils;

use ProgVars;
use MyApplication;

use Crystal::ForceField;
#use Crystal::LAMMPS;
#use Crystal::GULP;

#===============================================
# 大域変数
#===============================================
my $HOME = $ENV{'HOME'};
my $ProgramDir        = ProgVars::ProgramDir();
my $LAMMPSPerlDir     = ProgVars::LAMMPSPerlDir();
my $TemplateDir       = Deps::MakePath($LAMMPSPerlDir, 'Template', 0);
my $LAMMPSDatabaseDir = Deps::MakePath($LAMMPSPerlDir, 'Potentials', 0);
my $GULPDir           = ProgVars::GULPDir(); #Deps::MakePath($HOME, "gulp", 0);
my $GULPLibDir        = Deps::MakePath($GULPDir, "Libraries", 0);
my @PotentialLibs     = ForceField::PotentialLibs();

my $KC = $e / 4.0 / $pi / $e0 * 1.0e10;

#===============================
# Potential definition
#===============================
# Coulomb potential parameters
# $SplitEwaldSum: 1にするとPotential.csvにはCoulomb項を入れず、
#                 LAMMPS .inにewald項を入れる ($fCoulomb = 1)
my $SplitEwaldSum = 1;
my $alpha = 0.3;

# $fCoulomb: 最終的なポテンシャル中のCoulomb項の比率
my $fCoulomb = 1.0;

# Maximum connection radius allowed
my $rcmax = 1.4;

# R mesh for output
my $nr = 1500;
my ($rmin, $rstep) = (0.1, 0.01);
#my ($rmin, $rstep) = (0.01, 0.01);
my $rmax = $rmin + ($nr - 1) * $rstep;

# Numrical differentiation
my $h          = 0.0001; # delta r
my $nDiffOrder = 7;      # 2, 3, 5, 7, -1=Analytical

# Binary search criteria
my $nMaxIter = 100;
my $f2EPS = 1.0e-6;

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

#==========================================
# コマンドラインオプション読み込み
#==========================================
#$App->AddArgument("--Action",        "--Action=[ShowDatabases|FindDatabases|MakeInput|ConvCustomDump]",       '');
#$App->AddArgument("--Function",      "--Function=[Optimize|Fit|MD-NVT|MD-NPT] (Def:Optimize)", 'Optimize');
$App->AddArgument("--DebugMode", "--DebugMode: Set DebugMode", '');
exit 1 if($App->ReadArgs(0) != 1);
my $Args = $App->Args();

$App->SetAppName($ProgramName);
$App->SetUsage("$ProgramName [options] LibraryPath PotentialTablePath");

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

my $LibraryPath   = $Args->GetGetArg(0);
my $PotentialFile = $Args->GetGetArg(1);
my $DiffFile      = 'diff.csv';

if(!$PotentialFile) {
	print "\n";
	$App->Usage();
	exit;
}

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

my $ff = new ForceField;

print "\n";
print "Make LAMMPS connected potential file [$PotentialFile] for potential [$LibraryPath]\n";
print "\n";

print "r range: $rmin - $rmax, $rstep step for $nr mesh\n";
print "Differentiation: h=$h  nOrder=$nDiffOrder\n";
print "\n";

my $pPot = $ff->ReadPotentialLibrary($LibraryPath);
my $type = $pPot->{type};
$ff->PrintParameters();
print "  fCoulomb=$fCoulomb\n";

$ff->SavePotentials("Potential.csv", $rmin, $rstep, $nr);

my $name1 = 'O';
my $name2 = 'O';
my $pij = $pPot->{"$name1-core-$name2-core"};

#===
# U, 微分(1階〜3階)の計算
#===
my ($pr, $pf, $pf1, $pf2, $pf3);
($pr, $pf) = $ff->CalPotentialArray($name1, $name2, $rmin, $rstep, $nr);
($pf1, $pf2, $pf3) = $ff->CalDifferentials($pr, $pf, $nDiffOrder, $h);

#===
# U, 微分(1階〜3階)のをdiff.csvに保存
#===
open(OUT, ">$DiffFile") or die "Can not write to [$DiffFile].\n";
print OUT "i,r,f,f',f'',f'''\n";
for(my $k = 0 ; $k < $nr ; $k++) {
	printf OUT "%6d,%g,%g,%g,%g,%g\n", $k+1, $pr->[$k], $pf->[$k], $pf1->[$k], $pf2->[$k], $pf3->[$k];
}
close(OUT);

#===
# f'''=0となるrcを探す
#===
my ($rc, $fc, $f1c, $f2c, $f3c) = $ff->Find3rdDerivativeZero($nDiffOrder, $pr, $pf, $pf1, $pf2, $pf2);

#===
# rc > rcmaxのときはrc=rcmaxとし、rcmaxにおけるU,微分を計算しなおす
#===
if(defined $rc and $rc > $rcmax) {
	print "Warning: Calculated rc=$rc exceeds the allowed rcmax=$rcmax\n";
	$rc = $rcmax;
	$fc  = Algorism::Interpolate($pr, $pf, $rc);
	$f1c = Algorism::Interpolate($pr, $pf1, $rc);
	$f2c = Algorism::Interpolate($pr, $pf2, $rc);
	$f3c = Algorism::Interpolate($pr, $pf3, $rc);
	print "         Use rc=$rcmax (f=$fc  f'=$f1c  f''=$f2c  f'''=$f3c\n";
	print "\n";
}

#===
# 解析解で計算。うまくいかないので使わない
#===
#USR(r) = B / r^n + Dr^2
#my ($Bp, $Dp, $np, $Bm, $Dm, $nm) = $ff->CalBDnAnalytical($rc, $fc, $f1c, $f2c);

#===
# rcが見つかった場合のみ、Binary searchでdf''=0になるnを探す
#===
# Binary searchのnの初期範囲
my ($n0, $n1) = (0.01, 20.0);

my ($nm, $Bm, $Dm, $df2m);
if($rc) {
print "Binary search for df'' = 0\n";
my ($B0, $D0, $df20) = $ff->CalBD($rc, $fc, $f1c, $f2c, $n0, 0);
my ($B1, $D1, $df21) = $ff->CalBD($rc, $fc, $f1c, $f2c, $n1, 0);
print "Initial range:\n";
printf "  df''=%12.4g for n=$n0\n", $df20;
printf "  df''=%12.4g for n=$n1\n", $df21;
if($df20 * $df21 > 0.0) {
	print "Warning: Initial values must have different signs:\n";
	print "         rc is not found.\n";
	print "\n";
}

print "Binary search... (nMaxIter=$nMaxIter  EPS(f'')=$f2EPS)\n";
for(my $i = 0 ; $i < $nMaxIter ; $i++) {
	$nm = ($n0 + $n1) / 2.0;
	($Bm, $Dm, $df2m) = $ff->CalBD($rc, $fc, $f1c, $f2c, $nm, 0);
	printf "  df''=%12.4g for n=(n0+n1)/2=$nm\n", $df2m;

	if(abs($df2m) < $f2EPS) {
		print "** Final solution: n=$nm  B=$Bm  D=$Dm: df''=$df2m\n";
		last;
	}
	elsif($df2m * $df20 < 0.0) {
		$n1 = $nm;
		next;
	}
	elsif($df2m * $df21 < 0.0) {
		$n0 = $nm;
		next;
	}
}
print "\n";
}

#===
# Uを計算。$SplitEwaldSum=1の場合、逆空間和 $KC*Zi*Zj/r*erf(ar)を引く
#===
my @U;
for(my $k = 0 ; $k < $nr ; $k++) {
	my $r = $pr->[$k];
	if($r < $rc) {
		$U[$k] = $Bm / $r**$nm + $Dm * $r*$r;
	}
	else {
		$U[$k] = $ff->U($name1, $name2, $r);
	}
	$U[$k] -= $fCoulomb * $ff->Coulomb($name1, $name2, $r, 0.0) * erf($alpha * $r) if($SplitEwaldSum);
}

#===
# Potential.tableへ保存
# 微分で$nDiffOrder点公式を使った場合、最初と最後の$offset点ずつは計算できないので保存しない
#===
my $offset = int($nDiffOrder / 2);

print "Save potential table to [$PotentialFile]\n";
open(OUT, ">$PotentialFile") or die "Can not write to [$PotentialFile].\n";

print OUT "#\n";
print OUT "# Potential type: $pPot->{type}\n";
print     "# Potential type: $pPot->{type}\n";
print OUT "#  fCoulomb=$fCoulomb\n";
print     "#  fCoulomb=$fCoulomb\n";
if($SplitEwaldSum) {
	print OUT "#  SplitEwaldSum=$SplitEwaldSum (Short range Coulomb term included in this table)\n";
	print     "#  SplitEwaldSum=$SplitEwaldSum (Short range Coulomb term included in this table)\n";
	print OUT "#  Ewald parameter (G vector)=$alpha\n";
	print     "#  Ewald parameter (G vector)=$alpha\n";
}
else {
	print OUT "#  SplitEwaldSum=$SplitEwaldSum (Full Coulomb term included in this table)\n";
	print     "#  SplitEwaldSum=$SplitEwaldSum (Full Coulomb term included in this table)\n";
}
print OUT "#\n";
print OUT "#  connected at rc = $rc (df''=$df2m)\n";
if($pPot->{type} eq 'buckingham') {
	my $Aij = $pij->{A};
	my $bij = $pij->{rho};
	my $Cij = $pij->{C};
	print OUT "#  r > rc: $name1 - $name2: A=$Aij  rho=$bij  C=$Cij\n";
	print     "#  r > rc: $name1 - $name2: A=$Aij  rho=$bij  C=$Cij\n";
}
elsif($pPot->{type} eq 'morse') {
	my $De = $pij->{De};
	my $a0 = $pij->{a0};
	my $r0 = $pij->{r0};
	print OUT "#  r > rc: $name1 - $name2: De=$De  a0=$a0  r0=$r0\n";
	print     "#  r > rc: $name1 - $name2: De=$De  a0=$a0  r0=$r0\n";
}
print OUT "#  r < rc: n=$nm  B=$Bm  D=$Dm\n";
print OUT "# Differentiation: h=$h  nOrder=$nDiffOrder\n";
print     "# Differentiation: h=$h  nOrder=$nDiffOrder\n";
print OUT "#\n";
print OUT "\n";
print OUT "$name1-$name2\n";
printf OUT "N\t\t%d\tR\t%g\t%g\n", $nr - 6, $rmin + $rstep*3, $rmax - $rstep*3;
printf     "N\t\t%d\tR\t%g\t%g\n", $nr - 6, $rmin + $rstep*3, $rmax - $rstep*3;
print OUT "\n";
my $c = 1;
for(my $k = $offset ; $k < $nr-$offset ; $k++) {
	my $r = $pr->[$k];
	my $F;
	if($k < $offset or $k >= $nr-$offset) {
		$F = 0.0;
	}
	else {
		$F = -Algorism::Differentiate7WithIndex(\@U, $k, $rstep);
	}

	printf OUT "%6d\t\t%8.6f\t\t%12.6g\t\t%12.6g\n", $c, $r, $U[$k], $F;
#	printf     "%6d\t\t%8.6f\t\t%12.6g\t\t%12.6g\n", $c, $r, $U[$k], $F;
	$c++;
}
print OUT "\n";
close(OUT);

exit;


#===========================
# Subroutines
#===========================