#!/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 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]", '');
$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();
}
else {
	$App->print("Error: Invald Action: $Action\n");
}

#Utils::EndHTML();

exit $ret;

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

#==========================================
# &Subroutines
#==========================================
sub MakeCIFFile
{
	$App->print("\n\nMake CIF Files from MOPAC 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("MOPAC:  VERSION");
	$line = $in->SkipTo("\\*\\*\\*\\*\\*");
	$line = $in->ReadLine();
	my $IsOptimized = 1;
	$IsOptimized = 0 unless($line =~ /AM1/i);
	$App->print("  IsOptimized: $IsOptimized\n");

	$in->ReadLine();
	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("CARTESIAN COORDINATES");
		last if(!$line);
		$cycle++;
$App->print("\n  *** Read structure at cycle $cycle ***\n");
#print "line: $line\n";
		$in->ReadLine();
		$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, $AtName, $x, $y, $z) = Utils::Split("\\s+", $line);
			last if($i !~ /^\d+$/);
			$AtomCount{$AtName}++;
			$Mol0{"Name$count"} = $AtName;
			my $label = $Mol0{"Label$count"} = "$AtName$AtomCount{$AtName}";
$App->print("  $count: $AtName $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;

		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 $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);
		}
		$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("\nMake CIF Files from MOPAC out File: Finished\n");

	return 1;
}