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

use Utils;
use JFile;
use Crystal::CIF;
use Crystal::Crystal;
use Sci::Optimize;

my $nCase = 48;
my $EFormPath = "summary-FromAllBonds.csv";
my $OutCSV    = "fit-FromAllBonds.csv";
my $MaxR  = 2.5; # A

my $Opt = new Optimize;

#my $Method = "PDL::Simplex";
my $Method = "Amoeba::Simplex";
#my $Method = "ModifiedNewton";
#my $Method = "LinearOptimization"; # 線形パラメータのみ。

my $EPS         = 1.0e-7;
my $nMaxIter    = 3000;
my $iPrintLevel = 0;


#===================================
#===================================
my @EForm;
my @CoordinationArray;
for(my $im = 0 ; $im <= $nCase ; $im++) {
my $CIFPath    = "aIGZO-Dens608-PBE/VCRelax/a-IGZOm184-Dens608-MDRelaxed-final.cif";
my $OUTCARPath = "aIGZO-Dens608-PBE/SCF/OUTCAR";
if($im >= 1) {
	$CIFPath    = sprintf("aIGZO-Dens608-PBE-VO%02d/VCRelax/a-IGZOm184-Dens608-MDRelaxed-final.cif", $im);
	$OUTCARPath = sprintf("aIGZO-Dens608-PBE-VO%02d/SCF/OUTCAR", $im);
}

$EForm[$im] = &GetTotalEnergy($OUTCARPath);
print "$im: $EForm[$im] eV\n";

my $cif = new CIF;
$cif->Read($CIFPath) or die "Error: can not read [$CIFPath].\n";

my $Crystal = $cif->GetCCrystal();
my @latt = $Crystal->LatticeParameters();
print("Latt: ", join(', ', @latt), "\n");

$Crystal->ExpandCoordinates();
my @AtomTypeList = $Crystal->GetCAtomTypeList();
my @AtomSiteList = $Crystal->GetCAsymmetricAtomSiteList();
my @ExpandedAtomSiteList = $Crystal->GetCExpandedAtomSiteList();

my $site1 = $ExpandedAtomSiteList[0];
my $label1       = $site1->Label();
my $type1        = $site1->AtomName();
my ($x1,$y1,$z1) = $site1->Position(1);
my $site2 = $ExpandedAtomSiteList[74];
my $label2       = $site2->Label();
my $type2        = $site2->AtomName();
my ($x2,$y2,$z2) = $site2->Position(1);
#my $d = $Crystal->GetNearestInterAtomicDistance(0.313, 0.618, 0.877, 0.391, 0.461, 0.971);
my $d = $Crystal->GetNearestInterAtomicDistance($x1, $y1, $z1, $x2, $y2, $z2);
#print "d=$d: ($x1,$y1,$z1)-($x2,$y2,$z2): $label1-$label2\n";

my $c = 1;
my %Coordination;
for(my $i = 0 ; $i < @ExpandedAtomSiteList ; $i++) {
	my @List;

	my $site1 = $ExpandedAtomSiteList[$i];
	my $label1       = $site1->Label();
	my $type1        = $site1->AtomName();
	my ($x1,$y1,$z1) = $site1->Position(1);
	next if($type1 !~ /^O/);

	for(my $j = 0 ; $j < @ExpandedAtomSiteList ; $j++) {
		next if($i == $j);

		my $site2 = $ExpandedAtomSiteList[$j];
		my $label2       = $site2->Label();
		my $type2        = $site2->AtomName();
		my ($x2,$y2,$z2) = $site2->Position(1);

		my $d = $Crystal->GetNearestInterAtomicDistance($x1, $y1, $z1, $x2, $y2, $z2);
		if($d < $MaxR) {
			my %hash = (
				i      => $i,
				j      => $j,
				c      => $c,
				label1 => $label1,
				type1  => $type1,
				x1     => $x1,
				y1     => $y1,
				z1     => $z1,
				label2 => $label2,
				type2  => $type2,
				x2     => $x2,
				y2     => $y2,
				z2     => $z2,
				d      => $d,
			);
#			@List = (@List, \%hash);
			push(@List, \%hash);
			$c++;
		}
	}

	@List = sort { $a->{d} <=> $b->{d}; } @List;
	my ($nIn, $nGa, $nZn, $nO) = (0, 0, 0, 0);
	foreach my $l (@List) {
		if($l->{d} < $MaxR) {
			$nIn++ if($l->{label2} =~ /In/);
			$nGa++ if($l->{label2} =~ /Ga/);
			$nZn++ if($l->{label2} =~ /Zn/);
			$nO++  if($l->{label2} =~ /O/);
		}
	}
	my $nM = $nIn + $nGa + $nZn;

	$Coordination{$label1} = \@List;
#print "$label1: nIn,Ga,Zn,M = $nIn, $nGa, $nZn, $nM\n";
}
$CoordinationArray[$im] = \%Coordination;
}
#exit;
#===================================
#===================================

#my @Guess       = (  1,   1,   1,   0,   0,   0,   0,   0,   0,   0,   0,   0);
#my @Guess       = (-2.3, -3.0, -2.3,   0,   0,   0,   0,   0,   0,   0,   0,   0);
#my @Guess       = (-26.7, -20.23, -23.96, 0, 10.0, 7.5, 9.1, 0, 0, 0, 0, 0);
#my @Guess       = (-20.9, -32.8, -21.8, 0, 7.5, 13.4, 8.4, 0, 0, 0, 0, 0);
my @Guess       = (-10.5, -16.4, -10.9, 0, 3.8, 6.7, 4.2, 0, 0, 0, 0, 0);
#my @OptId       = (  1,   1,   1,   1,   1,   1,   1,   1,   1,   1,   1,   1);
#my @OptId       = (  1,   1,   1,   0,   1,   1,   1,   0,   1,   1,   1,   0);
my @OptId       = (  1,   1,   1,   0,   1,   1,   1,   0,   0,   0,   0,   0);
my @Scale       = (0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1);
#必要なければ、PDLDisplayFunc以降は指定する必要はない
my ($OptVars, $MinVal) = $Opt->Optimize($Method, \@Guess, \@OptId, \@Scale, 
				$EPS, $nMaxIter, $iPrintLevel,
				\&CalS2, sub {},
				sub { Optimize::BuildDifferentialMatrixes(@_); },
				);
print "Optimized at (",join(',',@{$OptVars}),")=$MinVal\n";

my (%EOM, %kOM, %kOM2);
($EOM{In},  $EOM{Ga},  $EOM{Zn},  $EOM{O}, 
 $kOM{In},  $kOM{Ga},  $kOM{Zn},  $kOM{O},
 $kOM2{In}, $kOM2{Ga}, $kOM2{Zn}, $kOM2{O}) = @$OptVars;

foreach my $atom qw(In Ga Zn O) {
	my $d;
	if($atom eq 'In') {
		$d = 2.217946;
	}
	elsif($atom eq 'Ga') {
		$d = 2.022151;
	}
	elsif($atom eq 'Zn') {
		$d = 2.049264;
	}
	printf "O-$atom: %6.3f + %6.3f*d + %6.3f*d^2 (%8.5f) eV\n", 
			$EOM{$atom}, $kOM{$atom}, $kOM2{$atom}, &CalBondEnergy($d, $EOM{$atom}, $kOM{$atom}, $kOM2{$atom});
#	printf "   %6.3f *exp(- %6.3f*d) (%8.5f) eV\n", 
#			$EOM{$atom}, abs($kOM{$atom}), &CalBondEnergy($d, $EOM{$atom}, $kOM{$atom}, $kOM2{$atom});
}

print "S2=$MinVal eV^2\n";
$MinVal = sqrt($MinVal / ($nCase+1));
print "S =$MinVal eV\n";

#rO 1.40
#rZn .6
#rIn .8
#rGa .47

#dO-In: 2.2
#dO-Ga: 1.87
#dO-Zn: 2.0

exit;

sub GetTotalEnergy
{
	my ($path) = @_;
	my $in = new JFile;
	$in->Open($path, "r") or die "Error: can not read [$path].\n";
	
	my $Etot;
	while(1) {
		my $line = $in->SkipTo("TOTEN\\s*=");
#print "l: $line";
		last if(!defined $line);

		($Etot) = ($line =~ /=\s*([\-\.\deEdD\+]+)/);
	}
	$in->Close();
	if(!defined $Etot) {
		die "Error: can not find Etot in [$path].\n";
	}
#print "Etot: $Etot\n";
#exit;
	return $Etot;
}

sub CalBondEnergy
{
	my ($d, $e0, $e1, $e2) = @_;
	return $e0 + $e1 * $d + $e2 * $d*$d;
#	return $e0 * exp(-abs($e1) * $d);
}

sub CalS2 {
	my ($pVars, $iPrintLevel) = @_;

	my (%EOM, %kOM, %kOM2);
	($EOM{In},  $EOM{Ga},  $EOM{Zn},  $EOM{O}, 
	 $kOM{In},  $kOM{Ga},  $kOM{Zn},  $kOM{O},
	 $kOM2{In}, $kOM2{Ga}, $kOM2{Zn}, $kOM2{O}) = @$pVars;

	my $S2 = 0.0;
	for(my $im = 0 ; $im < $nCase ; $im++) {
		my $pCoordination = $CoordinationArray[$im];

	my $Etot = 0.0;
	foreach my $label1 (sort keys %$pCoordination) {
		next if($label1 !~ /^O/);

		my $pList = $pCoordination->{$label1};
		foreach my $l (@$pList) {
			my ($atom) = ($l->{label2} =~ /(\D+)/);
			my $d = $l->{d};
			$Etot += &CalBondEnergy($d, $EOM{$atom}, $kOM{$atom}, $kOM2{$atom});
		}
	}
# 陽イオンだけから結合を数えているので、二重には数えていない
#	my $dE = 2*$EForm[$im] - $Etot;
	my $dE = $EForm[$im] - $Etot;
	$S2 += $dE*$dE;
printf("%02d: dE=%6.4g = %8.4g-%8.4g (%8.4g)\n", $im, $dE, $EForm[$im], $Etot, $dE);
	}

print "S2=$S2\n";


#printf "p: %4.3f, %4.3f, %4.3f, %4.3f - %4.3f, %4.3f, %4.3f, %4.3f - %4.3f, %4.3f, %4.3f, %4.3f: %4.2f\n",
#	 $EOM{In},  $EOM{Ga},  $EOM{Zn},  $EOM{O}, 
#	 $kOM{In},  $kOM{Ga},  $kOM{Zn},  $kOM{O},
#	 $kOM2{In}, $kOM2{Ga}, $kOM2{Zn}, $kOM2{O},$S2 ;

	return $S2;
}