#!/usr/bin/perl

use strict;
#use warnings;

#===============================
# 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 = 7;      # 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} =  1.8; #3.0;
$Zi{Ga} =  1.8; #3.0;
$Zi{Zn} =  1.2; #2.0;
$Zi{O}  = -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;
			if($nDiffOrder == -1) {
				$F = -&DiffPotential($r, @args);
			}
			elsif($nDiffOrder == 2) {
				my $U1 = &Potential($r-$h, @args);
				my $U2 = &Potential($r,    @args);
				$F = -($U2 - $U1) / $h;
			}
			elsif($nDiffOrder == 3) {
				my $U1 = &Potential($r-$h, @args);
#				my $U2 = &Potential($r,    @args);
				my $U3 = &Potential($r+$h, @args);
				$F = -($U3 - $U1) / $h / 2.;
			}
			elsif($nDiffOrder == 5) {
				my $U1 = &Potential($r-2*$h, @args);
				my $U2 = &Potential($r-$h,   @args);
#				my $U3 = &Potential($r,      @args);
				my $U4 = &Potential($r+$h,   @args);
				my $U5 = &Potential($r+2*$h, @args);
				$F = -(-$U5 + 8.0*$U4 - 8.0*$U2 + $U1) / 12.0 / $h;
			}
			elsif($nDiffOrder == 7) {
				my $U1 = &Potential($r-3*$h, @args);
				my $U2 = &Potential($r-2*$h, @args);
				my $U3 = &Potential($r-$h,   @args);
#				my $U4 = &Potential($r,      @args);
				my $U5 = &Potential($r+$h  , @args);
				my $U6 = &Potential($r+2*$h, @args);
				my $U7 = &Potential($r+3*$h, @args);
				$F = -($U7 - 9.0*$U6 + 45.0*$U5 - 45.0*$U3 + 9.0*$U2 - $U1)/ 60.0 / $h;
			}
			else {
				print "Error: Invalid nOrder=$nDiffOrder for differentiation "
					 ."(Choose from 2, 3, 5, or 7)\n";
				exit;
			}

#			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, $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 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;
}