#!/usr/bin/perl

use strict;

use lib 'd:/Programs/Perl/lib';
use JFile;
use Sci qw($e0 $e $pi);
use Crystal::MXD;

#F0...Born-Mayer Potential CONSTANT (J/A)
#  R(r) = F0 * (bi+bj) * exp( (ai+aj - rij) / (bi+bj) )
#my $F0 = 6.947700141e-21;
#K0...Coulomb Potential Coefficient (J.A)  e^2/4/PI/E0*1.0E+10
#K0=2.307100407D-18

my ($rmin, $rmax, $rstep) = (0.2, 5.0, 0.05);
my $nr = int( ($rmax - $rmin) / $rstep + 0.1);
my $nMaxIter = 100;
my $rEPS     = 0.001;
my $F05File = "FILE05.DAT";
my $CSVFile = "out.csv";

my $mxd = new MXD;

my $in = JFile->new($F05File, 'r') or die "$!: Can not read [$F05File].\n";
$in->SkipTo('BUSING');
my @Atoms;
my $count = 0;
while(1) {
	my $line = $in->ReadLine();
	$line =~ s/:\s*$//;
	Utils::DelSpace($line);
	last if($line eq '');

	my ($idx, $AtomName, $nAtom, $Z, $M, $a, $b, $c) = Utils::Split("\\s+", $line);
	$nAtom = int($nAtom + 0.1);
print "$idx, $AtomName, $nAtom, $Z, $a, $b, $c\n";
	$Atoms[$count] = {
		idx      => $idx,
		AtomName => $AtomName,
		nAtom    => $nAtom,
		Z        => $Z,
		Mass     => $M,
		a        => $a,
		b        => $b,
		c        => $c,
		};
	$count++;
#print "c=$count\n";
}

$in->Close();

print "Save potential files\n";
my @P;
for(my $i = 0 ; $i < $count ; $i++) {
#print "i=$i\n";
	$P[$i] = [];
	for(my $j = 0 ; $j <= $i ; $j++) {
		$P[$i][$j] = [];

		my $CSVFile = "$Atoms[$i]->{AtomName}-$Atoms[$j]->{AtomName}.csv";
print "Save to [$CSVFile]\n";
		my $out = JFile->new($CSVFile, 'w') or die "$!: Can not write to [$CSVFile].\n";
		$out->print("r(A),Uij(J)\n");

		for(my $k = 0 ; $k < $nr ; $k++) {
#print "    k=$k\n";
			my $r = $rmin + $k * $rstep;

			my $Uij = $mxd->UijWithParams($r, $Atoms[$i]->{Z}, $Atoms[$i]->{a}, $Atoms[$i]->{b}, $Atoms[$i]->{c},
											  $Atoms[$j]->{Z}, $Atoms[$j]->{a}, $Atoms[$j]->{b}, $Atoms[$j]->{c});
			$out->print("$r,$Uij\n");
		}
		$out->Close();
	}
}

print "r at minimum potential\n";
for(my $i = 0 ; $i < $count ; $i++) {
	for(my $j = 0 ; $j <= $i ; $j++) {
print "($i,$j): $Atoms[$i]->{AtomName}-$Atoms[$j]->{AtomName}\n";
		my $r1 = $rmin;
		my $r2 = $rmax;
		my $dU1 = $mxd->dUijdrWithParams($r1, $Atoms[$i]->{Z}, $Atoms[$i]->{a}, $Atoms[$i]->{b}, $Atoms[$i]->{c},
										  $Atoms[$j]->{Z}, $Atoms[$j]->{a}, $Atoms[$j]->{b}, $Atoms[$j]->{c});
		my $dU2 = $mxd->dUijdrWithParams($r2, $Atoms[$i]->{Z}, $Atoms[$i]->{a}, $Atoms[$i]->{b}, $Atoms[$i]->{c},
										  $Atoms[$j]->{Z}, $Atoms[$j]->{a}, $Atoms[$j]->{b}, $Atoms[$j]->{c});
		if($dU1 == 0.0) {
			print "  rmin = $r1 (dUmin=$dU1)\n";
			next;
		}
		if($dU2 == 0.0) {
			print "  rmin = $r2 (dUmin=$dU2)\n";
			next;
		}
		if($dU1 * $dU2 > 0.0) {
			print "  Error: Bad initial values (dU1*dU2 = ", $dU1*$dU2 > 0, ")\n";
			next;
		}
		my $IsConverged = 0;
		my ($rmin, $dUmin);
		for(my $k = 0 ; $k < $nMaxIter ; $k++) {
			my $r3 = ($r1 + $r2) / 2.0;
			my $dU3 = $mxd->dUijdrWithParams($r3, $Atoms[$i]->{Z}, $Atoms[$i]->{a}, $Atoms[$i]->{b}, $Atoms[$i]->{c},
											  $Atoms[$j]->{Z}, $Atoms[$j]->{a}, $Atoms[$j]->{b}, $Atoms[$j]->{c});
			if($dU3 == 0.0) {
				$IsConverged = 1;
				$rmin  = $r3;
				$dUmin = $dU3;
				last;
			}
			if($dU1 * $dU3 < 0.0) {
				$r2 = $r3;
			}
			else {
				$r1 = $r3;
			}
			if(abs($r1 - $r2) < $rEPS) {
				$IsConverged = 1;
				$rmin  = ($r1 + $r2) / 2.0;
				$dUmin = (abs($dU1) < abs($dU2))? $dU2 : $dU1;
				last;
			}
		}
		print "  rmin = $rmin (dUmin=$dUmin)\n";
	}
}



exit;