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

use strict;
#use warnings;
use File::Path;
use File::Basename;
use File::Find;

use Deps;
use Utils;
use JFile;

use MyApplication;
use Sci qw($pi $e $e0);

use Crystal::CIF;
use Crystal::Crystal;
use Crystal::Madel;

#===============================================
# デバッグ関係変数
#===============================================
#$PrintLevelが大きいほど、情報が詳しくなる
my $PrintLevel = 0;

#===============================================
# 文字コード関係変数
#===============================================
# sjis, euc, jis, noconv
my $PrintCharCode      = Deps::PrintCharCode();
my $OSCharCode         = Deps::OSCharCode();
my $FileSystemCharCode = Deps::FileSystemCharCode();

my $DirSep    = Deps::DirSep();
my $RegDirSep = Deps::RegDirSep();

#===============================================
# Applicationオブジェクト作成
#===============================================
my $App = new MyApplication;
exit if($App->Initialize() < 0);

#$App->SetLF("<br>\n");
#$App->SetPrintCharCode("sjis");
#$App->SetDebug($Debug);
$App->SetDeleteHTMLFlag(1);

#===============================================
# スクリプト大域変数
#===============================================
my $InitialDirectory = Deps::GetWorkingDirectory();
my %ParamHash;

my $Ke = $e / (4.0 * $pi * $e0); # eV

#==========================================
# コマンドラインオプション読み込み
#==========================================
$App->AddArgument("--Action",    "--Action=[Calc]", '');
$App->AddArgument("--Algo",      "--Algo=[Simple|Evjen] [Def: Simple]", "Simple");
$App->AddArgument("--Rmax",      "--Rmax=val [Def: 10.0 (A)]", 10.0);
$App->AddArgument("--Rmin",      "--Rmin=val [Def:  0.1 (A)]",  0.1);
$App->AddArgument("--DebugMode", "--DebugMode: Set DebugMode", '');
exit 1 if($App->ReadArgs(1, "sjis", 0) != 1);
my $Args = $App->Args();
#my $form = new CGI;
#$Args->SetCGIForm($form);
#$Args->parseInput($WebCharCode);

my %ArgHash = $Args->GetArgHash();
foreach my $key (keys %ArgHash) {
#print "key: $key: $ArgHash{$key}\n";
	if($key =~ /^Param:(.*?)$/i) {
		$ParamHash{$1} = $ArgHash{$key};
#print "$1:  $ArgHash{$key}\n";
	}
}
$App->{pParamHash} = \%ParamHash;

#==========================================
# メイン関数スタート
#==========================================

#Utils::InitHTML("Research", $WebCharSet, "_self");

my $Debug = $Args->GetGetArg("DebugMode");
$App->SetDebug($Debug);
my $Action = $Args->GetGetArg("Action");
$Action = 'Calc' if($Action eq '');

my $CIFFile = $Args->GetGetArg(0);
my $Algo    = $Args->GetGetArg("Algo");
$Algo       = 'Simple' if($Algo eq '');
my $Rmax    = $Args->GetGetArg('Rmax');
$Rmax       = 10.0 if($Rmax eq '');
my $Rmin    = $Args->GetGetArg('Rmin');
$Rmin       =  0.1 if($Rmin eq '');

$App->print("<b>CIFFile:</b> $CIFFile\n");
$App->print("<b>Algo   :</b> $Algo\n");
$App->print("<b>Rmax   :</b> $Rmax A\n");
$App->print("<b>Rmin   :</b> $Rmin A\n");

my $CIF = new CIF();
unless($CIF->Read($CIFFile)) {
	$App->print("  Error: Can not read [$CIFFile].\n\n");
	return -1;
}
my $Crystal = $CIF->GetCCrystal();
$Crystal->ExpandCoordinates();

my $ret = 0;
if($Action =~ /^Calc/i) {
	if($Algo eq 'Simple') {
		&CalMadelungPotentialSimple($App, $Crystal, $Rmax);
	}
	elsif($Algo eq 'Evjen') {
		&CalMadelungPotentialEvjen($App, $Crystal, $Rmax);
	}
	else {
		$App->print("Error: Invalid Algo: [$Algo]\n");
	}
}
else {
	$App->print("Error: Invalid Action: [$Action]\n");
}

#Utils::EndHTML();

exit $ret;

#===============================================
# スクリプト終了
#===============================================

#==========================================
# &Subroutines
#==========================================
sub PrepareCrystal
{
	my ($App, $Crystal, $Rmax) = @_;
	
print "Crystal structure\n";
	my ($a, $b, $c, $alpha, $beta, $gamma) = $Crystal->LatticeParametersByOutputMode(0);
print "  Lattice parameters: $a, $b, $c     $alpha, $beta, $gamma\n";
	my ($mina, $maxa) = Utils::CalMinMax([$a, $b, $c]);
print "    Minimum a: $mina\n";
	my ($nx, $ny, $nz) = (int($Rmax/$a + 1.0), int($Rmax/$b + 1.0), int($Rmax/$c + 1.0));
print "    N range: ($nx, $ny, $nz)\n";

	my @AtomType  = $Crystal->GetCAtomTypeList();
	my $nAtomType = @AtomType;
print "nAtomType: $nAtomType\n";
	my @AsymmetricAtomSite  = $Crystal->GetCAsymmetricAtomSiteList();
	my $nAsymmetricAtomSite = @AsymmetricAtomSite;
print "nAsymmetricAtomSite: $nAsymmetricAtomSite\n";
	my @AtomSite  = $Crystal->GetCExpandedAtomSiteList();
	my $nAtomSite = @AtomSite;
print "nAtomSite: $nAtomSite\n";
print "\n";

print "Asymmetric sites\n";
	my @Z0;
	for(my $is = 0 ; $is < $nAsymmetricAtomSite ; $is++) {
		my $atom0          = $AsymmetricAtomSite[$is];
		my $atomname0      = $atom0->AtomNameOnly();
		my $Z              = $atom0->Charge();
		my $occ0           = $atom0->Occupancy();
		my ($x0, $y0, $z0) = $atom0->Position(1);
		my ($charge, $sign) = ($Z =~ /(\d+)(\D*)/);
		if($sign eq '-') {
			$Z0[$is] = -$charge;
		}
		else {
			$Z0[$is] = $charge;
		}
print "  $is: $atomname0 ($x0, $y0, $z0) Z=$Z0[$is] occ=$occ0\n";
	}
print "\n";

print "All sites\n";
	my @Z1;
	for(my $it = 0 ; $it < $nAtomSite ; $it++) {
		my $atom1          = $AtomSite[$it];
		my $atomname1      = $atom1->AtomNameOnly();
		my $Z              = $atom1->Charge();
		my $occ1           = $atom1->Occupancy();
		my ($x1, $y1, $z1) = $atom1->Position(1);
		my ($charge, $sign) = ($Z =~ /(\d+)(\D*)/);
		if($sign eq '-') {
			$Z1[$it] = -$charge;
		}
		else {
			$Z1[$it] = $charge;
		}
print "  $it: $atomname1 ($x1, $y1, $z1) Z=$Z1[$it] occ=$occ1\n";
	}
print "\n";

	return ($mina, $nx, $ny, $nz, 
		$nAtomType, \@AtomType, $nAsymmetricAtomSite, \@AsymmetricAtomSite, $nAtomSite, \@AtomSite,
		\@Z0, \@Z1);
}

sub CalMadelungPotentialSimple
{
	my ($App, $Crystal, $Rmax) = @_;

	$App->print("<H2>Calculate Madelung potentials by simple summation.</H2>\n");
	
	my ($mina, $nx, $ny, $nz, 
		$nAtomType, $pAtomType, $nAsymmetricAtomSite, $pAsymmetricAtomSite, $nAtomSite, $pAtomSite,
		$pZ0, $pZ1)= &PrepareCrystal($App, $Crystal, $Rmax);

	my @AtomType           = @$pAtomType;
	my @AsymmetricAtomSite = @$pAsymmetricAtomSite;
	my @AtomSite           = @$pAtomSite;
	my @Z0 = @$pZ0;
	my @Z1 = @$pZ1;

print "Madelung potentials\n";
	my @MP;
	for(my $is = 0 ; $is < $nAsymmetricAtomSite ; $is++) {
		my $atom0          = $AsymmetricAtomSite[$is];
#		my $atomname0      = $atom0->AtomName();
		my $atomname0      = $atom0->AtomNameOnly();
#		my $iAtom          = $Crystal->FindiAtomType($atomname);
#		my $Z0             = $atom0->Charge();
		my $Z0             = $Z0[$is];
		my $occ0           = $atom0->Occupancy();
		my ($x0, $y0, $z0) = $atom0->Position();
#		my $mult      = $atom->Multiplicity();
#		$mult = 1 if($IsExpandCoordinate or $IsChooseRandomly);

print "  $is: $atomname0 ($x0, $y0, $z0) Z=$Z0 occ=$occ0: ";
		$MP[$is] = 0.0;
		for(my $ix = -$nx ; $ix <= $nx ; $ix++) {
		for(my $iy = -$ny ; $iy <= $ny ; $iy++) {
		for(my $iz = -$nz ; $iz <= $nz ; $iz++) {
#print "     i: ($ix, $iy, $iz)\n";
			my $sum = 0.0;
			for(my $it = 0 ; $it < $nAtomSite ; $it++) {
				my $atom1          = $AtomSite[$it];
				my $occ1           = $atom1->Occupancy();
#				my $Z1             = $atom1->Charge();
				my $Z1             = $Z1[$it];
				my ($x1, $y1, $z1) = $atom1->Position();

				my $r = $Crystal->GetInterAtomicDistance($x0, $y0, $z0, $x1+$ix, $y1+$iy, $z1+$iz);
				next if($r < $Rmin);
				next if($r > $Rmax);

				$sum += $occ1 * $Z1 * $Ke / ($r * 1.0e-10);
#print "r=$r ($Z1, $occ0, $occ1)\n";
			}
			$MP[$is] += $sum;
		}
		}
		}
print "MP=$MP[$is] eV\n";
	}

	return 1;
}

sub CalMadelungPotentialEvjen
{
	my ($App, $Crystal, $Rmax) = @_;

	$App->print("<H2>Calculate Madelung potentials by simple summation.</H2>\n");
	
	my ($mina, $nx, $ny, $nz, 
		$nAtomType, $pAtomType, $nAsymmetricAtomSite, $pAsymmetricAtomSite, $nAtomSite, $pAtomSite,
		$pZ0, $pZ1)= &PrepareCrystal($App, $Crystal, $Rmax);

	my @AtomType           = @$pAtomType;
	my @AsymmetricAtomSite = @$pAsymmetricAtomSite;
	my @AtomSite           = @$pAtomSite;
	my @Z0 = @$pZ0;
	my @Z1 = @$pZ1;

print "Madelung potentials\n";
	my @MP;
	for(my $is = 0 ; $is < $nAsymmetricAtomSite ; $is++) {
		my $atom0          = $AsymmetricAtomSite[$is];
#		my $atomname0      = $atom0->AtomName();
		my $atomname0      = $atom0->AtomNameOnly();
#		my $iAtom          = $Crystal->FindiAtomType($atomname);
#		my $Z0             = $atom0->Charge();
		my $Z0             = $Z0[$is];
		my $occ0           = $atom0->Occupancy();
		my ($x0, $y0, $z0) = $atom0->Position();
#		my $mult      = $atom->Multiplicity();
#		$mult = 1 if($IsExpandCoordinate or $IsChooseRandomly);

print "  $is: $atomname0 ($x0, $y0, $z0) Z=$Z0 occ=$occ0: ";
		$MP[$is] = 0.0;
		for(my $ix = -$nx ; $ix <= $nx ; $ix++) {
		for(my $iy = -$ny ; $iy <= $ny ; $iy++) {
		for(my $iz = -$nz ; $iz <= $nz ; $iz++) {
#print "     i: ($ix, $iy, $iz)\n";
			my $sum = 0.0;
			for(my $it = 0 ; $it < $nAtomSite ; $it++) {
				my $atom1          = $AtomSite[$it];
				my $occ1           = $atom1->Occupancy();
#				my $Z1             = $atom1->Charge();
				my $Z1             = $Z1[$it];
				my ($x1, $y1, $z1) = $atom1->Position();

				my $r = $Crystal->GetInterAtomicDistance($x0, $y0, $z0, $x1+$ix, $y1+$iy, $z1+$iz);
				next if($r < $Rmin);
#				next if($r > $Rmax);

				$sum += $occ1 * $Z1 * $Ke / ($r * 1.0e-10);
#print "r=$r ($Z1, $occ0, $occ1)\n";
			}
			$MP[$is] += $sum;
		}
		}
		}
print "MP=$MP[$is] eV\n";
	}

	return 1;
}