#!/usr/bin/perl

use lib 'd:/Programs/Perl/lib';

BEGIN {
#use lib 'd:/Programs/Perl/lib';
#use lib '/home/tkamiya/bin/lib';
my $BaseDir = $ENV{'PerlDir'};
$BaseDir = $ENV{'TkPerlDir'} unless(defined $BaseDir);
#print "BaseDir: $BaseDir\n";
@INC = ("$BaseDir/lib", "$BaseDir/VNL", @INC);
}

use strict;
#use warnings;

use Math::Matrix;
use Math::MatrixReal;

use Deps;
use Utils;
use ProgVars;

use Crystal::CIF;
use Crystal::ATOMS;

my $DBPath  = Deps::MakePath(ProgVars::DBDir(), "bvparm2016.cif", 0);
#my $DBPath  = Deps::MakePath(ProgVars::DBDir(), "bvparm2020.cif", 0);
my $CIFPath = ($ARGV[0])? $ARGV[0] : "LaCuOSe.cif";
my $MaxR    = ($ARGV[1])? $ARGV[1] : 2.7;

if(!-e $CIFPath) {
	my $path = Deps::MakePath("D:\\MyWebs\\Research\\CrystalStructure", $CIFPath, 0);
	my @path = glob($path);
	if(@path == 0) {
		print "Error: Can not find [$CIFPath].\n";
		exit;
	}
	$CIFPath = $path[0];
}

print "DBPath: $DBPath\n";
print "CIFPath: $CIFPath\n";

my $db = new CIF;
$db->Read($DBPath) or die "$!: Can not read [$DBPath].\n";

my %BV;
#foreach my $key (sort keys %$cif) {
#	print "$key: $cif->{$key}\n";
#}
for(my $i = 1 ; ; $i++) {
	last if(!defined $db->{"_valence_param_atom_1[$i]"});
	my $key = $db->{"_valence_param_atom_1[$i]"} . ":" . $db->{"_valence_param_atom_2[$i]"};
#	my $key = $db->{"_valence_param_atom_1[$i]"} . $db->{"_valence_param_atom_1_valence[$i]"}
#	   .":" . $db->{"_valence_param_atom_2[$i]"} . $db->{"_valence_param_atom_2_valence[$i]"};
	next if(defined $BV{$key});

	$BV{$key}->{atom1}   = $db->{"_valence_param_atom_1[$i]"};
	$BV{$key}->{atom2}   = $db->{"_valence_param_atom_2[$i]"};
	$BV{$key}->{Z1}      = $db->{"_valence_param_atom_1_valence[$i]"};
	$BV{$key}->{Z2}      = $db->{"_valence_param_atom_2_valence[$i]"};
	$BV{$key}->{r0}      = $db->{"_valence_param_Ro[$i]"};
	$BV{$key}->{B}       = $db->{"_valence_param_B[$i]"};
	$BV{$key}->{ref}     = $db->{"_valence_param_ref_id[$i]"};
	$BV{$key}->{details} = $db->{"_valence_param_details[$i]"};
#print "$i: $key: $BV{$key}->{r0}: $BV{$key}->{B}: $BV{$key}->{ref}: $BV{$key}->{details}\n";
}

my $cif = new CIF;
$cif->Read($CIFPath) or die "$!: Can not read [$CIFPath].\n";
my $Crystal = $cif->GetCCrystal();
$Crystal->ExpandCoordinates();

my @AtomSiteList           = $Crystal->GetCAsymmetricAtomSiteList();
my $nAsymmetricAtomSite    = @AtomSiteList;
my @ExpandedAtomSiteList   = $Crystal->GetCExpandedAtomSiteList();
my $nTotalExpandedAtomSite = $Crystal->nTotalExpandedAtomSite();
for(my $i = 0 ; $i < $nAsymmetricAtomSite ; $i++) {
	my $atom1   = $AtomSiteList[$i];
	my $name1   = $atom1->AtomNameOnly();
	my $charge1 = $atom1->Charge();
	my ($x1, $y1, $z1) = $atom1->Position(1);

print "$i: $name1$charge1\n";
	my $BVS = 0.0;
	for(my $j = 0 ; $j < $nTotalExpandedAtomSite ; $j++) {
		my $atom2   = $ExpandedAtomSiteList[$j];
		my $name2   = $atom2->AtomNameOnly();
		my $charge2 = $atom2->Charge();
		my ($x2, $y2, $z2) = $atom2->Position(1);

		for(my $ix = -1 ; $ix <= 1 ; $ix++) {
		for(my $iy = -1 ; $iy <= 1 ; $iy++) {
		for(my $iz = -1 ; $iz <= 1 ; $iz++) {

		my $d = $Crystal->GetInterAtomicDistance($x1, $y1, $z1, $x2+$ix, $y2+$iy, $z2+$iz);
#		my $d = $Crystal->GetNearestInterAtomicDistance($x1, $y1, $z1, $x2, $y2, $z2);
		next if($d == 0.0 or $d > $MaxR);

		my $bvkey = "$name1:$name2";
#print "key=$bvkey\n";
		if(!defined $BV{$bvkey}) {
			$bvkey = "$name2:$name1";
			if(!defined $BV{$bvkey}) {
				print "  Error: Can not find BV parameters for [$bvkey].\n";
				print("\nPress ENTER to terminate>> ");
				<STDIN>;
				next;
#				last;
			}
		}
		my $r0       = $BV{$bvkey}->{r0};
		my $B        = $BV{$bvkey}->{B};
		my $aname1   = $BV{$bvkey}->{atom1};
		my $aname2   = $BV{$bvkey}->{atom2};
		my $valence1 = $BV{$bvkey}->{Z1};
		my $valence2 = $BV{$bvkey}->{Z2};
		my $ref      = $BV{$bvkey}->{ref};
		my $bv = exp( ($r0 - $d) / $B);
		$BVS += $bv;
printf "  %2d: %-4s r0=%5.3f B=%4.2f  ri=%5.3f bv=%5.3f [%s for %s]\n", 
			$j, "$name2$charge2", $r0, $B, $d, $bv, $ref, 
			"$aname1 $valence1 - $aname2 $valence2";

		}
		}
		}
	}
printf "    Bond Valence Sum: %6.3f\n", $BVS;
}

print("\nPress ENTER to terminate>> ");
<STDIN>;
exit;