package ChemicalReaction;
use Common;
@ISA = qw(Common);

use strict;
#use warnings;

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

use CSV;
use Crystal::AtomType;

#===============================================
# 変数取得関数
#===============================================
sub SetMode                    { my ($this,$m)=@_; return $this->{Mode} = $m; }
sub GetMode                    { return shift->{Mode}; }
sub SetDebugMode               { my ($this,$f)=@_; return $this->{DebugMode} = $f; }
sub GetDebugMode               { return shift->{DebugMode}; }

sub GetpCSVDBs                 { my ($this)=@_; return $this->{pCVDDBs}; };
sub GetCSVDB                   { my ($this, $index)=@_; return $this->{pCVDDBs}->[$index]; };
sub GetAllDB                   { my ($this)=@_; return $this->{pAllDB}; }
sub GetElementaryDB            { return shift->{pElementaryDB}; }
sub GetCompoundDB              { return shift->{pCompoundDB}; }
sub GetCompositionAllDB        { return shift->{pCompositionAllDB}; }
sub GetCompositionAtomDB       { return shift->{pCompositionAtomDB}; }
sub GetCompositionElementaryDB { return shift->{pCompositionElementaryDB}; }
sub GetCompositionCompoundDB   { return shift->{pCompositionCompoundDB}; }

sub GetDB
{
	my ($this) = @_;
	my $m = $this->GetMode();

	if($m eq 'CSVDB') {
		return $this->GetAllDB();
	}
	if($m eq 'Composition') {
		return $this->GetCompositionCompoundDB();
	}
	return $this->GetCompoundDB();
}
sub GetAtomDB
{
	my ($this) = @_;
	my $m = $this->GetMode();

	if($m eq 'CSVDB') {
		return $this->{pAtomDB};
	}
	if($m eq 'Composition') {
		return $this->GetCompositionAtomDB();
	}
	return $this->{pAtomDB};
}
sub GetDBListHash
{
	my ($this) = @_;
	my %db = (
		all                   => $this->{pAllDB},
		elementary            => $this->{pElementaryDB},
		atom                  => $this->{pAtomDB},
		compound              => $this->{pCompoundDB},
		compositionall        => $this->{pCompositionAllDB},
		compositionelementary => $this->{pCompositionElementaryDB},
		compositionatom       => $this->{pCompositionAtomDB},
		compositioncompound   => $this->{pCompositionCompoundDB},
		);
	for(my $i = 0 ; ; $i++) {
		my $db = $this->GetCSVDB($i);
		last if(!defined $db);
		my $i1 = $i+1;
		$db{"csv$i1"} = $db;
	}
	return %db;
}

sub GetDBHash
{
	my ($this, $Compound, $pDB) = @_;
	my $m = $this->GetMode();

	my $c = $Compound;
	if($m eq 'Composition') {
		$c =~ s/^([a-zA-Z])+\-//;
		$c = $this->SortChemicalFormula($Compound);
	}
#print "GetDBHash for [$c][mode=$m]\n";

	if($m eq 'CSVDB') {
		for(my $i = 0 ; ; $i++) {
			my $pDB = $this->GetCSVDB($i);
			last if(!$pDB);
			my $pHash = $pDB->HashByXDataKey($c);
			return $pHash if($pHash and $pHash->{Status} eq '');
		}
		return undef;
	}

	if(!$pDB) {
		if($m eq 'Composition') {
			$pDB = $this->GetCompositionCompoundDB();
		}
		else {
			$pDB = $this->GetDB();
		}
	}
#foreach my $key (sort keys %$pDB) {
#print "key: $key ($pDB->{$key}{Type}): $pDB->{$key}{Emol}: $pDB->{$key}{Status}\n";
#}
#print "pDB=$pDB->{$c} for [$c]\n";
	return $pDB->{$c} if($pDB->{$c} and $pDB->{$c}{Status} eq '');
	return undef;
}

sub GetEmol
{
	my ($this, $Compound, $pDB) = @_;
#print "GetEmol for [$Compound] [mode=", $this->GetMode(), "]\n";
	$pDB = $this->GetDBHash($Compound, $pDB) if(!defined $pDB);
	return undef if(!defined $pDB);
	return $pDB->{Emol};
}

sub InitialReaction { my ($this)=@_; return $this->{Initial} ; };
sub FinalReaction   { my ($this)=@_; return $this->{Final} ; };
sub Reagents        { my ($this)=@_; return $this->{Initial} ; };
sub Products        { my ($this)=@_; return $this->{Final} ; };
sub ReagentsHash    { my ($this)=@_; return $this->{ReagentsHash} ; };
sub ProductsHash    { my ($this)=@_; return $this->{ProductsHash} ; };
sub pCompounds {
	my ($this, $key) = @_;
	if($key =~ /^P/i) {
		return $this->{ProductsHash}->{pCompounds};
	}
	return $this->{ReagentsHash}->{pCompounds};
}
sub nCompound {
	my ($this, $key, $compound) = @_;

	my $pHash;
	if($key =~ /^P/i) {
		$pHash = $this->ProductsHash();
	}
	else {
		$pHash = $this->ReagentsHash();
	}
	return $pHash->{$compound}{n};
}
sub ElementsArray {
	my ($this, $key) = @_;
	return $this->pElements($key);
}
sub pElements {
	my ($this, $key) = @_;
	if($key =~ /^P/i) {
		return $this->{ProductsHash}->{pElements};
	}
	return $this->{ReagentsHash}->{pElements};
}
sub nElement {
	my ($this, $key, $element) = @_;

	my $pHash;
	if($key =~ /^P/i) {
		$pHash = $this->ProductsHash();
	}
	else {
		$pHash = $this->ReagentsHash();
	}
	return $pHash->{$element}{nElement};
}

sub ListCompounds
{
	my ($this) = @_;
	my @a = $this->SearchCSVDBHit(".*");
	for(my $i = 0 ; $i < @a ; $i++) {
		print "$i: $a[$i]\n";
	}
}

sub ShowSpeculatedSpecies
{
	my ($this, $Compound) = @_;
	my @Hit = $this->SearchCSVDBHit("^${Compound}[\\-:]");
	for(my $i = 0 ; $i < @Hit ; $i++) {
		print "   $Hit[$i]\n";
	}
}

sub nElementByArray
{
	my ($this, $pCompounds, $pComposition, $IsPrint) = @_;

	my %nFinalElement;
	for(my $i = 0 ; $i < @$pComposition ; $i++) {
		my $c = $pCompounds->[$i];
#print "  $c: $pMinComposition->[$i]\n" if($IsPrint);
		my ($pElement, $pnElement) = $this->CompoundToElements($c);
		for(my $j = 0 ; $j < @$pElement ; $j++) {
			my $e = $pElement->[$j];
			my $n = $pnElement->[$j];
			$nFinalElement{$e} += $n * $pComposition->[$i];
		}
	}
	return \%nFinalElement;
}


sub SearchCSVDBHit
{
	my ($this, $RegExp) = @_;

	my @Hit;
	for(my $i = 0 ; ; $i++) {
		my $DB = $this->GetCSVDB($i);
		last if(!$DB);

		my $pX = $DB->{pX};
		$pX = $DB->GetData(0) if(!defined $pX);
		my $pS = $DB->GetData("Status");
#print "ps=$pS\n";
		my $nData = @$pX;
		for(my $i = 0 ; $i < $nData ; $i++) {
			if($pX->[$i] =~ /$RegExp/) {
				next if($pS and $pS->[$i] ne '');
				push(@Hit, $pX->[$i]);
			}
		}
	}
	return @Hit;
}

#===============================================
# コンストラクタ、デストラクタ
#===============================================
sub new
{
	my ($module, $pDBFiles, $BaseDBDir, $mode) = @_;
	$mode = "Compound" if(defined $mode and $mode eq '');
	my $this = {};
	bless $this;

	$this->Initialize($mode);
	$this->SetMode($mode);

	$this->ReadDBFiles($pDBFiles, $BaseDBDir) if(defined $pDBFiles and $pDBFiles ne '');
	return $this;
}

sub DESTROY
{
	my ($this) = @_;
	$this->SUPER::DESTROY(@_);
}

sub Initialize
{
	my ($this, $InitializeDB) = @_;
	$InitializeDB = 1 if(!defined $InitializeDB);

	$this->ClearHashes();
	if($InitializeDB) {
		delete $this->{pCVDDBs};
	}
	$this->SetDebugMode(0);
}

sub ClearHashes
{
	my ($this) = @_;
	$this->{pReagents}    = [];
	$this->{Reagents}     = [];
	$this->{ReagentsHash} = {};
	$this->{pProducts}    = [];
	$this->{Products}     = [];
	$this->{ProductsHash} = {};
}

#===============================================
# 一般関数
#===============================================
sub MolecularWeight
{
	my ($this, $Compound) = @_;

	my $R = new ChemicalReaction;
	$R->AnalyzeAReaction('Reagents', $Compound);
	my $pHash = $R->ReagentsHash();
	my $pElements = $R->ElementsArray("Reagents");
	my $w = 0.0;
	for(my $i = 0 ; $i < @$pElements ; $i++) {
		my $e = $pElements->[$i];
		my $n = $R->nElement("Reagents", $e);
#print "$e: $n\n";
		my $pAtomDB = AtomType::GetAtomDBHash($e);
#print "  Z=$pAtomDB->{Z}: M=$pAtomDB->{mass}\n";
		$w += $n * $pAtomDB->{mass};
	}
	return $w;
}

sub IsIncluded
{
	my ($this, $key, $pArray, $CaseSensitive) = @_;
	return Utils::IsIncludedInArray($key, $pArray, $CaseSensitive);
}

sub PrintDBData
{
	my ($this, %hash) = @_;
	my $StopByError = $hash{StopByError};
	my $db = $this->{pAllDB};
	foreach my $key (sort keys %$db) {
		if($StopByError and !defined $db->{$key}{SortedComposition}) {
			print "Error: Can not find SortedComposition: $key: $db->{$key}{Compound}\n";
			exit;
		}
		if($StopByError and !defined $db->{$key}{Emol}) {
			print "Error: Can not find Emol: $key: $db->{$key}{Compound}\n";
			exit;
		}
		print "$key [$db->{$key}{Type}]: $db->{$key}{Compound}: $db->{$key}{SortedComposition}: "
		     ."$db->{$key}{Etot} / $db->{$key}{Z} = $db->{$key}{Emol} eV\n";
	}
}

sub PrintElements
{
	my ($this, $key) = @_;

	my ($pReaction, $pList, $pHash);
	if($key =~ /^P/i or $key =~ /^F/i) {
		$pReaction = $this->{pProducts};
		$pList     = $this->{Products};
		$pHash     = $this->{ProductsHash};
	}
	else {
		$pReaction = $this->{pReagents};
		$pList     = $this->{Reagents};
		$pHash     = $this->{ReagentsHash};
	}

	print "$key:\n";
	for(my $i = 0 ; $i < @$pReaction ; $i++) {
		my $n        = $pList->[$i]{n};
		my $compound = $pList->[$i]{Compound};
		print "  $n * $compound: ";

		my $IsFirst = 1;
		my $pElements = $pHash->{pElements};
		foreach my $element (sort @$pElements) {
#			my $ne        = $pHash->{$element}{ne};
			my $ne        = $pHash->{$compound}{$element}{ne};;
			next if($ne eq '');
			my $nElement  = $pHash->{$element}{nElement};

			print " + " if(!$IsFirst);
			print "$ne * $element";
			$IsFirst = 0;
		}
		print "\n";
	}
}

sub PrintnElements
{
	my ($this, $key) = @_;

	my $pnElementHash;
	if($key =~ /^P/i or $key =~ /^F/i) {
		$pnElementHash = $this->{ProductsHash}{nElement};
	}
	else {
		$pnElementHash = $this->{ReagentsHash}{nElement};
	}

	my $IsFirst = 1;
	print "nElements in $key: ";
	foreach my $e (sort keys %$pnElementHash) {
		print " + " if(!$IsFirst);
		print "$pnElementHash->{$e} * $e";
		$IsFirst = 0;
	}
	print "\n";
}

sub ReadDBFiles
{
	my ($this, $pDBFiles, $BaseDBDir) = @_;

	my $pDBs = $this->{pCVDDBs};
	my @DBs = (defined $pDBs)? @$pDBs : ();
	my $nDBs = @DBs;
	for(my $i = 0 ; $i < @$pDBFiles ; $i++) {
#print "f=$pDBFiles->[$i]\n";
		$DBs[$nDBs+$i] = $this->ReadDBFile($pDBFiles->[$i], $BaseDBDir);
		if(!$DBs[$nDBs+$i]) {
#print "Warning: Can not read $pDBFiles->[$i]\n";
			return undef;
		}
#print "Read $pDBFiles->[$i]\n";
	}
	$this->{pCVDDBs} = \@DBs;
	return $this->{pCVDDBs};
}

sub ReadDBFile
{
	my ($this, $DBFile, $BaseDBDir, $IsPrint) = @_;
print "ReadDBFile [$DBFile]\n";
	$IsPrint = 0 if(!defined $IsPrint);
	$this->{BaseDBDir} = $BaseDBDir if(defined $BaseDBDir);

	my $csv = new CSV;
	$DBFile = Deps::MakePath($BaseDBDir, $DBFile, 0) if(defined $BaseDBDir);
	if(!$csv->Read($DBFile, 0, 1)) {
print "Error: $!: Can not read [$DBFile].\n";
		return undef;
	}
	my $pX = $csv->GetXData("Compound");
	if(!$pX) {
		print "Error: Can not find 'Compound' tag in [$DBFile].\n";
	}
	for(my $i = 0 ; $i < @$pX ; $i++) {
		$pX->[$i] =~ s/\s+//g;
		$pX->[$i] =~ s/Atom/atom/gi;
		$pX->[$i] =~ s/\((AFM|FM|NM|NP|SP|NM|atom)\)/:$1/gi;
		$pX->[$i] =~ s/([A-Z].*)-/$1:/g;
	}

	$this->{pAllDB}                   = {} if(!defined $this->{pAllDB});
	$this->{pAtomDB}                  = {} if(!defined $this->{pAtomDB});
	$this->{pElementaryDB}            = {} if(!defined $this->{pElementaryDB});
	$this->{pCompoundDB}              = {} if(!defined $this->{pCompoundDB});
	$this->{pCompositionAllDB}        = {} if(!defined $this->{pCompositionAllDB});
	$this->{pCompositionAtomDB}       = {} if(!defined $this->{pCompositionAtomDB});
	$this->{pCompositionElementaryDB} = {} if(!defined $this->{pCompositionElementaryDB});
	$this->{pCompositionCompoundDB}   = {} if(!defined $this->{pCompositionCompoundDB});
	my $pAll        = $this->{pAllDB};
	my $pAtom       = $this->{pAtomDB};
	my $pElementary = $this->{pElementaryDB};
	my $pCompound   = $this->{pCompoundDB};
	my $pCompositionAll        = $this->{pCompositionAllDB};
	my $pCompositionAtom       = $this->{pCompositionAtomDB};
	my $pCompositionElementary = $this->{pCompositionElementaryDB};
	my $pCompositionCompound   = $this->{pCompositionCompoundDB};
	for(my $i = 0 ; $i < @$pX ; $i++) {
		my $pHash = $csv->GetHashData($i);
		if($pHash->{Status} ne '') {
print "Error: Bad status for [$pHash->{Compound}] in [$pHash->{Ref}]: $pHash->{Status}\n";
			next;
		}
		$pHash->{Z} = 1 if(!defined $pHash->{Z} or $pHash->{Z} == 0.0);
		if(!defined $pHash->{Emol}) {
			$pHash->{Emol} = $pHash->{Etot} / $pHash->{Z};
		}
		$pHash->{SortedComposition} = $pHash->{"SortedComposition/Z"};
		if($pHash->{SortedComposition} eq '') {
			my $composition = $pHash->{"Composition/Z"};
			$pHash->{SortedComposition}
				= $this->SortChemicalFormula($composition);
		}
		$pHash->{Compound}        =~ s/\s//g;
		$pHash->{Composition}     =~ s/\s//g;
		$pHash->{"Composition/Z"} =~ s/\s//g;

		my $composition = $pHash->{SortedComposition};
		my $name        = $pHash->{Compound};
#		$name =~ s/[:\-\s].*$//;
#print "$composition => $name\n";
#print "$pHash->{Compound}: $pHash->{Etot} eV (Z=$pHash->{Z}) (Type=$pHash->{Type})\n";
		if($pHash->{Type} eq 'Elementary') {
			$name =~ s/[\d\s].*$//;
			if(!defined $pCompound->{$name} or $pCompound->{$name}->{Emol} > $pHash->{Emol}) {
				$pCompound->{$name} = $pHash;
			}
			if(!defined $pCompositionCompound->{$composition} 
			   or $pCompositionCompound->{$composition}->{Emol} > $pHash->{Emol}) {
				$pCompositionCompound->{$composition} = $pHash;
			}
			if(!defined $pElementary->{$name} or $pElementary->{$name}->{Emol} > $pHash->{Emol}) {
				$pElementary->{$name} = $pHash 
			}
			if(!defined $pCompositionElementary->{$composition} 
			   or $pCompositionElementary->{$name}->{Emol} > $pHash->{Emol}) {
				$pCompositionElementary->{$composition} = $pHash;
			}
print "  Elementary: $name: $pHash->{Emol} eV\n" if($IsPrint);
		}
		elsif($pHash->{Type} eq 'Atom') {
			$name =~ s/[\d\s].*$//;
			if(!defined $pAtom->{$name} or $pAtom->{$name}->{Emol} > $pHash->{Emol}) {
				$pAtom->{$name} = $pHash;
			}
			if(!defined $pCompositionAtom->{$name} or $pCompositionAtom->{$name}->{Emol} > $pHash->{Emol}) {
				$pCompositionAtom->{$composition} = $pHash;
			}
print "  Atom: $name: $pHash->{Emol} eV\n" if($IsPrint);
		}
		else {
			if($pHash->{Type} ne 'Compound') {
print "  Warning: Invalid Type [$pHash->{Type}] (should be Atom/Elementary/Compound): $name: $pHash->{Emol} eV\n";# if($IsPrint);
print "           Assume Compound\n";# if($IsPrint);
			}
			$pHash->{Type} = 'Compound';
#		elsif($pHash->{Type} eq 'Compound') {
			if(!defined $pCompound->{$name} or $pCompound->{$name}->{Emol} > $pHash->{Emol}) {
				$pCompound->{$name} = $pHash;
			}
			if(!defined $pCompositionCompound->{$composition} 
			   or $pCompositionCompound->{$composition}->{Emol} > $pHash->{Emol}) {
				$pCompositionCompound->{$composition} = $pHash;
			}
print "  Compound: $name: $pHash->{Emol} eV\n" if($IsPrint);
		}

		if(!defined $pAll->{$name} or $pAll->{$name}->{Emol} > $pHash->{Emol}) {
			$pAll->{$name} = $pHash;
		}
		if(!defined $pCompositionAll->{$name} or $pCompositionAll->{$name}->{Emol} > $pHash->{Emol}) {
			$pCompositionAll->{$composition} = $pHash;
		}
	}

	return $csv;
}

sub SplitReaction
{
	my ($this, $Reaction) = @_;
	return Utils::Split("=>?", $Reaction);
}

sub BuildReactionFromArrays
{
	my ($this, $pCompounds, $pComposition, $format, $min, $Sort) = @_;
	$format = "%0.3f" if($format eq '');
	$min = 1.0e-3 if(!defined $min);
	$Sort = 0 if(!defined $Sort);

	my $s = '';
	for(my $i = 0 ; $i < @$pCompounds ; $i++) {
		next if($pComposition->[$i] < $min);
		my $n = sprintf($format, $pComposition->[$i]);
		$n =~ s/[\.0]+$//;
		$s .= $n . $pCompounds->[$i] . "+";
	}
	if($Sort) {
		my @a = sort { return $a cmp $b; } Utils::Split("\\+", $s);
		$s = join('+', @a);
	}
	$s =~ s/\+$//;
	return $s;
}


sub ReactionToHTML
{
	my ($this, $reaction) = @_;

	my ($r, $p) = $this->SplitReaction($reaction);
	my $rh = $this->AReactionToHTML($r);
	return $rh if(!defined $p);;

	my $ph = $this->AReactionToHTML($p);
	return "$rh=>$ph";
}

sub AReactionToHTML
{
	my ($this, $reaction) = @_;

	my ($pChemicals, $pCompound, $pnCompound) = $this->ReactionToCompounds($reaction);
	my $r = '';
	for(my $i = 0 ; $i < @$pChemicals ; $i++) {
		my $c = $pCompound->[$i];
		my $n = $pnCompound->[$i];
		$c =~ s/([\d\.\/]+)/\<sub>$1\<\/sub\>/g;
		$n = '' if($n == 1);
		if($r eq '') {
			$r = "$n$c";
		}
		elsif($n == '') {
			$r .= "+$c";
		}
		elsif($n == -1) {
			$r .= "-$c";
		}
		elsif($n < -1) {
			$r .= "$n$c";
		}
		else {
			$r .= "+$n$c";
		}
	}
	$r =~ s/\+\-/\-/g;
	return $r;
}

sub SortChemicalFormula
{
	my ($this, $formula) = @_;

	my ($pe, $pn) = $this->CompoundToElements($formula);
	my %n;
	for(my $i = 0 ; $i < @$pe ; $i++) {
		$n{$pe->[$i]} += $pn->[$i];
	}
	my @a = sort { $a cmp $b; } keys %n;
	my $s = '';
	for(my $i = 0 ; $i < @a ; $i++) {
		$s .= "$a[$i]$n{$a[$i]}";
	}
	return $s;
#	return join('', @a);
}

sub NormalizeChemicalComposition
{
	my ($this, $pnAtom) = @_;

	my %AtomCount = %$pnAtom;

	my @Names = keys %AtomCount;
	my $Composition = '';
	my $MinN = 1000000;
	for(my $i = 0 ; $i < @Names ; $i++) {
		my $name = $Names[$i];
		next if($name =~ /^(O|S|Se|Te|F|Cl|I|N|As)$/i);
		my $n = $AtomCount{$name};
		$MinN = $n if($MinN > $n);
		$n = '' if($n == 1);
		$Composition = "$Composition $name$n";
#print "i=$i: [$n] [$MinN] [$name]\n";
	}

	for(my $i = 0 ; $i < @Names ; $i++) {
		my $name = $Names[$i];
		next unless($name =~ /^(O|S|Se|Te|F|Cl|I|N|As)$/i);
		my $n = $AtomCount{$name};
		$MinN = $n if($MinN > $n);
		$n = '' if($n == 1);
		$Composition = "$Composition $name$n";
	}
	$Composition =~ s/^\s*//;
	my $OriginalComposition = $Composition;

	my @elements = Algorism::Factorize($MinN);
	my $fac = 1;
#print "$MinN = ";
	for(my $i = 0 ; $i < @elements ; $i++) {
		my $num = $elements[$i];
#print " * "  if($i > 0);
#print "$num";

		my $IsCommonFactor = 1;
		for(my $j = 0 ; $j < @Names ; $j++) {
			my $name = $Names[$j];
			my $n = $AtomCount{$name};
			unless(Algorism::IsFactor($n, $num*$fac)) {
				$IsCommonFactor = 0;
			}
		}
		$fac *= $num if($IsCommonFactor);
#print("fac[$i]=$fac\n");
	}

	$Composition = '';
	for(my $i = 0 ; $i < @Names ; $i++) {
		my $name = $Names[$i];
		next if($name =~ /^(O|S|Se|Te|F|Cl|I|N|As|)$/i);
		my $n = $AtomCount{$name} / $fac;
		$n = '' if($n == 1);
		$Composition = "$Composition $name$n";
	}
	for(my $i = 0 ; $i < @Names ; $i++) {
		my $name = $Names[$i];
		next unless($name =~ /^(O|S|Se|Te|F|Cl|I|N|As|)$/i);
		my $n = $AtomCount{$name} / $fac;
		$n = '' if($n == 1);
		$Composition = "$Composition $name$n";
	}
	$Composition =~ s/^\s*//;
#print("C: [$Composition]\n");

	return ($OriginalComposition, $Composition);
}

sub ReactionToCompounds
{
	my ($this, $reaction) = @_;

	if($this->GetMode() eq 'Composition') {
		$reaction =~ s/\-/+-/g;
	}
	else {
		$reaction =~ s/([^\+])\-(\d)/$1+-$2/g;
	}
	my @Chemicals = Utils::Split("[\\+]", $reaction);
	my (@Compounds, @n);
	for(my $i = 0 ; $i < @Chemicals ; $i++) {
		($n[$i], $Compounds[$i]) = ($Chemicals[$i] =~ /^([\-\d\.\/]*)(\D.*)$/);
		$n[$i] =  1 if($n[$i] eq '');
		$n[$i] = -1 if($n[$i] eq '-');
		$n[$i] = eval($n[$i]);
#print "[$reaction]$i: $Compounds[$i]: $n[$i]<br>\n";
	}
	return (\@Chemicals, \@Compounds, \@n);
}

sub CompoundToElements
{
	my ($this, $Compound, $Debug) = @_;
	$Debug = $this->GetDebugMode() if(!defined $Debug);
#$Debug=1;

	$Compound =~ s/:.*$//;
	$Compound =~ s/\-.*$//;
	$Compound =~ s/\(.*(AFM|FM|NM|NP|SP|NM|atom).*\)$//i;
	while(1) {
#print "C: $Compound<br>\n" if($Debug);
		my ($head, $par, $n, $tail) = ($Compound =~ /^(.*?)[\(\[]([^\)\(\]\[]*)[\)\]](\d*)(.*)$/);
#		my ($head, $par, $n, $tail) = ($Compound =~ /^(.*?)[\(\[]([^\)\(\]\[]*)[\)\]]([\d\.\/]*)(.*)$/);
		last if(!defined $n);
#print "  [$head][$par][$n][$tail]<br>\n" if($Debug);
		my (@a) = Utils::FindAll("[A-Z][a-z]*\\d*", $par);
		my ($pe, $pn) = $this->CompoundToElements($par);
		my $s = '';
		for(my $i = 0 ; $i < @$pe ; $i++) {
			my $n2 = $pn->[$i] * $n;
			$s .= "$pe->[$i]$n2";
		}
		$Compound = "$head$s$tail";
	}

	my (@a) = Utils::FindAll("[A-Z][a-z]*[\\d\\\\/.]*", $Compound);
#print "a: ", join(',', @a), "<br>\n";

#	my (@Element, @n);
	my %nElement;
	for(my $i = 0 ; $i < @a ; $i++) {
#		($Element[$i], $n[$i]) = ($a[$i] =~ /^([A-Z][a-z]*)(\d*)$/);
#		$n[$i] = 1 if(!defined $n[$i] or $n[$i] eq '');
		my ($e, $n) = ($a[$i] =~ /^([A-Z][a-z]*)([\d\.\/]*)$/);
		$n = 1 if(!defined $n or $n eq '');
		$nElement{$e} += eval($n);
	}
	return ([keys %nElement], [values %nElement]);
#	return (\@Element, \@n);
}

sub RemoveCompounds
{
	my ($this, $pSource, $pRemove) = @_;
#for(my $i = 0 ; $i < @$pSource ; $i++) {
#print "$i: $pSource->[$i]\n";
#}
#exit;
	for(my $i = 0 ; $i < @$pRemove ; $i++) {
		$pRemove->[$i] = $this->SortChemicalFormula($pRemove->[$i]);
#print "R: $pRemove->[$i]\n";
	}
	my %s;
	for(my $i = 0 ; $i < @$pSource ; $i++) {
		my $IsHit = 0;
		for(my $j = 0 ; $j < @$pRemove ; $j++) {
			if($pSource->[$i] eq $pRemove->[$j]) {
print "RemoveCompounds: [$pRemove->[$j]] is removed.\n";
				$IsHit = 1;
				last;
			}
		}
		if(!$IsHit) {
			$s{$pSource->[$i]} = 1;
		}
	}
	return sort keys %s;
}

sub ExtractPossibleCompounds
{
	my ($this, $pElements, $pCompounds, $pCompoundDB) = @_;
	$pElements   = $this->pElements("Reagents") if(!defined $pElements);
	$pCompoundDB = $this->GetCompoundDB()       if(!defined $pCompoundDB);
	$pCompounds  = [keys %$pCompoundDB]         if(!defined $pCompounds);

	my @PC;
	foreach my $c (@$pCompounds) {
		my ($pe, $pn) = $this->CompoundToElements($c);
		my $passed = 1;
		for(my $i = 0 ; $i < @$pe ; $i++) {
			if($this->IsIncluded($pe->[$i], $pElements)) {
			}
			else {
				$passed = 0;
				last;
			}
		}
		if($passed) {
			push(@PC, $c);
		}
	}

	$this->{pPossibleCompounds} = \@PC;

	return @PC;
}

sub AnalyzeAReaction
{
	my ($this, $key, $reaction) = @_;

#print("Analyaze: $reaction<br>\n");
	my ($pChemicals, $pCompound, $pnCompound) = $this->ReactionToCompounds($reaction);
#for(my $i = 0 ; $i < @$pChemicals ; $i++) {
#	print("$i: $pChemicals->[$i]: $pCompound->[$i]: $pnCompound->[$i]<br>\n");
#}

	my ($pReaction, $pList, $pHash);
	if($key =~ /^P/i or $key =~ /^F/i) {
		$this->{pProducts} = $pChemicals;
		$pList     = $this->{Products};
		$pHash     = $this->{ProductsHash};
	}
	else {
		$this->{pReagents} = $pChemicals;
		$pList     = $this->{Reagents};
		$pHash     = $this->{ReagentsHash};
	}

#print "[$key]:<br>\n";
	my %ElementsHash;
	for(my $i = 0 ; $i < @$pCompound ; $i++) {
#print "R[$i]=$pChemicals->[$i]<br>\n";
		my $compound = $pCompound->[$i];
		my $n        = $pnCompound->[$i];
#print ">>$i: $compound: $n<br>\n";

		$pList->[$i]{n}        = $n;
		$pList->[$i]{Compound} = $compound;
		$pHash->{$compound}{n} = $n;
#print "  ${n}[$compound]: ";

		my ($pElement, $pnElement) = $this->CompoundToElements($compound);
		for(my $i = 0 ; $i < @$pElement ; $i++) {
			my $element = $pElement->[$i];
			my $ne      = $pnElement->[$i];
#print "  $ne * $element + ";
			$ElementsHash{$element}++;
			$pHash->{$element}{ne} = $ne;
			my $N = $n * $ne;
			$pHash->{$compound}{$element}{ne} += $ne;
#print "0:e=$element: n=$n: ne=$ne: N=$N  Ntot=$pHash->{$element}{nElement}<br>\n";
			$pHash->{$element}{nElement}      += $N;
			$pHash->{nElement}{$element}      += $N;
#print "1:e=$element: n=$n: ne=$ne: N=$N  Ntot=$pHash->{$element}{nElement}<br>\n";
		}
#print "<br>\n";
	}
#print "E: ", join(', ', keys %ElementsHash), "\n";
	$pHash->{pCompounds} = $pCompound;
	$pHash->{pElements} = [keys %ElementsHash];
}

sub Analyze
{
	my ($this, $Reaction, $IsPrint) = @_;
	$IsPrint = 1 if(!defined $IsPrint);

	if($Reaction eq 'list') {
		$this->ListCompounds();
		return ();
	}
	elsif($Reaction eq 'exit' or $Reaction eq 'quit' ) {
		exit;
	}

	$this->{Reaction} = $Reaction;
	($this->{Initial}, $this->{Final}) = $this->SplitReaction($Reaction);

#	$this->{Initial} =~ s/\-/+\-/g;
#	$this->{Final}   =~ s/\-/+\-/g;
	$this->AnalyzeAReaction("Reagents", $this->{Initial});
	$this->PrintElements("Reagents") if($IsPrint);
	$this->AnalyzeAReaction("Products", $this->{Final});
	$this->PrintElements("Products") if($IsPrint);
}

sub CheckMassBalance
{
	my ($this, $IsPrint) = @_;
	$IsPrint = 0 if(!defined $IsPrint);

	$this->PrintnElements("Reagents") if($IsPrint);
	$this->PrintnElements("Products") if($IsPrint);

	my $pnInitialElementHash = $this->{ReagentsHash}{nElement};
	my $pnFinalElementHash   = $this->{ProductsHash}{nElement};
	foreach my $e (sort keys %$pnInitialElementHash) {
		if($pnInitialElementHash->{$e} != $pnFinalElementHash->{$e}) {
print("Error: Mass balance violated for [$e]: "
	."$pnInitialElementHash->{$e} != $pnFinalElementHash->{$e}\n") if($IsPrint);
			return 0;
		}
	}
	foreach my $e (sort keys %$pnFinalElementHash) {
		if($pnInitialElementHash->{$e} != $pnFinalElementHash->{$e}) {
print("Error: Mass balance violated for [$e]: "
	."$pnInitialElementHash->{$e} != $pnFinalElementHash->{$e}\n") if($IsPrint);
			return 0;
		}
	}
	return 1;
}

sub MakePossibleComposition
{
	my ($this, $pReagentElements, $pCompounds, $pCompoundDB, $MaxNInReaction, $IsPrint) = @_;
	$pCompoundDB = $this->GetCompoundDB() if(!defined $pCompoundDB);

	my @Composition;
	my %nPElements;
#
# 試行的に組成比をつくる
#
	for(my $i = 0 ; $i < @$pCompounds ; $i++) {
#		my $c    = $pCompounds->[$i];
		my $c    = $this->SortChemicalFormula($pCompounds->[$i]);
		my $Type = $pCompoundDB->{$c}->{Type};
#print "$c: $Type\n";
		my $cp = 0.0;
		if($Type ne 'Elementary') {
			$cp = rand($MaxNInReaction);
		}
		$Composition[$i] = $cp;
		my ($pe, $pn) = $this->CompoundToElements($c);
		for(my $j = 0 ; $j < @$pe ; $j++) {
			$nPElements{$pe->[$j]} += $pn->[$j] * $Composition[$i];
		}
	}

	return $this->RepairComposition(\@Composition, $pReagentElements, $pCompounds, $pCompoundDB, $IsPrint);
}

sub NormalizeComopositionByArray
{
	my ($this, $pCompounds, $pComposition, $format) = @_;
	$format = "%f" if($format eq '');

	my %n;
	for(my $i = 0 ; $i < @$pCompounds ; $i++) {
		my ($pElement, $pnElement) = $this->CompoundToElements($pCompounds->[$i]);
		for(my $j = 0 ; $j < @$pElement ; $j++) {
			$n{$pElement->[$j]} += $pComposition->[$i] * $pnElement->[$j];
		}
	}
	my $s = '';
	foreach my $e (sort keys %n) {
		$s .= sprintf("$e$format", $n{$e});
	}
	return $s;
}

sub RepairComposition
{
	my ($this, $pComposition, $pReagentElements, $pCompounds, $pCompoundDB, $IsPrint) = @_;
my $Debug = 0;
	$pCompoundDB = $this->GetCompoundDB() if(!defined $pCompoundDB);

	my %nPElements;
	for(my $i = 0 ; $i < @$pCompounds ; $i++) {
#		my $c    = $pCompounds->[$i];
		my $c    = $this->SortChemicalFormula($pCompounds->[$i]);
		my $Type = $pCompoundDB->{$c}->{Type};

		if($pComposition->[$i] < 0.0) {
			$pComposition->[$i] = 0.0;
		}
		if($Type eq 'Elementary') {
			$pComposition->[$i] = 0.0;
		}
		my ($pe, $pn) = $this->CompoundToElements($c);
		for(my $j = 0 ; $j < @$pe ; $j++) {
			$nPElements{$pe->[$j]} += $pn->[$j] * $pComposition->[$i];
		}
	}
#print("rev: ", join(',', @$pComposition), "\n") if($Debug);

#
# 一番多い元素の組成が反応物の組成を越えないように係数 $K を決める
#
	my $K = 1.0;
	for(my $i = 0 ; $i < @$pReagentElements ; $i++) {
		my $e = $pReagentElements->[$i];
		my $Type = $pCompoundDB->{$e}->{Type};
		my $n = $nPElements{$e};
		my $nElementInReaction = $this->nElement("Reagents", $e);
print "$e: $nElementInReaction\n" if($Debug);
		my $k = ($n < 1e-10)? 1.0 : $nElementInReaction / $n;
		$K = $k if($K > $k); 
	}
#print "K=$K\n";
	for(my $i = 0 ; $i <  @$pCompounds ; $i++) {
		$pComposition->[$i] *= $K;
	}
	foreach my $e (keys %nPElements) {
		$nPElements{$e} *= $K;
	}
#print("rev3: ", join(',', @$pComposition), "\n") if($Debug);
#	return $pComposition;
#
# 不足元素分をElementaryで補う
#
	for(my $i = 0 ; $i <  @$pCompounds ; $i++) {
#		my $c    = $pCompounds->[$i];
		my $c    = $this->SortChemicalFormula($pCompounds->[$i]);
		my $e = $c;
		$e =~ s/\d+$//;
		my $Type = $pCompoundDB->{$c}{Type};
#print "$i: $Type: $c\n";
		next if($Type ne 'Elementary');
		my $nElementInReaction = $this->nElement("Reagents", $e);
		my $nDiff = $nElementInReaction - $nPElements{$e};
if(!defined $nElementInReaction) {
	print "Error in RepairComposition: Invalid nElementInReaction [$nElementInReaction]\n";
	exit;
}
if(!defined $nPElements{$e}) {
	print "Error in RepairComposition: Can not get nElement [ele=$e][n=$nPElements{$e}]\n";
	exit;
}
#print "$i: $Type: $e: diff=$nDiff ($nElementInReaction - $nPElements{$e})\n";
		$pComposition->[$i] += $nDiff;
	}
	my $TotalComposition = 0.0;
	for(my $i = 0 ; $i <  @$pCompounds ; $i++) {
		$TotalComposition += $pComposition->[$i];
	}
#print "C3: ", $this->NormalizeComopositionByArray($pCompounds, $pComposition, "%4.2f"), "\n";
print "TotalComposition=$TotalComposition (K=$K)\n" if($IsPrint);

	return $pComposition
}

sub CalCohesiveEnergy
{
	my ($this, $Compound, $IsPrint, $ErrorStop, $Debug) = @_;
$Debug = $this->GetDebugMode() if(!defined $Debug);
print "Ecoh for [$Compound]\n" if($Debug);
	$IsPrint   = 0 if(!defined $IsPrint);
	$ErrorStop = 0 if(!defined $ErrorStop);

	my $Etot = $this->CalTotalEnergy($Compound, $IsPrint, $ErrorStop, $Debug);

	my $pAtomDB = $this->GetAtomDB();
#foreach my $key (sort keys %$pAtomDB) {
#	print "AtomDB: $key: $pAtomDB->{$key}{Emol}: $pAtomDB->{$key}{Ref}\n";
#}
	my ($Nc, $c) = ($Compound =~ /^(\d+)(\D.*)$/);
	if(!defined $c) {
		$Nc = 1;
		$c = $Compound;
	}
	my ($pe, $pn) = $this->CompoundToElements($c);

	my $AllDefined = 1;
	my $Eele = 0.0;
	for(my $i = 0 ; $i < @$pe ; $i++) {
		my $e = $pe->[$i];
		my $n = $pn->[$i];
#print "**$i: $e: $n\n";

		my $Emol = $pAtomDB->{$e}{Emol};
		if(!defined $Emol) {
print "Error: Can not find [$e(Atom)] in DB\n";
			if($ErrorStop) {
				exit;
			}
			$AllDefined = 0;
			next;
		}
		my $E = $n * $Emol;
#print "  $E = $n * $Emol\n" if($Debug);
		$Eele += $E * $Nc;
printf "  A[%d]: %s%d * %d (E=%10.4g * %d * %d) [$pAtomDB->{$e}{Type}][$pAtomDB->{$e}{Compound}][$pAtomDB->{$e}{Source}]\n",
	  $i, $e, $n, $Nc, $Emol, $n, $Nc, if($IsPrint);
	}
	my $Ecoh = ($AllDefined)? $Etot - $Eele : undef;
print "  Ecoh for [$Compound]: $Ecoh eV ($Etot - $Eele)\n" if($IsPrint);
	return $Ecoh;
	return undef;
}

sub CalTotalEnergyByArray
{
	my ($this, $pCompounds, $pComposition, $IsPrint) = @_;

	my $Etot = 0.0;
	for(my $i = 0 ; $i < @$pCompounds ; $i++) {
		my $c    = $pCompounds->[$i];
		my $Emol = $this->GetEmol($c);
if(!defined $Emol) {
	print "Error in CalTotalEnergyByArray: Can not get Emol for [$c] [mode=", $this->GetMode(), "]\n";
	exit;
}
print "  $pCompounds->[$i]: $pComposition->[$i] (E=$Emol)\n" if($IsPrint);
		$Etot += $Emol * $pComposition->[$i];
	}
#exit;
	return $Etot;
}

sub CalTotalEnergy
{
	my ($this, $Reaction, $IsPrint, $ErrorStop, $Debug) = @_;
$Debug = $this->GetDebugMode() if(!defined $Debug);
#my $Debug = 1;
#print "Etot for [$Reaction]\n" if($Debug);
	$IsPrint   = 0 if(!defined $IsPrint);
	$ErrorStop = 0 if(!defined $ErrorStop);

	my ($pRChemicals, $pRCompounds, $pRnCompounds) = $this->ReactionToCompounds($Reaction);

	my $AllDefined = 1;
	my $EReaction = 0.0;
	for(my $i = 0 ; $i < @$pRCompounds ; $i++) {
		my $c = $pRCompounds->[$i];
		my $n = $pRnCompounds->[$i];
#print "CalTotalEnergy: $i: $c * $n\n" if($Debug);

# $c=F2などの場合
		my ($Atom, $nAtom) = ($c =~ /^([A-Z][a-z]?)(\d*)$/);
		$c = $Atom if($Atom ne '');
		$nAtom = 1 if(!defined $nAtom or $nAtom <= 0);
		my $pDB  = $this->GetDBHash($c);
		my $Emol = $this->GetEmol($c);
		if(!defined $Emol) {
print "Error in CalTotalEnergy: Can not find [$c] in DBs [mode=", $this->GetMode(), "]\n";
			if($ErrorStop) {
				exit;
			}
			$this->ShowSpeculatedSpecies($c);
			$AllDefined = 0;
			next;
		}
		my $E = $n * $nAtom * $Emol;
#print "  $E = $n * $nAtom * $Emol\n" if($Debug);
		$EReaction += $E;
		my $nAtomStr = ($nAtom == 1)? '' : $nAtom;
printf "  C[%d]: %s%s * %d (E=%10.4g * %d * %d) [$pDB->{Type}][$pDB->{Compound}][$pDB->{Source}]\n",
		$i, $c, $nAtomStr, $n, $Emol, $n, $nAtom if($IsPrint);
	}
#print "exit\n";
	$EReaction = undef if(!$AllDefined);
print "  Etot for [$Reaction]: $EReaction eV\n" if($IsPrint);
	return $EReaction;
}

sub CalFormationEnergy
{
	my ($this, $Reaction, $IsPrint, $ErrorStop) = @_;

	my ($Reagents, $Products) = $this->SplitReaction($Reaction);
	my $EFormInitial = $this->CalTotalEnergy($Reagents, $IsPrint, $ErrorStop);
	my $EFormFinal   = $this->CalTotalEnergy($Products, $IsPrint, $ErrorStop);
	my $EForm = (defined $EFormFinal and defined $EFormInitial)? $EFormFinal - $EFormInitial : undef;
print "EForm for [$Reaction]: $EForm eV ($EFormFinal - $EFormInitial)\n" if($IsPrint);
	return $EForm;
}

1;