#!/usr/bin/perl

BEGIN {
use lib 'c:/Programs/Perl/lib';
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 Utils;
use JFile;

use MyApplication;
use Sci qw($pi $e2_4pie0 $torad $todeg $e $e0 erfc erf log10);

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

#===============================================
# デバッグ関係変数
#===============================================
#$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
my $pi2  = 2.0 * $pi;

#==========================================
# コマンドラインオプション読み込み
#==========================================
$App->AddArgument("--Action",    "--Action=[Calc]", '');
$App->AddArgument("--Algo",      "--Algo=[Simple|Evjen|Evjen2|EWALD] [Def: EWALD]", "EWALD");
$App->AddArgument("--wBorder",   "--wBorder=val [Def: 0.1]",       0.1);
$App->AddArgument("--Shift",     "--Shift=val(dx,dy,dz) [Def: ]",  "");
$App->AddArgument("--Rmax",      "--Rmax=val [Def: 10.0 (A)]",    10.0);
$App->AddArgument("--Rmin",      "--Rmin=val [Def:  0.1 (A)]",     0.1);
$App->AddArgument("--EWALD",     "--EWALD=val [Def: 0.4]",         0.4);
$App->AddArgument("--EPS",       "--EPS=val [Def: 1.0e-4]",     1.0e-4);
$App->AddArgument("--ZEPS",      "--ZEPS=val [Def: 1.0e-3]",    1.0e-3);
$App->AddArgument("--UseP1",     "--UseP1=[0|1] [Def: 0]",           0);
$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       = 'EWALD' if($Algo eq '');
my $Rmax    = $Args->GetGetArg('Rmax');
$Rmax       = 10.0 if($Rmax eq '');
my $ZEPS    = $Args->GetGetArg('ZEPS');
$ZEPS       =  1.0e-3 if($ZEPS eq '');
my $Rmin    = $Args->GetGetArg('Rmin');
$Rmin       =  0.1 if($Rmin eq '');
my $wBorder = $Args->GetGetArg('wBorder');
$wBorder = 0.1 if(!defined $wBorder);
my $Shift   = $Args->GetGetArg('Shift');
my $UseP1   = $Args->GetGetArg('UseP1');
$UseP1 = 0 if(!defined $UseP1);

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

if(!-f $CIFFile) {
	$App->print("\n\nError: Can not read [$CIFFile]\n");
	$App->Usage();
	exit -1;
}

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

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

print("Press ENTER to terminate >>");
<STDIN>;
#Utils::EndHTML();

exit $ret;

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

#==========================================
# &Subroutines
#==========================================
sub ExpandBorderSite
{
	my ($atom, $wb, $Algo, $Shift) = @_;

	my ($dx, $dy, $dz) = Utils::Split("\\s*,\\s*", $Shift);
	$dx += 0.0;
	$dy += 0.0;
	$dz += 0.0;

	if($Algo ne 'Evjen2' or $wb eq '') {
		my ($x, $y, $z) = $atom->Position(1);
		$atom->SetPosition($x+$dx, $y+$dy, $z+$dz);
		return ($atom) 
	}

	my $atomname    = $atom->AtomName();
	my $label       = $atom->Label();
	my $Z           = $atom->Charge();
	my $occ         = $atom->Occupancy();
	my ($x, $y, $z) = $atom->Position(1);
	$x += $dx;
	$y += $dy;
	$z += $dz;
print "$label: ($x, $y, $z) occ=$occ\n";

	my (@pos);
	my $c    = 0;
	my $mult = 1;
	$pos[$c++] = [$x, $y, $z];
#print "xyz=$pos[0][0]. $pos[0][1]. $pos[0][2]\n";
#print "wb=$wb\n";

	if(abs($x) < $wb) {
#print "Boarder: $label: ($x, $y, $z)\n";
		my $n = @pos;
		for(my $i = 0 ; $i < $n ; $i++) {
			$pos[$c++] = [$pos[$i][0]+1.0, $pos[$i][1], $pos[$i][2]];
		}
		$occ  /= 2;
		$mult *= 2;
	}
	if(abs(1.0 - $x) < $wb) {
#print "Boarder: $label: ($x, $y, $z)\n";
		my $n = @pos;
		for(my $i = 0 ; $i < $n ; $i++) {
			$pos[$c++] = [$pos[$i][0]-1.0, $pos[$i][1], $pos[$i][2]];
		}
		$occ  /= 2;
		$mult *= 2;
	}
	if(abs($y) < $wb) {
#print "Boarder: $label: ($x, $y, $z)\n";
		my $n = @pos;
		for(my $i = 0 ; $i < $n ; $i++) {
			$pos[$c++] = [$pos[$i][0], $pos[$i][1]+1.0, $pos[$i][2]];
		}
		$occ  /= 2;
		$mult *= 2;
	}
	if(abs(1.0 - $y) < $wb) {
#print "Boarder: $label: ($x, $y, $z)\n";
		my $n = @pos;
		for(my $i = 0 ; $i < $n ; $i++) {
			$pos[$c++] = [$pos[$i][0], $pos[$i][1]-1.0, $pos[$i][2]];
		}
		$occ  /= 2;
		$mult *= 2;
	}
	if(abs($z) < $wb) {
#print "Boarder: $label: ($x, $y, $z)\n";
		my $n = @pos;
		for(my $i = 0 ; $i < $n ; $i++) {
			$pos[$c++] = [$pos[$i][0], $pos[$i][1], $pos[$i][2]+1.0];
		}
		$occ  /= 2;
		$mult *= 2;
	}
	if(abs(1.0 - $z) < $wb) {
#print "Boarder: $label: ($x, $y, $z)\n";
		my $n = @pos;
		for(my $i = 0 ; $i < $n ; $i++) {
			$pos[$c++] = [$pos[$i][0], $pos[$i][1], $pos[$i][2]-1.0];
		}
		$occ  /= 2;
		$mult *= 2;
	}

print "  mult=$mult  occ=$occ\n";
	my @s;
	for(my $i = 0 ; $i < @pos ; $i++) {
print "    ($pos[$i][0], $pos[$i][1], $pos[$i][2])  occ=$occ\n";
			$s[$i] = new AtomSite;
			$s[$i]->SetAtomName($atomname);
			$s[$i]->SetLabel($label);
			$s[$i]->SetPosition($pos[$i][0], $pos[$i][1], $pos[$i][2]);
			$s[$i]->SetOccupancy($occ);
	}

	return @s;
}

sub PrepareCrystal
{
	my ($App, $Crystal, $Algo, $Rmax, $wBorder, $Shift) = @_;

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";

	my @a = @AsymmetricAtomSite;
	@AsymmetricAtomSite = ();
	my $c = 0;
	for(my $is = 0 ; $is < $nAsymmetricAtomSite ; $is++) {
		my $atom0          = $a[$is];
		my @s = &ExpandBorderSite($atom0, $wBorder, $Algo, $Shift);
		for(my $i = 0 ; $i < @s ; $i++) {
			$AsymmetricAtomSite[$c++] = $s[$i];
		}
	}
	$nAsymmetricAtomSite = @AsymmetricAtomSite;

	@a = @AtomSite;
	@AtomSite = ();
	$c = 0;
	for(my $it = 0 ; $it < $nAtomSite ; $it++) {
		my $atom0          = $a[$it];
		my @s = &ExpandBorderSite($atom0, $wBorder, $Algo, $Shift);
		for(my $i = 0 ; $i < @s ; $i++) {
			$AtomSite[$c++] = $s[$i];
		}
	}
	$nAtomSite = @AtomSite;

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();
		if(undef $Z or $Z == '') {
			my $iAtomType = $Crystal->FindiAtomType($atomname0);
			$Z = $AtomType[$iAtomType-1]->Charge();
			$atom0->SetCharge($Z);
		}
		my $occ0           = $atom0->Occupancy();
		my ($x0, $y0, $z0) = $atom0->Position();
#		my ($charge, $sign) = ($Z =~ /([\d\.]+)(\D*)/);
#		if($sign eq '-') {
#			$Z0[$is] = -$charge;
#		}
#		else {
#			$Z0[$is] = $charge;
#		}
		$Z0[$is] = $Z;
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 $Z1             = $atom1->Charge();
		if(undef $Z1 or $Z1 == '') {
			my $iAtomType = $Crystal->FindiAtomType($atomname1);
			$Z1 = $AtomType[$iAtomType-1]->Charge();
			$atom1->SetCharge($Z1);
		}
		my $occ1           = $atom1->Occupancy();
		my ($x1, $y1, $z1) = $atom1->Position();
#		my ($charge, $sign) = ($Z =~ /([\d\.]+)(\D*)/);
#		if($sign eq '-') {
#			$Z1[$it] = -$charge;
#		}
#		else {
#			$Z1[$it] = $charge;
#		}
		$Z1[$it] = $Z1;
print "  $it: $atomname1 ($x1, $y1, $z1) Z=$Z1[$it] occ=$occ1\n";
	}
print "\n";

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

sub TotalCharge
{
	my ($App, $pAtomSite, $pZ) = @_;

	my $Z = 0.0;
	for(my $is = 0 ; $is < @$pAtomSite ; $is++) {
		my $atom   = $pAtomSite->[$is];
		my $charge = $pZ->[$is];
		my $occ    = $atom->Occupancy();
		$Z += $occ * $charge;
	}
	return $Z;
}

sub CalMadelungPotentialEWALD
{
	my ($App, $Crystal, $Algo, $Rmax, $wBorder, $Shift) = @_;

	$App->print("<H2>Calculate Madelung potentials by EWALD method\.</H2>\n");

	my ($mina, $maxa, $nx, $ny, $nz, 
		$nAtomType, $pAtomType, $nAsymmetricAtomSite, $pAsymmetricAtomSite, $nAtomSite, $pAtomSite,
		$pZ0, $pZ1)= &PrepareCrystal($App, $Crystal, $Algo, $Rmax, $wBorder, $Shift);

	my $EWalpha = $Args->GetGetArg('EWALD');
	$EWalpha    = 0.4 if($EWalpha eq '');
	my $EPS     = $Args->GetGetArg('EPS');
	$EPS        = 1.0e-4 if($EPS eq '');
	$App->print("<b>EWALD alpha:</b> $EWalpha\n");
	$App->print("<b>EPS        :</b> $EPS\n");

	my $N = -log10($EPS);
#print "N=$N\n";
	my $RDmax = (2.26 + 0.26 * $N) / $EWalpha;
	my $ErfcRDmax = erfc($EWalpha*$RDmax);
print "RDmax = $RDmax A  erfc(alpha*RDmax) = $ErfcRDmax\n";

	my ($a, $b, $c, $lalpha, $lbeta, $lgamma) = $Crystal->LatticeParametersByOutputMode(0);
	my ($g11, $g12, $g13, $g21, $g22, $g23, $g31, $g32, $g33) = $Crystal->Metrics();
	my $vol = $Crystal->Volume();
	my ($Ra, $Rb, $Rc, $Ralpha, $Rbeta, $Rgamma) = $Crystal->ReciprocalLatticeParameters();
	my ($Rg11, $Rg12, $Rg13, $Rg21, $Rg22, $Rg23, $Rg31, $Rg32, $Rg33) = $Crystal->RMetrics();
print "g11=$g11  Rg11=$Rg11   g22=$g22  Rg22=$Rg22   g33=$g33  Rg33=$Rg33\n";
print "g12=$g12  Rg12=$Rg12   g23=$g23  Rg23=$Rg23   g31=$g31  Rg31=$Rg31\n";
print "g21=$g21  Rg21=$Rg21   g32=$g32  Rg23=$Rg32   g13=$g13  Rg13=$Rg13\n";
		my $Rvol = $Crystal->CalculateReciprocalVolume();

	my (@lsin, @NRmax);
	$lsin[0] = sin($torad * $lalpha);
	$lsin[1] = sin($torad * $lbeta);
	$lsin[2] = sin($torad * $lgamma);
	$NRmax[0] = int($RDmax / sqrt($g11 * $lsin[1] * $lsin[2])) + 1;
	$NRmax[1] = int($RDmax / sqrt($g22 * $lsin[2] * $lsin[0])) + 1;
	$NRmax[2] = int($RDmax / sqrt($g33 * $lsin[0] * $lsin[1])) + 1;
print "Nmax = ($NRmax[0], $NRmax[1], $NRmax[2])\n";

#print "V=$vol nAtom=$nAtomSite\n";
	my $CalN = (int(2.0 / 3.0 * $pi2 * $RDmax*$RDmax*$RDmax / $vol) + 1.0) * $nAtomSite;
print "CalN = $CalN\n";

	my ($max, $max2, $max3);
	$max  = sqrt($g11 * $lsin[1] * $lsin[2]) * $NRmax[0];
	$max2 = sqrt($g22 * $lsin[0] * $lsin[2]) * $NRmax[1];
	$max = $max2 if($max2 > $max);
	$max2 = sqrt($g33 * $lsin[0] * $lsin[1]) * $NRmax[2];
	$max = $max2 if($max2 > $max);
	$max3 = sqrt($g11 * $lsin[1] * $lsin[2]) * ($NRmax[0] - 1);
	$max2 = sqrt($g22 * $lsin[0] * $lsin[2]) * ($NRmax[1] - 1);
	$max3 = $max2 if($max2 > $max3);
	$max2 = sqrt($g33 * $lsin[0] * $lsin[1]) * ($NRmax[2] - 1);
	$max3 = $max2 if($max2 > $max3);
print "Rmax(Calculation) = $max3 - $max\n";

#print "Ra = ($Ralpha, $Rbeta, $Rgamma)\n";
#print "log10=", log(10.0), "\n";
	my $G2max = 4.0 * $EWalpha* $EWalpha * $N * log(10.0) + 1.0;
#	my $G2max = $EWalpha * $EWalpha / ($pi*$pi) * (($N+10.0) * log(10.0));
#	my $G2max = $EWalpha * $EWalpha / ($pi*$pi) * (-log($pi * $vol / ($maxa*$maxa)) + ($N+10.0) * log(10.0));
print "G2max = $G2max\n";
#print "exp(-G2max / 4.0 / alpha^2) = ", exp(-$G2max / 4.0 / $EWalpha**2), "\n";
print "exp(-pi^2*G2max / alpha^2) = ", exp(-$G2max * $pi * $pi / $EWalpha**2), "\n";
	$lsin[0] = sin($torad * $Ralpha);
	$lsin[1] = sin($torad * $Rbeta);
	$lsin[2] = sin($torad * $Rgamma);
	my @HGmax;
	$HGmax[0] = int(sqrt($G2max / ($Rg11 * $lsin[1] * $lsin[2])) / $pi2) + 1;
	$HGmax[1] = int(sqrt($G2max / ($Rg22 * $lsin[0] * $lsin[2])) / $pi2) + 1;
	$HGmax[2] = int(sqrt($G2max / ($Rg33 * $lsin[0] * $lsin[1])) / $pi2) + 1;
print "HGmax = ($HGmax[0], $HGmax[1], $HGmax[2])\n";

	$CalN = int(2.0 / 3.0 * $pi2 * $G2max**1.5 / $Rvol / $pi2**3) + 1;
print "CalN = $CalN\n";
	$max  = sqrt($Rg11 * $lsin[1] * $lsin[2]) * $HGmax[0];
	$max2 = sqrt($Rg22 * $lsin[0] * $lsin[2]) * $HGmax[1];
	$max = $max2 if($max2 > $max);
	$max2 = sqrt($Rg33 * $lsin[0] * $lsin[1]) * $HGmax[2];
	$max = $max2 if($max2 > $max);
	$max3 = sqrt($Rg11 * $lsin[1] * $lsin[2]) * ($HGmax[0] - 1);
	$max2 = sqrt($Rg22 * $lsin[0] * $lsin[2]) * ($HGmax[1] - 1);
	$max3 = $max2 if($max2 > $max3);
	$max2 = sqrt($Rg33 * $lsin[0] * $lsin[1]) * ($HGmax[2] - 1);
	$max3 = $max2 if($max2 > $max3);
printf "G2max(Calculation) = %g - %g\n", ($max3*$pi2)**2, ($max*$pi2)**2;

#for(my $x = 0.0 ; $x <= 6.1 ; $x += 2.0) {
#print "erfc($x)=", erfc($x), ": erf($x)=", erf($x), "\n";
#}

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

	my $TotalCharge = &TotalCharge($App, \@AtomSite, \@Z1);
	if(abs($TotalCharge) > $ZEPS) {
		print "Error: Total charge ($TotalCharge) is greater than EPS=$ZEPS.\n";
		exit -1;
	}

print "Madelung potentials\n";
#print "n=$nAtomSite, $nx, $ny, $nz   Rmax=$Rmax\n";
	my $Kexp = $pi * $pi / $EWalpha / $EWalpha;
	my $Krec = $Ke / $pi / $vol / 1.0e-10;
	my (@UC1, @UC2, @UC3, @MP);
	my ($pA, $pZ);
	if($UseP1) {
		$pA = \@AtomSite;
		$pZ = \@Z1;
	}
	else {
		$pA = \@AsymmetricAtomSite;
		$pZ = \@Z0;
	}
	for(my $is = 0 ; $is < @$pA ; $is++) {
		my $atom0          = $pA->[$is];
#		my $atomname0      = $atom0->AtomName();
		my $atomname0      = $atom0->AtomNameOnly();
#		my $iAtom          = $Crystal->FindiAtomType($atomname);
#		my $Z0             = $atom0->Charge();
		my $Z0             = $pZ->[$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:\n";
		$UC1[$is] = 0.0;
		$UC2[$is] = 0.0;
		$UC3[$is] = 2.0 * $Ke * $occ0 * $Z0 * ($EWalpha*1.0e10) / sqrt($pi);
		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);

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

		for(my $h = -$HGmax[0] ; $h <= $HGmax[0] ; $h++) {
		for(my $k = -$HGmax[1] ; $k <= $HGmax[1] ; $k++) {
#		for(my $l = -$HGmax[2] ; $l <= $HGmax[2] ; $l++) {
		for(my $l = 0          ; $l <= $HGmax[2] ; $l++) {
			my $phi0  = $pi2 * ($h * $x0 + $k * $y0 + $l * $z0);
			my $cspi0 = cos($phi0);
			my $snpi0 = sin($phi0);

			my $G2 = $Rg11 * $h*$h + $Rg22 * $k*$k + $Rg33 * $l*$l
				   + 2.0*($Rg23 * $k*$l + $Rg13 * $h*$l + $Rg12 * $h*$k);
			next if($G2 > $G2max);
			next if($G2 == 0.0);

			my (@cspi, @snpi);
			my $cssum = 0.0;
			my $snsum = 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 $phi = $pi2 * ($h * $x1 + $k * $y1 + $l * $z1);
				$cspi[$it] = $occ1 * $Z1 * cos($phi);
				$snpi[$it] = $occ1 * $Z1 * sin($phi);
				$cssum += $cspi[$it];
				$snsum += $snpi[$it];
			}

			my $expg   = exp(-$Kexp * $G2) / $G2;
			my $fcal = $cspi0 * $cssum + $snpi0 * $snsum;
			if($l != 0) {
				$fcal *= 2.0;
			}
			$UC2[$is] += $Krec * $expg * $fcal;
		}
		}
		}

		$MP[$is] = $UC1[$is] + $UC2[$is] - $UC3[$is];
printf "    MP=%15.6g eV (= %15.6g + %15.6g - %15.6g)\n", $MP[$is], $UC1[$is], $UC2[$is], $UC3[$is];
		my $MPZ  = $MP[$is]  * $e * $Z0;
		my $UC1Z = $UC1[$is] * $e * $Z0;
		my $UC2Z = $UC2[$is] * $e * $Z0;
		my $UC3Z = $UC3[$is] * $e * $Z0;
printf "    MP=%15.6g J (= %15.6g + %15.6g - %15.6g)\n", $MPZ, $UC1Z, $UC2Z, $UC3Z;
	}

	return 1;
}

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

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

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

	my $TotalCharge = &TotalCharge($App, \@AtomSite, \@Z1);
	if(abs($TotalCharge) > $ZEPS) {
		print "Error: Total charge ($TotalCharge) is greater than EPS=$ZEPS.\n";
		exit -1;
	}

print "Madelung potentials\n";
	my @MP;
	my ($pA, $pZ);
	if($UseP1) {
		$pA = \@AtomSite;
		$pZ = \@Z1;
	}
	else {
		$pA = \@AsymmetricAtomSite;
		$pZ = \@Z0;
	}
	for(my $is = 0 ; $is < @$pA ; $is++) {
		my $atom0          = $pA->[$is];
#		my $atomname0      = $atom0->AtomName();
		my $atomname0      = $atom0->AtomNameOnly();
#		my $iAtom          = $Crystal->FindiAtomType($atomname);
#		my $Z0             = $atom0->Charge();
		my $Z0             = $pZ->[$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, $Algo, $Rmax, $wBorder, $Shift) = @_;

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

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

	my $TotalCharge = &TotalCharge($App, \@AtomSite, \@Z1);
	if(abs($TotalCharge) > $ZEPS) {
		print "Error: Total charge ($TotalCharge) is greater than EPS=$ZEPS.\n";
		exit -1;
	}

print "Madelung potentials\n";
	my @MP;
	my ($pA, $pZ);
	if($UseP1) {
		$pA = \@AtomSite;
		$pZ = \@Z1;
	}
	else {
		$pA = \@AsymmetricAtomSite;
		$pZ = \@Z0;
	}
	for(my $is = 0 ; $is < @$pA ; $is++) {
		my $atom0          = $pA->[$is];
#	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             = $pZ->[$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;
}