#!/usr/bin/perl

BEGIN {
#use lib 'd:/Programs/Perl/lib';
#use lib '/home/tkamiya/bin/lib';
my $BaseDir = $ENV{'TkPerlDir'};
print "\n\nBaseDir: $BaseDir\n";
@INC = ("$BaseDir/lib", "$BaseDir/VNL", "d:/Programs/Perl/lib", @INC);
}

use strict;
#use warnings;
use File::Path;
use File::Basename;
use File::Find;

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

use Deps;
use Utils;
use JFile;

use MyApplication;
use CSV;
use Sci::Algorism;
use Sci::GeneralFileFormat;

use Crystal::CIF;
use Crystal::Crystal;
use Crystal::SpaceGroup;
use Crystal::AtomType;
use Crystal::VASP;
use Crystal::XCrySDen;

#===============================================
# デバッグ関係変数
#===============================================
#$PrintLevelが大きいほど、情報が詳しくなる
my $PrintLevel = 0;

#===============================================
# 文字コード関係変数
#===============================================
# sjis, euc, jis, noconv
my $PrintCharCode      = Deps::PrintCharCode();
my $OSCharCode         = Deps::OSCharCode();
my $FileSystemCharCode = Deps::FileSystemCharCode();

my $LF        = Deps::LF();
my $DirSep    = Deps::DirSep();
my $RegDirSep = Deps::RegDirSep();

#===============================================
# Applicationオブジェクト作成
#===============================================
my $App = new MyApplication;
exit if($App->Initialize() < 0);

#$App->SetLF("<br>\n");
#$App->SetPrintCharCode("sjis");
#$App->SetDebug($Debug);
$App->SetDeleteHTMLFlag(1);

#===============================================
# スクリプト大域変数
#===============================================
my $InitialDirectory = Deps::GetWorkingDirectory();

#==========================================
# コマンドラインオプション読み込み
#==========================================
$App->AddArgument("--Action",
		"--Action=[MakeCIF|MakeLevelsCSV]", '');
$App->AddArgument("--LatticeParameterFactor", "--LatticeParameterFactor=[val] (Def: 1.2)", 1.2);
$App->AddArgument("--DebugMode", "--DebugMode: Set DebugMode", '');
exit 1 if($App->ReadArgs(0) != 1);
my $Args = $App->Args();
#my $form = new CGI;
#$Args->SetCGIForm($form);
#$Args->parseInput($WebCharCode);

#==========================================
# メイン関数スタート
#==========================================

my $Debug = $Args->GetGetArg("DebugMode");
$App->SetDebug($Debug);
my $Action = $Args->GetGetArg("Action");

my $ret = 0;
if($Action =~ /MakeCIF/i) {
	&MakeCIFFile();
}
elsif($Action =~ /MakeLevelsCSV/i) {
	&MakeLevelsCSV();
}
else {
	$App->print("Error: Invald Action: $Action\n");
}

#Utils::EndHTML();

exit $ret;

#===============================================
# スクリプト終了
#===============================================

#==========================================
# &Subroutines
#==========================================
sub MakeLevelsCSV
{
	$App->print("\n\nMake Levels CSV File from Gaussian03 out File:\n");
	my $toeV = Sci::HartreeToeV();

	my $OUTFile = $Args->GetGetArg(0);
	my $CSVFile = $Args->GetGetArg(1);
	unless($CSVFile) {
		$App->print("File names should be specified.\n");
		$App->print("   Usage: Gaussian.pl --Action=MakeLevelsCSV OUTFile CSVFile\n");
		return 0;
	}
	$App->print("  Read from [$OUTFile].\n");
	$App->print("  CSV File to be written: $CSVFile.\n");

	my $in = new JFile($OUTFile, "r");
	if(!$in) {
		$App->print("  Error!: Can not read [$OUTFile].\n");
		return -1;
	}

	my (@occ, @virt);
	for(my $iter = 0 ; ; $iter++) {
		my $line = $in->SkipTo("occ\\. eigenvalues");
		last if($in->eof());
#		$App->print("  [$iter] $line\n");
		@occ = ();
		@virt = ();
		my ($a) = ($line =~ /--\s*(.+)$/);
		my @a = Utils::Split("\\s+", $a);
#		print @a, "\n";
		if($line =~ /occ\./i) {
			@occ = (@occ, @a);
		}
		else {
			@virt = (@virt, @a);
		}
		while(1) {
			my $line = $in->ReadLine();
			last if($in->eof());
			last if($line !~ /eigenvalues/i);
#			$App->print("  [$iter] $line\n");
			($a) = ($line =~ /--\s*(.+)$/);
			@a = Utils::Split("\\s+", $a);
#			print @a, "\n";
			if($line =~ /occ\./i) {
				@occ = (@occ, @a);
			}
			else {
				@virt = (@virt, @a);
			}
		}
	}
	$in->Close();

	my @Ne;
	my @Energy;
	my @ELabel;
	my $count = 0;
	for(my $i = 0 ; $i < @occ ; $i++) {
		$Ne[$count] = 1.0;
		$Energy[$count] = $occ[$i] * $toeV;
		$App->print("  [$count] $Energy[$count] eV  $Ne[$count]\n");
		$count++;
		$ELabel[$count-1] = "O$count";
	}
	$ELabel[$count-1] = "HOMO";
	for(my $i = 0 ; $i < @virt ; $i++) {
		$Ne[$count] = 0.0;
		$Energy[$count] = $virt[$i] * $toeV;
		$App->print("  [$count] $Energy[$count] eV  $Ne[$count]\n");
		$count++;
		$ELabel[$count-1] = "V$count";
		$ELabel[$count-1] = "LUMO" if($i == 0);
	}

	my $csv = new CSV;
	my @LabelList = ("Label", "Energy(eV)", "Ne");
	my @DataArray = (\@ELabel, \@Energy, \@Ne);
#	$csv->SetFileType("CSV");
#	$csv->SetLabelArray(\@LabelList);
#	$csv->SetDataArray(\@DataArray);
	if($csv->Save($CSVFile, \@LabelList, \@DataArray)) {
		$App->print("  Save energy levels to [$CSVFile].\n");
	}
	else {
		$App->print("  Error!: Can not write to [$CSVFile].\n");
		return -2;
	}

	$App->print("\n\nMake Levels CSV File from Gaussian03 out File: Finished\n");
}

sub MakeCIFFile
{
	$App->print("\n\nMake CIF Files from Gaussian03 out File:\n");

	my $OUTFile = $Args->GetGetArg(0);
	my $CIFFile = $Args->GetGetArg(1);
	my $OptCIFFile = $Args->GetGetArg(2);
	my $XSFFile = $Args->GetGetArg(3);
	unless($CIFFile) {
		$App->print("File names should be specified.\n");
		$App->print("   Usage: Gaussian.pl --Action=MakeCIF CAR_Dir CIFFile\n");
		return 0;
	}
	$App->print("  Read from [$OUTFile].\n");
	$App->print("  CIF File for Initial Structure: $CIFFile.\n");
	$App->print("  CIF File for Optimized Structure: $OptCIFFile.\n");
	$App->print("  XSF File: $XSFFile.\n");
	my $LatticeParameterFactor = $Args->GetGetArg('LatticeParameterFactor');
	$App->print("  LatticeParameterFactor=$LatticeParameterFactor\n");

	my $in = new JFile($OUTFile, "r");
	if(!$in) {
		$App->print("  Error!: Can not read [$OUTFile].\n");
		return -1;
	}

	my $line = $in->SkipTo("# ");
#print "line: $line\n";	
	my $IsOptimized = 0;
	$IsOptimized = 1 if($line =~ /opt/i);
	$IsOptimized = 0 if(!$OptCIFFile);
	$App->print("  IsOptimized: $IsOptimized\n");
	$in->SkipTo("-----");

	$in->SkipTo("-----");
	my $title = $in->ReadLine();
	$App->print("  Title: $title\n");

	my $cycle = 0;
	my @Crystal;
	my ($xcenter,$ycenter,$zcenter);
	my ($a,$b,$c);
	my (@pfx,@pfy,@pfz);
	while(!$in->eof()) {
		$line = $in->SkipTo(" Center     Atomic     Atomic              Coordinates");
		last if(!$line);
		$cycle++;
$App->print("\n  *** Read structure at cycle $cycle ***\n");
#print "line: $line\n";
		$in->ReadLine();
		$in->ReadLine();

		my ($xmin,$xmax,$ymin,$ymax,$zmin,$zmax)
			= (1.0e10,-1.0e10,1.0e10,-1.0e10,1.0e10,-1.0e10);
		my ($xsum,$ysum,$zsum) = (0.0, 0.0, 0.0);
		my %Mol0;
		my %AtomCount;
		my $count = 0;
		while(!$in->eof()) {
			$line = $in->ReadLine();
			last if(!$line);

			my ($i, $AtNum, $AtType, $x, $y, $z) = Utils::Split("\\s+", $line);
			last if($i !~ /^\d+$/);
			my $AtName = AtomType::GetAtomName($AtNum);
			$AtomCount{$AtName}++;
			$Mol0{"Name$count"} = $AtName;
			my $label = $Mol0{"Label$count"} = "$AtName$AtomCount{$AtName}";
#$App->print("  $count: $AtName [$AtNum] $label ($x, $y, $z)\n");

			$Mol0{"x$count"} = $x;
			$Mol0{"y$count"} = $y;
			$Mol0{"z$count"} = $z;
			$xmin = $x if($x < $xmin);
			$xmax = $x if($x > $xmax);
			$ymin = $y if($y < $ymin);
			$ymax = $y if($y > $ymax);
			$zmin = $z if($z < $zmin);
			$zmax = $z if($z > $zmax);
			$xsum += $x;
			$ysum += $y;
			$zsum += $z;

			$count++;
		}
		$App->print("  nSites: $count\n");
		$App->print("  X range: $xmin - $xmax\n");
		$App->print("  Y range: $ymin - $ymax\n");
		$App->print("  Z range: $zmin - $zmax\n");
		if($cycle == 1) {
			$xcenter = ($xmin+$xmax)/2.0; #$xsum / $count;
			$ycenter = ($ymin+$ymax)/2.0; #$ysum / $count;
			$zcenter = ($zmin+$zmax)/2.0; #$zsum / $count;
			$App->print("  Center: ($xcenter, $ycenter, $zcenter)\n");
			$a = Utils::Round($LatticeParameterFactor * ($xmax - $xmin), 6);
			$b = Utils::Round($LatticeParameterFactor * ($ymax - $ymin), 6);
			$c = Utils::Round($LatticeParameterFactor * ($zmax - $zmin), 6);
			$App->print("  Lattice Parameters: a=$a  b=$b  c=$c\n");
		}
		$Mol0{"nSites"} = $count;

#Read Force
		$line = $in->SkipTo(" Center     Atomic                   Forces");
		if($line) {
			$in->ReadLine();
			$in->ReadLine();

			for(my $i = 0 ; $i < $count ; $i++) {
				$line = $in->ReadLine();
				last if(!$line);

				my ($i0, $AtNum, $fx, $fy, $fz) = Utils::Split("\\s+", $line);
				last if($i !~ /^\d+$/);
#$App->print("  $i: F=($fx, $fy, $fz)\n");

				$Mol0{"fx$i"} = $fx;
				$Mol0{"fy$i"} = $fy;
				$Mol0{"fz$i"} = $fz;
			}
		}

		my $ic = $cycle-1;
		$Crystal[$ic] = new Crystal();
		$Crystal[$ic]->SetSampleName($title);
		$Crystal[$ic]->SetCrystalName($title);
#print "latt: $a,$b,$c\n";
		$Crystal[$ic]->SetLatticeParameters($a, $b, $c, 90.0, 90.0, 90.0);

		for(my $i = 0 ; $i < $Mol0{'nSites'} ; $i++) {
			my $xc = ($Mol0{"x$i"}- $xcenter);
			my $yc = ($Mol0{"y$i"}- $ycenter);
			my $zc = ($Mol0{"z$i"}- $zcenter);
			my $fx = $Mol0{"fx$i"};
			my $fy = $Mol0{"fy$i"};
			my $fz = $Mol0{"fz$i"};
			if(!defined $fz) {
				$fx = $pfx[$i];
				$fy = $pfy[$i];
				$fz = $pfz[$i];
			}
			$pfx[$i] = $fx;
			$pfy[$i] = $fy;
			$pfz[$i] = $fz;
#print "*F[$i]=$fx,$fy,$fz\n";
			my $AtName = $Mol0{"Name$i"};
			my $Label = $Mol0{"Label$i"};
			my ($x,$y,$z) = $Crystal[$ic]->CartesianToFractional($xc,$yc,$zc);
#$App->print("  $i: $AtName $Label ($x, $y, $z)\n");

			$Crystal[$ic]->AddAtomType($AtName);
			$Crystal[$ic]->AddAtomSite($Label, $AtName, $x+0.5, $y+0.5, $z+0.5, 1.0,
					$fx, $fy, $fz);
		}
		$Crystal[$ic]->ExpandCoordinates();
	}
	$in->Close();

	$App->print("\n");
	my $CIF = new CIF;
	if($CIF->CreateCIFFileFromCCrystal($Crystal[0], $CIFFile, 0, 'unix')) {
		$App->print("  Optimized structure is written to [$CIFFile].\n");
	}
	else {
		$App->print("  Error!: Can not write to [$OptCIFFile].\n");
	}

	if($IsOptimized) {
		my $OptCIF = new CIF;
		if($OptCIF->CreateCIFFileFromCCrystal($Crystal[@Crystal-1], $OptCIFFile, 0, 'unix')) {
			$App->print("  Optimized structure is written to [$OptCIFFile].\n");
		}
		else {
			$App->print("  Error!: Can not write to [$OptCIFFile].\n");
		}

		$App->print("  XSF file is written to [$XSFFile].\n");
		my $XS = new XCrySDen;
		$XS->WriteXSFAnimationFileHeader($XSFFile, scalar @Crystal);
		for(my $i = 0 ; $i < @Crystal ; $i++) {
			$XS->WriteXSFFileFromCrystal($XSFFile, $i+1, $Crystal[$i]);
next;
			if($i == 0) {
				$XS->WriteXSFFileFromCrystal($XSFFile, $i+1, $Crystal[0], $Crystal[0]);
			}
			else {
				$XS->WriteXSFFileFromCrystal($XSFFile, $i+1, $Crystal[$i], $Crystal[$i-1]);
			}
		}
	}

	$App->print("\n");
	my ($drive, $dir, $filename, $ext, $lastdir, $filebody)
		= Deps::SplitFilePath($CIFFile);
	my $GJFFile    = Deps::MakePath("$drive$dir", "$filebody.gjf", 0);
	my $NewGJFFile = Deps::MakePath("$drive$dir", "$filebody-new-Cartesian.gjf", 0);
	$App->print("  Original input file: $GJFFile\n");
	$App->print("  New input file     : $NewGJFFile\n");

	my $in = new JFile($GJFFile, "r");
	if($in) {
		$App->print("  Read [$GJFFile].\n");
	}
	else {
		$App->print("  Error!!: Can not read [$GJFFile].\n");
		return -2;
	}
	my $out = new JFile($NewGJFFile, "w");
	if($out) {
		$App->print("  Writing to [$NewGJFFile].\n");
	}
	else {
		$App->print("  Error!!: Can not write to [$NewGJFFile].\n");
		$in->Close();
		return -3;
	}

	for(my $i = 0 ; ; $i++) {
		my $line = $in->ReadLine();
		last if($in->eof());
		$line =~ s/^(#.*)geom=connectivity(.*)$/$1$2/i;
		$out->print("$line");
		last if($line !~ /^[%#]/);
	}
	for(my $i = 0 ; $i < 3 ; $i++) {
		my $line = $in->ReadLine();
		last if($in->eof());
		$out->print("$line");
	}
	my $Crystal = $Crystal[@Crystal - 1];
	my @AtomList = $Crystal->GetCExpandedAtomSiteList();
	my $nAtom = @AtomList;
	for(my $i = 0 ; $i < $nAtom ; $i++) {
		my $atom  = $AtomList[$i];
		my $atname = $atom->AtomNameOnly();
		my ($x,$y,$z) = $atom->Position(0);
		my ($xc,$yc,$zc) = $Crystal->FractionalToCartesian($x,$y,$z);
		$out->printf("%2s     %10.6f  %10.6f  %10.6f\n", $atname, $xc, $yc, $zc);
	}

	$in->Close();
	$out->Close();

	$App->print("\nMake CIF Files from Gaussian03 out File: Finished\n");

	return 1;
}