#!/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 Utils;
use JFile;

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;
print "KC = $e / 4.0 / $pi / $e0 * 1.0e10 = $KC\n";
print "\n";

#===============================
# File configuration
#===============================
print "ARGV[0]: $ARGV[0]\n";

if($ARGV[0] eq '') {
	print "\nUsage: perl ComparePotential.pl InputFile\n\n";
	exit;
}

my $InFile = $ARGV[0];
my $Dir = $InFile;
$Dir =~ s/\..*?$/_lmp_tmp/;

$InFile = "$Dir/lmp_tmp.in";

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

print "\n";
print "Make LAMMPS potential csv for comparison\n";
print "\n";

print "Read [$InFile]\n";
print "\n";

my $pHash = &ReadPotentialFromInput($InFile);
my $pAtoms = $pHash->{pAtoms};

for(my $i = 0 ; $i < @$pAtoms ; $i++) {
	for(my $j = $i ; $j < @$pAtoms ; $j++) {
print "\n";
print "pair $i,$j: style=$pHash->{pair_style}\n";
		my $ip = $i+1;
		my $jp = $j+1;
		my $key = "$ip-$jp";
		my $Z1 = $pHash->{"Z$ip"};
		my $Z2 = $pHash->{"Z$jp"};
		my $pParams  = $pHash->{"PotentialParameters$key"};
		my $FileName = $pHash->{"PotentialFile$key"};

		my ($De, $a0, $r0);
		my ($A, $rho, $sigma, $C);
		my ($file, $kye, $rcut);
		if($pHash->{pair_style} =~ /coul/i) {
			print "Z1=$Z1  Z2=$Z2\n";
		}
		if($pHash->{pair_style}  =~ /morse/i) {
			($De, $a0, $r0) = @$pParams;
			print "De=$De  a0=$a0  r0=$r0\n";
		}
		if($pHash->{pair_style}  =~ /born/i) {
			($A, $rho, $sigma, $C) = @$pParams;
			print "A=$A  rho=$rho  sigma=$sigma  C=$C\n";
		}
		if($pHash->{pair_style}  =~ /buck/i) {
			($A, $rho, $C) = @$pParams;
			print "A=$A  rho=$rho  C=$C\n";
		}
		if($pHash->{pair_style}  =~ /table/i) {
			($file, $key, $rcut) = @$pParams;
			print "file=$file  key=$key  rcut=$rcut\n";
		}
#print "$FileName = pHash->{PotentialFile$key}\n";

		my $WriteFile = "$Dir/$FileName";
		print "Read potential from [$WriteFile]\n";
		my $in = JFile->new($WriteFile, "r") or die "Can not read [$WriteFile].\n";
		my (@r, @U);
		if($pHash->{pair_style}  =~ /table/i) {
			my $tbl = JFile->new($file, 'r') or die "Can not read [$file]\n";

			my $line = $tbl->SkipTo("^$key");
#print "l: $line\n";
			$line = $tbl->ReadLine();
			my ($key, $nr, $type, $rmin, $rmax) = Utils::Split("\\s+", $line);
print "  r = ($rmin, $rmax) nr=$nr [type=$type]\n";
			for(my $i = 0 ; $i < $nr ; ) {
				$line = $tbl->ReadLine();
				last if(!defined $line);
				my ($idx, $r, $U, $F) = Utils::Split("\\s+", $line);
				next if(!defined $U);
				
				$r[$i] = $r + 0.0;
				$U[$i] = $U + 0.0;
#print "$i: r=$r  U=$U\n";
				$i++;
			}

			$tbl->Close();
		}
#exit;

		my $OutFile = "Comp$pAtoms->[$i]-$pAtoms->[$j].csv";
		print "Write potentials to [$OutFile]\n";

		my $out = JFile->new($OutFile, "w") or die "Can not write to [$OutFile].\n";
		$out->print("r,Write,Calc\n");
		for(my $i = 0 ; $i < @r ; $i++) {
			my $r = $r[$i];
			my $Ucal = 0.0;
			if($pHash->{pair_style}  =~ /coul/i) {
				$Ucal += &Coulomb($r, $Z1, $Z2);
			}
			if($pHash->{pair_style}  =~ /morse/i) {
				my $UM = &Morse($r, $De, $a0, $r0);
#				my $UR = &Rn($r, 12, $CR12);
				$Ucal += $UM;
			}
			if($pHash->{pair_style}  =~ /born/i) {
				my $UB = &Born($r - $sigma, $A, $rho);
				my $UD = &Dispersion($r, $C);
				$Ucal += $UB + $UD;
			}
			if($pHash->{pair_style}  =~ /buck/i) {
				my $UB = &Born($r, $A, $rho);
				my $UD = &Dispersion($r, $C);
				$Ucal += $UB + $UD;
			}
			if($pHash->{pair_style}  =~ /table/i) {
				$Ucal += Algorism::Interpolate(\@r, \@U, $r);
			}
			$out->print("$r[$i],$U[$i],$Ucal\n");
		}

		$out->Close();
	}
}
exit;


#===============================
# Functions
#===============================
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 = exp(-$a0 * ($r - $r0));
#	my $U = $De * ($v*$v - 2.0*$v);
	my $v = 1.0 - exp(-$a0 * ($r - $r0));
	my $U = $De * ($v*$v - 1.0);
#print "UM=$U ($r, $De, $a0, $r0)\n";
	return $U;
}

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

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

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

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

sub ReadPotentialFromInput
{
	my ($InFile) = @_;

	my %hash;
	my $in = JFile->new($InFile, "r") or die "Can not read [$InFile]\n";

	my $line = $in->SkipTo("pair_style");
	my @a;
	my $keyword;
	($keyword, $hash{pair_style}, $hash{global_cutoff}, @a) = Utils::Split("\\s+", $line);
print "pair_style: [$hash{pair_style}][$hash{global_cutoff}][$a[0]][$a[1]]\n";
	if($hash{pair_style} eq 'table') {
		$hash{table_type} = $hash{global_cutoff};
		$hash{global_cutoff} = '';
		$hash{nr} = $a[0];
		$hash{ewald} = 1 if($a[1] eq 'ewald');
	}

	$line = $in->ReadLine();
	$line =~ s/^#//;
	my @atoms = Utils::Split("\\s+", $line);
	$hash{pAtoms} = \@atoms;
print "atoms: ", join(', ', @atoms), "\n";
print "\n";

	while(1) {
		$line = $in->ReadLine();
		last if(!defined $line);
		next if($line =~ /^#/);
		last if($line !~ /pair_coeff/i);

		my ($keyword, $idx1, $idx2, @params) = Utils::Split("\\s+", $line);
		$hash{"PotentialParameters$idx1-$idx2"} 
			= $hash{"PotentialParameters$idx2-$idx1"} 
			= \@params;
		my ($idx1m, $idx2m) = ($idx1-1, $idx2-1);
print "$atoms[$idx1m] - $atoms[$idx2m]: ", join(', ', @params), "\n";
	}

	$in->rewind();
	while(1) {
		$line = $in->SkipTo("pair_write");
		last if(!defined $line);
		next if($line =~ /^#/);

		my ($keyword, $idx1, $idx2, $nr, $meshtype, $rmin, $rmax, 
			$filename, $key, $Z1, $Z2)
			= Utils::Split("\\s+", $line);
		$hash{"PotentialFile$idx1-$idx2"} 
			= $hash{"PotentialFile$idx2-$idx1"} = $filename;
		$hash{"Z$idx1"} = $Z1;
		$hash{"Z$idx2"} = $Z2;
		my ($idx1m, $idx2m) = ($idx1-1, $idx2-1);
print "$atoms[$idx1m] - $atoms[$idx2m]: PotentialFile$idx1-$idx2: $filename\n";
	}

	$in->Close();

	return \%hash;
}