#!/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", "c:/Programs/Perl/lib", "d:/Programs/Perl/lib", @INC);
}

use strict;
#use warnings;

use Sci::Algorism;

#===============================
# Fundamental parameters
#===============================
our $pi          = 3.14159265358979323846;
our $e           = 1.60218e-19;      # C";
our $e0          = 8.854418782e-12;  # C<sup>2</sup>N<sup>-1</sup>m<sup>-2</sup>";

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


#===============================
# File configuration
#===============================
my $OutFile = 'Potential.table';

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

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


#===============================
# Potential definition
#===============================
my @atoms = ('O', 'Zn', 'Ga', 'In');
my $type = 'buckingham';

# Coulomb potential parameters

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

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

my %Zi;
$Zi{In} =  3.0; #1.8; #3.0;
$Zi{Ga} =  3.0; #1.8; #3.0;
$Zi{Zn} =  2.0; #1.2; #2.0;
$Zi{O}  = -2.0; #-1.2; #-2.0;

# Born repulsion parameters
my (%Aij, %rhoij);

$Aij{In}{In} = $Aij{In}{In} =    0.0;
$Aij{In}{Ga} = $Aij{Ga}{In} =    0.0;
$Aij{In}{Zn} = $Aij{Zn}{In} =    0.0;
$Aij{In}{O}  = $Aij{O}{In}  = 1293.600;
$Aij{Ga}{Ga} = $Aij{Ga}{Ga} =    0.0;
$Aij{Ga}{Zn} = $Aij{Zn}{Ga} =    0.0;
$Aij{Ga}{O}  = $Aij{O}{Ga}  = 2339.766;
$Aij{Zn}{Zn} = $Aij{Zn}{Zn} =    0.0;
$Aij{Zn}{O}  = $Aij{O}{Zn}  =  600.300;
$Aij{O}{O}   = $Aij{O}{O}   =   25.410;

$rhoij{In}{In} = $rhoij{In}{In} = 1.0;
$rhoij{In}{Ga} = $rhoij{Ga}{In} = 1.0;
$rhoij{In}{Zn} = $rhoij{Zn}{In} = 1.0;
$rhoij{In}{O}  = $rhoij{O}{In}  = 0.3310;
$rhoij{Ga}{Ga} = $rhoij{Ga}{Ga} = 1.0;
$rhoij{Ga}{Zn} = $rhoij{Zn}{Ga} = 1.0;
$rhoij{Ga}{O}  = $rhoij{O}{Ga}  = 0.2740;
$rhoij{Zn}{Zn} = $rhoij{Zn}{Zn} = 1.0;
$rhoij{Zn}{O}  = $rhoij{O}{Zn}  = 0.3370;
$rhoij{O}{O}   = $rhoij{O}{O}   = 0.6940;

# Dispersion parameters
my %Cij;

$Cij{In}{In} = $Cij{In}{In} = 0.0;
$Cij{In}{Ga} = $Cij{Ga}{In} = 0.0;
$Cij{In}{Zn} = $Cij{Zn}{In} = 0.0;
$Cij{In}{O}  = $Cij{O}{In}  = 4.325;
$Cij{Ga}{Ga} = $Cij{Ga}{Ga} = 0.0;
$Cij{Ga}{Zn} = $Cij{Zn}{Ga} = 0.0;
$Cij{Ga}{O}  = $Cij{O}{Ga}  = 0.0;
$Cij{Zn}{Zn} = $Cij{Zn}{Zn} = 0.0;
$Cij{Zn}{O}  = $Cij{O}{Zn}  = 0.0;
$Cij{O}{O}   = $Cij{O}{O}   = 32.32;

# Morse potential parameters
my %CR12;
$CR12{In}{In} = 0.0;

my (%Deij, %a0ij, %r0ij);

$Deij{In}{In} = $Deij{In}{In} = 0.0;
$Deij{In}{Ga} = $Deij{Ga}{In} = 0.0;
$Deij{In}{Zn} = $Deij{Zn}{In} = 0.0;
$Deij{In}{O}  = $Deij{O}{In}  = 0.072975;
$Deij{Ga}{Ga} = $Deij{Ga}{Ga} = 0.0;
$Deij{Ga}{Zn} = $Deij{Zn}{Ga} = 0.0;
$Deij{Ga}{O}  = $Deij{O}{Ga}  = 0.038263;
$Deij{Zn}{Zn} = $Deij{Zn}{Zn} = 0.0;
$Deij{Zn}{O}  = $Deij{O}{Zn}  = 0.001221;
$Deij{O}{O}   = $Deij{O}{O}   = 0.042395;

$a0ij{In}{In} = $a0ij{In}{In} = 0.0;
$a0ij{In}{Ga} = $a0ij{Ga}{In} = 0.0;
$a0ij{In}{Zn} = $a0ij{Zn}{In} = 0.0;
$a0ij{In}{O}  = $a0ij{O}{In}  = 1.9342;
$a0ij{Ga}{Ga} = $a0ij{Ga}{Ga} = 0.0;
$a0ij{Ga}{Zn} = $a0ij{Zn}{Ga} = 0.0;
$a0ij{Ga}{O}  = $a0ij{O}{Ga}  = 2.1052;
$a0ij{Zn}{Zn} = $a0ij{Zn}{Zn} = 0.0;
$a0ij{Zn}{O}  = $a0ij{O}{Zn}  = 2.347100;
$a0ij{O}{O}   = $a0ij{O}{O}   = 1.379316;

$r0ij{In}{In} = $r0ij{In}{In} = 0.0;
$r0ij{In}{Ga} = $r0ij{Ga}{In} = 0.0;
$r0ij{In}{Zn} = $r0ij{Zn}{In} = 0.0;
$r0ij{In}{O}  = $r0ij{O}{In}  = 2.82499;
$r0ij{Ga}{Ga} = $r0ij{Ga}{Ga} = 0.0;
$r0ij{Ga}{Zn} = $r0ij{Zn}{Ga} = 0.0;
$r0ij{Ga}{O}  = $r0ij{O}{Ga}  = 2.668100;
$r0ij{Zn}{Zn} = $r0ij{Zn}{Zn} = 0.0;
$r0ij{Zn}{O}  = $r0ij{O}{Zn}  = 3.167030;
$r0ij{O}{O}   = $r0ij{O}{O}   = 3.618701;

# 1/r^n potential parameters
my $nRn = 12;
my %CRnij;

$CRnij{In}{In} = $CRnij{In}{In} = 0.0;
$CRnij{In}{Ga} = $CRnij{Ga}{In} = 0.0;
$CRnij{In}{Zn} = $CRnij{Zn}{In} = 0.0;
$CRnij{In}{O}  = $CRnij{O}{In}  = 0.0;
$CRnij{Ga}{Ga} = $CRnij{Ga}{Ga} = 0.0;
$CRnij{Ga}{Zn} = $CRnij{Zn}{Ga} = 0.0;
$CRnij{Ga}{O}  = $CRnij{O}{Ga}  = 0.0;
$CRnij{Zn}{Zn} = $CRnij{Zn}{Zn} = 0.0;
$CRnij{Zn}{O}  = $CRnij{O}{Zn}  = 0.0;
$CRnij{O}{O}   = $CRnij{O}{O}   = 0.0;


#===============================
# Main routine
#===============================

print "\n";
print "Make LAMMPS potential table [$OutFile] for $type potential\n";
print "\n";
print "r range: $rmin - $rmax, $rstep step for $N mesh\n";
print "\n";

open(OUT, ">$OutFile") or die "Can not write to [$OutFile].\n";

print OUT "#\n";
print OUT "# Potential type: $type\n";
print     "# Potential type: $type\n";
print OUT "#  fCoulomb=$fCoulomb\n";
print     "#  fCoulomb=$fCoulomb\n";
print OUT "#  SplitEwaldSum=$SplitEwaldSum (Coulomb term included in this table)\n";
print     "#  SplitEwaldSum=$SplitEwaldSum (Coulomb term included in this table)\n";

for(my $i = 0 ; $i < @atoms ; $i++) {
	my $n1 = $atoms[$i];
	print OUT "#  $n1: Z=$Zi{$n1}\n";
	for(my $j = $i ; $j < @atoms ; $j++) {
		my $n2 = $atoms[$j];
		print OUT "#  $n1 - $n2: A=$Aij{$n1}{$n2}  rho=$rhoij{$n1}{$n2}  C=$Cij{$n1}{$n2}\n";
#		print     "#  $n1 - $n2: De=$Deij{$n1}{$n2}  a0=$a0ij{$n1}{$n2}  r0=$r0ij{$n1}{$n2}\n";
	}
}
print OUT "#\n";
print OUT "# Differentiation: h=$h  nOrder=$nDiffOrder\n";
print     "# Differentiation: h=$h  nOrder=$nDiffOrder\n";
print OUT "#\n";
print OUT "\n";
print "\n";

for(my $i = 0 ; $i < @atoms ; $i++) {
	my $n1 = $atoms[$i];
	my $Z1 = $Zi{$n1};
	for(my $j = $i ; $j < @atoms ; $j++) {
		my $n2 = $atoms[$j];
		my $Z2 = $Zi{$n2};
#		my $De = $Deij{$n1}{$n2};
#		my $a0 = $a0ij{$n1}{$n2};
#		my $r0 = $r0ij{$n1}{$n2};
#		my @args = ($fCoulomb, $Z1, $Z2, $De, $a0, $r0);
		my $A   = $Aij{$n1}{$n2};
		my $rho = $rhoij{$n1}{$n2};
		my $C   = $Cij{$n1}{$n2};
		my @args = ($fCoulomb, $Z1, $Z2, $A, $rho, $C);
		print "for $n1-$n2: ", join(', ', @args), "\n";

		print OUT "$n1-$n2\n";
		print OUT "N\t\t$N\tR\t$rmin\t$rmax\n";
		print OUT "\n";
		for(my $k = 0 ; $k < $N ; $k++) {
			my $r = $rmin + $k * $rstep;
			my $U = &Potential($r, @args);
			my $F = -&DifferentiatePotential($nDiffOrder, $r, @args);

#			printf     "%6d\t\t%8.6f\t\t%12.6g\t\t%12.6g\n", $k+1, $r, $U, $F;
			printf OUT "%6d\t\t%8.6f\t\t%12.6g\t\t%12.6g\n", $k+1, $r, $U, $F;
		}
		print OUT "\n";
	}
}


close(OUT);
exit;


#===============================
# Functions
#===============================
sub Potential
{
	my ($r, $fCoulomb, $Z1, $Z2, $A, $rho, $C) = @_;
#	my ($r, $fCoulomb, $Z1, $Z2, $De, $a0, $r0) = @_;

	return $fCoulomb * &Coulomb($r, $Z1, $Z2) 
#		 + &Morse($r, $De, $a0, $r0);
#		 + &Rn($r, 12, $CR12);
		 + &Born($r, $A, $rho)
		 - &Dispersion($r, $C);
}

sub DiffPotential
{
	my ($r, $fCoulomb, $Z1, $Z2, $De, $a0, $r0) = @_;

	my $v = 1.0 - exp(-$a0 * ($r - $r0));
	return -$fCoulomb * $KC * $Z1 * $Z2 / $r / $r
		   -2.0*$De * $v * $a0 * exp(-$a0 * ($r - $r0));
#		   -$A / $rho * exp(-$r / $rho)
#		   +6.0 * $C / $r**7;
}

sub DifferentiatePotential
{
	my ($nDiffOrder, $r, @args) = @_;

		if($nDiffOrder == -1) {
			return &DiffPotential($r, @args);
		}
		elsif($nDiffOrder == 2) {
			return Algorism::Differentiate2(
					sub { return &Potential($_[0], @args); },
					$r, $h);
		}
		elsif($nDiffOrder == 3) {
			return Algorism::Differentiate3(
					sub { return &Potential($_[0], @args); },
					$r, $h);
		}
		elsif($nDiffOrder == 5) {
			return Algorism::Differentiate5(
					sub { return &Potential($_[0], @args); },
					$r, $h);
		}
		elsif($nDiffOrder == 7) {
			return Algorism::Differentiate7(
					sub { return &Potential($_[0], @args); },
					$r, $h);
		}

		print "Error: Invalid nOrder=$nDiffOrder for differentiation "
			 ."(Choose from 2, 3, 5, or 7)\n";
		exit;
}

sub Coulomb
{
	my ($r, $Z1, $Z2) = @_;
	my $U = $KC * $Z1 * $Z2 / $r;
#print "UC=$U (KC=$KC)\n";
	return $U;
}

sub Born
{
	my ($r, $A, $rho) = @_;
	return $A * exp(-$r / $rho);
}

sub Morse
{
	my ($r, $De, $a0, $r0) = @_;

	my $v = 1.0 - exp(-$a0 * ($r - $r0));
	my $U = $De * ($v*$v - 1.0);
#print "UM=$U\n";
	return $U;
}

sub Dispersion
{
	my ($r, $C) = @_;

	return  &Rn($r, 6, $C);
}

sub Rn
{
	my ($r, $nRn, $CR) = @_;

	return $CR / $r**$nRn;
}

sub Differentiate7
{
	my ($pFunc, $x, $h) = @_;
	return (&$pFunc($x + 3*$h) - 9*&$pFunc($x + 2*$h) + 45*&$pFunc($x + $h) 
			- 45*&$pFunc($x - $h) + 9*&$pFunc($x - 2*$h) - &$pFunc($x - 3*$h)) / 60.0 / $h;
}