package Gaussian;
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;

#===============================================
# スクリプト大域変数
#===============================================

#============================================================
#　変数等取得、再定義関数
#============================================================
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 =~ /\.(log|out)$/i) {
		my $in = new JFile($path, "r");
		return undef unless($in);
		for(my $i = 0 ; $i < 10 ; $i++) {
			my $line1 = $in->ReadLine();
			if($line1 =~ /Copyright\s.*Gaussian,\s*Inc\./) {
				$in->Close();
				return "Gaussian03 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'} = Gaussian::CheckFileType($filename);
	$this->SetFileName($filename);
#print "aFileType: $FileType\n";
	if($FileType eq "Gaussian03 Output file") {
		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 $HrToeV = Sci::HartreeToeV();

	my $MO = new MolecularOrbital;

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

	$MO->SetSourcePath($path);

	my $line = $in->SkipTo("Gaussian 03:");
	$line = $in->SkipTo(" ---------------------------------");
	$line = $in->ReadLine();
	my $IsUHF = 0;
	$IsUHF = 1 if($line =~ /uhf/i);
	$MO->SetUHF($IsUHF);
print "IsUHF: $IsUHF\n";

	$line = $in->SkipTo(" ---------------------------------");
	$line = $in->SkipTo(" ---------");
	my $SampleName = $in->ReadLine();
	Utils::DelSpace($SampleName);
	if($SampleName eq '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();
	$line = $in->SkipTo("NBasis= +\\d+");
	my ($nBasisFunctions) = ($line =~ /NBasis= +(\d+)/);
print("nBasisFunctions: $nBasisFunctions\n");
 
	my ($nAlphaElectrons, $nBetaElectrons);
	$in->rewind();
	my $nAO = 0;
	my ($previatom, $prevatom) = (0, '');
	for(my $ispin = 0 ; $ispin <= $IsUHF ; $ispin++) {
		if($IsUHF) {
			if($ispin == 0) {
				$line = $in->SkipTo(" Alpha Molecular Orbital Coefficients");
			}
			else {
				$line = $in->SkipTo(" Beta Molecular Orbital Coefficients");
			}
		}
		else {
			$line = $in->SkipTo(" Molecular Orbital Coefficients");
		}
		while(1) {
			$line = $in->ReadLine();
			last if(!defined $line);
			my @RN = Utils::Split("\\s+", $line);
			$line = $in->ReadLine();
			my @Orbs = Utils::Split("\\s+", $line);
			$line = $in->ReadLine();
			my ($eigstr, $minusstr, @energies) = Utils::Split("\\s+", $line);
			for(my $i = 0 ; $i < @energies ; $i++) {
				$energies[$i] *= $HrToeV;
#print "  i=$i: e=$energies[$i]\n";
			}

			my $nAOBase = $nAO;
			my (@AO, @IsOccupied);
			for(my $i = 0 ; $i < @Orbs ; $i++) {
				($AO[$i], $IsOccupied[$i]) = ($Orbs[$i] =~ /\((.*)\)--(.?)/);
				$AO[$i] = lc $AO[$i];
				$IsOccupied[$i] = $IsOccupied[$i-1] if($IsOccupied[$i] eq '' and $i > 0);
				$IsOccupied[$i] = ($IsOccupied[$i] eq 'O')? 1 : 0;
#print "ispin=$ispin: iAO=$nAO: i=$i: $AO[$i], $IsOccupied[$i]\n";
			}

			for(my $j = 0 ; $j < @Orbs ; $j++) {
				my $ne;
				if($IsOccupied[$j]) {
					if($ispin == 0) {
						$nAlphaElectrons++;
					}
					else {
						$nBetaElectrons++;
					}
					$ne = 1.0;
				}
				else {
					$ne = 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) {
				my $pos = $in->tell();
				$line = $in->ReadLine();
				last if(!defined $line);
#print "l[nO=$nOrb]: $line";
				if($line =~ /^\s+(\d+\s+)+\s*$/) {
					$in->seek($pos, 0);
					last;
				}

				my ($idx, $iatom, $atom, $orb, @coeff);
				if($line =~ /^\s+\d+\s+\d+\s+\w+\s/) {
					($idx, $iatom, $atom, $orb, @coeff) = Utils::Split("\\s+", $line);
					 $previatom = $iatom;
					 $prevatom  = $atom;
				}
				else {
					($idx,$orb, @coeff) = Utils::Split("\\s+", $line);
					 $iatom = $previatom;
					 $atom  = $prevatom;
				}
#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));
#last if($nAO >= 5);
		}
print "nAtomicOrbitals(spin=$ispin): $nAO\n";
	}
	$in->Close();

	$MO->SetnAlphaElectrons($nAlphaElectrons);
	$MO->SetnBetaElectrons($nBetaElectrons);
print "# of alpha electrons: $nAlphaElectrons\n";
print "# of beta  electrons: $nBetaElectrons\n";
	$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;