package GULP;
use Exporter;
use Crystal::MD;
use Crystal::ForceField;
@ISA = qw(Exporter ForceField MD);

#公開したいサブルーチン
@EXPORT = qw();

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

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::Rietan;
use Crystal::XCrySDen;
use Sci qw($torad $e0 $a0 $e $pi log10 erfc);
use Sci::Science;
use Sci::EnergyBandArray;

#===============================================
# 大域変数
#===============================================
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();
my $GULPLibDir        = ProgVars::GULPLibDir();
my @PotentialLibs     = ForceField::PotentialLibs();

#===============================================
# コンストラクタ、デストラクタ
#===============================================
sub new
{
	my ($module) = @_;
	my $this = {};
	bless $this;
	$this->{pLibraryLines} = [];

	return $this;
}

sub DESTROY
{
	my $this = shift;
}

#============================================================
#　継承クラスで定義しなおす関数
#============================================================
sub CheckFileType
{
	my ($path) = @_;
	my ($drive, $dir, $filename, $ext, $lastdir, $filebody)
		= Deps::SplitFilePath($path);

	if($filename =~ /\.dens$/i) {
		return "GULP DENS file";
	}
	if($filename =~ /\.disp$/i) {
		return "GULP DISP file";
	}
	return undef;
}

sub ReadFiles
{
	my ($this, $filename, $TargetData) = @_;
	$this->ClearAll();
	my $FileType = $this->{'FileType'} = GULP::CheckFileType($filename);
	$this->SetFileName($filename);
	if($FileType eq "GULP DENS file") {
		return $this->ReadDENS($filename, $TargetData);
	}
	if($FileType eq "GULP DISP file") {
		return $this->ReadDISP($filename);
	}
	return undef;
}

#===============================================
# 変数取得関数
#===============================================
sub SetSampleName
{
	my ($this, $name) = @_;
	return $this->{'SampleName'} = $name;
}
sub SampleName
{
	my ($this) = @_;
	return $this->{'SampleName'};
}

#===============================================
# 一般関数
#===============================================
sub GetnAtomFromMDHistory
{
	my ($this, $fname) = @_;

	return 0 unless(open(IN,"<$fname"));
	while(<IN>) {
		last if($_ =~ /^\s*timestep\s/i);
	}
	<IN>;
	<IN>;
	<IN>;
	my $count = 0;
	while(<IN>) {
		last if($_ =~ /^\s*timestep\s/i);
		<IN>;
		<IN>;
		$count++;
	}
	close(IN);
	return $count;
}

sub ReadDENS
{
	my ($this, $path) = @_;

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

	$this->SetDataArray(new GraphDataArray);
	my $Data0 = new GraphData;

	my ($drive, $directory, $filename, $ext, $lastdir, $filebody)
		= Deps::SplitFilePath($path);
	my $SampleName = $filename;
	Utils::DelSpace($SampleName);
	$Data0->SetTitle($SampleName);

	my $line = $in->ReadLine();
	$line = $in->ReadLine();

	my @wn;
	my @DOS;
	my $i = 0;
	while(!$in->eof()) {
		$line = $in->ReadLine();
		last if(not defined $line);
		my ($w,$d) = Utils::Split("\\s+", $line);
		last if(not defined $d);
		$wn[$i]  = $w;
		$DOS[$i] = $d;
print "$i: $w $d\n";
		last if($in->eof());
		$i++;
	}

	$Data0->{'x0'}      = \@wn;
	$Data0->{'x0_Name'} = "Wave number Energy / eV";
	$Data0->{'y0'}      = \@DOS;
	$Data0->{'y0_Name'} = "DOS";
	$Data0->CalMinMax();
	$this->{'DataArray'}->AddGraphData($Data0);

	$in->Close();
	return $path;
}

sub ReadDISP
{
	my ($this, $path) = @_;

	my $outfile = Deps::ReplaceExtension($path, ".out");
print "Read .out [$outfile]\n";
	my $CIF = $this->ExtractCrystalStructuresFromOptimiseOutput($outfile, 1);
	if(!$CIF) {
		print "Warning!!: Can not read [$outfile].\n";
	}
	my $Crystal;
	$Crystal = $CIF->GetCCrystal() if(!$CIF);

	if(!defined $Crystal) {
		my $glpfile = Deps::ReplaceExtension($path, ".glp");
print "Read .glp [$glpfile]\n";
		$Crystal = $this->GetCCrystalFromInput($glpfile);
	}

#print "Crystal: $Crystal\n";
	my ($a,$b,$c,$alpha,$beta,$gamma) = $Crystal->LatticeParameters();
print "Latt: $a $b $c  $alpha $beta $gamma\n";

	my ($drive, $directory, $filename, $ext, $lastdir, $filebody)
		= Deps::SplitFilePath($path);
	my $SampleName = $filename;
	Utils::DelSpace($SampleName);

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

	my $nK = 0;
	while(!$in->eof()) {
		$line = $in->SkipTo("#  Section number = ");
		last if(!defined $line);
		$nK++;
	}
print "nK: $nK\n";

	$in->rewind();
	$line = $in->SkipTo("#  Final K point = ");
	my $nData = 0;
	my $nKMesh = 0;
	while(!$in->eof()) {
		my $line = $in->ReadLine();
		last if(!defined $line);
		last if($line =~ /^#\s/);
		$nData++;
		($nKMesh) = ($line =~ /(\d+)/);
	}
	my $nLevels = $nData / $nKMesh;
print "nLevels: $nLevels\n";
print "nKMesh: $nKMesh\n";

	$in->rewind();

	my $EnergyBandArray = new EnergyBandArray;
	$EnergyBandArray->SetTitle($SampleName);
	$EnergyBandArray->SetCrystal($Crystal);
	for(my $iL = 0 ; $iL < $nLevels ; $iL++) {
		$EnergyBandArray->AddEnergyBand();
	}

	my $Distance = 0.0;
#	$nK = 1;
#	$nKMesh = 15;
	for(my $iK = 0 ; $iK < $nK ; $iK++) {
##  Section number =   1  Configuration =   1
	$line = $in->ReadLine();
##  Start K point =   0.500000  0.000000  0.000000
		$line = $in->ReadLine();
		($line) = ($line =~ /=\s*(\w.*)$/);
		my ($kx0, $ky0, $kz0) = Utils::Split("\\s+", $line);
##  Final K point =   0.000000  0.000000  0.000000
		$line = $in->ReadLine();
		($line) = ($line =~ /=\s*(\w.*)$/);
		my ($kx1, $ky1, $kz1) = Utils::Split("\\s+", $line);
print "k: $kx0,$ky0,$kz0 - $kx1,$ky1,$kz1\n";

		for(my $im = 0 ; $im < $nKMesh ; $im++) {
			my $kx = $kx0 + ($kx1-$kx0)/($nKMesh-1) * $im;
			my $ky = $ky0 + ($ky1-$ky0)/($nKMesh-1) * $im;
			my $kz = $kz0 + ($kz1-$kz0)/($nKMesh-1) * $im;
			for(my $iL = 0 ; $iL < $nLevels ; $iL++) {
				$line = $in->ReadLine();
				my ($ib,$wn) = Utils::Split("\\s+", $line);
				$EnergyBandArray->AddByKE($iL, $kx, $ky, $kz, $wn);
#print "  $kx,$ky,$kz: $wn\n";
			}
		}
	}
	$EnergyBandArray->CalMinMax();
	$EnergyBandArray->GetBandBoundary();

	my $DataArray = new GraphDataArray;
	$this->SetDataArray($DataArray);

	$EnergyBandArray->SetGraphDataArray($DataArray);

	return $path;

}

sub GetMDStepFromMDHistory
{
	my ($this, $fname) = @_;

	return 0 unless(open(IN,"<$fname"));
	my $count = 0;
	while(<IN>) {
		$count++ if($_ =~ /^\s*timestep\s/i);
	}
	close(IN);
	return $count;
}

sub GetNextCrystalStructureFromMDHistory
{
	my ($this, $in, $nAtom) = @_;

	my ($a,$b,$c,$alpha,$beta,$gamma);
	my $TotalTime;
#	my $CIF = new CIF();

	my $SPGName = "P 1";
	my $Rietan = new Rietan();
	my $SPG = $Rietan->ReadSpaceGroupFromSPGName($SPGName);
	if($SPG == 0) {
		return -2;
	}
	my $iSPG = $SPG->iSPG();
	my $iSet = $SPG->iSet();
	my $CIF = new CIF();
#	$CIF->SetCrystalName($Title);
	$CIF->SetSpaceGroupInformation($SPG);

	while(<$in>) {
		last if($_ =~ /^timestep\s/);
	}
	$_ =~ /(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)/;
	my $iTotalStep = $2;
	my $TimeStep   = $6;
	$TotalTime = $iTotalStep * $TimeStep;
#print "TotalTime: $TotalTime = $iTotalStep * $TimeStep<br>";

	my $line = <$in>;
	my ($a11,$a12,$a13) = ($line =~ /(\S+)\s+(\S+)\s+(\S+)/);
	$line = <$in>;
#print "line: $line<br>\n";
	my ($a21,$a22,$a23) = ($line =~ /(\S+)\s+(\S+)\s+(\S+)/);
	$line = <$in>;
#print "line: $line<br>\n";
	my ($a31,$a32,$a33) = ($line =~ /(\S+)\s+(\S+)\s+(\S+)/);
	my $Crystal = new Crystal();
	$Crystal->SetLatticeVector($a11, $a12, $a13,
				   $a21, $a22, $a23,
				   $a31, $a32, $a33);
	($a,$b,$c,$alpha,$beta,$gamma) = $Crystal->CalculateLatticeConstantFromVector();
#print "Latt: $a $b $c $alpha $beta $gamma<br>";

	my %nAtomName;
	for(my $i = 0 ; $i < $nAtom ; $i++) {
		my $IsShell = 0;
		my $line1 = <$in>;
		my $line2 = <$in>;
#print "line2: $line2<br>\n" if($i == 0);
		my $line3 = <$in>;
#print "line1: $line1<br>\n";
		$line1 =~ /(\S+)/;
		my $name = $1;
		if($name =~ /_s$/) {
			$IsShell = 1;
			$name =~ s/_s$//;
		}
		$nAtomName{$name}++;
		my $label = sprintf "%s%d", $name, $nAtomName{$name};

		$line2 =~ /(\S+)\s+(\S+)\s+(\S+)/;
		my $xc = $1;
		my $yc = $2;
		my $zc = $3;
		my ($x,$y,$z) = $Crystal->CartesianToFractional($xc,$yc,$zc);
#print "$i: x: ($x,$y,$z)<br>\n";
		my $occ = 1.0;

		$line3 =~ /(\S+)\s+(\S+)\s+(\S+)/;
		my $vx = $1;
		my $vy = $2;
		my $vz = $3;
		unless($IsShell) {
			my $nSite = $CIF->AddAsymmetricAtomSiteWithVelocity(
					$label, $name, $x, $y, $z, $occ, $vx, $vy, $vz);
			$CIF->SetAsymmetricAtomSiteVelocity($nSite, $vx, $vy, $vz);
		}
	}
#print "a3<br>\n";

	return (1,1) if($a eq '');

	$CIF->SetLatticeParameters($a, $b, $c, $alpha, $beta, $gamma);
	$CIF->FillCIFData();

	return ($CIF,$TotalTime);
}

sub ExtractCrystalStructureFromMDHistory
{
	my ($this, $fname, $istep) = @_;
	my $Title;
	my $Formula;
	my $nAtom;
	my $Dimensionality;
	my $LatticeSystem;
	my $SPGName;
	my ($a, $b, $c, $alpha, $beta, $gamma);
	my $TotalTime;

	return (0,0) unless(open(IN,"<$fname"));
	my $line = <IN>;
	$line =~ s/^\s*//;
	$line =~ s/\s*$//;
	$Title = $line;
	$Title = 'undefined' if($Title eq '');
	<IN>;

	$SPGName = "P 1";
	my $Rietan = new Rietan();
	my $SPG = $Rietan->ReadSpaceGroupFromSPGName($SPGName);
	if($SPG == 0) {
		return -2;
	}
	my $iSPG = $SPG->iSPG();
	my $iSet = $SPG->iSet();
	my $CIF = new CIF();
	$CIF->SetCrystalName($Title);
	$CIF->SetSpaceGroupInformation($SPG);
	my $count = 0;
	while(<IN>) {
		$count++ if($_ =~ /^timestep\s/);
		last if($count >= $istep);
	}
	$_ =~ /(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)/;
	my $iTotalStep = $2;
	my $TimeStep   = $6;
	$TotalTime = $iTotalStep * $TimeStep;
#print "TotalTime: $TotalTime = $iTotalStep * $TimeStep<br>";

	while(<IN>) {
		my $line = $_;
		my ($a11,$a12,$a13) = ($line =~ /(\S+)\s+(\S+)\s+(\S+)/);
		$line = <IN>;
		my ($a21,$a22,$a23) = ($line =~ /(\S+)\s+(\S+)\s+(\S+)/);
		$line = <IN>;
		my ($a31,$a32,$a33) = ($line =~ /(\S+)\s+(\S+)\s+(\S+)/);
		my $Crystal = new Crystal();
		$Crystal->SetLatticeVector($a11, $a12, $a13,
					   $a21, $a22, $a23,
					   $a31, $a32, $a33);
		($a,$b,$c,$alpha,$beta,$gamma) = $Crystal->CalculateLatticeConstantFromVector();
#print "Latt: $a $b $c $alpha $beta $gamma<br>";
		my $PrevPos = -1;
		$nAtom = 0;
		my %nAtomName;
		last unless($line);
		for(my $i = 0 ; $i < 10000 ; $i++) {
			my $IsShell = 0;
			$line = <IN>;
			last unless($line);
			if($line =~ /^timestep\s/) {
				seek(IN, $PrevPos, 0) if($PrevPos > 0);
				last;
			}
			$line =~ /(\S+)/;
			my $name = $1;
			if($name =~ /_s$/) {
				$IsShell = 1;
				$name =~ s/_s$//;
			}
			$nAtomName{$name}++;
			my $label = sprintf "%s%d", $name, $nAtomName{$name};
			$nAtom++;
			my $line1 = <IN>;
			$line1 =~ /(\S+)\s+(\S+)\s+(\S+)/;
			my $xc = $1;
			my $yc = $2;
			my $zc = $3;
			my ($x,$y,$z) = $Crystal->CartesianToFractional($xc,$yc,$zc);
			my $occ = 1.0;
			my $line2 = <IN>;
			$line2 =~ /(\S+)\s+(\S+)\s+(\S+)/;
			my $vx = $1;
			my $vy = $2;
			my $vz = $3;
#print "$i: $label : $name : ($x,$y,$z) ($vx,$vy,$vz) [$occ]<br>\n";
			unless($IsShell) {
				my $nSite = $CIF->AddAsymmetricAtomSiteWithVelocity(
						$label, $name, $x, $y, $z, $occ, $vx, $vy, $vz);
#				$CIF->SetAsymmetricAtomSiteVelocity($nSite, $vx, $vy, $vz);
			}
			$PrevPos = tell(IN);
		}
		last;
	}

	close(IN);
	return (1,1) if($a eq '');

	$CIF->SetLatticeParameters($a, $b, $c, $alpha, $beta, $gamma);
	$CIF->FillCIFData();

	return ($CIF,$TotalTime);
}

sub ExtractCrystalStructuresFromOptimiseOutput
{
	my ($this, $fname, $iFileType) = @_;
	return $this->ExtractCrystalStructureFromOptimiseOutput2($fname) if($iFileType == 2);
	return $this->ExtractCrystalStructureFromOptimiseOutput1($fname) if($iFileType == 1);
}

sub ExtractCrystalStructureFromOptimiseOutput2
{
	my ($this, $fname) = @_;

	my $Title;
	my $Formula;
	my $nAtom;
	my $nAsymmetricAtomsShells;
	my $nAtomsShells;
	my $Dimensionality;
	my $LatticeSystem;
	my $SPGName;
	my ($a, $b, $c, $alpha, $beta, $gamma);

	return 0 unless(open(IN,"<$fname"));
	my $sepcount = 0;
	while(<IN>) {
		$sepcount++ if($_ =~ /^\*{80}/);
		if($sepcount == 3) {
			my $line = <IN>;
			($Title) = ($line =~ /([^\*\s]+)/);
			last;
		}
	}
	$Title = 'undefined' if($Title eq '');
	while(<IN>) {
		last if($_ =~ /^\*\s+Input for Configuration/i);
	}
	while(<IN>) {
#print "line1: $_<br>\n";
		if($_ =~ /^\s*Formula\s+=\s+(\S*)\s/i) {
			$Formula = $1;
		}
		elsif($_ =~ /^\s*Number of irreducible atoms\/shells\s*=\s*(\S*)\s/i) {
			$nAsymmetricAtomsShells = $1;
		}
		elsif($_ =~ /^\s*Total number atoms\/shells\s+=\s+(\S*)\s/i) {
			$nAtomsShells = $1;
		}
		elsif($_ =~ /^\s*Dimensionality\s+=\s+(\S*)\s/i) {
			$Dimensionality = $1;
			if($Dimensionality != 3) {
				close(IN);
				return -1;
			}
		}
		elsif($_ =~ /^\s*Crystal\s+family\s*:\s*(\S*)\s/i) {
			$LatticeSystem = $1;
		}
		elsif($_ =~ /^\s*Space\s+group.*?:\s*(.*)$/i) {
			$SPGName = $1;
			$SPGName =~ s/\s*$//;
		}
		if($_ =~ /^\s*Final asymmetric unit/i or $_ =~ /^\s*Final fractional/i) {
			last;
		}
	}
	while(<IN>) {
		if($_ =~ /^\s*Label\s+\(Frac\)/i) {
			last;
		}
	}
	<IN>;

	my $Rietan = new Rietan();
	Utils::DelSpace($SPGName);
	$SPGName = "P 1" if(not defined $SPGName or $SPGName eq '');
	my $SPG = $Rietan->ReadSpaceGroupFromSPGName($SPGName);
	if($SPG == 0) {
		print "Error!!: Can not find the space group [$SPGName]\n";
		return -2;
	}
	my $iSPG = $SPG->iSPG();
	my $iSet = $SPG->iSet();
	my $CIF = new CIF();
#	$CIF->SetFileName($FileName);
	$CIF->SetCrystalName($Title);
	$CIF->SetFormula($Formula);
	$CIF->SetSpaceGroupInformation($SPG);

	$nAtom = 0;
	my %nAtomName;
#print "nAsymmetricAtomsShells: $nAsymmetricAtomsShells<br>\n";
	for(my $i = 0 ; $i < $nAsymmetricAtomsShells ; $i++) {
		my $line = <IN>;
		$line =~ /(\d+)[\s\*]+(\S+)[\s\*]+(\S+)[\s\*]+(\S+)[\s\*]+(\S+)[\s\*]+(\S+)/i;
		my $name = $2;
		my $type = $3;
		my $x = $4;
		my $y = $5;
		my $z = $6;
		$name =~ s/\d+//;
		$nAtomName{$name}++;
		my $label = sprintf "%s%d", $name, $nAtomName{$name};
		my $occ = 1.0;
		next if($type =~ /s/i);
#print "$i: $label: $name [t=$type] ($x,$y,$z) [$occ]<br>\n";
		$CIF->AddAsymmetricAtomSite($label, $name, $x, $y, $z, $occ);
	}

	while(<IN>) {
#print "line: $_<br>\n";
		if($_ =~ /\s*Final cell parameters/i) {
			last;
		}
	}
	while(<IN>) {
		if($_ =~ /-{10}/) {
			last;
		}
	}
	my $line = <IN>;
#print "line: $line\n";
	($a) = ($line =~ /\S+\s+([\d\.\+-eE]+)/);
	$line = <IN>;
	($b) = ($line =~ /\S+\s+([\d\.\+-eE]+)/);
	$line = <IN>;
	($c) = ($line =~ /\S+\s+([\d\.\+-eE]+)/);
	$line = <IN>;
	($alpha) = ($line =~ /\S+\s+([\d\.\+-eE]+)/);
	$line = <IN>;
	($beta)  = ($line =~ /\S+\s+([\d\.\+-eE]+)/);
	$line = <IN>;
	($gamma) = ($line =~ /\S+\s+([\d\.\+-eE]+)/);

	while(<IN>) {
		if($_ =~ /^\s*Non-primitive lattice parameters\s*:/i) {
			<IN>;
			my $line = <IN>;
#print "l:$line<br>\n";
			($a,$b,$c) = ($line =~ /a\s*=\s*(\S+)\s+b\s*=\s*(\S+)\s+c\s*=\s*(\S+)/i);
			$line = <IN>;
#print "l:$line<br>\n";
			($alpha,$beta,$gamma) = ($line =~ /alpha\s*=\s*(\S+)\s+beta\s*=\s*(\S+)\s+gamma\s*=\s*(\S+)/i);
		}
	}
#print "Lattice(GULP.pm): $a $b $c $alpha $beta $gamma<br>\n";

	$CIF->SetLatticeParameters($a, $b, $c, $alpha, $beta, $gamma);

	$CIF->FillCIFData();

	close(IN);
	return $CIF;
}

sub ExtractCrystalStructureFromOptimiseOutput1
{
	my ($this, $fname) = @_;
my $Debug=1;
	my $Title;
	my $Formula;
	my $nAtom;
	my $nAtomsShells;
	my $Dimensionality;
	my $LatticeSystem;
	my $SPGName;
	my ($a, $b, $c, $alpha, $beta, $gamma);

	return 0 unless(open(IN,"<$fname"));
	my $sepcount = 0;
	while(<IN>) {
		$sepcount++ if($_ =~ /^\*{80}/);
#print "line: $_ : $sepcount<br>\n";
		if($sepcount == 3) {
			my $line = <IN>;
#print "line: $line<br>\n";
			($Title) = ($line =~ /([^\*\s]+)/);
print "Title: $Title<br>\n" if($Debug);
			last;
		}
	}
	$Title = 'undefined' if($Title eq '');
	while(<IN>) {
		last if($_ =~ /^\*\s+Input for Configuration/i);
	}
	while(<IN>) {
#print "line: $_<br>\n";
		if($_ =~ /^\s*Formula\s+=\s+(\S*)\s/i) {
			$Formula = $1;
print "Formula: $Formula<br>\n" if($Debug);
		}
		elsif($_ =~ /^\s*Total number atoms\/shells\s+=\s+(\S*)\s/i) {
			$nAtomsShells = $1;
print "nAtomsShells: $nAtomsShells<br>\n" if($Debug);
		}
		elsif($_ =~ /^\s*Dimensionality\s+=\s+(\S*)\s/i) {
			$Dimensionality = $1;
print "Dimensionality: $Dimensionality<br>\n" if($Debug);
			if($Dimensionality != 3) {
				close(IN);
				return -1;
			}
		}
		elsif($_ =~ /^\s*Crystal\s+family\s*:\s*(\S*)\s/i) {
			$LatticeSystem = $1;
print "LatticeSystem: $LatticeSystem<br>\n" if($Debug);
		}
		elsif($_ =~ /^\s*Space\s+group.*?:\s*(.*)$/i) {
			$SPGName = $1;
			$SPGName =~ s/\s*$//;
print "SPGName: $SPGName<br>\n" if($Debug);
		}
		elsif($_ =~ /^\s*Cell\s+parameters.*?:/i) {
			<IN>;
			my $line = <IN>;
			($a,$alpha) 
				= ($line =~ /a\s*=\s*([\d\.eE\+-]+)\s+alpha\s*=\s*([\d\.eE\+-]+)/);
			$line = <IN>;
			($b,$beta ) 
				= ($line =~ /b\s*=\s*([\d\.eE\+-]+)\s+beta\s*=\s*([\d\.eE\+-]+)/);
			$line = <IN>;
			($c,$gamma) 
				= ($line =~ /c\s*=\s*([\d\.eE\+-]+)\s+gamma\s*=\s*([\d\.eE\+-]+)/);
print "Lattice: $a $b $c $alpha $beta $gamma<br>\n" if($Debug);
		}
		elsif($_ =~ /^\s*Primitive\s+cell\s+parameters\s*:/i) {
			<IN>;
			my $line = <IN>;
			($a,$alpha) 
				= ($line =~ /alpha\s*=.*?a\s*=\s*([\d\.eE\+-]+)\s+alpha\s*=\s*([\d\.eE\+-]+)/);
			$line = <IN>;
			($b,$beta ) 
				= ($line =~ /beta\s*=.*?b\s*=\s*([\d\.eE\+-]+)\s+beta\s*=\s*([\d\.eE\+-]+)/);
			$line = <IN>;
			($c,$gamma) 
				= ($line =~ /gamma\s*=.*?c\s*=\s*([\d\.eE\+-]+)\s+gamma\s*=\s*([\d\.eE\+	-]+)/);
print "Lattice: $a $b $c $alpha $beta $gamma<br>\n" if($Debug);
		}
		elsif($_ =~ /^\s*Label\s+\(Frac\)/i) {
			last;
		}
	}
	<IN>;

	my $Rietan = new Rietan();
	Utils::DelSpace($SPGName);
	$SPGName = "P 1" if(not defined $SPGName or $SPGName eq '');
print "SPGName: [$SPGName]\n";
	my $SPG = $Rietan->ReadSpaceGroupFromSPGName($SPGName);
	if($SPG == 0) {
		print "Error!!: Can not find the space group [$SPGName]\n";
		return -2;
	}
	my $iSPG = $SPG->iSPG();
	my $iSet = $SPG->iSet();
print "iSPG, iSet: $iSPG-$iSet<br>\n" if($Debug);
	my $CIF = new CIF();
#	$CIF->SetFileName($FileName);
	$CIF->SetCrystalName($Title);
	$CIF->SetFormula($Formula);
print "Lattice: $a $b $c $alpha $beta $gamma<br>\n" if($Debug);
	$CIF->SetLatticeParameters($a, $b, $c, $alpha, $beta, $gamma);
#	$CIF->SetSpaceGroup($SPGName, $SPG->iSPG());
	$CIF->SetSpaceGroupInformation($SPG);

	$nAtom = 0;
	my %nAtomName;
	for(my $i = 0 ; $i < $nAtomsShells ; $i++) {
		my $line = <IN>;
		$line =~ /(\d+)[\s\*]+(\S+)[\s\*]+(\S+)[\s\*]+(\S+)[\s\*]+(\S+)[\s\*]+(\S+)[\s\*]+(\S+)[\s\*]+(\S+)[\s\*]+/i;
		my $name = $2;
		my $type = $3;
		my $x = $4;
		my $y = $5;
		my $z = $6;
		my $occ = $8;
		$name =~ s/\d+//g;
		$nAtomName{$name}++;
		my $label = sprintf "%s%d", $name, $nAtomName{$name};
		$occ = 1.0 if($occ eq '');
		next if($type =~ /s/i);
print "$i: $label: $name ($x,$y,$z) [$occ]<br>\n" if($Debug);
		$CIF->AddAsymmetricAtomSite($label, $name, $x, $y, $z, $occ);
	}

	$CIF->FillCIFData();

	close(IN);
	return $CIF;
}

sub GetCCrystalFromInput
{
	my ($this, $fname) = @_;

	my $Crystal = new Crystal();

	my $in = new JFile($fname, "r");
	return 0 unless($in);

	my $line;
	my $Title;
	my $iSPG;
	my ($a, $b, $c, $alpha, $beta, $gamma);
	my %AtomCount;
	while(!$in->eof()) {
		$line = $in->ReadLine();
		last if(!defined $line);

#print "line: $line\n";
		if($line =~ /^\s*title\s/i) {
			$Title = $in->ReadLine();
			Utils::DelSpace($Title);
			$Crystal->SetCrystalName($Title);
			$Crystal->SetSampleName($Title);
		}
		elsif($line =~ /^\s*cell\s/i) {
			$line = $in->ReadLine();
			($a, $b, $c, $alpha, $beta, $gamma) = Utils::Split("\\s+", $line);
#print "Lattice: $a $b $c $alpha $beta $gamma<br>\n";
			$Crystal->SetLatticeParameters($a, $b, $c, $alpha, $beta, $gamma);
		}
		elsif($line =~ /^\s*frac\s/i) {
			while(!$in->eof()) {
				$line = $in->ReadLine();
#print "line: $line\n";
				my ($atom, $pot, $x, $y, $z, $charge, $occ) = Utils::Split("\\s+", $line);
				last if($pot !~ /(core|shel)/i);
				last if(!defined $occ);

				$atom = ucfirst lc $atom;
				$AtomCount{$atom}++;
				$Crystal->AddAtomType($atom);
				my $Label = $atom . $AtomCount{$atom};
				$Crystal->AddAtomSite($Label, $atom, $x, $y, $z, $occ);
			}
		}
		elsif($line =~ /^\s*space\s/i) {
			$iSPG = $in->ReadLine();
			my $iSet = 1;
			my $Rietan = new Rietan();
			my $SPG = $Rietan->ReadSpaceGroup($iSPG, $iSet);
			if($SPG == 0) {
				print "Error!!: Can not find iSPG=$iSPG and iSet=$iSet\n";
				return -2;
			}
#print "iSPG, iSet: $iSPG-$iSet<br>\n";
		}
	}
	$in->Close();

	$Crystal->ExpandCoordinates();

	return $Crystal;
}

#Save GULP md history file
sub SaveMDLatticeHistoryToCSV
{
	my ($this, $infile, $outfile) = @_;

	return -1 unless(open(IN, "<$infile"));
	unless(open(OUT, ">$outfile")) {
		close(IN);
		return -2;
	}

	print OUT "sn,time,a,b,c,alpha,beta,gamma\n";
	my $count = 0;
	while(<IN>) {
		my $line = $_;
		next unless($line =~ /^\s*timestep\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)/);

		$count++;
#print "c: $count ($istep): $line<br>";
		last if($count > 10000);

		my $iTotalStep = $1;
		my $TimeStep   = $5;
		my $TotalTime = $iTotalStep * $TimeStep;
#print "TotalTime: $TotalTime = $iTotalStep * $TimeStep<br>";

		$line = <IN>;
		my ($a11,$a12,$a13) = ($line =~ /(\S+)\s+(\S+)\s+(\S+)/);
		$line = <IN>;
		my ($a21,$a22,$a23) = ($line =~ /(\S+)\s+(\S+)\s+(\S+)/);
		$line = <IN>;
		my ($a31,$a32,$a33) = ($line =~ /(\S+)\s+(\S+)\s+(\S+)/);
		my $Crystal = new Crystal();
		$Crystal->SetLatticeVector($a11, $a12, $a13,
					   $a21, $a22, $a23,
					   $a31, $a32, $a33);
		my ($a,$b,$c,$alpha,$beta,$gamma) = $Crystal->CalculateLatticeConstantFromVector();
#print "Latt: $a $b $c $alpha $beta $gamma<br>";

		print OUT "$count,$TotalTime,$a,$b,$c,$alpha,$beta,$gamma\n";
	}

	close(OUT);
	close(IN);

	return 1;
}

#Read GULP md output file
sub SaveMDPropertyHistoryToCSV
{
	my ($this, $infile, $outfile) = @_;

	return -1 unless(open(IN, "<$infile"));
#出力内容をリストする
	my @OutputKeys;
	while(<IN>) {
		my $line = $_;
		last if($line =~ /^\s*Molecular dynamics equilibration\s*:/i);
	}
	my $timeunit;
	while(<IN>) {
		my $line = $_;
		if($line =~ /Time\s*:\s*(\S+)\s*(\S+)/i) {
			$timeunit = $2;
			last;
		}
	}
	<IN>;
	while(<IN>) {
		my $line = $_;
		last if($line =~ /Time\s*:/i);
		$line =~ /\s*?(\S.*)\s*?=\s*(\S+)\s*(\S+)/;
		my $l = $1;
		$l =~ s/\s+/ /g;
		push(@OutputKeys, $l);
#print "line: $line<br>";
#print "val: $1, $2, $3<br>";
	}

	seek(IN,0,0);
	unless(open(OUT, ">$outfile")) {
		close(IN);
		return -2;
	}
	
	while(<IN>) {
		my $line = $_;
		last if($line =~ /^\s*Molecular dynamics equilibration\s*:/i);
	}

	print OUT "sn,time($timeunit)";
	for(my $i = 0 ; $i < @OutputKeys ; $i++) {
		print OUT ",$OutputKeys[$i](Instantaneous),$OutputKeys[$i](Averaged)";
	}
	print OUT "\n";
	my $count = 0;
	my $timeoffset = 0;
	my $time = 0;
	while(<IN>) {
		if($_ =~ /^\s*Molecular dynamics production\s*:/) {
			$timeoffset = $time;
		}
		next unless($_ =~ /Time\s*:\s*([\d\.\+-eE]+)\s*(\S+)/);
		$count++;
		<IN>;
		$time = $1 + $timeoffset;

		print OUT "$count,$time";
		for(my $i = 0 ; $i < @OutputKeys ; $i++) {
			my $line = <IN>;
			$line =~ /\s*?(\S.*)\s*?=\s*(\S+)\s*(\S+)/;
			print OUT ",$2,$3";
		}
		print OUT "\n";
	}
	
	close(OUT);
	close(IN);
	return 1;
}

sub SaveHistoryFromOutput
{
	my ($this, $Crystal, $OutFile, $HistoryCSVIns, $HistoryCSVAvr) = @_;

print "\n";
print "Save hisotry to [$HistoryCSVIns] and [$HistoryCSVAvr] from [$OutFile].\n";

	my $in   = JFile->new($OutFile, 'r') or return undef;
	my $out1 = JFile->new($HistoryCSVIns, 'w') or return undef;
	my $out2 = JFile->new($HistoryCSVAvr, 'w') or return undef;
	$out1->print("step,time");
	$out2->print("step,time");

	my $line = $in->SkipTo("\\*\\* Time\\s*:");
	my ($timeunit) = ($line =~ /(\S+)\s*:\s*$/);
print "Output properties:\n";	
print "  Time unit: $timeunit\n";
#print "l:$line";
	my @Labels;
	$in->ReadLine();
	my $nProperties = 0;
	my $idxl = -1;
	while(1) {
		$line = $in->ReadLine();
		last if(!defined $line);
		last if($line =~ /\*\*/);
		
		my ($tag, $unit) = ($line =~ /(\S.*)\s+\((.*?)\)\s*=/);
print "  $nProperties: $tag $unit\n";
		Utils::DelSpace($tag);
		Utils::DelSpace($unit);
		$Labels[$nProperties] = "$tag ($unit)";
		$idxl = $nProperties if($idxl < 0 and $tag =~ /Cell\s+parameter/i);
		$out1->print(",$tag ($unit)");
		$out2->print(",$tag ($unit)");
		$nProperties++;
	}
	if($idxl >= 0) {
		$out1->print(",Volume (A3),Density (g/cm3)");
		$out2->print(",Volume (A3),Density (g/cm3)");
	}
	$out1->print("\n");
	$out2->print("\n");
	print "\n";

	my $C1 = $Crystal->Duplicate();
	my $C2 = $Crystal->Duplicate();

	$in->rewind();
	my $istep = 0;
	$line = $in->SkipTo("\\*\\* Time\\s*:");
	my ($time) = ($line =~ /Time\s*:\s*(\S+)/);
	if(!defined $time or $time eq '') {
		$in->Close();
		$out1->Close();
		$out2->Close();
print "\n";
		return ();
	}

	while(1) {
print "  time: $time $timeunit\n" if($istep % 100 == 0);
		$out1->print("$istep,$time");
		$out2->print("$istep,$time");
		$in->ReadLine();
		my (@val1, @val2);
		for(my $i = 0 ; $i < $nProperties ; $i++) {
			$line = $in->ReadLine();
			my ($val1, $val2) = ($line =~ /=\s*(\S+)\s+(\S+)/);
			$val1[$i] = $val1;
			$val2[$i] = $val2;
			$out1->print(",$val1");
			$out2->print(",$val2");
		}
#print "idxl=$idxl\n";
#print "l: $val1[$idxl+0], $val1[$idxl+1], $val1[$idxl+2]\n";
#print "   $val1[$idxl+3], $val1[$idxl+4], $val1[$idxl+5]\n";
		if($idxl >= 0) {
			$C1->SetLatticeParameters(
				$val1[$idxl+0], $val1[$idxl+1], $val1[$idxl+2],
				$val1[$idxl+3], $val1[$idxl+4], $val1[$idxl+5],
				);
			$C2->SetLatticeParameters(
				$val2[$idxl+0], $val2[$idxl+1], $val2[$idxl+2],
				$val2[$idxl+3], $val2[$idxl+4], $val2[$idxl+5],
				);
			$out1->print(",", $C1->Volume(), ",", $C1->CalculateDensity());
			$out2->print(",", $C2->Volume(), ",", $C2->CalculateDensity());
		}
		$out1->print("\n");
		$out2->print("\n");
		$istep++;

		$line = $in->SkipTo("\\*\\* Time\\s*:");
		last if(!defined $line);
		($time) = ($line =~ /Time\s*:\s*(\S+)/);
	}

	$in->Close();
	$out1->Close();
	$out2->Close();
print "\n";

	return (\@Labels);
}

sub GetnStepInHistoryFile
{
	my ($this, $HisFile) = @_;
	
	my $in = JFile->new($HisFile, 'r') or die "Can not read [$HisFile].\n";
	my $n = 0;
	while(1) {
		my $line = $in->SkipTo("timestep");
		last if(!defined $line);
		$n++;
	}
	$in->Close();
	return $n;
}

sub ReadNextCrystalFromHistoryFile
{
	my ($this, $idx, $in, $cry, $IsPrint) = @_;
	$IsPrint = 1 if(!defined $IsPrint);

	my $k = 1.0;

	my @ExpandedAtomSiteList   = $cry->GetCExpandedAtomSiteListByOutputMode();
	my $nSite = @ExpandedAtomSiteList;
#print "nSite=$nSite\n";

	my $Crystal = new Crystal;

	my $line;
	$line = $in->SkipTo("timestep");
	my ($header, $step, $natom, $x, $y, $timestep) = Utils::Split("\\s+", $line);
	$Crystal->{iStep} = $step;
	$Crystal->{Time} = $step * $timestep;

	my ($a11, $a12, $a13) = Utils::Split("\\s+", $in->ReadLine());
	my ($a21, $a22, $a23) = Utils::Split("\\s+", $in->ReadLine());
	my ($a31, $a32, $a33) = Utils::Split("\\s+", $in->ReadLine());
printf "step %6d:  Lattice vector: ($a11, $a12, $a13)  ($a21, $a22, $a23)  "
      ."($a31, $a32, $a33)\n", $idx  if($IsPrint);
	$Crystal->SetLatticeVector(
			$a11, $a12, $a13,
			$a21, $a22, $a23,
			$a31, $a32, $a33);
	my ($a, $b, $c, $alpha, $beta, $gamma) = $Crystal->LatticeParameters();
print "      Lattice parameters: $a, $b, $c, $alpha, $beta, $gamma\n" if($IsPrint);

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

		my ($name, $idx, $mass, $charge) = Utils::Split("\\s+", $line);
		my ($xr, $yr, $zr)  = Utils::Split("\\s+", $in->ReadLine());
		my ($vx, $vy, $vz) = Utils::Split("\\s+", $in->ReadLine());
		next if(!defined $zr);

		my ($x, $y, $z);
#		if($coordmode eq 'frac') {
#			($x, $y, $z) = ($xr, $yr, $zr);
#		}
#		else {
			($x, $y, $z) = $Crystal->CalFractionalCoordinate($xr, $yr, $zr);
#		}

		$AtomCount{$name}++;
		my $Label = $name . $AtomCount{$name};
printf " %06d: %6s: %6s (%8.4g, %8.4g, %8.4g) => (%8.5g, %8.5g, %8.5g)\n", 
	$i, $name, $Label, $xr, $yr, $zr, $x, $y, $z if($IsPrint);
		$Crystal->AddAtomSite($Label, $name, $x, $y, $z, 1.0);
	}
	$Crystal->ExpandCoordinates();

	return ($Crystal, 1);
}

sub MakeXSFFromHistoryFile
{
	my ($this, $HisFile, $XSFFile, $IniCrystal, $nOutputInterval) = @_;

	print("\n\n<b>Make XSF File from GULP history file:</b>\n");
	print("  Read from [$HisFile].\n");
	print("  Save to   [$XSFFile].\n");

#my $nStep = 1462;
	my $nStep = $this->GetnStepInHistoryFile($HisFile);
	my $nOutput = int($nStep / $nOutputInterval);
print "nStep: Total $nStep steps, $nOutput output every $nOutputInterval steps will be saved to [$HisFile] and [$XSFFile].\n";

	my $in = JFile->new($HisFile, 'r') or die "Can not read [$HisFile].\n";

	my $xc = new XCrySDen;
	$xc->WriteXSFAnimationFileHeader($XSFFile, $nOutput);

	my (@pos0, @prevpos, @pos, @offset);
	my ($Crystal);
	my $ret;
	my $iStep = 0;
	my $iOutput = 1;
	my $msd = 0.0;
	for(my $i = 0 ; $i < $nStep ; $i++) {
#$site->SetVelocity($dx, $dy, $dz);
#$site->SetForce(0, 0, $tot);
		if($i % $nOutputInterval == 0) {
			($Crystal, $ret) = $this->ReadNextCrystalFromHistoryFile($iOutput, $in, $IniCrystal);
			last if(!$ret or !$Crystal);
print "\n";
print "iStep=$iStep / $nStep\n";
			$xc->WriteXSFFileFromCrystal($XSFFile, $iStep, $Crystal);
			$iOutput++;
		}
		else {
			$in->SkipTo("timestep");
		}
		$iStep++;
	}
	
	print("\nMake XSF File from history file: Finished\n");

	return $Crystal;
}

sub ReadFinalCrystalStructureFromOutput
{
	my ($this, $OutFile) = @_;

	my $k = 1.0;
#	my $k = $a0*1.0e10;

	my $Crystal = new Crystal();
	
	my $in = new JFile($OutFile, "r");
	return 0 unless($in);

	my $line;
	$line = $in->SkipTo("Formula\\s*=");
	my ($name) = ($line =~ /Formula\s*=\s*(\S.*)\s*$/i);
print "Formula: $name\n";
	$Crystal->SetCrystalName($name);
	$Crystal->SetSampleName($name);

	$line = $in->SkipTo("Fractional coordinates of");
	$in->SkipTo("---");
	$in->SkipTo("---");
	my @occ;
	my $i = 0;
	while(1) {
		$line = $in->ReadLine();
		last if(!defined $line);
		
		$line =~ s/\*/ /g;
		my ($idx, $name, $type, $x, $y, $z, $charge, $occ) = Utils::Split("\\s+", $line);
		next if(!defined $z);
		last if($idx !~ /^\d/);
		$occ[$i] = $occ;
		$i++;
	}

	$in->rewind();
	my ($a, $b, $c, $alpha, $beta, $gamma);
	$line = $in->SkipTo("Final cell parameters and derivatives");
	$line = $in->SkipTo("---");
	my ($s1, $a)     = Utils::Split("\\s+", $in->ReadLine());
	my ($s2, $b)     = Utils::Split("\\s+", $in->ReadLine());
	my ($s3, $c)     = Utils::Split("\\s+", $in->ReadLine());
	my ($s4, $alpha) = Utils::Split("\\s+", $in->ReadLine());
	my ($s5, $beta)  = Utils::Split("\\s+", $in->ReadLine());
	my ($s6, $gamma) = Utils::Split("\\s+", $in->ReadLine());
print "  Lattice parameters: $a, $b, $c, $alpha, $beta, $gamma\n";
	$Crystal->SetLatticeParameters($a, $b, $c, $alpha, $beta, $gamma);

	$in->rewind();
	$line = $in->SkipTo("Final fractional coordinates of");
	$line = $in->SkipTo("---");
	$line = $in->SkipTo("---");

	my $nAtoms = 0;
	my %AtomCount;
	while(1) {
		$line = $in->ReadLine();
#print "l:$line\n";
		last if(!defined $line);

		my ($idx, $name, $type, $x, $y, $z, $r) = Utils::Split("\\s+", $line);
		next if(!defined $z);
		last if($idx !~ /^\d/);

		$AtomCount{$name}++;
		my $Label = $name . $AtomCount{$name};
printf " %06d: %6s: %6s (%8.5g, %8.5g, %8.5g)\n", 
	$nAtoms, $name, $Label, $x, $y, $z;
		$Crystal->AddAtomSite($Label, $name, $x, $y, $z, $occ[$nAtoms]);
		$nAtoms++;
	}
print "  nAtoms = $nAtoms\n";
	$in->Close();

	$Crystal->ExpandCoordinates();
my $V = $Crystal->Volume();
my $d = $Crystal->Density();
print "  Volume  = $V A3\n";
print "  Density = $d g/cm3\n";

	return $Crystal;
}

sub ReadCrystalStructureFromOutput
{
	my ($this, $OutFile) = @_;

	my $k = 1.0;
#	my $k = $a0*1.0e10;

	my $Crystal = new Crystal();
	
	my $in = new JFile($OutFile, "r");
	return 0 unless($in);

	my $line;
	$line = $in->SkipTo("Formula\\s*=");
	my ($name) = ($line =~ /Formula\s*=\s*(\S.*)\s*$/i);
print "Formula: $name\n";
	$Crystal->SetCrystalName($name);
	$Crystal->SetSampleName($name);

	my ($a, $b, $c, $alpha, $beta, $gamma);
	$line = $in->SkipTo("Cell parameters");
	$in->ReadLine();
	$line = $in->ReadLine();
	my ($a, $alpha) = ($line =~ /a\s*=\s*(\S+).*alpha\s*=\s*(\S+)/i);
	$line = $in->ReadLine();
	my ($b, $beta)  = ($line =~ /b\s*=\s*(\S+).*beta\s*=\s*(\S+)/i);
	$line = $in->ReadLine();
	my ($c, $gamma) = ($line =~ /c\s*=\s*(\S+).*gamma\s*=\s*(\S+)/i);
print "  Lattice parameters: $a, $b, $c, $alpha, $beta, $gamma\n";
	$Crystal->SetLatticeParameters($a, $b, $c, $alpha, $beta, $gamma);

	$line = $in->SkipTo("Fractional coordinates of");
	$in->ReadLine();
	$in->ReadLine();
	$in->ReadLine();
	$in->ReadLine();
	$in->ReadLine();

	my $nAtoms = 0;
	my %AtomCount;
	while(1) {
		$line = $in->ReadLine();
#print "l:$line\n";
		last if(!defined $line);

		$line =~ s/\*/ /g;
		my ($idx, $name, $type, $x, $y, $z, $charge, $occ) = Utils::Split("\\s+", $line);
		next if(!defined $z);
		last if($idx !~ /^\d/);

		$AtomCount{$name}++;
		my $Label = $name . $AtomCount{$name};
printf " %06d: %6s: %6s (%8.5g, %8.5g, %8.5g)\n", 
	$nAtoms, $name, $Label, $x, $y, $z;
		$Crystal->AddAtomSite($Label, $name, $x, $y, $z, $occ);
		$nAtoms++;
	}
print "  nAtoms = $nAtoms\n";
	$in->Close();

	$Crystal->ExpandCoordinates();
my $V = $Crystal->Volume();
my $d = $Crystal->Density();
print "  Volume  = $V A3\n";
print "  Density = $d g/cm3\n";

	return $Crystal;
}

sub ReadCrystalStructureFromInput
{
	my ($this, $InFile) = @_;

	my $k = 1.0;
#	my $k = $a0*1.0e10;

	my $Crystal = new Crystal();
	
	my $in = new JFile($InFile, "r");
	return 0 unless($in);

	my $line = $in->SkipTo("name ");
	my ($name) = ($line =~ /name\s+(\S.*)\s*$/i);
	$Crystal->SetCrystalName($name);
	$Crystal->SetSampleName($name);

	$in->rewind();
	my ($a, $b, $c, $alpha, $beta, $gamma);
	$line = $in->SkipTo("cell\\s*\$");
#print "cell=$line\n";
	if(defined $line) {
		$line = $in->ReadLine();
		($a, $b, $c, $alpha, $beta, $gamma) = Utils::Split("\\s+", $line);
		$Crystal->SetLatticeParameters($a, $b, $c, $alpha, $beta, $gamma);
	}
	else {
		$line = $in->SkipTo("vectors\\s*\$");
		my ($a11, $a12, $a13) = Utils::Split("\\s+", $in->ReadLine());
		my ($a21, $a22, $a23) = Utils::Split("\\s+", $in->ReadLine());
		my ($a31, $a32, $a33) = Utils::Split("\\s+", $in->ReadLine());
print "  Lattice vector: ($a11, $a12, $a13)  ($a21, $a22, $a23)  ($a31, $a32, $a33)\n";
		$Crystal->SetLatticeVector(
				$a11, $a12, $a13,
				$a21, $a22, $a23,
				$a31, $a32, $a33);
		($a, $b, $c, $alpha, $beta, $gamma) = $Crystal->LatticeParameters();
	}
print "  Lattice parameters: $a, $b, $c, $alpha, $beta, $gamma\n";

	my $coordmode = '';
	while(1) {
		$line = $in->ReadLine();
		last if(!defined $line);
		if($line =~ /frac/i) {
			$coordmode = 'frac';
			last;
		}
		if($line =~ /cart/i) {
			$coordmode = 'cart';
			last;
		}
	}
print "  Coordinate mode: $coordmode\n";
#exit;
	Utils::DelSpace($coordmode);
	my $nAtoms = 0;
	my %AtomCount;
	while(1) {
		$line = $in->ReadLine();
#print "l:$line\n";
		last if(!defined $line);

		my ($name, $type, $xr, $yr, $zr, $a, $occ, $c, $idx, $idy, $idz)
			= Utils::Split("\\s+", $line);
		next if(!defined $zr);
		last if($type !~ /core/i and $type !~ /shel/i);
		last if($name =~ /space/i or $name =~ /^\d/);

		my ($x, $y, $z);
		if($coordmode eq 'frac') {
			($x, $y, $z) = ($xr, $yr, $zr);
		}
		else {
			($x, $y, $z) = $Crystal->CalFractionalCoordinate($xr, $yr, $zr);
		}

		$AtomCount{$name}++;
		my $Label = $name . $AtomCount{$name};
printf " %06d: %6s: %6s (%8.4g, %8.4g, %8.4g) => (%8.5g, %8.5g, %8.5g)\n", 
	$nAtoms, $name, $Label, $xr, $yr, $zr, $x, $y, $z;
		$Crystal->AddAtomSite($Label, $name, $x, $y, $z, $occ);
		$nAtoms++;
	}
print "  nAtoms = $nAtoms\n";
	$in->Close();

	$Crystal->ExpandCoordinates();
my $V = $Crystal->Volume();
my $d = $Crystal->Density();
print "  Volume  = $V A3\n";
print "  Density = $d g/cm3\n";

	return $Crystal;
}

sub ReadInputToHash
{
	my ($this, $InFile) = @_;

	my %hash;
	my $in = JFile->new($InFile, 'r') or return undef;
	while(1) {
		my $line = $in->ReadLine();
		last if(!defined $line);
		
		$line =~ s/#.*$//;
		my ($key, $dummy, $val) = ($line =~ /^\s*(\S+)(\s+(\S.*)\s*)?$/);
		next if(!defined $key or $key =~ /^[A-Z][a-z]?[_\-\d]*$/);
		next if(lc $key eq 'end');
		
		$val = 1 if($key =~ /(fit|optimize|md)/ and !defined $val);
		if($key =~ /[a-zA-Z]+/) {
			if($val eq '') {
				$val = $in->ReadLine();
				Utils::DelSpace($val);
			}
			if($key =~ /^opti/i) {
				$key = 'optimize';
			}

			$hash{$key} = $val;
#print "hash: $key: $hash{$key}\n";
		}
	}	
	$in->Close();

	if(defined $hash{fit}) {
		$hash{Task} = 'fit';
	}
	elsif(defined $hash{optimize}) {
		$hash{Task} = 'optimize';
	}
	elsif(defined $hash{md} and $hash{md} =~ /conp/) {
		$hash{Task} = 'MD-NPT';
	}
	elsif(defined $hash{md}) {
		$hash{Task} = 'MD-NVT';
	}

	return \%hash;
}

sub ReadOutputToHash
{
	my ($this, $OutFile) = @_;

	my %hash;
	my $in = JFile->new($OutFile, 'r') or return undef;
	my $line = $in->SkipTo("Version\\s*=");
#print "l$line\n";
	($hash{version}) = ($line =~ /Version\s*=\s*(\S+)/);
#print "Version: $hash{version}\n";
	$in->ReadLine();
	while(1) {
		$line = $in->ReadLine();
#print "l:$line";
		last if(!defined $line);
		last if($line =~ /^\*\*/);

		$line =~ s/^\*\s*//;
		my ($key) = Utils::Split("\\s+", $line);
		$key = 'optimize' if($key =~ /^opti/i);
		$hash{$key} = 1;
#print "key: $key\n";
	}
	$in->ReadLine();	
	$in->ReadLine();	
	for(my $i = 0 ; $i < 100 ; $i++) {
		$line = $in->ReadLine();
#print "l2:$line";
		last if(!defined $line);
		next if($line =~ /^\*/);
		my ($key, $val) = ($line =~ /^\s*(\S.*)\s*=\s*(\S.*)\s*$/);
		$hash{$key} = $val if(defined $val); 
	}

	$in->Close();

	if(defined $hash{fit}) {
		$hash{Task} = 'fit';
	}
	elsif(defined $hash{optimize}) {
		$hash{Task} = 'optimize';
	}
	elsif(defined $hash{md} and $hash{md} =~ /conp/) {
		$hash{Task} = 'MD-NPT';
	}
	elsif(defined $hash{md}) {
		$hash{Task} = 'MD-NVT';
	}

	return \%hash;
}

sub WriteSpecies
{
	my ($this, $out, $pPot, $pAtomTypeList) = @_;

	$out->print("species ", scalar @$pAtomTypeList, "\n");
	for(my $i = 0 ; $i < @$pAtomTypeList ; $i++) {
		my $type1 = $pAtomTypeList->[$i];
		my $name1 = $type1->AtomNameOnly();
		printf       "%3s %4s %7.4f\n", $name1, 'core', $pPot->{"$name1-core"}{charge};
		$out->printf("%3s %4s %7.4f\n", $name1, 'core', $pPot->{"$name1-core"}{charge});
	}

	return 1;
}

sub WritePotential
{
	my ($this, $out, $pPot, $pAtomTypeList) = @_;

	if($pPot->{type} eq 'buckingham') {
		$out->print("buckingham\n");
		for(my $i = 0 ; $i < @$pAtomTypeList ; $i++) {
			my $type1 = $pAtomTypeList->[$i];
			my $name1 = $type1->AtomNameOnly();
			for(my $j = $i ; $j < @$pAtomTypeList ; $j++) {
				my $type2 = $pAtomTypeList->[$j];
				my $name2 = $type2->AtomNameOnly();
#print "i,j: $i,$j: $name1 $name2$LF";
				my $p = $pPot->{"$name1-core-$name2-core"};
				my $Aij = $p->{A};
				my $bij = $p->{rho};
				my $Cij = $p->{C};
				$Aij = 0.0 if(!defined $Aij);
				$bij = 1.0 if(!defined $bij);
				$Cij = 0.0 if(!defined $Cij);
				my ($idA, $idb, $idC) = (1, 0, 0);
				$idA = 0 if($Aij == 0.0);
#print "  $name1 - $name2: $Aij  $bij  $Cij\n";
				printf       "%3s %4s %3s %4s %11.6g %11.6g %11.6g %3.1f %4.1f  %d %d %d\n",
					$name1, 'core', $name2, 'core', $Aij, $bij, $Cij, 0.0, 10.0, $idA, $idb, $idC;
				$out->printf("%3s %4s %3s %4s %11.6g %11.6g %11.6g %3.1f %4.1f  %d %d %d\n",
					$name1, 'core', $name2, 'core', $Aij, $bij, $Cij, 0.0, 10.0, $idA, $idb, $idC);
			}
		}
	}
	elsif($pPot->{type} eq 'morse') {
		$out->print("morse\n");
		for(my $i = 0 ; $i < @$pAtomTypeList ; $i++) {
			my $type1 = $pAtomTypeList->[$i];
			my $name1 = $type1->AtomNameOnly();
			for(my $j = $i ; $j < @$pAtomTypeList ; $j++) {
				my $type2 = $pAtomTypeList->[$j];
				my $name2 = $type2->AtomNameOnly();
#print "i,j: $i,$j: $name1 $name2$LF";
				my $p = $pPot->{"$name1-core-$name2-core"};
				my $De        = $p->{De};
				my $a0        = $p->{a0};
				my $r0        = $p->{r0};
				my $fsCoulomb = $p->{fsCoulomb};
#print "  $name1 - $name2: $De  $a0  $r0  $fsCoulomb\n";
				my ($idDe, $ida0, $idr0) = (1, 0, 0);
				$idDe = 0 if($De == 0.0);
				printf       "%3s %4s %3s %4s %11.6g %11.6g %11.6g %6.3f %3.1f %4.1f  %d %d %d\n",
					$name1, 'core', $name2, 'core', $De, $a0, $r0, $fsCoulomb, 0.0, 10.0, $idDe, $ida0, $idr0;
				$out->printf("%3s %4s %3s %4s %11.6g %11.6g %11.6g %6.3f %3.1f %4.1f  %d %d %d\n",
					$name1, 'core', $name2, 'core', $De, $a0, $r0, $fsCoulomb, 0.0, 10.0, $idDe, $ida0, $idr0);
			}
		}
	}
#	if($nShellAtoms > 0) {
#		print OUT "spring\n";
#		for(my $i = 0 ; $i < @Spring ; $i++) {
#			print OUT "$Spring[$i]\n";
#		}
#	}

	return 1;
}

sub ReadPotentialFromInput
{
	my ($this, $InFile, $pPot, $SpeciesOnly) = @_;
	$SpeciesOnly = 0 if(!defined $SpeciesOnly);

	print "Potential parameters from [$InFile]:\n";

	my %hash = ($pPot)? %$pPot : ();

	my $in = JFile->new($InFile, 'r') or return undef;
	my $line = $in->SkipTo("species\\s+(\\d+)");
#print "line: $line";
	my ($nAtom) = ($line =~ /species\s+(\d+)/);
print "  nAtom in species: $nAtom\n";
	for(my $i = 0 ; $i < $nAtom ; $i++) {
		$line = $in->ReadLine();
		last if(!defined $line);
		
		my ($name, $type, $charge) = Utils::Split("\\s+", $line);
print "  $name $type: $charge\n";
		$hash{"$name-$type"}{charge} = $charge;
	}
	if($SpeciesOnly) {
		$in->Close();
		return \%hash;
	}

	$in->rewind();
	$line = $in->SkipTo("^\\s*buck");
#print "l:$line";
	if($line) {
		$hash{buckingham} = 1;
		$hash{type}       = 'buckingham';
		for(my $i = 0 ; $i < $nAtom ; $i++) {
			for(my $j = $i ; $j < $nAtom ; $j++) {
				$line = $in->ReadLine();
				my ($name1, $type1, $name2, $type2, $A, $rho, $C, 
					$rmin, $rmax, $idA, $idrho, $idC)
						= Utils::Split("\\s+", $line);

print "  Buchkingham for $name1 $type1 - $name2 $type2: $A, $rho, $C, $rmin, $rmax, $idA, $idrho, $idC\n";
				$hash{"$name1-$type1-$name2-$type2"} 
					= $hash{"$name2-$type2-$name1-$type1"} 
						= {
						A     => $A,
						rho   => $rho,
						C     => $C,
						rmin  => $rmin,
						rmax  => $rmax,
						idA   => $idA,
						idrho => $idrho, 
						idC   => $idC
						
						};
			}
		}
		return \%hash;
	}

	$in->rewind();
	$line = $in->SkipTo("^\\s*mors");
#print "l:$line";
	if($line) {
		$hash{morse} = 1;
		$hash{type}  = 'morse';
		for(my $i = 0 ; $i < $nAtom ; $i++) {
			for(my $j = $i ; $j < $nAtom ; $j++) {
				$line = $in->ReadLine();
				my ($name1, $type1, $name2, $type2, $De, $a0, $r0, $fsCoulomb, 
					$rmin, $rmax, $idDe, $ida0, $idr0)
						= Utils::Split("\\s+", $line);

print "  Morse for $name1 $type1 - $name2 $type2: $De, $a0, $r0, $fsCoulomb, $rmin, $rmax, $idDe, $ida0, $idr0\n";
				$hash{"$name1-$type1-$name2-$type2"} 
					= $hash{"$name2-$type2-$name1-$type1"} 
						= {
						De        => $De,
						a0        => $a0,
						r0        => $r0,
						fsCoulomb => $fsCoulomb,
						rmin      => $rmin,
						rmax      => $rmax,
						idDe      => $idDe,
						ida0      => $ida0,
						idr0      => $idr0,
						};
			}
		}
		return \%hash;
	}

	return undef;
}

sub ReadNewPotentialFromOutput
{
	my ($this, $OutFile, $pPot) = @_;

	print "New potential parameters from [$OutFile]:\n";

	my $in = JFile->new($OutFile, 'r') or return undef;
	my $line = $in->SkipTo("Observable\\s+no");
	   $line = $in->SkipTo("Comparison.+initial.*final");
	return undef if(!defined $line);
	$line = $in->SkipTo("General interatomic potentials");
	$line = $in->SkipTo("Atom\\s+Types");
	$line = $in->SkipTo("---");

	$this->ClearPotentialLines();

	my %hash = ($pPot)? %$pPot : ();
	my $IsEnd = 0;
	while(1) {
		while(1) {
			$line = $in->ReadLine();
			last if($line =~ /---/);

			$line =~ s/(\d)eV/$1 eV/;
			my ($name1, $type1, $name2, $type2, $pottype, $pottype1, 
				$para, $val, $unit, $rmin, $rmax) 
				= Utils::Split("\\s+", $line);
			if($type1 ne 'c' and $type1 ne 's') {
				$IsEnd = 1;
				last;
			}
			$type1 = 'core'  if($type1 eq 'c');
			$type1 = 'shell' if($type1 eq 's');
			$type2 = 'core'  if($type2 eq 'c');
			$type2 = 'shell' if($type2 eq 's');

#			my $key = "$name1-$type1-$name2-$type2";
			if($pottype =~ /buck/i) {
				$this->AddPotentialLines($line);
				my $A = $val;
print "  A for $name1-$name2: $A\n";

				$line = $in->ReadLine();
				$this->AddPotentialLines($line);
				($pottype1, $para, $val, $unit, $rmin, $rmax) = Utils::Split("\\s+", $line);
				my $rho = $val;
print "  rho for $name1-$name2: $rho\n";

				$line = $in->ReadLine();
				$this->AddPotentialLines($line);
				($pottype1, $para, $val, $unit, $rmin, $rmax) = Utils::Split("\\s+", $line);
				my $C = $val;
print "  C for for $name1-$name2: $C\n";

				$hash{"$name1-$type1-$name2-$type2"} 
					= $hash{"$name2-$type2-$name1-$type1"} 
						= {
						A    => $A,
						rho  => $rho,
						C    => $C,
						rmin => $rmin,
						rmax => $rmax,
						};
			}
			elsif($pottype =~ /morse/i) {
				$this->AddPotentialLines($line);
				my $De = $val;
print "  De for $name1-$name2: $De\n";

				$line = $in->ReadLine();
				$this->AddPotentialLines($line);
				($pottype1, $para, $val, $unit, $rmin, $rmax) = Utils::Split("\\s+", $line);
				my $a0 = $val;
print "  a0 for $name1-$name2: $a0\n";

				$line = $in->ReadLine();
				$this->AddPotentialLines($line);
				($pottype1, $para, $val, $unit, $rmin, $rmax) = Utils::Split("\\s+", $line);
				my $r0 = $val;
print "  r0 for for $name1-$name2: $r0\n";

				$hash{"$name1-$type1-$name2-$type2"} 
					= $hash{"$name2-$type2-$name1-$type1"} 
						= {
						De   => $De,
						a0   => $a0,
						r0   => $r0,
						rmin => $rmin,
						rmax => $rmax,
						};
			}
			
		}
		last if($IsEnd);
	}
	$in->Close();

	return \%hash;
}

sub UpdateInputForFit
{
	my ($this, $InpFile, $OutFile, $NewFile, %args) = @_;

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

	print "\n";
	my $pPot = $this->ReadPotentialFromInput($InpFile);
	print "\n";
	   $pPot = $this->ReadNewPotentialFromOutput($OutFile, $pPot);
	
	print "\n";
	print "Update [$InpFile] to [$NewFile].\n";
	my $in = JFile->new($InpFile, 'r') or return undef;
	my $out = JFile->new($NewFile, 'w') or return undef;
	my $line;
	while(1) {
		$line = $in->ReadLine();
		last if(!defined $line);
		last if($line =~ /spec/i);

		$out->print($line);
	}

	for(my $i = 0 ; $i < @AtomTypeList ; $i++) {
		$in->ReadLine();
	}
	$this->WriteSpecies($out, $pPot, \@AtomTypeList);

	while(1) {
		$line = $in->ReadLine();
		last if(!defined $line);
		last if($line =~ /buck/i or $line =~ /morse/i);

		$out->print($line);
	}

	for(my $i = 0 ; $i < @AtomTypeList ; $i++) {
		for(my $j = $i ; $j < @AtomTypeList ; $j++) {
			$in->ReadLine();
		}
	}
	$this->WritePotential($out, $pPot, \@AtomTypeList);

	while(1) {
		$line = $in->ReadLine();
		last if(!defined $line);
		$out->print($line);
	}

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

	my ($rmin, $rstep, $nr) = (0.01, 0.01, 400);
	$this->SavePotentials('Potentials.csv', $rmin, $rstep, $nr, \@AtomTypeList, $pPot);

	return 1;
}

sub UpdateInputForOptimize
{
	my ($this, $InpFile, $OutFile, $NewFile, %args) = @_;

	my $PrevFile = "$InpFile.prev";

	my $Function    = $args{Function};
	my $production  = $args{production};
	my $timestep    = $args{timestep};
	my $nstep       = (defined $timestep)? int($production / $timestep) + 1 : 0;
	my $tscale      = $args{tscale};
	my $temperature = $args{temperature};
	my $pressure    = $args{pressure};
	my $sample      = $args{sample};
	my $write       = $args{write};

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

	print("  Update $InpFile, backup to $PrevFile.\n");
	print("  New input file: $NewFile\n");
	print("  Function      : $Function\n");
#	print("  equilibration : $equilibration\n");
	print("  production    : $production\n");
	print("  timestep      : $timestep\n");
	print("    production interval at $timestep ps with $nstep steps\n");
	print("  tscale     : $tscale\n");
	print("  temperature: $temperature\n");
	print("  pressure   : $pressure\n");
	print("  sample     : $sample\n");
	print("  write      : $write\n");

	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;
	}

	print "Read final crystal structure from [$OutFile].\n";
	my $Crystal = $this->ReadFinalCrystalStructureFromOutput($OutFile);

	my $SourceFile = $InpFile;

	my $pPotHash = {};
	$this->ModifyInput($SourceFile, $NewFile, $Function, $Crystal, $pPotHash, 
		%args,
		nstep => $nstep,
		T0    => $T0,
		T1    => $T1,
		Tstep => $Tstep,
		);
}

sub UpdateInputForMD
{
	my ($this, $InpFile, $OutFile, $ResFile, $NewFile, %args) = @_;

	my $PrevFile = "$InpFile.prev";
	if(-f $ResFile) {
		print("  Read restart file [$ResFile].\n");
	}
	else {
		print("  Restart file [$ResFile] is not found.\n");
	}

	my $Function    = $args{Function};
	my $production  = $args{production};
	my $timestep    = $args{timestep};
	my $nstep       = (defined $timestep)? int($production / $timestep) + 1 : 0;
	my $tscale      = $args{tscale};
	my $temperature = $args{temperature};
	my $pressure    = $args{pressure};
	my $sample      = $args{sample};
	my $write       = $args{write};

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

	print("  Update $InpFile, backup to $PrevFile.\n");
	print("  New input file: $NewFile\n");
	print("  Function      : $Function\n");
#	print("  equilibration : $equilibration\n");
	print("  production    : $production\n");
	print("  timestep      : $timestep\n");
	print("    production interval at $timestep ps with $nstep steps\n");
	print("  tscale     : $tscale\n");
	print("  temperature: $temperature\n");
	print("  pressure   : $pressure\n");
	print("  sample     : $sample\n");
	print("  write      : $write\n");

	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;
	}

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

	my $SourceFile = $InpFile;
	if(-f $ResFile) {
		$SourceFile = $ResFile;
#		print("Copy [$ResFile] to [$InpFile].\n");
#		if(!Utils::CopyFile($ResFile, $InpFile)) {
#			print("  Error in GULP.pl::UpdateInputFile: Can not copy [$InpFile] to [$ResFile].\n");
#			exit;
#		}
	}

	my $pPotHash = {};
	$this->ModifyInput($SourceFile, $NewFile, $Function, $Crystal, $pPotHash, 
		%args,
		nstep => $nstep,
		T0    => $T0,
		T1    => $T1,
		Tstep => $Tstep,
		);

#	print("Delete [$InpFile].\n");
#	if(!Utils::DeleteFile($InpFile)) {
#		print("  Error in GULP.pl::UpdateInputFile: Can not delete [$InpFile].\n");
#		exit;
#	}

#	print("Copy [$NewFile] to [$InpFile].\n");
#	if(!Utils::CopyFile($NewFile, $InpFile)) {
#		print("  Error in GULP.pl::UpdateInputFile: Can not copy [$InpFile] to [$PrevFile].\n");
#		exit;
#	}
}

sub ModifyInput
{
	my ($this, $SourceFile, $NewFile, $Function, $Crystal, $pPotHash, %args) = @_;
	$args{ResetTime} = 0 if(!defined $args{ResetTime});

	print "\n";
	print "Modify input file [$SourceFile] to new file [$NewFile].\n";
	
	my $in  = JFile->new($SourceFile, 'r') or die "Can not read [$SourceFile].\n";
	my $out = JFile->new($NewFile, 'w') or die "Can not write to [$NewFile].\n";
	while(1) {
		my $line = $in->ReadLine();
		last if(!defined $line);
		
		$line =~ s/[\r\n]+$//;

		if($Function =~ /Optimize/i) {
			if($line =~ /^(\s*)#(\s*opti.*)$/) {
				$out->print("$1$2\n");
				next;
			}
			elsif($line =~ /^(\s*)(#?)(\s*(md|fit|phonon).*)$/) {
				$out->print("$1#$3\n");
				next;
			}
			elsif($line =~ /^(\s*)(#?)(\s*(ftol|gtol|xtol|maxcyc).*)$/) {
				$out->print("$1$3\n");
				next;
			}
			elsif($line =~ /^(\s*)(#?)(\s*omega.*)$/) {
				$out->print("$1#$3\n");
				next;
			}
			elsif($line =~ /^(\s*)(#?)(\s*dispersion.*)$/) {
				$out->print("$1#$3\n");
				while(1) {
					$line = $in->ReadLine();
					last if(!defined $line);
					Utils::DelSpace($line);
					if($line =~ /^(\s*)(#?)(.*\sto\s.*)$/) {
						$out->print("$1#$3\n");
						next;
					}
					else {
						$out->print("$line\n");
						last;
					}
				}
				next;
			}
			elsif($line =~ /^(\s*)(#?)(\s*(ensemble).*)$/) {
				$out->print("$1#$3\n");
				next;
			}
		}
		elsif($Function =~ /Fit/i) {
			if($line =~ /^(\s*)#(\s*fit.*)$/) {
				$out->print("$1$2\n");
				next;
			}
			elsif($line =~ /^(\s*)(#?)(\s*(md|opti|phonon).*)$/) {
				$out->print("$1#$3\n");
				next;
			}
			elsif($line =~ /^(\s*)(#?)(\s*(ftol|gtol|xtol|maxcyc).*)$/) {
				$out->print("$1#$3\n");
				next;
			}
			elsif($line =~ /^(\s*)(#?)(\s*omega.*)$/) {
				$out->print("$1#$3\n");
				next;
			}
			elsif($line =~ /^(\s*)(#?)(\s*dispersion.*)$/) {
				$out->print("$1#$3\n");
				while(1) {
					$line = $in->ReadLine();
					last if(!defined $line);
					Utils::DelSpace($line);
					if($line =~ /^(\s*)(#?)(.*\sto\s.*)$/) {
						$out->print("$1#$3\n");
						next;
					}
					else {
						$out->print("$line\n");
						last;
					}
				}
				next;
			}
			elsif($line =~ /^(\s*)(#?)(\s*ensemble.*)$/) {
				$out->print("$1#$3\n");
				next;
			}
		}
		elsif($Function =~ /MD-NVT/i) {
			if($line =~ /^(\s*)(#?)(\s*md.+conv.*)$/) {
				$out->print("$1$3\n");
				next;
			}
			elsif($line =~ /^(\s*)(#?)(\s*(md|fit|opti|phonon).*)$/) {
				$out->print("$1#$3\n");
				next;
			}
			elsif($line =~ /^(\s*)(#?)(\s*(ftol|gtol|xtol|maxcyc).*)$/) {
				$out->print("$1#$3\n");
				next;
			}
			elsif($line =~ /^(\s*)(#?)(\s*omega.*)$/) {
				$out->print("$1#$3\n");
				next;
			}
			elsif($line =~ /^(\s*)(#?)(\s*dispersion.*)$/) {
				$out->print("$1#$3\n");
				while(1) {
					$line = $in->ReadLine();
					last if(!defined $line);
					Utils::DelSpace($line);
					if($line =~ /^(\s*)(#?)(.*\sto\s.*)$/) {
						$out->print("$1#$3\n");
						next;
					}
					else {
						$out->print("$line\n");
						last;
					}
				}
				next;
			}
			elsif($line =~ /^(\s*)(#?)(\s*ensemble\s+nvt.*)$/) {
				$out->print("$1$3\n");
				next;
			}
			elsif($line =~ /^(\s*)(#?)(\s*ensemble\s+npt.*)$/) {
				$out->print("$1#$3\n");
				next;
			}
		}
		elsif($Function =~ /MD-NPT/i) {
			if($line =~ /^(\s*)(#?)(\s*md.+conp.*)$/) {
				$out->print("$1$3\n");
				next;
			}
			elsif($line =~ /^(\s*)(#?)(\s*(md|fit|opti|phonon).*)$/) {
				$out->print("$1#$3\n");
				next;
			}
			elsif($line =~ /^(\s*)(#?)(\s*(ftol|gtol|xtol|maxcyc).*)$/) {
				$out->print("$1#$3\n");
				next;
			}
			elsif($line =~ /^(\s*)(#?)(\s*omega.*)$/) {
				$out->print("$1#$3\n");
				next;
			}
			elsif($line =~ /^(\s*)(#?)(\s*dispersion.*)$/) {
				$out->print("$1#$3\n");
				while(1) {
					$line = $in->ReadLine();
					last if(!defined $line);
					Utils::DelSpace($line);
					if($line =~ /^(\s*)(#?)(.*\sto\s.*)$/) {
						$out->print("$1#$3\n");
						next;
					}
					else {
						$out->print("$line\n");
						last;
					}
				}
				next;
			}
			elsif($line =~ /^(\s*)(#?)(\s*ensemble\s+npt.*)$/) {
				$out->print("$1$3\n");
				next;
			}
			elsif($line =~ /^(\s*)(#?)(\s*ensemble\s+nvt.*)$/) {
				$out->print("$1#$3\n");
				next;
			}
		}
		elsif($Function =~ /phonon/i) {
			if($line =~ /^(\s*)#(\s*phon.*)$/) {
				$out->print("$1$2\n");
				next;
			}
			elsif($line =~ /^(\s*)(#?)(\s*(md|fit|opti).*)$/) {
				$out->print("$1#$3\n");
				next;
			}
			elsif($line =~ /^(\s*)(#?)(\s*(ftol|gtol|xtol|maxcyc).*)$/) {
				$out->print("$1#$3\n");
				next;
			}
			elsif($line =~ /^(\s*)(#?)(\s*omega.*)$/) {
				$out->print("$1$3\n");
				next;
			}
			elsif($line =~ /^(\s*)(#?)(\s*dispersion.*)$/) {
				$out->print("$1$3\n");
				while(1) {
					$line = $in->ReadLine();
					last if(!defined $line);
					Utils::DelSpace($line);
					if($line =~ /^(\s*)(#?)(.*\sto\s.*)$/) {
						$out->print("$1$3\n");
						next;
					}
					else {
						$out->print("$line\n");
						last;
					}
				}
				next;
			}
			elsif($line =~ /^(\s*)(#?)(\s*ensemble\s+npt.*)$/) {
				$out->print("$1$3\n");
				next;
			}
			elsif($line =~ /^(\s*)(#?)(\s*ensemble\s+nvt.*)$/) {
				$out->print("$1#$3\n");
				next;
			}
		}

		if($Function =~ /MD/i) {
			if($line =~ /^(\s*)(#?)(temp\S*\s+)(\d.*)$/i) {
#print "/$1/$2/$3/$4 [$args{T0}]\n";
#exit;
				my $head = "$1$3";
				my $val  = $4;
				if(defined $args{T0} and $args{T0} =~ /^[\d+\-\.eE]+$/) {
					$out->print("$head$args{T0} $args{Tstep} $args{nstep}\n");
				}
				else {
					$out->print("$head$val\n");
				}
				next;
			}
			elsif($line =~ /^(\s*)(#?)(equi\S*\s+)(\S*)(\s*)(\S+\s*)?$/i) {
				if(defined $args{equilibration}) {
					$out->print("$1$3$args{equilibration}$5$6\n");
				}
				else {
					$out->print("$1$3$4$5$6\n");
				}
				next;
			}
			elsif($line =~ /^(\s*)(#?)(prod\S*\s+)(\S*)(\s*)(\S+\s*)?$/i) {
				if(defined $args{production}) {
					$out->print("$1$3$args{production}$5$6\n");
				}
				else {
					$out->print("$1$3$4$5$6\n");
				}
				next;
			}
			elsif($line =~ /^(\s*)(#?)(timestep\S*\s+)(\S*)(\s*)(\S+\s*)?$/i) {
				if(defined $args{timestep}) {
#print "/$1/$2/$3/$4/$5/$6 [$args{timestep}]\n";
#exit;
					$out->print("$1$3$args{timestep}$5$6\n");
				}
				else {
					$out->print("$1$3$4$5$6\n");
				}
				next;
			}
			elsif($line =~ /^(\s*)(#?)(tsca\S*\s+)(\S*)(\s*)(\S+\s*)?$/i) {
				if(defined $args{tscale}) {
					$out->print("$1$3$args{tscale}$5$6\n");
				}
				else {
					$out->print("$1$3$4$5$6\n");
				}
				next;
			}
			elsif($line =~ /^(\s*)(#?)(tau_t\S*\s+)(\S*)(\s*)(\S+\s*)?$/i) {
				if(defined $args{tau_thermostat}) {
					$out->print("$1$3$args{tau_thermostat}$5$6\n");
				}
				else {
					$out->print("$1$3$4$5$6\n");
				}
				next;
			}
			elsif($line =~ /^(\s*)(#?)(tau_b\S*\s+)(\S*)(\s*)(\S+\s*)?$/i) {
				if(defined $args{tau_barostat}) {
					$out->print("$1$3$args{tau_barostat}$5$6\n");
				}
				else {
					$out->print("$1$3$4$5$6\n");
				}
				next;
			}
			elsif($line =~ /^(\s*)(#?)(samp\S*\s+)(\S*)(\s*)(\S+\s*)?$/i) {
				if(defined $args{sample}) {
					$out->print("$1$3$args{sample}$5$6\n");
				}
				else {
					$out->print("$1$3$4$5$6\n");
				}
				next;
			}
			elsif($line =~ /^(\s*)(#?)(write\S*\s+)(\S*)(\s*)(\S+\s*)?$/i) {
				if(defined $args{write}) {
					$out->print("$1$3$args{write}$5$6\n");
				}
				else {
					$out->print("$1$3$4$5$6\n");
				}
				next;
			}
			elsif($line =~ /^(\s*)(#?)(dump\S*\s.*)$/i) {
				$out->print("$1$3\n");
				next;
			}
			elsif($line =~ /^(\s*)(#?)(outp\S*\s.*)$/i) {
				$out->print("$1$3\n");
				next;
			}
			elsif($line =~ /^current_time.*(\S+)\s*?$/) {
				if($args{ResetTime}) {
					$out->print("current_time     0.000000 $1\n");
					$in->ReadLine();
					$in->ReadLine();
					$in->ReadLine();
					$in->ReadLine();
				}
				else {
					$out->print("$line\n");
				}
				next;
			}
		}
		else {
			if($line =~ /^(\s*)(#?)(pres\S*\s+)([\-+\d\.]\S+)/i) {
				if(defined $args{pressure}) {
					$out->print("$1$3$args{pressure}\n");
				}
				else {
					$out->print("$1$3$4$5$6\n");
				}
				next;
			}
			elsif($line =~ /^(\s*)(#?)(equi.*)$/i) {
				$out->print("$1#$3\n");
				next;
			}
			elsif($line =~ /^(\s*)(#?)(prod.*)(\S*)(\s*)(\S+\s*)?$/i) {
				$out->print("$1#$3$4$5$6\n");
				next;
			}
			elsif($line =~ /^(\s*)(#?)(timestep.*)(\S*)(\s*)(\S+\s*)?$/i) {
				$out->print("$1#$3$4$5$6\n");
				next;
			}
			elsif($line =~ /^(\s*)(#?)(tsca.*)(\S*)(\s*)(\S+\s*)?$/i) {
				$out->print("$1#$3$4$5$6\n");
				next;
			}
			elsif($line =~ /^(\s*)(#?)(tau_.*)(\S*)(\s*)(\S+\s*)?$/i) {
				$out->print("$1#$3$4$5$6\n");
				next;
			}
			elsif($line =~ /^(\s*)(#?)(samp.*)(\S*)(\s*)(\S+\s*)?$/i) {
				$out->print("$1#$3$4$5$6\n");
				next;
			}
			elsif($line =~ /^(\s*)(#?)(write.*)(\S*)(\s*)(\S+\s*)?$/i) {
				$out->print("$1#$3$4$5$6\n");
				next;
			}
			elsif($line =~ /^(\s*)(#?)(dump\S*\s.*)$/i) {
				$out->print("$1#$3\n");
				next;
			}
			elsif($line =~ /^(\s*)(#?)(outp\S*\s.*)$/i) {
				$out->print("$1#$3\n");
				next;
			}
			elsif($line =~ /^(\s*)(#?)(mdar.*)$/i) {
				$out->print("$1#$3\n");
				next;
			}
		}

# 結晶構造出力部分
		my $EndLoop = 0;
		my $EnterCrystal = 0;
		if(($line =~ /^\s*##Crystal\s+structure/ or $line =~ /^\s*cell/i or $line =~ /^\s*(frac|cart)/i)
		    and defined $Crystal) {
			$EnterCrystal = 1;
			
			my @ExpandedAtomSiteList   = $Crystal->GetCExpandedAtomSiteListByOutputMode();
			my $nExpandedAtomSiteList  = @ExpandedAtomSiteList;
			my ($a,$b,$c,$alpha,$beta,$gamma) = $Crystal->LatticeParametersByOutputMode(0);

			while(1) {
				$line = $in->ReadLine();
				if(!defined $line) {
					$EndLoop = 1;
					last;
				}

				$line =~ s/[\r\n]+//;

				if($line =~ /^\s*(#|name)/) {
					$out->print("$line\n");
					next;
				}
				elsif($line =~ /^\s*cell/) {
					$in->ReadLine();
					$out->print("cell\n");
					$out->printf("%10.6f %10.6f %10.6f %12.6f %12.6f %12.6f\n",
						$a, $b, $c, $alpha, $beta, $gamma);
					next;
				}
				elsif($line =~ /^\s*space/) {
					$line = $in->ReadLine();
#print "l2:$line";
					my $SPGName = $Crystal->SPGNameByOutputMode();
					my $iSPG    = $Crystal->iSPGByOutputMode();
					$out->print("space\n");
					$out->print("$iSPG\n");
					$EndLoop = 1;
					last;
				}
				elsif($line =~ /^\s*(frac|cart)/) {
					$out->print("frac\n");
					my @Site;
					for(my $i = 0 ; $i < @ExpandedAtomSiteList ; $i++) {
						$line = $in->ReadLine();
#print "l:$line";
						my ($name0, $type0, $x0, $y0, $z0, $charge0, $occ, $a, $idx, $idy, $idz)
							= Utils::Split("\\s+", $line);
						$Site[$i] = {
							name      => $name0,
							type      => $type0,
							x         => $x0,
							y         => $y0,
							z         => $z0,
							charge    => $charge0,
							occupancy => $occ,
							a         => $a,
							idx       => $idx,
							idy       => $idy,
							idz       => $idz,
							};
					}

#Occupancyが1のサイト
					for(my $i = 0 ; $i < @ExpandedAtomSiteList ; $i++) {
						my $atom      = $ExpandedAtomSiteList[$i];
						my $atomname  = $atom->AtomNameOnly();
						my $charge    = $Site[$i]{charge};
#						my $charge    = ($args{UseSpeciesCharges} eq "no" )?
#								$pPot->{"$atomname-core"}{charge} : 0.0;
						my ($x,$y,$z) = $atom->Position(1);
						my $occupancy = $Site[$i]{occupancy}; #$atom->Occupancy();
#						my $mult      = $atom->Multiplicity();
#						my $id        = $atom->Id();
						next if($occupancy < 0.9999);

						$out->printf("%2s core %10.6f %10.6f %10.6f  "
							    ."%7.4f %7.4f %7.4f  %d %d %d\n", 
							$atomname, $x, $y, $z, $charge, $occupancy, 0.0,
							$Site[$i]{idx}, $Site[$i]{idy}, $Site[$i]{idz});
					}
#Occupancyが1未満のサイト
					for(my $i = 0 ; $i < @ExpandedAtomSiteList ; $i++) {
						my $atom      = $ExpandedAtomSiteList[$i];
						my $atomname  = $atom->AtomNameOnly();
						my $charge    = $Site[$i]{charge};
#						my $charge    = ($args{UseSpeciesCharges} eq "no" )?
#								$pPot->{"$atomname-core"}{charge} : 0.0;
						my ($x,$y,$z) = $atom->Position(1);
						my $occupancy = $Site[$i]{occupancy}; #$atom->Occupancy();
#						my $mult      = $atom->Multiplicity();
#						my $id        = $atom->Id();
						next if($occupancy >= 0.9999);

						$out->printf("%2s core %10.6f %10.6f %10.6f  %7.4f %7.4f %7.4f  "
							    ."%d %d %d\n", 
							$atomname, $x, $y, $z, $charge, $occupancy, 0.0, 
							$Site[$i]{idx}, $Site[$i]{idy}, $Site[$i]{idz});
					}
				}
			}

#print "loop exited\n";
#			last if($EndLoop);
			next;
		}
		$out->print("$line\n");
#print "loop2 exited\n";
	}

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

#Function: "Fit", "Opti", "MD-NVT", "MD-NPT"
sub SaveInputFile
{
	my ($this, $Crystal, $Function, $filename, $IsChooseRandomly, %args) = @_;
#	my $LF = "<br>\n";
#	$CopyScript = 0 if(!defined $CopyScript);

	$args{UseSpeciesCharges} = "no" if(!defined $args{UseSpeciesCharges});
	$args{temperature}   = "300,3000" if(!defined $args{temperature});
	$args{pressure}      = 0.1   if(!defined $args{pressure});
	$args{timestep}      = 0.001 if(!defined $args{timestep});
	$args{tscale}        = 0.10  if(!defined $args{tscale});
	$args{equilibration} = 0.10  if(!defined $args{equilibration});
	$args{production}    = 10.0  if(!defined $args{production});
	$args{sample}        = 0.001 if(!defined $args{sample});
	$args{write}         = 0.005 if(!defined $args{write});
	$args{CopyScript}    = 1     if(!defined $args{CopyScript});

#ファイル名を（ベース名, ディレクトリ名, 拡張子）に分解
	my @filenames = fileparse($filename, "\.[^\.]+");

	my $LibraryPath = $this->LibraryPath();

	my $SampleName = $this->SampleName();
	my ($a,$b,$c,$alpha,$beta,$gamma) 
		= $Crystal->LatticeParametersByOutputMode(0);
	my ($a11, $a12, $a13) = ($a,                       0.0,                     0.0);
	my ($a21, $a22, $a23) = ($b * cos($torad*$gamma),  $b * sin($torad*$gamma), 0.0);
	my $a31 = $c * cos($torad * $beta);
	my $a32 = $c * (cos($torad * $alpha) - cos($torad * $beta) * cos($torad * $gamma)) / sin($torad * $gamma);
	my $a33 = sqrt($c*$c - $a31*$a31 - $a32*$32);
	$a31 = 0.0 if(abs($a31) < 1.0e-6);
	$a32 = 0.0 if(abs($a31) < 1.0e-6);
#print "a: $a31, $a32, $a33\n";
	my $V  = $Crystal->Volume();
	my $RV = 1.0 / $V;

	my @AtomTypeList           = $Crystal->GetCAtomTypeList();
	my $nAtomType              = @AtomTypeList;
	my @AtomSiteList           = $Crystal->GetCAsymmetricAtomSiteList();
	my $nAtomSite              = @AtomSiteList;
	my @ExpandedAtomSiteList   = $Crystal->GetCExpandedAtomSiteListByOutputMode();
	my $nExpandedAtomSiteList  = @ExpandedAtomSiteList;
	my $UseShellModelForCation = $this->UseShellModelForCation();
	my $UseShellModelForAnion  = $this->UseShellModelForAnion();
print "\n";
print "Atom types:\n";
for(my $i = 0 ; $i < @AtomTypeList ; $i++) {
	print "  $i: ", $AtomTypeList[$i]->AtomNameOnly(), "\n";
}

print "\n";
print "Potentials:\n";
	my %CoreShellType;
	my $nShellAtoms = 0;
	for(my $i = 0 ; $i < @AtomTypeList ; $i++) {
		my $type           = $AtomTypeList[$i];
		my $name           = $type->AtomNameOnly();
		my $charge         = $type->Charge();
		my ($core, $shell) = $this->ReadSpecies($name);
#print "Core: $core<br>\n";
#print "Shell: $shell<br>\n";
		if($charge == 0) {
			my ($CoreCharge)  = ($core  =~ /([\+-\.\d]+)\s*$/);
			my ($ShellCharge) = ($shell =~ /([\+-\.\d]+)\s*$/);
			$charge = $CoreCharge + $ShellCharge;
#print "CoreCharge : $CoreCharge<br>\n";
#print "ShellCharge: $ShellCharge<br>\n";
#print "Charge: $charge<br>\n";
		}
		if($core eq '') {
			$CoreShellType{$name} = 'undefined';
		}
		elsif($shell eq '') {
			$CoreShellType{$name} = 'core';
		}
		else {
			$CoreShellType{$name} = 'core-shell';
		}
#print "$name: $CoreShellType{$name}<br>\n";
#print "UseAnion: $UseShellModelForAnion<br>\n";
		if(!$UseShellModelForCation and $charge > 0 and $CoreShellType{$name} eq 'core-shell') {
			$CoreShellType{$name} = 'core';
		}
		if(!$UseShellModelForAnion and $charge < 0 and $CoreShellType{$name} eq 'core-shell') {
			$CoreShellType{$name} = 'core';
		}
		$nShellAtoms++ if($CoreShellType{$name} eq 'core-shell');
	}

	print "<b>Core-Shell type:</b>\n";
	my @list = keys %CoreShellType;
	for(my $i = 0 ; $i < @list ; $i++) {
		my $key = $list[$i];
		my $val = $CoreShellType{$key};
		print "  $key: $val\n";
	}

	print "\n";
	printf "Read potential parameters from [%s]\n", $LibraryPath;
	$this->ClearPotentialLines();
	my $pPot = $this->ReadPotentialLibrary($LibraryPath);
print "  Potential type: $pPot->{type}\n";
foreach my $key (sort keys %$pPot) {
	my $p = $pPot->{$key};
	if($p =~ /hash/i) {
		foreach my $key2 (sort keys %$p) {
			print "   $key:$key2:$p->{$key2}\n";
		}
	}
	else {
		print "   $key: $p\n";
	}
}

	my ($rmin, $rstep, $nr) = (0.01, 0.01, 400);
	$this->SavePotentials('Potentials.csv', $rmin, $rstep, $nr, \@AtomTypeList, $pPot);

#ファイル作製開始
	my $out = JFile->new($filename, 'w');
	unless($out) {
		print "Can not write to [$filename].\n\n";
		return;
	}

	$out->print("##Tasks to be calculated\n");
	$out->print("##  see Kywords in help\n");
	$out->print("##  e.g.: opti md fit\n");
	$out->print("##        prop(elastic,dielectric,piezoe,wave v,bulk/shear/Yongs modulus)\n");
	$out->print("##        phonon raman\n");
	$out->print("##        md conp conse qok nomod pres spat pot efg\n");
	$out->print("##        fit libff libd conp simu conse qok nomod pres spat pot efg\n");
	if($Function =~ /Fit/i) {
		$out->print("#opti prop conp\n");
		$out->print("#md conv\n");
		$out->print("#md conp\n");
		$out->print("fit prop conp\n");
		$out->print("#phonon raman\n");
	}
	elsif($Function =~ /Optimize/i) {
		$out->print("opti prop conp\n");
		$out->print("#md conv\n");
		$out->print("#md conp\n");
		$out->print("#fit prop conp\n");
		$out->print("#phonon raman\n");
	}
	elsif($Function =~ /MD-NVT/i) {
		$out->print("#opti prop conp\n");
		$out->print("md conv\n");
		$out->print("#md conp\n");
		$out->print("#fit prop conp\n");
		$out->print("#phonon raman\n");
	}
	elsif($Function =~ /MD-NPT/i) {
		$out->print("#opti prop conp\n");
		$out->print("#md conv\n");
		$out->print("md conp\n");
		$out->print("#fit prop conp\n");
		$out->print("#phonon raman\n");
	}
	elsif($Function =~ /Phonon/i) {
		$out->print("#opti prop conp\n");
		$out->print("#md conv\n");
		$out->print("#md conp\n");
		$out->print("#fit prop conp\n");
		$out->print("phonon raman\n");
	}
	else {
		$out->print("opti prop conp\n");
		$out->print("#md conv\n");
		$out->print("#md conp\n");
		$out->print("#fit prop conp\n");
	}
	$out->print("\n");

	my $SampleName = $this->SampleName();
	$out->print("title\n");
	$out->print("$SampleName\n");
	$out->print("end\n");
	$out->print("\n");

	$out->print("##Crystal structure\n");
	$out->print("name $SampleName\n");
	my ($a,$b,$c,$alpha,$beta,$gamma) 
		= $Crystal->LatticeParametersByOutputMode(0);
	$out->print("cell\n");
	$out->printf("%10.6f %10.6f %10.6f %12.6f %12.6f %12.6f\n",
		$a, $b, $c, $alpha, $beta, $gamma);
	$out->print("frac\n");

#Occupancyが1のサイト
	for(my $i = 0 ; $i < @ExpandedAtomSiteList ; $i++) {
		my $atom      = $ExpandedAtomSiteList[$i];
		my $atomname  = $atom->AtomNameOnly();
		my $charge    = ($args{UseSpeciesCharges} eq "no" )?
					$pPot->{"$atomname-core"}{charge} : 0.0;
		my ($x,$y,$z) = $atom->Position(1);
		my $occupancy = $atom->Occupancy();
		my $mult      = $atom->Multiplicity();
		my $id        = $atom->Id();
		next if($occupancy < 0.9999);

		$out->printf("%2s core %10.6f %10.6f %10.6f  %7.4f %7.4f %7.4f  %d %d %d\n", 
			$atomname, $x, $y, $z, $charge, $occupancy, 0.0, 1,1,1);
	}

#Occupancyが1未満のサイト
	for(my $i = 0 ; $i < @ExpandedAtomSiteList ; $i++) {
		my $atom      = $ExpandedAtomSiteList[$i];
		my $atomname  = $atom->AtomNameOnly();
		my $charge    = ($args{UseSpeciesCharges} eq "no" )?
					$pPot->{"$atomname-core"}{charge} : 0.0;
		my ($x,$y,$z) = $atom->Position(1);
		my $occupancy = $atom->Occupancy();
		$occupancy = 1 if($IsChooseRandomly);
		my $mult      = $atom->Multiplicity();
		my $id        = $atom->Id();
		next if($occupancy >= 0.9999);

		$out->printf("%2s core %10.6f %10.6f %10.6f  %7.4f %7.4f %7.4f  %d %d %d\n", 
			$atomname, $x, $y, $z, $charge, $occupancy, 0.0, 1,1,1);
	}

	my $SPGName = $Crystal->SPGNameByOutputMode();
	my $iSPG    = $Crystal->iSPGByOutputMode();
	$out->print("space\n");
	$out->print("$iSPG\n");
	$out->print("\n");

#原子種をコメントで出力
	$out->print("##");
	for(my $i = 0 ; $i < $nAtomType ; $i++) {
		my $name1 = $AtomTypeList[$i]->AtomNameOnly();
		$out->print("$name1 ");
	}
	$out->print("\n");

	$out->print("##Original potentials taken from [$LibraryPath]\n");
	my @lines = $this->GetPotentialLines();
	for(my $i = 0 ; $i < @lines ; $i++) {
		$out->print("##  $lines[$i]\n");
	}

if(0) {
foreach my $key (sort keys %$pPot) {
	my $p = $pPot->{$key};
	print "key=$key: $p\n";
	if($p =~ /hash/i) {
		foreach my $key2 (sort keys %$p) {
			print "key2=$key2: $p->{$key2}\n";
		}
	}
}
}

	$this->WriteSpecies($out, $pPot, \@AtomTypeList);
	$out->print("\n");
	$this->WritePotential($out, $pPot, \@AtomTypeList);
	$out->print("\n");

	$out->print("##Short range potential cut off (taper radius can be added)\n");
	$out->print("#cutp 12.0\n");
	$out->print("##pressure P <GPa/kPa/MPa/Pa/atm/Nm-2/kbar (unit, Def=GPa)>\n");
	$out->print("pressure    $args{pressure}\n");
	$out->print("\n");

	my $nstep = int($args{production} / $args{timestep}) + 1;
	print "  Production tine: $args{production} ps at $args{timestep} ps interval "
	     ."with $nstep steps\n";
	my ($T0, $T1) = Utils::Split("\\s*,\\s*", $args{temperature});
	my $Tstep;
	if(defined $T1) {
		$Tstep = ($T1 - $T0) / ($nstep - 1);
	}
	if($Function =~ /MD/i) {
		$out->print("##optimize options\n");
		$out->print("#ftol   1e-005\n");
		$out->print("#gtol   0.0001\n");
		$out->print("#xtol   1e-005\n");
		$out->print("#maxcyc 1000\n");
		$out->print("\n");
		$out->print("##phonon options (omega in cm-1)\n");
		$out->print("#omega 1000 50 180\n");
		$out->print("#dispersion 4 10\n");
		$out->print("#0 0 0 to 0 0.5 0\n");
		$out->print("#0 0.5 0 to 0 0.5 0.5\n");
		$out->print("#0 0.5 0.5 to 0 0 0.5\n");
		$out->print("#0 0 0.5 to 0 0 0\n");
		$out->print("\n");
		$out->print("##MD options\n");
		$out->print("supercell 1 1 1\n");
		$out->print("shellmass 0.1\n");
		$out->print("##integrator <gear/velocity verlet/leapfrog verlet/stochastic> <iter>\n");
		$out->print("integrator leapfrog verlet\n");
		$out->print("##ensemble <NVE/NVT qnose/NPT qnose qpress/NPH>\n");
		$out->print("tau_thermostat  0.1 ps\n");
		$out->print("tau_barostat    0.1 ps\n");
		if($Function =~ /NVT/i) {
			$out->print("ensemble nvt 0.1\n");
			$out->print("#ensemble npt 0.05 0.10\n");
		}
		else {
			$out->print("#ensemble nvt 0.1\n");
			$out->print("ensemble npt 0.05 0.10\n");
		}
		$out->print("##temperature T <C/F/K (unit, Def=K)> <Tstep> <# of steps> <1st step of MD>\n");
		$out->print("temperature $T0 $Tstep $nstep\n");
		$out->print("##Simulation time for equilibration prior to MD production phase\n");
		$out->print("equilibration $args{equilibration} ps\n");
		$out->print("##Temperature scaling time: Def = equilibration\n");
		$out->print("tscale        $args{tscale} ps\n");
		$out->print("##Simulation time to collect production data\n");
		$out->print("production    $args{production} ps\n");
		$out->print("timestep      $args{timestep} ps\n");
		$out->print("sample        $args{sample} ps\n");
		$out->print("write         $args{write} ps\n");
		$out->print("dump every        $filenames[0].res\n");
#		$out->print("dump              $filenames[0].grs\n");
		$out->print("output trajectory $filenames[0].trg\n");
		$out->print("output history    $filenames[0].his\n");
		$out->print("output pressure   $filenames[0].pre\n");
		$out->print("output movie arc  $filenames[0]\n");
#		$out->print("mdarchive         $filenames[0].arc\n");
	}
	elsif($Function =~ /opti/i) {
		$out->print("##optimize options\n");
		$out->print("ftol   1e-005\n");
		$out->print("gtol   0.0001\n");
		$out->print("xtol   1e-005\n");
		$out->print("maxcyc 1000\n");
		$out->print("\n");
		$out->print("##phonon options\n");
		$out->print("#omega 1000 50 180\n");
		$out->print("#dispersion 4 10\n");
		$out->print("#0 0 0 to 0 0.5 0\n");
		$out->print("#0 0.5 0 to 0 0.5 0.5\n");
		$out->print("#0 0.5 0.5 to 0 0 0.5\n");
		$out->print("#0 0 0.5 to 0 0 0\n");
		$out->print("\n");
		$out->print("##MD options\n");
		$out->print("#supercell 1 1 1\n");
		$out->print("#shellmass 0.1\n");
		$out->print("##integrator <gear/velocity verlet/leapfrog verlet/stochastic> <iter>\n");
		$out->print("#integrator leapfrog verlet\n");
		$out->print("##ensemble <NVE/NVT qnose/NPT qnose qpress/NPH>\n");
		$out->print("#ensemble nvt 0.1\n");
		$out->print("#ensemble npt 0.05 0.10\n");
		$out->print("#tau_thermostat  0.1 ps\n");
		$out->print("#tau_barostat    0.1 ps\n");
		$out->print("##temperature T <C/F/K (unit, Def=K)> <step> <# of steps> <1st step of MD>\n");
		$out->print("#temperature $T0 $Tstep $nstep\n");
		$out->print("##Simulation time for equilibration prior to MD production phase\n");
		$out->print("#equilibration $args{equilibration} ps\n");
		$out->print("##Temperature scaling time: Def = equilibration\n");
		$out->print("#tscale        $args{tscale} ps\n");
		$out->print("##Simulation time to collect production data\n");
		$out->print("#production    $args{production} ps\n");
		$out->print("#timestep      $args{timestep} ps\n");
		$out->print("#sample        $args{sample} ps\n");
		$out->print("#write         $args{write} ps\n");
		$out->print("dump every        $filenames[0].res\n");
#		$out->print("dump              $filenames[0].grs\n");
		$out->print("##output filtype filename\n");
		$out->print("##  filetype: marvin, thbrel, xtl, xr, cssr, arc, xyz, trajectory, history, fdf\n");
		$out->print("##            drv, frc, cif, str, pressure, osc, phonon, pdf, qbo, cosmo, lammps\n");
		$out->print("##            lammps_pots, dcd, eig, cml, inertia, kpt, castep_phonon, shengBTE\n");
		$out->print("output trajectory $filenames[0].trg\n");
		$out->print("output history    $filenames[0].his\n");
		$out->print("output pressure   $filenames[0].pre\n");
		$out->print("output movie arc  $filenames[0]\n");
#		$out->print("mdarchive         $filenames[0].arc\n");
	}
	elsif($Function =~ /phon/i) {
		$out->print("##optimize options\n");
		$out->print("#ftol   1e-005\n");
		$out->print("#gtol   0.0001\n");
		$out->print("#xtol   1e-005\n");
		$out->print("#maxcyc 1000\n");
		$out->print("\n");
		$out->print("##phonon options\n");
		$out->print("omega 1000 50 180\n");
		$out->print("dispersion 4\n");
		$out->print("0 0 0 to 0 0.5 0\n");
		$out->print("0 0.5 0 to 0 0.5 0.5\n");
		$out->print("0 0.5 0.5 to 0 0 0.5\n");
		$out->print("0 0 0.5 to 0 0 0\n");
		$out->print("\n");
		$out->print("##MD options\n");
		$out->print("#supercell 1 1 1\n");
		$out->print("#shellmass 0.1\n");
		$out->print("##integrator <gear/velocity verlet/leapfrog verlet/stochastic> <iter>\n");
		$out->print("#integrator leapfrog verlet\n");
		$out->print("##ensemble <NVE/NVT qnose/NPT qnose qpress/NPH>\n");
		$out->print("#ensemble nvt 0.1\n");
		$out->print("#ensemble npt 0.05 0.10\n");
		$out->print("#tau_thermostat  0.1 ps\n");
		$out->print("#tau_barostat    0.1 ps\n");
		$out->print("##temperature T <C/F/K (unit, Def=K)> <step> <# of steps> <1st step of MD>\n");
		$out->print("#temperature $T0 $Tstep $nstep\n");
		$out->print("##Simulation time for equilibration prior to MD production phase\n");
		$out->print("#equilibration $args{equilibration} ps\n");
		$out->print("##Temperature scaling time: Def = equilibration\n");
		$out->print("#tscale        $args{tscale} ps\n");
		$out->print("##Simulation time to collect production data\n");
		$out->print("#production    $args{production} ps\n");
		$out->print("#timestep      $args{timestep} ps\n");
		$out->print("#sample        $args{sample} ps\n");
		$out->print("#write         $args{write} ps\n");
		$out->print("#dump every        $filenames[0].res\n");
#		$out->print("#dump              $filenames[0].grs\n");
		$out->print("##output filtype filename\n");
		$out->print("##  filetype: marvin, thbrel, xtl, xr, cssr, arc, xyz, trajectory, history, fdf\n");
		$out->print("##            drv, frc, cif, str, pressure, osc, phonon, pdf, qbo, cosmo, lammps\n");
		$out->print("##            lammps_pots, dcd, eig, cml, inertia, kpt, castep_phonon, shengBTE\n");
		$out->print("output trajectory $filenames[0].trg\n");
		$out->print("#output history    $filenames[0].his\n");
		$out->print("#output pressure   $filenames[0].pre\n");
		$out->print("output movie arc  $filenames[0]\n");
#		$out->print("#mdarchive         $filenames[0].arc\n");
	}
	else {
		$out->print("##optimize options\n");
		$out->print("#ftol   1e-005\n");
		$out->print("#gtol   0.0001\n");
		$out->print("#xtol   1e-005\n");
		$out->print("#maxcyc 1000\n");
		$out->print("\n");
		$out->print("##phonon options\n");
		$out->print("#omega 1000 50 180\n");
		$out->print("#dispersion 4\n");
		$out->print("#0 0 0 to 0 0.5 0\n");
		$out->print("#0 0.5 0 to 0 0.5 0.5\n");
		$out->print("#0 0.5 0.5 to 0 0 0.5\n");
		$out->print("#0 0 0.5 to 0 0 0\n");
		$out->print("\n");
		$out->print("##MD options\n");
		$out->print("#supercell 1 1 1\n");
		$out->print("#shellmass 0.1\n");
		$out->print("##integrator <gear/velocity verlet/leapfrog verlet/stochastic> <iter>\n");
		$out->print("#integrator leapfrog verlet\n");
		$out->print("##ensemble <NVE/NVT qnose/NPT qnose qpress/NPH>\n");
		$out->print("#ensemble nvt 0.1\n");
		$out->print("#ensemble npt 0.05 0.10\n");
		$out->print("#tau_thermostat  0.1 ps\n");
		$out->print("#tau_barostat    0.1 ps\n");
		$out->print("##temperature T <C/F/K (unit, Def=K)> <step> <# of steps> <1st step of MD>\n");
		$out->print("#temperature $T0 $Tstep $nstep\n");
		$out->print("##Simulation time for equilibration prior to MD production phase\n");
		$out->print("#equilibration $args{equilibration} ps\n");
		$out->print("##Temperature scaling time: Def = equilibration\n");
		$out->print("#tscale        $args{tscale} ps\n");
		$out->print("##Simulation time to collect production data\n");
		$out->print("#production    $args{production} ps\n");
		$out->print("#timestep      $args{timestep} ps\n");
		$out->print("#sample        $args{sample} ps\n");
		$out->print("#write         $args{write} ps\n");
		$out->print("#dump every        $filenames[0].res\n");
#		$out->print("#dump              $filenames[0].grs\n");
		$out->print("##output filtype filename\n");
		$out->print("##  filetype: marvin, thbrel, xtl, xr, cssr, arc, xyz, trajectory, history, fdf\n");
		$out->print("##            drv, frc, cif, str, pressure, osc, phonon, pdf, qbo, cosmo, lammps\n");
		$out->print("##            lammps_pots, dcd, eig, cml, inertia, kpt, castep_phonon, shengBTE\n");
		$out->print("#output trajectory $filenames[0].trg\n");
		$out->print("#output history    $filenames[0].his\n");
		$out->print("#output pressure   $filenames[0].pre\n");
		$out->print("#output movie arc  $filenames[0]\n");
#		$out->print("#mdarchive         $filenames[0].arc\n");
	}

#	$out->print("output xr $filenames[0]\n");
#	$out->print("output marvin $filenames[0].mvn\n");

	$out->print("\n\n");
	
	close OUT;

	return 1;
}


1;