package GAMESS;
use Exporter;
@ISA = qw(Exporter);

#公開したいサブルーチン
@EXPORT = qw();

use strict;
use File::Basename;

use Utils;
use GraphData;
use Sci::Science;
use Sci::MolecularOrbital;
use MyTk::GraphFrameArray;

use Crystal::Crystal;
use Crystal::AtomType;
use Crystal::AtomSite;
use Crystal::XCrySDen;
use Crystal::CIF;

#===============================================
# スクリプト大域変数
#===============================================
my $LF = "<br>\n";
my $DirectorySeparator = "\\";

#===============================================
# パス（読み込みDB）
# Web関係変数
# CGI の仮想パス
# プログラム名
#===============================================
my $Program         = basename($0);
my $ProgramDir      = "d:\\Programs";

my $WebRootDir      = "d:\\MyWebs";
my $KListDBDir      = Deps::MakePath($WebRootDir, "Research", 1);

#============================================================
#　変数等取得、再定義関数
#============================================================
sub ClearAll      { my $this=shift; undef $this->{'FileType'}; undef $this->{'DataArray'}; }
sub FileType      { return shift->{'FileType'}; }
sub FileName      { return shift->{'FileName'}; }
sub SetFileName   { my ($this,$f)=@_; return $this->{'FileName'} = $f; }
#sub DataArray     { return shift->{'DataArray'}; }
#sub SetDataArray  { my ($this, $da) = @_; return $this->{'DataArray'} = $da; }
sub SetSampleName { my ($this,$n)=@_; return $this->{'SampleName'} = $n; }
sub SampleName    { return shift->{'SampleName'}; }

#============================================================
#　継承クラスで定義しなおす関数
#============================================================
sub CheckFileType
{
	my ($path) = @_;
#print "Gaussian::CheckFileType\n";
	my ($drive, $dir, $filename, $ext, $lastdir, $filebody) = Deps::SplitFilePath($path);
#print "f: $filename\n";
	if($filename =~ /\.out$/i) {
		my $in = new JFile($path, "r");
		return undef unless($in);
		for(my $i = 0 ; $i < 10 ; $i++) {
			my $line1 = $in->ReadLine();
#print "l: $line1";
			if($line1 =~ /\*(.*GAMESS.*?),/i) {
				my $type = $1;
				$in->Close();
				return "$type Output file";
			}
		}
		$in->Close();
	}
	return undef;
}

sub Read
{
	my ($this, $filename) = @_;
#print "this: $this\n";
#print "Gaussian::Read\n";
	$this->ClearAll();
	my $FileType = $this->{FileType} = GAMESS::CheckFileType($filename);
	$this->SetFileName($filename);
#print "aFileType: $FileType\n";
	if($FileType =~ /GAMESS.*Output/i) {
		my $MO = $this->ReadOutputFile($filename);
		return $this->MoToGraphData($MO);
	}
}

#============================================================
#　コンストラクタ、デストラクタ
#============================================================
sub new
{
	my ($module) = @_;
	my $this = {};
	bless $this;
	return $this;
}

sub DESTROY { my $this = shift; }

#============================================================
#　ファイル読み込み関数
#============================================================
sub ReadOutputFile
{
	my ($this, $path) = @_;

	my $MO = new MolecularOrbital;

	my $in = new JFile($path, "rb");
	return undef unless($in);

	$MO->SetSourcePath($path);

	my $line = $in->SkipTo(" INPUT CARD\\> \\\$DATA");
	my $SampleName = $in->ReadLine();
	Utils::DelSpace($SampleName);
	if($SampleName =~ /\>Winmostar/) {
		my ($drive, $directory, $filename, $ext, $lastdir, $filebody) = Deps::SplitFilePath($path);
		$SampleName = $filebody;
		Utils::DelSpace($SampleName);
	}
	$MO->SetTitle($SampleName);
print("Sample: $SampleName\n");

	$in->rewind();
	my $IsUHF = 0;
	$IsUHF = 1 if($in->SkipTo("SCFTYP\\s*=\\s*UHF "));
	$MO->SetUHF($IsUHF);
print "IsUHF: $IsUHF\n";

	my ($nBasisFunctions, $nAlphaElectrons, $nBetaElectrons);
	$in->rewind();
 	$line = $in->SkipTo(" TOTAL NUMBER OF BASIS FUNCTIONS *= *(\\d+)");
	($nBasisFunctions) = ($line =~ / TOTAL NUMBER OF BASIS FUNCTIONS *= *(\d+)/);
	$line = $in->SkipTo(" NUMBER OF OCCUPIED ORBITALS \\(ALPHA\\) *= *(\\d+)");
	($nAlphaElectrons) = ($line =~ / NUMBER OF OCCUPIED ORBITALS \(ALPHA\) *= *(\d+)/);
	$MO->SetnAlphaElectrons($nAlphaElectrons);
	$line = $in->SkipTo(" NUMBER OF OCCUPIED ORBITALS \\(BETA \\) *= *(\\d+)");
	($nBetaElectrons)  = ($line =~ / NUMBER OF OCCUPIED ORBITALS \(BETA \) *= *(\d+)/);
	$MO->SetnBetaElectrons($nBetaElectrons);
print "# of basis functions: $nBasisFunctions\n";
print "# of alpha electrons: $nAlphaElectrons\n";
print "# of beta  electrons: $nBetaElectrons\n";

	$in->rewind();
	my $nAO = 0;
	for(my $ispin = 0 ; $ispin <= $IsUHF ; $ispin++) {
		if($IsUHF) {
			if($ispin == 0) {
				$line = $in->SkipTo("---\\s*ALPHA SET\\s+--");
			}
			else {
				$line = $in->SkipTo("---\\s*BETA SET\\s+--");
			}
			$line = $in->SkipTo(" EIGENVECTORS");
		}
		else {
			$line = $in->SkipTo(" EIGENVECTORS");
		}
		while(1) {
#read:                      1          2          3          4          5
			while(1) {
				$line = $in->ReadLine();
				last if(!defined $line);
				last if($line =~ /^\s*(\d+\s+)+$/);
			}
			my @RN = Utils::Split("\\s+", $line);
#read:                  -20.2436    -1.2675    -0.6142    -0.4549    -0.3917
			$line = $in->ReadLine();
			last if(!defined $line);
			my @energies = Utils::Split("\\s+", $line);
			$line = $in->ReadLine();
			my @AO = Utils::Split("\\s+", $line);

			my $nAOBase = $nAO;
			for(my $j = 0 ; $j < @energies ; $j++) {
				my $ne;
				if($ispin == 0) {
					$ne = ($RN[$j] <= $nAlphaElectrons)? 1.0 : 0.0;
				}
				else {
					$ne = ($RN[$j] <= $nBetaElectrons)? 1.0 : 0.0;
				}
print "Add: $nAO, $RN[$j], $AO[$j], $energies[$j], $ne, $ispin\n";
				$MO->Add($nAO, $RN[$j], $AO[$j], $energies[$j], $ne, $ispin);
				$nAO++;
			}

			my $nOrb = 0;
			while(1) {
				$line = $in->ReadLine();
				last if(!defined $line);
				my ($idx, $atom, $iatom, $orb, @coeff) = Utils::Split("\\s+", $line);
				last if(@coeff == 0);
#print "$nOrb: $orb: $atom$iatom: ", join(', ', @coeff), "\n";
				$MO->SetBase($nOrb, "$atom$iatom $orb");
				for(my $j = 0 ; $j < @coeff ; $j++) {
					$MO->SetWFCoeff($nAOBase+$j, $nOrb, $coeff[$j]);
				}
				$nOrb++;
			}

			last if($nAO >= $nBasisFunctions*($ispin+1));
		}
print "nAtomicOrbitals(spin=$ispin): $nAO\n";
	}
	$in->Close();
	$MO->NormalizeWFCoeff();

	return $MO;
}

sub MoToGraphData
{
	my ($this, $MO, $target) = @_;
	$target = '' unless(defined $target);

	my $nAlphaOccupied  = 0;
	my @AlphaOccupiedEnergy;
	my @XAlphaOccupied;
	my @SymmetryAlphaOccupied;
	my @strOrbitalsAlphaOccupied;
	my $nAlphaVirtual   = 0;
	my @AlphaVirtualEnergy;
	my @XAlphaVirtual;
	my @strOrbitalsAlphaVirtual;
	my @SymmetryAlphaVirtual;
	my $nBetaOccupied   = 0;
	my @BetaOccupiedEnergy;
	my @XBetaOccupied;
	my @strOrbitalsBetaOccupied;
	my @SymmetryBetaOccupied;
	my $nBetaVirtual    = 0;
	my @BetaVirtualEnergy;
	my @XBetaVirtual;
	my @SymmetryBetaVirtual;
	my @strOrbitalsBetaVirtual;

	my $nAtomicOrbitals = $MO->nRedMOs();
print "nAtomicOrbitals(alpha+beta): $nAtomicOrbitals\n";

	for(my $i = 0 ; $i < $nAtomicOrbitals ; $i++) {
		my $RootNumber    = $MO->RootNumber($i);
		my $AtomicOrbital = $MO->Symmetry($i);
		my $Spin          = $MO->Spin($i);
		my $Energy        = $MO->Energy($i);
		my $Ne            = $MO->Ne($i);
print "$i: spin=$Spin ne=$Ne [RN:$RootNumber]: $AtomicOrbital: $Energy\n";
		my $orb = '';
		for(my $j = 0 ; $j < $MO->nBase() ; $j++) {
			my $c = $MO->WFCoeff($i, $j);
			next if(abs($c) < 0.01);
			$orb .= "\n" . $MO->Base($j) . ": " . sprintf("%3.2f", $c);
		}
		if($Spin == 0 and $Ne != 0) {
			push(@SymmetryAlphaOccupied, $AtomicOrbital);
			$XAlphaOccupied[$nAlphaOccupied]           = $Ne;
			$AlphaOccupiedEnergy[$nAlphaOccupied]      = $Energy;
			$strOrbitalsAlphaOccupied[$nAlphaOccupied] = " alpha occupied$orb";
			$nAlphaOccupied++;
		}
		elsif($Spin == 0 and $Ne == 0) {
			push(@SymmetryAlphaVirtual, $AtomicOrbital);
			$XAlphaVirtual[$nAlphaVirtual]           = $Ne;
			$AlphaVirtualEnergy[$nAlphaVirtual]      = $Energy;
			$strOrbitalsAlphaVirtual[$nAlphaVirtual] = " alpha virtual$orb";
			$nAlphaVirtual++;
		}
		elsif($Spin == 1 and $Ne != 0) {
			push(@SymmetryBetaOccupied, $AtomicOrbital);
			$XBetaOccupied[$nBetaOccupied]           = $Ne;
			$BetaOccupiedEnergy[$nBetaOccupied]      = $Energy;
			$strOrbitalsBetaOccupied[$nBetaOccupied] = " beta occupied$orb";
			$nBetaOccupied++;
		}
		elsif($Spin == 1 and $Ne == 0) {
			push(@SymmetryBetaVirtual, $AtomicOrbital);
			$XBetaVirtual[$nBetaVirtual]           = $Ne;
			$BetaVirtualEnergy[$nBetaVirtual]      = $Energy;
			$strOrbitalsBetaVirtual[$nBetaVirtual] = " beta virtual$orb";
			$nBetaVirtual++;
		}
	}

	$this->SetDataArray(new GraphDataArray);
	my $Data = new GraphData;
	$Data->SetTitle($MO->Title());
	$Data->{x0}          = \@XAlphaOccupied;
	$Data->{x1}          = \@XAlphaVirtual;
	$Data->{x2}          = \@XBetaOccupied;
	$Data->{x3}          = \@XBetaVirtual;
	$Data->{x0_Name}     = "alph occ";
	$Data->{x1_Name}     = "alph virt";
	$Data->{x2_Name}     = "beta occ";
	$Data->{x3_Name}     = "beta virt";
	$Data->{y0}          = \@AlphaOccupiedEnergy;
	$Data->{y1}          = \@AlphaVirtualEnergy;
	$Data->{y2}          = \@BetaOccupiedEnergy;
	$Data->{y3}          = \@BetaVirtualEnergy;
	$Data->{y0_Name}     = "Energy (Alpha occupied)";
	$Data->{y1_Name}     = "Energy (Alpha virtual)";
	$Data->{y2_Name}     = "Energy (Beta occupid)";
	$Data->{y3_Name}     = "Energy (Beta virtual)";
	$Data->{Symmetry0}   = \@SymmetryAlphaOccupied;
	$Data->{Symmetry1}   = \@SymmetryAlphaVirtual;
	$Data->{Symmetry2}   = \@SymmetryBetaOccupied;
	$Data->{Symmetry3}   = \@SymmetryBetaVirtual;
	$Data->{strOrbital0} = \@strOrbitalsAlphaOccupied;
	$Data->{strOrbital1} = \@strOrbitalsAlphaVirtual;
	$Data->{strOrbital2} = \@strOrbitalsBetaOccupied;
	$Data->{strOrbital3} = \@strOrbitalsBetaVirtual;

	$Data->CalMinMax();
	$this->{DataArray}->AddGraphData($Data);
	$this->{DataArray}->{TargetData} = $target;

	return $MO->SourcePath();
}

#============================================================
#　一般関数
#============================================================

1;