#!/usr/bin/perl

BEGIN {
use lib 'c:/Programs/Perl/lib';
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/GULP", @INC);
#use lib "$BaseDir/lib";
#use lib "$BaseDir/GULP";
}

use strict;
#use warnings;

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

use Deps;
use Utils;

use ProgVars;
use MyApplication;
use Sci::GeneralFileFormat;

use Crystal::CIF;
use Crystal::Crystal;
use Crystal::SpaceGroup;
use Crystal::AtomType;
use Crystal::MD;
use Crystal::GULP;

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

#===============================================
# 大域変数
#===============================================
my $HOME = $ENV{'HOME'};
my $ProgramDir        = ProgVars::ProgramDir();
my $LAMMPSPerlDir     = ProgVars::LAMMPSPerlDir();
my $TemplateDir       = Deps::MakePath($LAMMPSPerlDir, 'Template', 0);
my $LAMMPSDatabaseDir = Deps::MakePath($LAMMPSPerlDir, 'Potentials', 0);
my $GULPDir           = ProgVars::GULPDir(); #Deps::MakePath($HOME, "gulp", 0);
my $GULPLibDir        = Deps::MakePath($GULPDir, "Libraries", 0);
my @PotentialLibs     = ForceField::PotentialLibs();

#===============================================
# 文字コード関係変数
#===============================================
# 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=[ShowDatabases|FindDatabases"
							    ."|MakeInput|UpdateInputFile|ModifyInput|MakeOptCIF|MakeHistoryCSV|PostProcess"
							    ."|MergeCSV|CalMSD"
							    ."]",       '');
$App->AddArgument("--Function", "--Function=[Optimize|Fit|MD-NVT|MD-NPT|Phonon] (Def:Optimize)", 'Optimize');
$App->AddArgument("--PotentialType", "--PotentialType=[buck|morse] (Def:)", '');
$App->AddArgument("--CopyScript",    "--CopyScript=[0|1] (Def:0)",     0);
$App->AddArgument("--UseShellModelForCation", "--UseShellModelForCation=[0|1] (Def:0)",   0);
$App->AddArgument("--UseShellModelForAnion",  "--UseShellModelForAnion=[0|1] (Def:0)",    0);
$App->AddArgument("--LibraryFile",  "--LibraryFile=[Kamiya.lib|Kamiya-NoShell.lib|bush.lib|catlow.lib|"
				   ."lewis.lib|suttonchen.lib|glass.lib] (Def:0)",   "Kamiya-NoShell.lib");
$App->AddArgument("--UseSpeciesCharges",  "--UseSpeciesCharges=[yes|no] (Def:no)", 'no');
$App->AddArgument("--TargetFiles",        "--TargetFiles=[files] (Def:)", '');

$App->AddArgument("--temperature",     "--Temperature=start,end (Def:300.0,3000.0)", "300.0,3000.0");
$App->AddArgument("--pressure",        "--Pressure=val          (Def:0.1)",        0.1);
$App->AddArgument("--timestep",        "--timestep=val          (Def:0.001)",      0.001);
$App->AddArgument("--tscale",          "--tscale=val            (Def:0.10)",       0.10);
$App->AddArgument("--equilibration",   "--equilibration=val     (Def:0.10)",       0.10);
$App->AddArgument("--production",      "--production=val        (Def:10.0)",      10.001);
$App->AddArgument("--sample",          "--sample=val            (Def:0.1)",        0.001);
$App->AddArgument("--write",           "--write=val             (Def:0.1)",        0.005);
$App->AddArgument("--ResetTime",       "--ResetTime=[0|1]       (Def:0)",              0);

$App->AddArgument("--nOutputInterval", "--nOutputInterval       (Def:100)",          100);
$App->AddArgument("--DebugMode", "--DebugMode: Set DebugMode", '');
exit 1 if($App->ReadArgs(0) != 1);
my $Args = $App->Args();

$App->SetAppName("GULP.pl");
$App->SetUsage("GULP.pl [options] Files");

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

#Utils::InitHTML("Research", $WebCharSet, "_self");

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

if($Action =~ /MakeInput/i) {
	&MakeInputFile();
}
elsif($Action =~ /UpdateInputFile/i) {
	&UpdateInputFile();
}
elsif($Action =~ /ModifyInput/i) {
	&ModifyInput();
}
elsif($Action =~ /ShowDatabases/i) {
	&ShowDatabases();
}
elsif($Action =~ /FindDatabases/i) {
	&FindDatabases();
}
elsif($Action =~ /MakeOptCIF/i) {
	&MakeOptCIFFile();
}
elsif($Action =~ /PostProcess/i or $Action =~ /MakeHistoryCSV/i) {
	&PostProcess();
}
elsif($Action =~ /MergeCSV/i or $Action =~ /MergeCSV/i) {
	&MergeCSV();
}
elsif($Action =~ /PostProcess/i or $Action =~ /CalMSD/i) {
	&CalMSD();
}
else {
	$App->print("Error: Invald Action: $Action\n");
	$App->Usage();
}

#Utils::EndHTML();

exit;

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

#==========================================
# &Subroutines
#==========================================

# this is obsolete
sub MakeOptCIFFile
{
	$App->print("\n\nMake CIF File from Optimization output:\n");

	my $OutputFile = $Args->GetGetArg(0);
	my $CIFFile = $Args->GetGetArg(1);
	unless($CIFFile) {
		$App->print("File name should be specified.\n");
		$App->print("  usage: GULP.pl --Action=MakeOptCIF OutputFile CIFFile.\n");
		return 0;
	}
	
	my $InitialFile = $CIFFile;
	$InitialFile =~ s/{to_be_replaced}/-Initial/g;
	my $OptimizedFile = $CIFFile;
	$OptimizedFile =~ s/{to_be_replaced}/-Optimized/g;

	$App->print("  Read $OutputFile.\n");
	$App->print("  Save Initial structure   to $InitialFile.\n");
	$App->print("  Save Optimized structure to $OptimizedFile.\n");


	my $GULP = new GULP;
	my $InitialCrystalCIF
		= $GULP->ExtractCrystalStructuresFromOptimiseOutput($OutputFile, 1);
	my $ret;
	if(($ret = $InitialCrystalCIF) <= 0) {
		$App->print("  Error in GULP.pm::"
			."ExtractCrystalStructuresFromOptimiseOutput1 [Ret: $ret]\n");
		return -1;
	}
	my $OptimizedCrystalCIF
		= $GULP->ExtractCrystalStructuresFromOptimiseOutput($OutputFile, 2);
	if(($ret = $OptimizedCrystalCIF) <= 0) {
		$App->print("  Error in GULP.pm::"
			."ExtractCrystalStructuresFromOptimiseOutput2 [Ret: $ret]\n");
		return -2;
	}

	unless($InitialCrystalCIF->WriteSimpleCIFFile($InitialFile)) {
		$App->print("  Error: Could not write to ($InitialFile).\n");
		return -1;
	}
	unless($OptimizedCrystalCIF->WriteSimpleCIFFile($OptimizedFile)) {
		$App->print("  Error: Could not write to ($OptimizedFile).\n");
		return -1;
	}
}

sub MergeCSV
{
	$App->print("\n\nMerge hisotry CSV files in directories:\n");

	my $GULP = new GULP;

	my $InFile  = $Args->GetGetArg(0);
	my $OutFile = $Args->GetGetArg(1);
	
	print "  Merged csv will be saved to [$OutFile].\n";
	my $out = JFile->new($OutFile, 'w') or return undef;
	my $lasttime   = 0;
	my $offsettime = 0;
	my $ifc = 1;
	my $totalcount = 1;
	for(my $c = 2 ; ; $c++) {
		my $d = $Args->GetGetArg($c);
		last if(!defined $d);
		
		my @dirs = sort glob($d);
		for(my $d = 0 ; $d < @dirs ; $d++) {
			my $path = Deps::MakePath($dirs[$d], $InFile, 0);
			print "  $ifc: $path [offset time = $offsettime]\n";
			my $in = JFile->new($path, 'r') or last;
			for(my $i = 0 ; ; $i++) {
				my $line = $in->ReadLine();
				if(!defined $line) {
					$offsettime = $lasttime;
					last;
				}
			
				if($i == 0 and $ifc > 1) {
					next;
				}
				elsif($i == 0) {
					$out->print("step,itot,int. time,$line");
					next;
				}
				
				Utils::DelSpace($line);
				my @a = Utils::Split("\\s*,\\s*", $line);
				$lasttime = $offsettime + $a[1];
				$out->print("$ifc,$totalcount,$lasttime,$a[0],$a[1]");
				for(my $j = 2 ; $j < @a ; $j++) {
					$out->print(",$a[$j]");
				}
				$out->print("\n");
				$totalcount++;
			}
			$in->Close();
			$ifc++;
		}
	}
	
	$out->Close();
	
	return 1;
}

sub CalMSD
{
	$App->print("\n\nCalculate mean squared displacement from MD output/history:\n");

	my $GULP = new GULP;

	my $nOutputInterval = $Args->GetGetArg('nOutputInterval');
	$nOutputInterval = 1 if($nOutputInterval == 0);
	my $InFile  = $Args->GetGetArg(0);
	my $OutFile = $Args->GetGetArg(1);
	my $HisFile = $Args->GetGetArg(2);
	my $MSDCSV  = $Args->GetGetArg(3);

	$App->print("  Read initial crystal structure from [$InFile].\n");
	$App->print("  Read output                    from [$OutFile].\n");
	$App->print("  Read hisotry                   from [$HisFile].\n");
	$App->print("  Save MSD to                      to [$MSDCSV].\n");
	$App->print("  nOutputInterval: $nOutputInterval\n");

	unless($MSDCSV) {
		$App->print("File name should be specified.\n");
		$App->print("  usage: GULP.pl --Action=CalculateMSD --nOutputInterval=val "
				   ."InFile OutputFile HisFile CSVFile\n");
		return 0;
	}

	print "Read configuration from [$OutFile].\n";
	my $pHash = $GULP->ReadOutputToHash($OutFile);
#	my $pHash = $GULP->ReadInputToHash($InFile);

	if(defined $pHash->{md}) {
	}
	else {
		$App->print("  Nothing to do for postprocessing fit calculation\n");
		exit;
	}

	print "Delete old [$MSDCSV]..\n";
	unlink($MSDCSV);
	print "\n";


	my $pOutHashArray = $GULP->ReadOutFileToHashArray($OutFile);
	my (@time, @temperature, @pressure);
	for(my $i = 0 ; $i < @$pOutHashArray ; $i++) {
		my $p = $pOutHashArray->[$i];
		$time[$i] = $p->{time};
		$temperature[$i] = $p->{T}{Instantaneous};
		$pressure[$i]    = $p->{P}{Instantaneous};
	}
my $out = JFile->new("temp.csv", 'w');
$out->print("time,T,P\n");
for(my $i = 0 ; $i < @time ; $i++) {
$out->print("$time[$i],$temperature[$i],$pressure[$i]\n");
}
$out->Close();
#exit;

	print "Read initial crystal structure from [$OutFile].\n";
	my $IniCrystal = $GULP->ReadCrystalStructureFromOutput($OutFile);
#	print "Read initial crystal structure from [$InFile].\n";
#	my $IniCrystal = $GULP->ReadCrystalStructureFromInput($InFile);
	print "\n";

#	my $OutputMode = 'site';
	my $OutputMode = 'type';
	
	my $in = JFile->new($HisFile, 'r') or die "Can not read [$HisFile].\n";
	my $csv = JFile->new($MSDCSV, 'w') or die "Can not read [$MSDCSV].\n";

	my ($Crystal, $ret) = $GULP->ReadNextCrystalFromHistoryFile(0, $in, $IniCrystal, 0);
	my $pAtomType = $Crystal->GetpAtomType();
	my $pAtomSite = $Crystal->GetpAtomSite(1);
#	$GULP->InitializeMSD($Crystal);
	$GULP->InitializeMSD();

	$csv->print("iStep,iTotalStep,time,T(K),P");
	if($OutputMode eq 'type') {
		for(my $i = 0; $i < @$pAtomType ; $i++) {
			$csv->print(",msd(", $pAtomType->[$i]->AtomNameOnly(), ")");
		}
	}
	elsif($OutputMode eq 'site') {
		for(my $i = 0; $i < @$pAtomSite ; $i++) {
			$csv->print(",msd(", $pAtomSite->[$i]->Label(), ")");
		}
	}
	$csv->print(",msd(total)\n");

	my $iStep = 1;
	my $msd = 0.0;
#	for(my $i = 0 ; $i < 15 ; $i++) {
	for(my $i = 0 ; ; $i++) {
print "iStep=$iStep\n";
		if($i % $nOutputInterval == 0) {
#			($Crystal, $ret) = $GULP->ReadNextCrystalFromHistoryFile($iStep, $in, $IniCrystal, 1);
			($Crystal, $ret) = $GULP->ReadNextCrystalFromHistoryFile($iStep, $in, $IniCrystal, 0);
			if(!$ret or !$Crystal) {
				last;
			}

print "  iTotalStep=$Crystal->{iStep}  Time=$Crystal->{Time}\n";
			$GULP->AddMSD($Crystal);
			my ($pAtomSiteMSD, $pAtomTypeMSD, $TotalMSD) = $GULP->GetMSD($Crystal);
			my $temperature = Algorism::Interpolate(\@time, \@temperature, $Crystal->{Time});
			my $pressure    = Algorism::Interpolate(\@time, \@pressure,    $Crystal->{Time});
			$csv->print("$iStep,$Crystal->{iStep},$Crystal->{Time},$temperature,$pressure");
			if($OutputMode eq 'type') {
				for(my $j = 0 ; $j < @$pAtomTypeMSD ; $j++) {
					$csv->print(",$pAtomTypeMSD->[$j]");
				}
			}
			elsif($OutputMode eq 'site') {
				for(my $j = 0 ; $j < @$pAtomSiteMSD ; $j++) {
					$csv->print(",$pAtomSiteMSD->[$j]");
				}
			}
			$csv->print(",$TotalMSD\n");
		}
		else {
			my $line = $in->SkipTo("timestep");
			last if(!defined $line);
		}
		$iStep++;
	}

}

sub MakeHistoryCSV
{
	&PostProcess(@_);
}

sub PostProcess
{
	$App->print("\n\nMake Hisotry CSV/CIF/AXSF files from output:\n");

	my $GULP = new GULP;

	my $nOutputInterval = $Args->GetGetArg('nOutputInterval');

	my $InFile         = $Args->GetGetArg(0);
	my $OutFile        = $Args->GetGetArg(1);
	my $HisFile        = $Args->GetGetArg(2);
	my $InitialCIFFile = $Args->GetGetArg(3);
	my $FinalCIFFile   = $Args->GetGetArg(4);
	my $XSFFile        = $Args->GetGetArg(5);
	my $HistoryCSVIns  = $Args->GetGetArg(6);
	my $HistoryCSVAvr  = $Args->GetGetArg(7);

	print "Read configuration from [$OutFile].\n";
	my $pHash = $GULP->ReadOutputToHash($OutFile);
#	my $pHash = $GULP->ReadInputToHash($InFile);
#print "target:$pHash->{fit}/$pHash->{optimize}/$pHash->{md}\n";
#exit;
	if(defined $pHash->{fit}) {
		$App->print("  Nothing to do for postprocessing fit calculation\n");
		exit;
	}
	elsif(defined $pHash->{optimize}) {
		&PostProcessForOptimize($InFile, $OutFile, $InitialCIFFile, $FinalCIFFile);
	}
	elsif(defined $pHash->{md}) {
		&PostProcessForMD($InFile, $OutFile, $HisFile, $InitialCIFFile, $FinalCIFFile, $XSFFile,
			$HistoryCSVIns, $HistoryCSVAvr, $nOutputInterval);
	}
	else {
		print "  Error: Postprocess target not specified\n";
		exit;
	}
}

sub PostProcessForOptimize
{
	my ($InFile, $OutFile, $InitialCIFFile, $FinalCIFFile) = @_;

	my $GULP = new GULP;

	$App->print("  Read initial crystal structure from [$InFile].\n");
	$App->print("  Read output                    from [$OutFile].\n");
	$App->print("  Save initial structure         to $InitialCIFFile.\n");
	$App->print("  Save final structure           to $FinalCIFFile.\n");

	print "Delete old [$InitialCIFFile]..\n";
	print "Delete old [$FinalCIFFile]..\n";
	unlink($InitialCIFFile);
	unlink($FinalCIFFile);
	print "\n";

	my $CIF = new CIF;

	print "Read initial crystal structure from [$OutFile].\n";
	my $IniCrystal = $GULP->ReadCrystalStructureFromOutput($OutFile);
#	print "Read initial crystal structure from [$InFile].\n";
#	my $IniCrystal = $GULP->ReadCrystalStructureFromInput($InFile);
	print "Save initial crystal structure to [$InitialCIFFile].\n";
	my $ret = $CIF->CreateCIFFileFromCCrystal($IniCrystal, $InitialCIFFile, 0);
	if(!$ret) {
		print "Error: Can not write to [$InitialCIFFile].\n";
		exit;
	}
	print "\n";

	print "Read final crystal structure from [$OutFile].\n";
	my $FinalCrystal = $GULP->ReadFinalCrystalStructureFromOutput($OutFile);
	print "Save final crystal structure to [$FinalCIFFile].\n";
	$ret = $CIF->CreateCIFFileFromCCrystal($FinalCrystal, $FinalCIFFile, 0);
	if(!$ret) {
		print "Error: Can not write to [$FinalCIFFile].\n";
		exit;
	}
}

sub PostProcessForMD
{
	my ($InFile, $OutFile, $HisFile, $InitialCIFFile, $FinalCIFFile, $XSFFile,
		$HistoryCSVIns, $HistoryCSVAvr, $nOutputInterval) = @_;

	my $GULP = new GULP;

	$App->print("  Read initial crystal structure from [$InFile].\n");
	$App->print("  Read output                    from [$OutFile].\n");
	$App->print("  Read hisotry                   from [$HisFile].\n");
	$App->print("  Save initial structure         to $InitialCIFFile.\n");
	$App->print("  Save final structure           to $FinalCIFFile.\n");
	$App->print("  Save Relaxation/MD animation   to $XSFFile.\n");
	$App->print("  Save hisotry (instantaneous)   to $HistoryCSVIns.\n");
	$App->print("  Save hisotry (average)         to $HistoryCSVAvr.\n");

	unless($HistoryCSVAvr) {
		$App->print("File name should be specified.\n");
		$App->print("  usage: GULP.pl --Action=MakeHistoryCSV InFile OutputFile HisFile "
				   ."InitialCIFFile FinalCIFFile XSFFile HistoryCSVFile(instantaneous) HistoryCSVFile(average).\n");
		return 0;
	}

	print "Delete old [$InitialCIFFile]..\n";
	print "Delete old [$FinalCIFFile]..\n";
	print "Delete old [$XSFFile]..\n";
	print "Delete old [$HistoryCSVIns]..\n";
	print "Delete old [$HistoryCSVAvr]..\n";
	unlink($InitialCIFFile);
	unlink($FinalCIFFile);
	unlink($XSFFile);
	unlink($HistoryCSVIns);
	unlink($HistoryCSVAvr);
	print "\n";

	print "Read initial crystal structure from [$OutFile].\n";
	my $IniCrystal = $GULP->ReadCrystalStructureFromOutput($OutFile);
#	print "Read initial crystal structure from [$InFile].\n";
#	my $IniCrystal = $GULP->ReadCrystalStructureFromInput($InFile);
	print "\n";

	my $CIF = new CIF;
	print "Save initial crystal structure to [$InitialCIFFile].\n";
	my $ret = $CIF->CreateCIFFileFromCCrystal($IniCrystal, $InitialCIFFile, 0);
	if(!$ret) {
		print "Error: Can not write to [$InitialCIFFile].\n";
		exit;
	}

	print "Save hisotry to [$HistoryCSVIns] and [$HistoryCSVIns] from [$OutFile].\n";
	$GULP->SaveHistoryFromOutput($IniCrystal, $OutFile, $HistoryCSVIns, $HistoryCSVAvr);

	my $FinalCrystal = $GULP->MakeXSFFromHistoryFile($HisFile, $XSFFile, $IniCrystal, $nOutputInterval);

	print "Save final crystal structure to [$FinalCIFFile]..\n";
	my $ret = $CIF->CreateCIFFileFromCCrystal($FinalCrystal, $FinalCIFFile, 0);
	if(!$ret) {
		print "Error: Can not write to [$FinalCIFFile].\n";
		exit;
	}
}

sub ModifyInput
{
	$App->print("\n\Modify GULP input File:\n");

	my $GULP = new GULP;

	my $InpFile  = $Args->GetGetArg(0);
	my $NextFile = $Args->GetGetArg(1);
	my $PrevFile = "$InpFile.prev";

	my $Function       = $Args->GetGetArg("Function");
	my $temperature    = $Args->GetGetArg("temperature");
	my $pressure       = $Args->GetGetArg("pressure");
	my $equilibration  = $Args->GetGetArg("equilibration");
	my $production     = $Args->GetGetArg("production");
	my $timestep       = $Args->GetGetArg("timestep");
	my $tscale         = $Args->GetGetArg("tscale");
	my $sample         = $Args->GetGetArg("sample");
	my $write          = $Args->GetGetArg("write");
	my $ResetTime      = $Args->GetGetArg("ResetTime");

	if(!defined $Function or $Function eq '') {
		print "Read configuration from [$InpFile].\n";
#		my $pHash = $GULP->ReadOutputToHash($OutFile);
		my $pHash = $GULP->ReadInputToHash($InpFile);
		$Function = $pHash->{Task};
	}

	print "  Update for [$Function]\n";

	my $pHash = $GULP->ReadInputToHash($InpFile);
#foreach my $key (sort keys %$pHash) {
#print "$key: $pHash->{$key}\n";
#}
	my ($originalproduction, $punit) = Utils::Split("\\s+", $pHash->{production});
	my $newproduction;
	my $newtimerange;
print "p:$production $ResetTime\n";
	if($ResetTime) {
		$newproduction = $production + 0;
		$newtimerange  = $production + 0;
	}
	else {
		$newproduction = ($production =~ /^+/)? $production + $originalproduction : $production;
		$newtimerange = ($production =~ /^+/)? $production+0 : $production - $originalproduction;
	}
	my $nstep = int($newtimerange / $timestep) + 1;
#print "nstep = int($newtimerange / $timestep) + 1\n";

	print("  Input file     : $InpFile, backup to $PrevFile.\n");
	print("  Next input file: $NextFile\n");
	print("  Function       : $Function\n");
#	print("  equilibration  : $equilibration\n");
	print("  production     : $production\n");
	print("    original production is $originalproduction $punit, to be updated to $newproduction $punit\n");
	print("  timestep       : $timestep\n");
	print("    additional prduction time is $newtimerange $punit\n");
	print("    nstep = $nstep\n");
	print("  tscale         : $tscale\n");
	print("  temperature    : $temperature\n");
	print("  pressure       : $pressure\n");
	print("  sample         : $sample\n");
	print("  write          : $write\n");
	print("  Reset MD time  : $ResetTime\n");

	unless($NextFile) {
		$App->print("File name should be specified.\n");
		$App->print("  usage: GULP.pl --Action=ModifyInput InputFile NextInputFile parameters\n");
		return 0;
	}

	my ($T0, $T1) = Utils::Split("\\s*,\\s*", $temperature);
	my $Tstep;
	if(defined $T1) {
		$Tstep = ($T1 - $T0) / ($nstep - 1);
	}

	print("\n");
	print("Delete previous input [$PrevFile].\n");
	Utils::DeleteFile($PrevFile);
	print("Backup [$InpFile] to [$PrevFile].\n");
	if(!Utils::CopyFile($InpFile, $PrevFile)) {
#		print("  Error in GULP.pl::UpdateInputFile: Can not move [$InpFile] to [$PrevFile].\n");
#		exit;
	}

	my $Crystal  = undef;
	my $pPotHash = undef;

	$GULP->ModifyInput($InpFile, $NextFile, $Function, $Crystal, $pPotHash, 
		Function    => $Function,
		production  => $newproduction,
		timestep    => $timestep,
		tscale      => $tscale,
		temperature => $temperature,
		pressure    => $pressure,
		sample      => $sample,
		write       => $write,
		nstep       => $nstep,
		T0          => $T0,
		T1          => $T1,
		Tstep       => $Tstep,
		ResetTime => $ResetTime,
		);
}

sub UpdateInputFile
{
	$App->print("\n\nUpdate GULP input File:\n");

	my $GULP = new GULP;

	my $InpFile = $Args->GetGetArg(0);
	my $OutFile = $Args->GetGetArg(1);
	my $ResFile = $Args->GetGetArg(2);
	my $NewFile = $Args->GetGetArg(3);
	unless($NewFile) {
		$App->print("File name should be specified.\n");
		$App->print("  usage: GULP.pl --Action=UpdateInputFile InputFile OutputFile ResFile NewInputFile\n");
		return 0;
	}

	my $Function       = $Args->GetGetArg("Function");
	my $temperature    = $Args->GetGetArg("temperature");
	$temperature =~ s/,.*$//;
	my $pressure       = $Args->GetGetArg("pressure");
	my $timestep       = $Args->GetGetArg("timestep");
	my $tscale         = $Args->GetGetArg("tscale");
	my $equilibration  = $Args->GetGetArg("equilibration");
	my $production     = $Args->GetGetArg("production");
	my $sample         = $Args->GetGetArg("sample");
	my $write          = $Args->GetGetArg("write");

	my $CopyScript     = $Args->GetGetArg("CopyScript");

	print "Read initial crystal structure from [$OutFile].\n";
	my $Crystal = $GULP->ReadCrystalStructureFromOutput($OutFile);

	print "\n";
	print "Read configuration from [$InpFile].\n";
	my $pHash = $GULP->ReadOutputToHash($OutFile);
#	my $pHash = $GULP->ReadInputToHash($InpFile);
	if(defined $pHash->{fit}) {
		print "  Update for fit\n";
		$GULP->UpdateInputForFit($InpFile, $OutFile, $NewFile);
	}
	elsif(defined $pHash->{optimize}) {
		print "  Update for optimize\n";
		$GULP->UpdateInputForOptimize($InpFile, $OutFile, $NewFile,
			Function    => $Function,
			production  => $production,
			timestep    => $timestep,
			tscale      => $tscale,
			temperature => $temperature,
			pressure    => $pressure,
			sample      => $sample,
			write       => $write,
			);
	}
	elsif(defined $pHash->{md}) {
		print "  Update for md\n";
		$GULP->UpdateInputForMD($InpFile, $OutFile, $ResFile, $NewFile,
			Function    => $Function,
			production  => $production,
			timestep    => $timestep,
			tscale      => $tscale,
			temperature => $temperature,
			pressure    => $pressure,
			sample      => $sample,
			write       => $write,
			);
	}
	else {
		print "  Error: Update target not specified\n";
		exit;
	}
}

sub MakeInputFile
{
	$App->print("\n\nMake GULP .inp File:\n");

	my $GULP = new GULP;

	my $Action                 = $Args->GetGetArg("Action");
	my $Function               = $Args->GetGetArg("Function");
	$Function = "Optimize"      unless($Function);
	my $UseShellModelForCation = $Args->GetGetArg("UseShellModelForCation");
	$UseShellModelForCation = 0 unless($UseShellModelForCation);
	my $UseShellModelForAnion  = $Args->GetGetArg("UseShellModelForAnion");
	$UseShellModelForAnion  = 0 unless($UseShellModelForAnion);
	my $LibraryFile            = $Args->GetGetArg("LibraryFile");
	$LibraryFile = "Kamiya-NoShell.lib" unless($LibraryFile);
	for(my $i = @PotentialLibs - 1 ; $i >= 0  ; $i--) {
		my $path = Deps::MakePath($PotentialLibs[$i], $LibraryFile, 0);
		if(-f $path) {
			$LibraryFile = $path;
			last;
		}
	}
	if(-f $LibraryFile) {
		$GULP->SetLibraryFile($LibraryFile);
	}
	else {
		print "Error: Can not find library [$LibraryFile].\n";
		return;
	}
	my $TargetFiles = $Args->GetGetArg("TargetFiles");
	my $UseSpeciesCharges = $Args->GetGetArg("UseSpeciesCharges");
	   $UseSpeciesCharges = 'no' if(!defined $UseSpeciesCharges);

	my $timestep       = $Args->GetGetArg("timestep");
	my $tscale         = $Args->GetGetArg("tscale");
	my $equilibration  = $Args->GetGetArg("equilibration");
	my $production     = $Args->GetGetArg("production");
	my $sample         = $Args->GetGetArg("sample");
	my $write          = $Args->GetGetArg("write");
	my $CopyScript     = $Args->GetGetArg("CopyScript");

	my $CIFFile = $Args->GetGetArg(0);
	my $InpFile = $Args->GetGetArg(1);
	unless($InpFile) {
		$App->print("File name should be specified.\n");
		$App->print("  usage: GULP.pl --Action=MakeInp CIFFile StructFile.\n");
		return 0;
	}

#ファイル名を（ベース名, ディレクトリ名, 拡張子）に分解
	my ($drive, $dir, $filename, $ext, $lastdir, $filebody)
		= Deps::SplitFilePath($CIFFile);
	my $SampleName = $filebody;

	$App->print("  Read $CIFFile.\n");
	if($Function =~ /fit/i) {
		$App->print("  Target structures to be fitted: $TargetFiles\n");
	}
	$App->print("  Save to $InpFile.\n");
	$App->print("  Function: $Function\n");
	$App->print("  UseShellModelForCation: $UseShellModelForCation\n");
	$App->print("  UseShellModelForAnion : $UseShellModelForAnion\n");
	$App->print("  LibraryFile: $LibraryFile in [$GULPLibDir]\n");
	$App->print("  Use Species charges: $UseSpeciesCharges\n");

	$App->print("  timestep     : $timestep\n");
	$App->print("  tscale       : $tscale\n");
	$App->print("  equilibration: $equilibration\n");
	$App->print("  production   : $production\n");
	$App->print("  sample       : $sample\n");
	$App->print("  write        : $write\n");
	$App->print("  CopyScript   : $CopyScript\n");

	$App->print("\n");
	$App->print("  Read [$CIFFile].\n");

	my $CIF = new CIF();
	unless($CIF->Read($CIFFile)) {
		$App->print("Error: Can not read $CIFFile.\n\n");
		return 0;
	}

#CIFクラスの内容から、Crystalクラスを作成
	my $Crystal = $CIF->GetCCrystal();
	$Crystal->ExpandCoordinates();

	$App->print("  Save to [$InpFile].\n");
	$GULP->SetSampleName($SampleName);
	$GULP->SetUseShellModelForCation($UseShellModelForCation);
	$GULP->SetUseShellModelForAnion($UseShellModelForAnion);
	$GULP->SetLibraryDir($GULPLibDir);
	$GULP->SetLibraryFile($LibraryFile);
	$GULP->SaveInputFile($Crystal, $Function, $InpFile, 0,
				UseSpeciesCharges => $UseSpeciesCharges,
				timestep      => $timestep,
				tscale        => $tscale,
				equilibration => $equilibration,
				production    => $production,
				sample        => $sample,
				write         => $write,
				CopyScript    => $CopyScript,
				TargetFiles   => $TargetFiles,
				);

	$App->print("\n***Make GULP Input file Finished.\n\n");

	return;
}

sub FindDatabases
{
	my ($this, $IsPrint) = @_;
	$IsPrint = 0 if(!defined $IsPrint);

	my $PotentialType               = $Args->GetGetArg("PotentialTyp");
	my $CIFFile  = $Args->GetGetArg(0);
	unless($CIFFile) {
		$App->print("File name should be specified.\n");
		$App->print("  usage: GULP.pl --Action=FindDatabases --PotentialType=[buck|lennard] CIFFile\n");
		return 0;
	}

	$App->print("\nFind available potential databases from the chemical composition in [$CIFFile].\n");
	$App->print("  PotentialType: $PotentialType\n");
	$App->print("  Read $CIFFile...\n");

	my $CIF = new CIF();
	unless($CIF->Read($CIFFile)) {
		$App->print("Error: Can not read $CIFFile.\n\n");
		return 0;
	}

	my $GULP = new GULP;

#CIFクラスの内容から、Crystalクラスを作成
	my $Crystal = $CIF->GetCCrystal();
#	$Crystal->ExpandCoordinates();
	my @AtomTypeList           = $Crystal->GetCAtomTypeList();
	my $nAtomType              = @AtomTypeList;

	print "\n";
	print "Search databases for ";
	for(my $i = 0 ; $i < $nAtomType ; $i++) {
		my $name = $AtomTypeList[$i]->AtomNameOnly();
		print "$name ";
	}
	print "\n";

#	print "Databases in [$LAMMPSDatabaseDir] and [$GULPLibDir].\n" if($IsPrint);
#	my @files = sort glob(Deps::MakePath($LAMMPSDatabaseDir, '*', 0));
#	my @files = sort glob(Deps::MakePath($GULPLibDir, '*', 0));
	my @files;
	for(my $k = 0 ; $k < @PotentialLibs ; $k++) {
		my $dir = $PotentialLibs[$k];
		next if(!-d $dir);

		print "Search databases in [$dir]...\n";
		my @f = sort glob(Deps::MakePath($dir, '*', 0));
		for(my $i = 0 ; $i < @f ; $i++) {
			next if(!-f $f[$i]);
#			next if($f[$i] =~ /\.htm/i);
			next if($f[$i] !~ /\.lib/i);
			push(@files, $f[$i]);
		}
	}

	my $nHit = 0;
	for(my $i = 0 ; $i < @files ; $i++) {
		next if($files[$i] =~ /eledata$/i);
#next if($files[$i] !~ /kamiya/i);

#print "$files[$i]\n" if($IsPrint);
		my ($phash, $pAtomArray) = $GULP->ReadPotentialToHash($files[$i], $IsPrint);
		my $AllFound = 1;
		for(my $i = 0 ; $i < $nAtomType ; $i++) {
			my $name = $AtomTypeList[$i]->AtomNameOnly();
			my $IsFound = 0;
			for(my $j = 0 ; $j < @$pAtomArray ; $j++) {
				if($name eq $pAtomArray->[$j]) {
					$IsFound = 1;
					last;
				}
			}
			if(!$IsFound) {
				$AllFound = 0;
				last;
			}
		}
		if($AllFound) {
			print "All potentials are found in [$files[$i]].\n";
			print "   Available atoms:";
			for(my $i = 0 ; $i < @$pAtomArray ; $i++) {
				print " $pAtomArray->[$i]";
			}
			print "\n\n";
			$nHit++;
		}
	}
	print "\n$nHit databases are found.\n";
	print "\n";
}

sub ShowDatabases
{
	my ($this, $IsPrint) = @_;
	$IsPrint = 1 if(!defined $IsPrint);
	
	print "\n";
	print "Databases in [$LAMMPSDatabaseDir] and [$GULPLibDir].\n" if($IsPrint);
	my @files = sort glob(Deps::MakePath($LAMMPSDatabaseDir, '*', 0));
#	for(my $i = 0 ; $i < @files ; $i++) {
#		print "  $files[$i]\n";
#	}
	print "\n";

	my $GULP = new GULP;
	for(my $k = 0 ; $k < @PotentialLibs ; $k++) {
		my $dir = $PotentialLibs[$k];
		next if(!-d $dir);

		print "\n";
		print "Search [$dir]...\n";
		my @files = sort glob(Deps::MakePath($dir, '*', 0));
		for(my $i = 0 ; $i < @files ; $i++) {
			next if(!-f $files[$i]);
#			next if($files[$i] =~ /eledata$/i);
			next if($files[$i] =~ /\.lib/i);
#			next if($files[$i] =~ /\.htm/i);

			print "$files[$i]\n" if($IsPrint);
			my ($phash, $pAtomArray) = $GULP->ReadPotentialToHash($files[$i], $IsPrint);
			foreach my $key (sort keys %$phash) {
				my $p2 = $phash->{$key};
				print "    $key: $p2->{Type} $p2->{nAtom} ";
				for(my $j = 1 ; $j <= $p2->{nAtom} ; $j++) {
					print " " . $p2->{"AtomName$j"};
				}
				print "\n";
			}
			print "    atoms: ", join(', ', @$pAtomArray), "\n";
#die if($files[$i] =~ /morse/i);
		}
	}
}