package GULP;
use Exporter;
@ISA = qw(Exporter);

#公開したいサブルーチン
@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   = (
	$GULPLibDir,
	Deps::MakePath($ProgramDir, ["MD", "gulp-OtherLibraries",
		"www.ucl.ac.uk", "klmc", "Potentials", "Library"], 0),
	Deps::MakePath($ProgramDir, ["MD", "gulp-OtherLibraries"], 0),
	);

#===============================================
# コンストラクタ、デストラクタ
#===============================================
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 LibraryDir { 
	my $this = shift;
	return $this->{'LibraryDir'} if($this->{'LibraryDir'});
	return $this->{'LibraryDir'} = $GULPLibDir;
}
sub SetLibraryDir { my($this,$d)=@_; return $this->{'LibraryDir'} = $d; }

sub LibraryPath
{
	my ($this) = @_;
	return $this->{'LibraryPath'};
}
sub SetLibraryFile
{
	my ($this, $file) = @_;
	$this->{'LibraryFile'} = $file;
	$GULPLibDir = $this->LibraryDir();
	my $GULPLibaryPath = ($file =~ /[\\\/]/)? $file : Deps::MakePath($GULPLibDir, $file, 0);
	$this->{'LibraryPath'} = $GULPLibaryPath;
	return $this->{'LibraryPath'};
}

sub SetLibraryFile
{
	my ($this, $file) = @_;
	$this->{'LibraryFile'} = $file;
	$GULPLibDir = $this->LibraryDir();
	my $GULPLibaryPath = ($file =~ /[\\\/]/)? $file : Deps::MakePath($GULPLibDir, $file, 0);
	$this->{'LibraryPath'} = $GULPLibaryPath;
	return $this->{'LibraryPath'};
}
sub LibraryPath
{
	my ($this) = @_;
	return $this->{'LibraryPath'};
}

sub ClearPotentialLines
{
	my ($this) = @_;
	return $this->{pLibraryLines} = [];
}

sub AddPotentialLines
{
	my ($this, @a) = @_;

	$this->{pLibraryLines} = [] if(!defined $this->{pLibraryLines});

	for(my $i = 0 ; $i < @a ; $i++) {
		Utils::DelSpace($a[$i]);
	}
	my $p = $this->{pLibraryLines};
	for(my $j = 0 ; $j < @a ; $j++) {
		my $IsHit = 1;
		for(my $i = 0 ; $i < @$p ; $i++) {
			if($p->[$i] eq $a[$j]) {
				$IsHit = 0;
				last;
			}
		}
		push(@$p, $a[$j]) if($IsHit);
	}
	return @$p;
}

sub GetPotentialLines
{
	my ($this) = @_;

	$this->{pLibraryLines} = [] if(!defined $this->{pLibraryLines});
	my $p = $this->{pLibraryLines};
	return @$p;
}

sub SetSampleName
{
	my ($this, $name) = @_;
	return $this->{'SampleName'} = $name;
}
sub SampleName
{
	my ($this) = @_;
	return $this->{'SampleName'};
}

#===============================================
# 一般関数
#===============================================
sub SetUseShellModelForCation
{
	my ($this,$val) = @_;
	return $this->{'UseShellModelForCation'} = $val;
}

sub UseShellModelForCation
{
	my ($this) = @_;
	return $this->{'UseShellModelForCation'};
}

sub SetUseShellModelForAnion
{
	my ($this,$val) = @_;
	return $this->{'UseShellModelForAnion'} = $val;
}

sub UseShellModelForAnion
{
	my ($this) = @_;
	return $this->{'UseShellModelForAnion'};
}

sub ReadSpecies
{
	my ($this, $atom, $cstype) = @_;
#print "cs=$atom, $cstype\n";	
	return (undef, undef, undef) if($cstype eq 'undefined');

	my $GULPLibaryPath = $this->LibraryPath();
#print "lib: $GULPLibaryPath\n";

	if(!open(IN, "<$GULPLibaryPath")) {
		print "Error in LAMMPS::ReadBuckinghamPotential: Can not read [$GULPLibaryPath]\n";
		return ();
	}

	while(<IN>) {
		last if($_ =~ /^\s*species/i);
	}
#Read species (ionic charges)
	my ($CoreCharge, $ShellCharge);
	while(<IN>) {
		my $line = $_;
		chomp($line);
#print "l:$line\n";
		if($line =~ /^$atom\s+core\s+(\S+)/i) {
			$CoreCharge = $1;
			$this->AddPotentialLines($line);
		}
		if($line =~ /^$atom\s+shel\s+(\S+)/i) {
			$ShellCharge = $1;
			$this->AddPotentialLines($line);
		}
		last if($line =~ /^\s*buck/i);
	}
	close(IN);

#print "Z: $CoreCharge, $ShellCharge\n";
	return ($CoreCharge+$ShellCharge, $CoreCharge, $ShellCharge);
}

sub ReadSpring
{
	my ($this, $atom, $cstype) = @_;

	my $GULPLibaryPath = $this->LibraryPath();
	return ('', '') if($cstype eq 'core' or $cstype eq 'undefined');
	if(!open(IN, "<$GULPLibaryPath")) {
		print "Error in GULP::ReadBuckinghamPotential: Can not read [$GULPLibaryPath]\n";
		return ();
	}

	my $Spring;
	while(<IN>) {
		last if($_ =~ /^\s*spring/i);
	}
#Read spring parameters
	while(<IN>) {
		my $line = $_;
		chomp($line);
		if($line =~ /^$atom\s/i) {
			$Spring = $line;
			$this->AddPotentialLines($line);
		}
	}
	close(IN);
	unless($Spring =~ /^(\S+)\s+(\S+)\s+(\S+)/) {
		$Spring = "$Spring 0" if($Spring ne '');
	}

	return $Spring;
}

sub AddAllPotentialsToHash
{
	my ($this, $in, $type, $natom, $pHash, $IsPrint) = @_;

	$IsPrint = 1 if(!defined $IsPrint);
#print "type:$type\n";

	my %array;
	while(1) {
		my $pos = $in->tell();
		my $line = $in->ReadLine();
#print "l2:$line";
		last if(!defined $line);
#		Utils::DelSpace($line);
#		next if($line =~ /^#/ or $line eq '');
		next if($line =~ /^#/ or $line =~ /^\s/);
#		Utils::DelSpace($line);
		next if($line eq '');

		if($line !~ /^cova\s/ 
		   and $line !~ /^mass\s/i
		   and $line !~ /^[A-Z][a-z0-9\-+]*\s/ 
		   and $line !~ /^\s*[A-Z][a-z0-9\-]*_+[A-Za-z0-9+\-_]*\s/) {
			$in->seek($pos, 0);
			last;
		}

		my @a = Utils::Split("\\s+", $line);
		if($a[scalar @a - 1] eq '&') {
			pop @a;
			@a = (@a, Utils::Split($in->ReadLine()));
		}
		if($a[0] eq 'mass') {
			my $n = @a;
			@a = @a[1..$n];
		}
#print "   natom=$natom: $a[0] - $a[2] - $a[4] - $a[6]\n";
		if($type eq 'element') {
			my %hash = (
				Type      => $type,
				nAtom     => 1,
				AtomName1 => $a[1],
				);
			$array{"$a[1]:$type"} = \%hash;
		}
		elsif($type eq 'SX-1') {
# Born-Mayer Potential CONSTANT (J/A) $F0 = 6.947700141e-21
#   A(I,J)=F0*B(I,J)*DEXP((PA(I)+PA(J))/B(I,J))
#Zn	65.39	2	1.5145	0.1200	0	0
			if($a[4] != 0.0) {
				my %hash = (
					Type      => $type,
					nAtom     => 1,
					AtomName1 => $a[0],
					);
				$array{"$a[0]:$type"} = \%hash;
			}
		}
		elsif($natom == 1) {
			my %hash = (
				Type      => $type,
				nAtom     => 1,
				AtomName1 => $a[0],
				);
			$array{"$a[0]:$type"} = \%hash;
		}
		elsif($natom == 2) {
			my %hash = (
				Type      => $type,
				nAtom     => 2,
				AtomName1 => $a[0],
				AtomName2 => $a[2],
				);
			@a = sort ($a[0], $a[2]);	
			$array{"$a[0]-$a[1]:$type"} = \%hash;
		}
		elsif($natom == 3) {
			my %hash = (
				Type      => $type,
				nAtom     => 3,
				AtomName1 => $a[0],
				AtomName1 => $a[2],
				AtomName1 => $a[4],
				);
			@a = sort ($a[0], $a[2], $a[4]);
			$array{"$a[0]-$a[1]-$a[2]:$type"} = \%hash;
		}
		elsif($natom == 4) {
			my %hash = (
				Type      => $type,
				nAtom     => 4,
				AtomName1 => $a[0],
				AtomName1 => $a[2],
				AtomName1 => $a[4],
				);
			@a = sort ($a[0], $a[2], $a[4], $a[6]);
			$array{"$a[0]-$a[1]-$a[2]-$a[3]:$type"} = \%hash;
		}
	}
	foreach my $key (keys %array) {
#print "key:$key\n";		
		$pHash->{$key} = $array{$key};
	}
#foreach my $key (keys %$pHash) {
#print "key:$key\n";		
#}

	return $pHash;
}

sub ReadPotentialToHash
{
	my ($this, $db, $IsPrint) = @_;
	$IsPrint = 1 if(!defined $IsPrint);
#$IsPrint=1;

	my %array;
	my $c = 0;
	my $in = JFile->new($db, 'r') or return undef;
	while(1) {
		my $line = $in->ReadLine();
		last if(!defined $line);
		next if($line =~ /^#/ or $line =~ /^\s/);
		Utils::DelSpace($line);
		next if($line eq '');
#print "l0:$line\n";

		my $type = $line;
		Utils::DelSpace($type);
print "  potential type=$type\n" if($IsPrint);

		if($type =~ /^end$/i) {
			next;
#			last;
		}
		elsif($type =~ /^keyword/i) {
			next;
		}
		elsif($type =~ /^keywords/i) {
			next;
		}
		elsif($type =~ /^cutp/i) {
			next;
		}
		elsif($type =~ /^rtol/i) {
			next;
		}
		elsif($type =~ /^shellmass/i) {
			$in->ReadLine();
			next;
		}
		elsif($type =~ /^reaxFFvdwcutoff/) {
		}
		elsif($type =~ /^reaxFFqcutoff/) {
		}
		elsif($type =~ /^reaxFFtol/) {
		}
		elsif($type =~ /^reaxff0_bond/i) {
		}
		elsif($type =~ /^reaxff0_over/i) {
		}
		elsif($type =~ /^reaxff0_valence/i) {
		}
		elsif($type =~ /^reaxff0_penalty/i) {
		}
		elsif($type =~ /^reaxff0_torsion/i) {
		}
		elsif($type =~ /^reaxff0_vdw/i) {
		}
		elsif($type =~ /^reaxff0_lonepair/i) {
		}
		elsif($type =~ /^cuts\s/i) {
		}
		elsif($type =~ /^morse/i) {
			$this->AddAllPotentialsToHash($in, $type, 2, \%array);
		}
		elsif($type =~ /^buckingham/i or $type =~ /^buck/i) {
			$this->AddAllPotentialsToHash($in, $type, 2, \%array);
		}
		elsif($type =~ /^SX-1/i) {
			$type =~ s/\s.*$//;
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^reaxff1_radii/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^reaxff1_valence/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^reaxff1_over/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^reaxff1_under/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^reaxff1_lonepair/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^reaxff1_angle/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^reaxff1_lonepair/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^reaxff1_morse/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^reaxff_chi/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^reaxff_mu/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^reaxff_gamma/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^reaxff2_bo/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^reaxff2_bond/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^reaxff2_over/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^reaxff2_morse/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^reaxff2_over/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^reaxff2_penalty/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^reaxff2_pen/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^reaxff3_angle/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^reaxff3_penalty/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^reaxff3_conjugation/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^reaxff3_hbond/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^reaxff4_torsion/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^torsion/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^inversion/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^element/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^species/i) {
#print "hit\n";
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^spring/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^general/i) {
			$this->AddAllPotentialsToHash($in, $type, 2, \%array);
		}
		elsif($type =~ /^bcross/i) {
			$this->AddAllPotentialsToHash($in, $type, 3, \%array);
		}
		elsif($type =~ /^urey-bradley/i) {
			$this->AddAllPotentialsToHash($in, $type, 3, \%array);
		}
		elsif($type =~ /^many/i) {
			$this->AddAllPotentialsToHash($in, $type, 2, \%array);
		}
		elsif($type =~ /^meam_functional/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^meam_density/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^meam_screening/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^meam_screen/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^meam_rhotype/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^baskes/i) {
			$this->AddAllPotentialsToHash($in, $type, 2, \%array);
		}
		elsif($type =~ /^mei-davenport/i) {
			$this->AddAllPotentialsToHash($in, $type, 2, \%array);
		}
		elsif($type =~ /^eam_fun/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^eam_den/i) {
			$this->AddAllPotentialsToHash($in, $type, 2, \%array);
		}
		elsif($type =~ /^eam_alloy/i) {
			$this->AddAllPotentialsToHash($in, $type, 2, \%array);
		}
		elsif($type =~ /^cfm_harm/i) {
			$this->AddAllPotentialsToHash($in, $type, 2, \%array);
		}
		elsif($type =~ /^cfm_power/i) {
			$this->AddAllPotentialsToHash($in, $type, 2, \%array);
		}
		elsif($type =~ /^cfm_gauss/i) {
			$this->AddAllPotentialsToHash($in, $type, 2, \%array);
		}
		elsif($type =~ /^cfm_fermi/i) {
			$this->AddAllPotentialsToHash($in, $type, 2, \%array);
		}
		elsif($type =~ /^vbo_twobody/i) {
			$this->AddAllPotentialsToHash($in, $type, 2, \%array);
		}
		elsif($type =~ /^uff1/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^smelectronegativity/i) {
			$this->AddAllPotentialsToHash($in, $type, 1, \%array);
		}
		elsif($type =~ /^rydberg/i) {
			$this->AddAllPotentialsToHash($in, $type, 2, \%array);
		}
		elsif($type =~ /^qerfc/i) {
			$this->AddAllPotentialsToHash($in, $type, 2, \%array);
		}
		elsif($type =~ /^fermi/i) {
			$this->AddAllPotentialsToHash($in, $type, 2, \%array);
		}
		elsif($type =~ /^sw2/i) {
			$this->AddAllPotentialsToHash($in, $type, 2, \%array);
		}
		elsif($type =~ /^sw3/i) {
			$this->AddAllPotentialsToHash($in, $type, 3, \%array);
		}
		elsif($type =~ /^srglue/i) {
			$this->AddAllPotentialsToHash($in, $type, 2, \%array);
		}
		elsif($type =~ /^hydrogen-bond/i) {
			$this->AddAllPotentialsToHash($in, $type, 3, \%array);
		}
		elsif($type =~ /^three bond/i or $type =~ /^three/i) {
			$this->AddAllPotentialsToHash($in, $type, 3, \%array);
		}
		elsif($type =~ /^outofplane/i) {
			$this->AddAllPotentialsToHash($in, $type, 4, \%array);
		}
		elsif($type =~ /^polynomial/i) {
			my $npoly = $in->ReadLine() + 0;
			$this->AddAllPotentialsToHash($in, $type, 2, \%array);
		}
		elsif($type =~ /^harm/i) {
			$this->AddAllPotentialsToHash($in, $type, 2, \%array);
		}
		elsif($type =~ /^epsilon/i) {
			$this->AddAllPotentialsToHash($in, $type, 2, \%array);
		}
		elsif($type =~ /^edip_coordination/i) {
			$this->AddAllPotentialsToHash($in, $type, 2, \%array);
		}
		elsif($type =~ /^edip_twobody/i) {
			$this->AddAllPotentialsToHash($in, $type, 2, \%array);
		}
		elsif($type =~ /^edip_threebody/i) {
			$this->AddAllPotentialsToHash($in, $type, 3, \%array);
		}
		elsif($type =~ /^botwobody/i) {
			$this->AddAllPotentialsToHash($in, $type, 2, \%array);
		}
		elsif($type =~ /^borepulsive/i) {
			$this->AddAllPotentialsToHash($in, $type, 2, \%array);
		}
		elsif($type =~ /^boattractive/i) {
			$this->AddAllPotentialsToHash($in, $type, 2, \%array);
		}
		elsif($type =~ /^bocoordination/i) {
			$this->AddAllPotentialsToHash($in, $type, 2, \%array);
		}
		elsif($type =~ /^tsuneyuki/i) {
			$this->AddAllPotentialsToHash($in, $type, 2, \%array);
		}
		elsif($type =~ /^edip_zmax/i) {
		}
		elsif($type =~ /^lennard/i) {
			while(1) {
				my $pos = $in->tell();
				$line = $in->ReadLine();
				last if(!defined $line);
				if($line !~ /^[0-9\.eEdD+\-]+\s/ and $line !~ /^[A-Z][a-z]*\s/) {
					$in->seek($pos, 0);
					last;
				}
			}
		}
		elsif($type =~ /^damped_dispersion/i) {
			$this->AddAllPotentialsToHash($in, $type, 2, \%array);
		}
		elsif($type =~ /^lenn\s/i) {
			$this->AddAllPotentialsToHash($in, $type, 2, \%array);
		}
		elsif($type =~ /^aolm/i) {
			$this->AddAllPotentialsToHash($in, $type, 2, \%array);
		}
		elsif($type =~ /^tortaper/i) {
			$this->AddAllPotentialsToHash($in, $type, 4, \%array);
		}
		else {
			print "Error: Invalid type [$type].\n";
			exit;
		}
	}
	$in->Close();

	my %atomshash;
	foreach my $key (keys %array) {
#print "key:$key\n";
		$key =~ s/:.*$//;
		my @a = Utils::Split("-", $key);
		for(my $i = 0 ; $i < @a ; $i++) {
			$a[$i] =~ s/[\d_].*$//;
			$atomshash{$a[$i]}++;
		}
	}
	my @atoms = keys %atomshash;
#for(my $i = 0 ; $i < @atoms ; $i++) {
#print " $atoms[$i]";
#}
#print "\n";

	return (\%array, \@atoms);
}

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 SX1toBMH
{
	my ($this, $ai, $bi, $ci, $aj, $bj, $cj) = @_;

# Born-Mayer Potential CONSTANT (J/A)
	my $f0 = 6.947700141e-21;

#   A(I,J)=F0*B(I,J)*DEXP((PA(I)+PA(J))/B(I,J))
	return (0.0, 0.1, 0.0) if($bi == 0.0 or $bj == 0.0);

	my $bij = $bi + $bj;
	my $Cij = $ci + $cj;
	my $Aij = $f0 * $bij * exp(($ai + $aj) / $bij) / $e;
#print "aa: $Aij = $f0 * $bij * exp(($ai + $aj) / $bij) / $e\n";

	return ($Aij, $bij, $Cij);
}

sub ReadSX1
{
	my ($this, $atom, $cstype) = @_;

	my $GULPLibaryPath = $this->LibraryPath();
#	return ('', '') if($cstype eq 'shell' or $cstype eq 'undefined');
	
	my $in = JFile->new($GULPLibaryPath, 'r');
	if(!$in) {
		print "Error in LAMMPS::ReadBuckinghamPotential: Can not read [$GULPLibaryPath]\n";
		return {};
	}

	my %hash;
	$in->SkipTo("SX-1");
	while(1) {
		my $line = $in->ReadLine();
		last if(!defined $line);
#print "l:$line";
		my @a = Utils::Split("\\s+", $line);
		$a[0] =~ s/[0-9+\-_].*$//;
		if($a[0] eq $atom) {
			%hash = (
				AtomName => $a[0],
				Mass     => $a[1],
				Charge   => $a[2],
				ai       => $a[3],
				bi       => $a[4],
				ci       => $a[5],
				ARAD     => $a[6],
				line     => $line,
				);
			$this->AddPotentialLines($line);
			last;
		}
	}
	$in->Close();

	return \%hash;
}

sub PotentialLibrary
{
	my ($this, $LibraryPath) = @_;

	my %hash;
	my $in = JFile->new($LibraryPath, 'r') or return undef;
	my $pos;
	while(1) {
		my $line = $in->ReadLine();
		last if(!defined $line);
		$line =~ s/#.*$//;
		Utils::DelSpace($line);
		next if($line eq '');
#print "l:$line\n";

		if($line =~ /^spec/i) {
			my ($s, $n) = Utils::Split("\\s+", $line);
			$hash{species} = $n;
#print "l:$line, $n\n";
			while(1) {
				$pos = $in->tell();
				$line = $in->ReadLine();
				my ($name, $charge) = Utils::Split("\\s+", $line);
				if(!defined $charge or $name !~ /^[A-Z][a-z]?[+\-_\d]*$/) {
					$in->seek($pos);
					last;
				}
				$name =~ s/[+\-_\d]+//g;
				$hash{"$name-core"}{charge} = $charge;
#print "$name-$type:charge = $charge\n";
#exit;
			}
		}
		elsif($line =~ /^buck/i) {
			$hash{buckingham} = 1;
			$hash{type} = 'buckingham';
			while(1) {
				$pos = $in->tell();
				$line = $in->ReadLine();
				my ($name1, $type1, $name2, $type2, $A, $rho, $C, $rmin, $rmax) 
					= Utils::Split("\\s+", $line);
				if(!defined $C or $name1 !~ /^[A-Z][a-z]?[+\-_\d]*$/ 
				   or $name2 !~ /^[A-Z][a-z]?[+\-_\d]*$/) {
					$in->seek($pos);
					last;
				}
				$name1 =~ s/[+\-_\d]+//g;
				$name2 =~ s/[+\-_\d]+//g;
				$hash{"$name1-$type1-$name2-$type2"} 
					= $hash{"$name2-$type2-$name1-$type1"} 
						= {
						A    => $A,
						rho  => $rho,
						C    => $C,
						rmin => $rmin,
						rmax => $rmax,
						}
#print "$name1-$type1-$name2-$type2:A = $A\n";
#exit;
			}
		}
		elsif($line =~ /^mors/i) {
			$hash{morse} = 1;
			$hash{type} = 'morse';
			while(1) {
				$pos = $in->tell();
				$line = $in->ReadLine();
				my ($name1, $type1, $name2, $type2, $De, $a, $r0, $fsCoulomb, $rmin, $rmax) 
					= Utils::Split("\\s+", $line);
				if(!defined $r0 or $name1 !~ /^[A-Z][a-z]?[+\-_\d]*$/
				   or $name2 !~ /^[A-Z][a-z]?[+\-_\d]*$/) {
					$in->seek($pos);
					last;
				}
				$name1 =~ s/[+\-_\d]+//g;
				$name2 =~ s/[+\-_\d]+//g;
				$hash{"$name1-$type1-$name2-$type2"} 
					= $hash{"$name2-$type2-$name1-$type1"} 
						= {
						De        => $De,
						a         => $a,
						r0        => $r0,
						fsCoulomb => $fsCoulomb,
						rmin      => $rmin,
						rmax      => $rmax,
						}
#print "$name1-$type1-$name2-$type2:De = $De\n";
#exit;
			}
		}
	}
	$in->Close();
	return \%hash;
}

sub ReadBuckinghamPotential
{
	my ($this, $atom1, $atom2, $cstype1, $cstype2) = @_;

	my $GULPLibaryPath = $this->LibraryPath();
#print "read [$GULPLibaryPath] for $atom1 - $atom2\n";

	my $pSX1Hash1 = $this->ReadSX1($atom1, $cstype1);
	my $pSX1Hash2 = $this->ReadSX1($atom2, $cstype1);
#print "$pSX1Hash1, $pSX1Hash2\n";
#print "$atom1: $pSX1Hash1->{Charge}, $pSX1Hash1->{ai}, $pSX1Hash1->{bi}\n";
#print "$atom2: $pSX1Hash2->{Charge}, $pSX1Hash2->{ai}, $pSX1Hash2->{bi}\n";
 	if(defined $pSX1Hash1->{ai} and defined $pSX1Hash2->{ai}) {
		my ($Aij, $bij, $Cij) = $this->SX1toBMH(
			$pSX1Hash1->{ai}, $pSX1Hash1->{bi}, $pSX1Hash1->{ci}, 
			$pSX1Hash2->{ai}, $pSX1Hash2->{bi}, $pSX1Hash2->{ci}
			);
#print " return ($pSX1Hash1->{Charge}, $pSX1Hash2->{Charge}, $Aij, $bij, $Cij)\n";
		return ($pSX1Hash1->{Charge}, $pSX1Hash2->{Charge}, $Aij, $bij, $Cij);
	}

	my ($Charge1, $CoreCharge1, $ShellCharge1) = $this->ReadSpecies($atom1, $cstype1);
	my ($Charge2, $CoreCharge2, $ShellCharge2) = $this->ReadSpecies($atom2, $cstype2);
	my $Spring1 = $this->ReadSpring($atom1, $cstype1);
	my $Spring2 = $this->ReadSpring($atom2, $cstype2);
#print "A1: $atom1 = $Charge1, $CoreCharge1, $ShellCharge1\n";
#print "A2: $atom2 = $Charge2, $CoreCharge2, $ShellCharge2\n";

	my $in = JFile->new($GULPLibaryPath, 'r');
	if(!$in) {
		print "Error in LAMMPS::ReadBuckinghamPotential: Can not read [$GULPLibaryPath]\n";
		return ();
	}
	$in->SkipTo("\\s*buck");
#Read buckingham two body parameters
	while(1) {
		my $line = $in->ReadLine();
		last if(!defined $line);

		if($line =~ /^$atom1\s+core\s+$atom2\s+core\s/i) {
			my @a = Utils::Split("\\s+", $line);
			my ($Aij, $bij, $Cij) = ($a[4], $a[5], $a[6]);
#print "1 $atom1 - $atom2: ($Charge1, $Charge2, $Aij, $bij, $Cij)\n";
			$this->AddPotentialLines($line);
			return ($Charge1, $Charge2, $Aij, $bij, $Cij);
		}
		elsif($line =~ /^$atom2\s+core\s+$atom1\s+core\s/i) {
			my @a = Utils::Split("\\s+", $line);
			my ($Aij, $bij, $Cij) = ($a[4], $a[5], $a[6]);
#print "2 $atom2 - $atom1: ($Charge2, $Charge1, $Aij, $bij, $Cij)\n";
			$this->AddPotentialLines($line);
			return ($Charge1, $Charge2, $Aij, $bij, $Cij);
		}
		last if($line =~ /^\s*spring/i);
	}
	return ();
}

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($line =~ /\*\*/);
		
		my ($tag, $unit) = ($line =~ /(\S.*)\s+\((.*?)\)\s*=/);
print "  $nProperties: $tag $unit\n";
		$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+)/);
	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) = @_;

	my $k = 1.0;

	my @ExpandedAtomSiteList   = $cry->GetCExpandedAtomSiteListByOutputMode();
	my $nSite = @ExpandedAtomSiteList;

	my $Crystal = new Crystal;

	my $line;
	$in->SkipTo("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;
	$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";

	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 ($x,  $y,  $z)  = Utils::Split("\\s+", $in->ReadLine());
		my ($vx, $vy, $vz) = Utils::Split("\\s+", $in->ReadLine());

		$AtomCount{$name}++;
		my $Label = $name . $AtomCount{$name};
printf "   $i: $name: $Label (%12.4g, %12.4g, %12.4g)\n", $x, $y, $z;
		$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);
print "nStep: $nStep at each $nOutputInterval steps\n";

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

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

	my (@pos0, @prevpos, @pos, @offset);
	my ($Crystal);
	my $ret;
	my $iStep = 0;
	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($iStep, $in, $IniCrystal);
			last if(!$ret or !$Crystal);
print "\n";
print "iStep=$iStep / $nStep\n";
			$xc->WriteXSFFileFromCrystal($XSFFile, $iStep, $Crystal);
		}
		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');
		
		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();
	
	return \%hash;
}

sub SavePotentials
{
	my ($this, $CSVPath, $rmin, $rstep, $nr, $pAtomTypes, $pPotHash) = @_;

print "\n";
print "Save potentials for potential type [$pPotHash->{type}]\n";
#foreach my $key (sort keys %$pPot) {
#	my $p = $pPot->{$key};
#	foreach my $key2 (sort keys %$p) {
#		print "   $key:$key2:$p->{$key2} \n";
#	}
#}

	my $KC = $e / 4.0 / $pi / $e0 * 1.0e10;

	my $csv = JFile->new($CSVPath, 'w') or return 0;
	$csv->print("r(A)");
	for(my $i = 0 ; $i < @$pAtomTypes ; $i++) {
		my $atom1 = $pAtomTypes->[$i];
		my $name1 = $atom1->AtomNameOnly();
		for(my $j = $i ; $j < @$pAtomTypes ; $j++) {
			my $atom2 = $pAtomTypes->[$j];
			my $name2 = $atom2->AtomNameOnly();
			$csv->print(",$name1-$name2");
		}
	}
	$csv->print("\n");

	for(my $l = 0 ; $l < $nr ; $l++) {
		my $r = $rmin + $l * $rstep;
		$csv->print("$r");
		for(my $i = 0 ; $i < @$pAtomTypes ; $i++) {
			my $atom1 = $pAtomTypes->[$i];
			my $name1 = $atom1->AtomNameOnly();
			my $Z1    = $pPotHash->{"$name1-core"}{charge};
			for(my $j = $i ; $j < @$pAtomTypes ; $j++) {
				my $atom2 = $pAtomTypes->[$j];
				my $name2 = $atom2->AtomNameOnly();
				my $Z2    = $pPotHash->{"$name2-core"}{charge};

				if($pPotHash->{type} eq 'buckingham') {
					my $Aij   = $pPotHash->{"$name1-core-$name2-core"}{A};
					my $bij   = $pPotHash->{"$name1-core-$name2-core"}{rho};
					my $Cij   = $pPotHash->{"$name1-core-$name2-core"}{C};
#print "$name1-$name2: $Z1, $Z2, $Aij, $bij, $Cij\n";
					my $UC = $KC * $Z1 * $Z2 / $r;
					my $UA = $Aij * exp(-$r / $bij);
					my $UD = $Cij / $r**6;
					$csv->print(",", $UC + $UA + $UD);
				}
				elsif($pPotHash->{type} eq 'morse') {
					my $De = $pPotHash->{"$name1-core-$name2-core"}{De};
					my $a  = $pPotHash->{"$name1-core-$name2-core"}{a};
					my $r0 = $pPotHash->{"$name1-core-$name2-core"}{r0};
					my $fC = $pPotHash->{"$name1-core-$name2-core"}{fsCoulomb};
#print "$name1-$name2: $Z1, $Z2, $Aij, $bij, $Cij\n";
					my $UC = (1.0 - $fC) * $KC * $Z1 * $Z2 / $r;
					my $v  = 1.0 - exp(-$a * ($r - $r0));
					my $UM = $De * ($v*$v - 1.0);
					$csv->print(",", $UC + $UM);
				}
			}
		}
		$csv->print("\n");
	}
	$csv->Close();
#exit;

	return 1;
}

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

	my @AtomTypeList = $Crystal->GetCAtomTypeList();

	my $in = JFile->new($OutFile, 'r') or return undef;
	my $line = $in->SkipTo("Observable\\s+no");
	return undef if(!defined $line);
	$line = $in->SkipTo("General interatomic potentials");
	$line = $in->SkipTo("Atom\\s+Types");
	$line = $in->SkipTo("---");
	my (%BParams, %Charge);
	my $IsEnd = 0;
	print "New potential parameters:\n";
	while(1) {
		while(1) {
			$line = $in->ReadLine();
#print "l:$line";
			last if($line =~ /---/);
			my ($name1, $type1, $name2, $type2, $pottype, $pottype1, $para, $val, $unit) 
				= 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";
			my $A = $val;
print "  A for $name1-$name2: $A\n";

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

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

			$BParams{"$name1-$name2"} = $BParams{"$name2-$name1"} = {
				Aij => $A,
				bij => $rho,
				Cij => $C,
				}
		}
		last if($IsEnd);
	}
	$in->Close();
	
	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;
	while(1) {
		$line = $in->ReadLine();
		last if(!defined $line);

		$out->print($line);
		if($line =~ /spec/i) {
			last;
		}
	}
	for(my $i = 0 ; $i < @AtomTypeList ; $i++) {
		my ($name, $charge) = Utils::Split("\\s+", $in->ReadLine());
		$Charge{$name} = $charge;
print "  Charge for $name: $charge\n";
	}

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

		$out->print($line);
		if($line =~ /buck/i or $line =~ /morse/i) {
			last;
		}
	}
	while(1) {
		$line = $in->ReadLine();
		last if(!defined $line);
		
		my ($name1, $type1, $name2, $type2, $A, $rho, $C,$rmin, $rmax, $idA, $idrho, $idC)
			= Utils::Split("\\s+", $line);
		$A   = $BParams{"$name1-$name2"}{Aij};
		$rho = $BParams{"$name1-$name2"}{bij};
		$C   = $BParams{"$name1-$name2"}{Cij};
#print "$name1-$name2: A=$A  rho=$rho  C=$C\n";

		if(!defined $C) {
			$out->print($line);
			last;
		}
		$out->printf("%3s %4s %3s %4s %12.6g %12.6g % 12.6g  %3.1f %4.1f    %d %d %d\n",
				$name1, $type1, $name2, $type2, $A, $rho, $C, $rmin, $rmax, $idA, $idrho, $idC);
	}
	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-new.csv', $rmin, $rstep, $nr, \@AtomTypeList, \%BParams, \%Charge);

	return 1;
}

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

}

sub UpdateInputForMD
{
	my ($this, $InpFile, $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};
	
	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");

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

	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 $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);

		if($Function =~ /Optimize/i) {
			if($line =~ /^(\s*)#(\s*opti.*)$/) {
				$out->print("$1$2\n");
			}
			elsif($line =~ /^(\s*)(#?)(\s*(md|fit).*)$/) {
				$out->print("$1#$3\n");
			}
			elsif($line =~ /^(\s*)(#?)(\s*ensemble.*)$/) {
				$out->print("$1#$3\n");
			}
		}
		elsif($Function =~ /Fit/i) {
			if($line =~ /^(\s*)#(\s*fit.*)$/) {
				$out->print("$1$2\n");
			}
			elsif($line =~ /^(\s*)(#?)(\s*(md|opti).*)$/) {
				$out->print("$1#$3\n");
			}
			elsif($line =~ /^(\s*)(#?)(\s*ensemble.*)$/) {
				$out->print("$1#$3\n");
			}
		}
		elsif($Function =~ /MD-NVT/i) {
			if($line =~ /^(\s*)#(\s*md.+conv.*)$/) {
				$out->print("$1$2\n");
			}
			elsif($line =~ /^(\s*)(#?)(\s*(md|fit|opti).*)$/) {
				$out->print("$1#$3\n");
			}
			elsif($line =~ /^(\s*)(#?)(\s*ensemble\s+nvt.*)$/) {
				$out->print("$1$3\n");
			}
			elsif($line =~ /^(\s*)(#?)(\s*ensemble\s+npt.*)$/) {
				$out->print("$1#$3\n");
			}
		}
		elsif($Function =~ /MD-NPT/i) {
			if($line =~ /^(\s*)#(\s*md.+conp.*)$/) {
				$out->print("$1$2\n");
			}
			elsif($line =~ /^(\s*)(#?)(\s*(md|fit|opti).*)$/) {
				$out->print("$1#$3\n");
			}
			elsif($line =~ /^(\s*)(#?)(\s*ensemble\s+npt.*)$/) {
				$out->print("$1$3\n");
			}
			elsif($line =~ /^(\s*)(#?)(\s*ensemble\s+nvt.*)$/) {
				$out->print("$1#$3\n");
			}
		}

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

			if($line =~ /^(\s*)(#?)(temp\S*\s+)(\d.*)$/i) {
				my $header = "$1$3";
#				my $s = $4;
#				my ($Tstart, $Tstep0, $nTstep0) = Utils::Split("\\s+", $s);
				$out->print("$header$T0 $Tstep $nstep\n");
			}
			elsif($line =~ /^(\s*)(#?)(equi.*)$/i) {
				$out->print("$1$3\n");
			}
			elsif($line =~ /^(\s*)(#?)(prod.*)(\S*)(\s*)(\S+\s*)?$/i) {
				$out->print("$1$3$production$5$6\n");
			}
			elsif($line =~ /^(\s*)(#?)(timestep.*)(\S*)(\s*)(\S+\s*)?$/i) {
				$out->print("$1$3$timestep$5$6\n");
			}
			elsif($line =~ /^(\s*)(#?)(tsca.*)(\S*)(\s*)(\S+\s*)?$/i) {
				$out->print("$1$3$tscale$5$6\n");
			}
			elsif($line =~ /^(\s*)(#?)(samp.*)(\S*)(\s*)(\S+\s*)?$/i) {
				$out->print("$1$3$sample$5$6\n");
			}
			elsif($line =~ /^(\s*)(#?)(write.*)(\S*)(\s*)(\S+\s*)?$/i) {
				$out->print("$1$3$write$5$6\n");
			}
			elsif($line =~ /^(\s*)(#?)(dump.*)$/i) {
				$out->print("$1$3\n");
			}
		}
		else {
			if($line =~ /^(\s*)(#?)(pres\S*\s+)(\d\S+)/i) {
				$out->print("$1$pressure\n");
			}
			elsif($line =~ /^(\s*)(#?)(equi.*)$/i) {
				$out->print("$1#$3\n");
			}
			elsif($line =~ /^(\s*)(#?)(prod.*)(\S*)(\s*)(\S+\s*)?$/i) {
				$out->print("$1#$3$production$5$6\n");
			}
			elsif($line =~ /^(\s*)(#?)(timestep.*)(\S*)(\s*)(\S+\s*)?$/i) {
				$out->print("$1#$3$timestep$5$6\n");
			}
			elsif($line =~ /^(\s*)(#?)(tsca.*)(\S*)(\s*)(\S+\s*)?$/i) {
				$out->print("$1#$3$tscale$5$6\n");
			}
			elsif($line =~ /^(\s*)(#?)(samp.*)(\S*)(\s*)(\S+\s*)?$/i) {
				$out->print("$1#$3$sample$5$6\n");
			}
			elsif($line =~ /^(\s*)(#?)(write.*)(\S*)(\s*)(\S+\s*)?$/i) {
				$out->print("$1#$3$write$5$6\n");
			}
			elsif($line =~ /^(\s*)(#?)(dump.*)$/i) {
				$out->print("$1#$3\n");
			}
		}
	}

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

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

#Function: "Fit", "Opti", "MD-NVT", "MD-NPT"
sub SaveGULPInputFile
{
	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->PotentialLibrary($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);

#ファイル作製開始
	unless(open(OUT,">$filename")) {
		print "Can not write to [$filename].\n\n";
		return;
	}

	print OUT "#Tasks to be calculated\n";
	print OUT "#  see Kywords in help\n";
	print OUT "#  e.g.: opti md fit\n";
	print OUT "#        prop(elastic,dielectric,piezoe,wave v,bulk/shear/Yongs modulus)\n";
	print OUT "#        phonon raman\n";
	if($Function =~ /Fit/i) {
		print OUT "#opti prop conp\n";
		print OUT "#md conv\n";
		print OUT "#md conp\n";
		print OUT "fit prop conp\n";
	}
	elsif($Function =~ /MD-NVT/i) {
		print OUT "#opti prop conp\n";
		print OUT "md conv\n";
		print OUT "#md conp\n";
		print OUT "#fit prop conp\n";
	}
	elsif($Function =~ /MD-NPT/i) {
		print OUT "#opti prop conp\n";
		print OUT "#md conv\n";
		print OUT "md conp\n";
		print OUT "#fit prop conp\n";
	}
	else {
		print OUT "opti prop conp\n";
		print OUT "#md conv\n";
		print OUT "#md conp\n";
		print OUT "#fit prop conp\n";
	}
	print OUT "\n";

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

	print OUT "#Crystal structure\n";
	print OUT "name $SampleName\n";
	my ($a,$b,$c,$alpha,$beta,$gamma) 
		= $Crystal->LatticeParametersByOutputMode(0);
	print OUT "cell\n";
	printf OUT "%10.6f %10.6f %10.6f %12.6f %12.6f %12.6f\n",
		$a, $b, $c, $alpha, $beta, $gamma;
	print OUT "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);

		printf OUT "%2s core %12.8f  %12.8f  %12.8f  %9.6f %9.6f %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);

		printf OUT "%2s core %12.8f  %12.8f  %12.8f  %9.6f %9.6f %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();
	print OUT "space\n";
	print OUT "$iSPG\n";
	print OUT "\n";

#	print OUT "library bush\n";

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

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

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

	if($pPot->{type} eq 'buckingham') {
		print OUT "buckingham\n";
		for(my $i = 0 ; $i < @AtomTypeList ; $i++) {
			my $type1 = $AtomTypeList[$i];
			my $name1 = $type1->AtomNameOnly();
			for(my $j = $i ; $j < @AtomTypeList ; $j++) {
				my $type2 = $AtomTypeList[$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);
print "  $name1 - $name2: $Aij  $bij  $Cij\n";
				printf OUT "%3s %4s %3s %4s %12.6g %12.6g %12.6g  0.0 10.0    1 0 0\n",
					$name1, 'core', $name2, 'core', $Aij, $bij, $Cij;
			}
		}
	}
	elsif($pPot->{type} eq 'morse') {
		print OUT "morse\n";
		for(my $i = 0 ; $i < @AtomTypeList ; $i++) {
			my $type1 = $AtomTypeList[$i];
			my $name1 = $type1->AtomNameOnly();
			for(my $j = $i ; $j < @AtomTypeList ; $j++) {
				my $type2 = $AtomTypeList[$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 $a         = $p->{a};
				my $r0        = $p->{r0};
				my $fsCoulomb = $p->{fsCoulomb};
print "  $name1 - $name2: $De  $a  $r0  $fsCoulomb\n";
				printf OUT "%3s %4s %3s %4s %12.6g %12.6g %12.6g %12.6g  0.0 10.0    1 0 0\n",
					$name1, 'core', $name2, 'core', $De, $a, $r0, $fsCoulomb;
			}
		}
	}

#	if($nShellAtoms > 0) {
#		print OUT "spring\n";
#		for(my $i = 0 ; $i < @Spring ; $i++) {
#			print OUT "$Spring[$i]\n";
#		}
#	}
	print OUT "\n";

	print OUT "#Use 1 Angstrom taper on the short range potential for smoothness\n";
	print OUT "#cutp 12.0 1.0\n";
	print OUT "cutp 12.0 1.0\n";
	print OUT "#pressure P <GPa/kPa/MPa/Pa/atm/Nm-2/kbar (unit, Def=GPa)>\n";
	print OUT "pressure    $args{pressure}\n";
	print OUT "\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) {
		print OUT "#MD options\n";
		print OUT "supercell 1 1 1\n";
		print OUT "shellmass 0.1\n";
		print OUT "#integrator <gear/velocity verlet/leapfrog verlet/stochastic> <iter>\n";
		print OUT "integrator leapfrog verlet\n";
		print OUT "#ensemble <NVE/NVT qnose/NPT qnose qpress/NPH>\n";
		if($Function =~ /NVT/i) {
			print OUT "ensemble nvt 0.1\n";
			print OUT "#ensemble npt 0.05 0.10\n";
		}
		else {
			print OUT "#ensemble nvt 0.1\n";
			print OUT "ensemble npt 0.05 0.10\n";
		}
		print OUT "#temperature T <C/F/K (unit, Def=K)> <Tstep> <# of steps> <1st step of MD>\n";
		print OUT "temperature $T0 $Tstep $nstep\n";
		print OUT "#Simulation time for equilibration prior to MD production phase\n";
		print OUT "equilibration $args{equilibration} ps\n";
		print OUT "#Temperature scaling time: Def = equilibration\n";
		print OUT "tscale        $args{tscale} ps\n";
		print OUT "#Simulation time to collect production data\n";
		print OUT "production    $args{production} ps\n";
		print OUT "timestep      $args{timestep} ps\n";
		print OUT "sample        $args{sample} ps\n";
		print OUT "write         $args{write} ps\n";
		print OUT "dump every        $filenames[0].res\n";
		print OUT "output trajectory $filenames[0].trg\n";
		print OUT "output history    $filenames[0].his\n";
	}
	else {
		print OUT "#MD options\n";
		print OUT "#supercell 1 1 1\n";
		print OUT "#shellmass 0.1\n";
		print OUT "#integrator <gear/velocity verlet/leapfrog verlet/stochastic> <iter>\n";
		print OUT "#integrator leapfrog verlet\n";
		print OUT "#ensemble <NVE/NVT qnose/NPT qnose qpress/NPH>\n";
		print OUT "#ensemble nvt 0.1\n";
		print OUT "#ensemble npt 0.05 0.10\n";
		print OUT "#temperature T <C/F/K (unit, Def=K)> <step> <# of steps> <1st step of MD>\n";
		print OUT "#temperature $T0 $Tstep $nstep\n";
		print OUT "#Simulation time for equilibration prior to MD production phase\n";
		print OUT "#equilibration $args{equilibration} ps\n";
		print OUT "#Temperature scaling time: Def = equilibration\n";
		print OUT "#tscale        $args{tscale} ps\n";
		print OUT "#Simulation time to collect production data\n";
		print OUT "#production    $args{production} ps\n";
		print OUT "#timestep      $args{timestep} ps\n";
		print OUT "#sample        $args{sample} ps\n";
		print OUT "#write         $args{write} ps\n";
		print OUT "#dump every        $filenames[0].res\n";
		print OUT "#output trajectory $filenames[0].trg\n";
		print OUT "#output history    $filenames[0].his\n";
	}

#	print OUT "output xr $filenames[0]\n";
#	print OUT "output marvin $filenames[0].mvn\n";

	print OUT "\n\n";
	
	close OUT;

	return 1;
}


1;