package MOPAC;
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::Molecule;
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 =~ /\.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 =~ /\s(MOPAC\S+\s+(VERSION\s+[\d\.]+)?)/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} = MOPAC::CheckFileType($filename);
	$this->SetFileName($filename);
#print "aFileType: $FileType\n";
	if($FileType =~ /^MOPAC.*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 SaveDatFileFromMolecule
{
	my ($this, $path, $mol, $IsPrint) = @_;

	my $out = new JFile;
	$out->Open($path, "w") or return undef;

	$out->print("AM1 EF PRECISE GNORM=0.05 NOINTER GRAPHF\n");
	$out->print("\n");
print($mol->Title(), "\n") if($IsPrint);
	$out->print($mol->Title(), "\n");
	my @AtomSite = $mol->GetCAtomSiteList();
	for(my $i = 0 ; $i < @AtomSite ; $i++) {
		my $site = $AtomSite[$i];
		my $atomname = $site->AtomNameOnly();
		my ($x, $y, $z) = $site->Position();
		my ($R, $Angle1, $Angle2)        = (0, 0, 0);
		my ($idR, $idAngle1, $idAngle2) = (1, 1, 1);
		my ($iAtom2, $iAtom3, $iAtom4)  = (0, 0, 0);
		my ($site2, $site3, $site4);
		my ($x2, $y2, $z2);
		my ($x3, $y3, $z3);
		my ($x4, $y4, $z4);
		if($i == 0) {
		}
		if($i >= 1) {
			$iAtom2 = 0;
			$site2 = $AtomSite[$iAtom2];
			($x2, $y2, $z2) = $site2->Position();
			my ($dx, $dy, $dz) = ($x2-$x, $y2-$y, $z2-$z);
			my $R2 = Vector::InnerProduct($dx, $dy, $dz, $dx, $dy, $dz);
			$R   = sqrt($R2);
			$idR = 1;
		}
		if($i >= 2) {
			$iAtom3 = 1;
			$site3 = $AtomSite[$iAtom3];
			($x3, $y3, $z3) = $site3->Position();
			$Angle1 = Vector::Angle([$x-$x2,  $y-$y2,  $z-$z2], 
									[$x3-$x2, $y3-$y2, $z3-$z2],
									1);
			$idAngle1 = 1;
		}
		if($i >= 3) {
			$iAtom4 = 2;
			$site4 = $AtomSite[$iAtom4];
			($x4, $y4, $z4) = $site4->Position();
			my @r123 = Vector::OuterProduct($x-$x2,  $y-$y2,  $z-$z2,  $x3-$x2, $y3-$y2, $z3-$z2);
			my @r234 = Vector::OuterProduct($x2-$x3, $y2-$y3, $z2-$z3, $x4-$x3, $y4-$y3, $z4-$z3);
			$Angle2 = Vector::Angle(\@r123, \@r234, 1);
			$idAngle2 = 1;
			$Angle2 = -$Angle2 if($z < 0.0);
		}
		$out->printf("%2s  %9.5f    %d %10.5f  %d %10.5f %d %6d %6d %6d\n",
			$atomname, $R, $idR, $Angle1, $idAngle1, $Angle2, $idAngle2,
			$iAtom2+1, $iAtom3+1, $iAtom4+1);
printf("%2s  %9.5f    %d %10.5f  %d %10.5f %d %6d %6d %6d\n",
		$atomname, $R, $idR, $Angle1, $idAngle1, $Angle2, $idAngle2,
		$iAtom2+1, $iAtom3+1, $iAtom4+1) if($IsPrint);
	}
	$out->Close();
	return 1;
}

sub ReadDatFile
{
	my ($this, $path, $TranslateToCenter, $IsPrint) = @_;

	my $mol = new Molecule;

	my $in  = new JFile;
	if(!$in->Open($path,  "r")) {
		return undef;
	}

	$mol->SetCalculationCondition($in->ReadLine(1));
	$in->ReadLine(1);
	$mol->SetTitle($in->ReadLine(1));

	while(1) {
		my $line = $in->ReadLine(1);
		last if($line eq '');
		my $i = $mol->nZMatrix();
			my ($atom, $d, $idd, $angle1, $idangle1, $angle2, $idangle2, 
				$atom2, $atom3, $atom4) = Utils::Split("\\s+", $line);
			$mol->AddZMatrix($atom, $d, $angle1, $angle2, 
						 $atom2, $atom3, $atom4,
						 $idd, $idangle1, $idangle2);
	}
	$in->Close();

	$mol->ZMatrixToXYZ($TranslateToCenter, $IsPrint);

	return $mol;
}

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(" CALCULATION RESULTS");
	$line = $in->SkipTo("\\*\\*\\*\\*\\*");
	$line = $in->SkipTo("\\*\\*\\*\\*\\*");
	$line = $in->ReadLine();
	$line = $in->ReadLine();
	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();
	my $IsUHF = 0;
	$IsUHF = 1 if($in->SkipTo("\\* +UHF +\\-"));
	$MO->SetUHF($IsUHF);
print "IsUHF: $IsUHF\n";

	my ($nAlphaElectrons, $nBetaElectrons);
	$in->rewind();
	if($IsUHF) {
		$in->SkipTo(" FINAL HEAT OF FORMATION *=");
		$line = $in->SkipTo(" NO. OF ALPHA ELECTRONS *= *(\\d+)");
		($nAlphaElectrons) = ($line =~ / NO. OF ALPHA ELECTRONS *= *(\d+)/);
		$MO->SetnAlphaElectrons($nAlphaElectrons);
		$line = $in->SkipTo(" NO. OF BETA  ELECTRONS *= *(\\d+)");
		($nBetaElectrons) = ($line =~ / NO. OF BETA  ELECTRONS *= *(\d+)/);
		$MO->SetnBetaElectrons($nBetaElectrons);
print "# of alpha electrons: $nAlphaElectrons\n";
print "# of beta electrons: $nBetaElectrons\n";
	}
	else {
		$line = $in->SkipTo(" NO. OF DOUBLY OCCUPIED LEVELS *= *(\\d+)");
		($nAlphaElectrons) = ($line =~ / NO. OF DOUBLY OCCUPIED LEVELS *= *(\d+)/);
		$MO->SetnAlphaElectrons($nAlphaElectrons);
print "# of doubly occupied levels: $nAlphaElectrons\n";
	}

	$in->rewind();
	if($IsUHF) {
		$line = $in->SkipTo(" ALPHA EIGENVECTORS ");
	}
	else {
		$line = $in->SkipTo(" EIGENVECTORS ");
	}
	$line = $in->SkipTo(" Root\\s+No");
	my @RootNumbers = Utils::Split("\\s+", $line);
	my $nReadLines = 1;
	if($RootNumbers[@RootNumbers-1] < $nAlphaElectrons) {
		$nReadLines = 2;
	}

	$in->rewind();
	my $nAO = 0;
	for(my $ispin = 0 ; $ispin <= $IsUHF ; $ispin++) {
		for(my $i = 0 ; $i < $nReadLines ; $i++) {
			if($ispin == 0) {
				$line = $in->SkipTo(" Root\\s+No");
			}
			else {
				if($i == 0) {
					$in->SkipTo(" BETA EIGENVECTORS ");
				}
				$line = $in->SkipTo(" Root\\s+No");
			}
			$line =~ s/Root No\.//;
			my @RN = Utils::Split("\\s+", $line);
			$line = $in->ReadLine();
			$line = $in->ReadLine();
			$line =~ s/(\d)\s+(\w)/$1$2/g;
			$line = lc $line;
			my @AO = Utils::Split("\\s+", $line);
			$in->ReadLine();
			$line = $in->ReadLine();
			my @energies = 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;
				}
				$MO->Add($nAO, $RN[$j], $AO[$j], $energies[$j], $ne, $ispin);
				$nAO++;
			}

			$line = $in->ReadLine();
			my $nOrb = 0;
			while(1) {
				$line = $in->ReadLine();
				last if(!defined $line);
				my ($orb, $atom, $iatom, @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++;
			}
		}
	}
	$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;