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

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

use strict;
use File::Basename;

use ProgVars;
use Crystal::Crystal;
use Crystal::SpaceGroup;
use Crystal::Rietan;
use Sci::Science;
use Sci::EnergyBandArray;

#my $DirectorySeparator = "\\";

#my $WebRootDir      = "d:\\MyWebs";
#my $GULPDir         = "d:\\Programs\\gulp";
#my $GULPLibDir      = Deps::MakePath($GULPDir, "Libraries",     0);
#my $GULPLibaryPath = Deps::MakePath($GULPLibDir, "Kamiya.lib", 0);
my $ProgramDir      = ProgVars::ProgramDir();
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 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 SetSampleName
{
	my ($this, $name) = @_;
	return $this->{'SampleName'} = $name;
}
sub SampleName
{
	my ($this) = @_;
	return $this->{'SampleName'};
}

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

sub DESTROY
{
	my $this = shift;
}

#===============================================
# 一般関数
#===============================================
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 ReadSpecies
{
	my ($this, $atom, $cstype) = @_;

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

	my $SpeciesCore;
	my $SpeciesShell;
	while(<IN>) {
		last if($_ =~ /^\s*species/i);
	}
#Read species (ionic charges)
	while(<IN>) {
		my $line = $_;
		chomp($line);
		$SpeciesCore  = $line if($line =~ /^$atom\s+core\s/i);
		$SpeciesShell = $line if($line =~ /^$atom\s+shel\s/i);
		last if($line =~ /^\s*buck/i);
	}
	close(IN);

	if($cstype eq 'core') {
		my ($name1, $cs1, $charge1) = ($SpeciesCore  =~ /^(\S*)\s+(\S*)\s+(\S*)/);
		my ($name2, $cs2, $charge2) = ($SpeciesShell =~ /^(\S*)\s+(\S*)\s+(\S*)/);

		$SpeciesCore = sprintf "%2s %s %8.5f", 
				$name1, 'core', $charge1+$charge2;
		$SpeciesShell = '';
	}

	return ($SpeciesCore, $SpeciesShell);
}

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);
		$Spring = $line if($line =~ /^$atom\s/i);
	}
	close(IN);
	unless($Spring =~ /^(\S+)\s+(\S+)\s+(\S+)/) {
		$Spring = "$Spring 0" if($Spring ne '');
	}

	return $Spring;
}

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

#print "atoms: $atom1, $atom2<br>";
#print "read [$GULPLibaryPath]<br>";

	my ($SpeciesCore1, $SpeciesShell1) = $this->ReadSpecies($atom1, $cstype1);
	my ($SpeciesCore2, $SpeciesShell2) = $this->ReadSpecies($atom2, $cstype2);
	my $Spring1 = $this->ReadSpring($atom1, $cstype1);
	my $Spring2 = $this->ReadSpring($atom2, $cstype2);

	my $GULPLibaryPath = $this->LibraryPath();
	if(!open(IN, "<$GULPLibaryPath")) {
		print "Error in GULP::ReadBuckinghamPotential: Can not read [$GULPLibaryPath]\n";
		return ();
	}
	while(<IN>) {
		last if($_ =~ /^\s*buck/i);
	}
#Read buckingham two body parameters
	my $Buckingham;

	while(<IN>) {
		my $line = $_;
		chomp($line);
		if($line =~ /^$atom1\s+(\S+)\s+$atom2\s/i or 
		   $line =~ /^$atom2\s+(\S+)\s+$atom1\s/i) {
			$Buckingham = $line;
			last;
		}
		last if($line =~ /^\s*spring/i);
	}

	close(IN);

	if($cstype1 eq 'core' or $cstype1 eq 'undefined') {
		$Buckingham =~ s/$atom1\s+shell?/$atom1 core/gi;
	}
	if($cstype2 eq 'core' or $cstype2 eq 'undefined') {
		$Buckingham =~ s/$atom2\s+shell?/$atom2 core/gi;
	}
	unless($Buckingham =~ /^(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+/) {
		$Buckingham = "$Buckingham 0 0 0"
			if($Buckingham ne '');
	}

	return ($SpeciesCore1, $SpeciesShell1, 
		$SpeciesCore2, $SpeciesShell2,
		$Buckingham,
		$Spring1, $Spring2);
}

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);
		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 "   $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;
print "  type=$type\n" if($IsPrint);
		Utils::DelSpace($type);

		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 =~ /^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 =~ /^buckingham/i or $type =~ /^buck/i) {
			$this->AddAllPotentialsToHash($in, $type, 2, \%array);
		}
		elsif($type =~ /^morse/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 =~ /^buck/i) {
			$this->AddAllPotentialsToHash($in, $type, 2, \%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 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 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'};
}

#Function: "Fit", "Opti", "MD-NVT", "MD-NPT"
sub SaveGULPInputFile
{
	my ($this, $Crystal, $Function, $filename, $IsChooseRandomly) = @_;
	my $LF = "<br>\n";

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

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

	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>$LF";
	my @list = keys %CoreShellType;
	for(my $i = 0 ; $i < @list ; $i++) {
		my $key = $list[$i];
		my $val = $CoreShellType{$key};
		print "  $key: $val$LF";
	}

#ファイル作製開始
	unless(open(OUT,">$filename")) {
		print "Can not write to [$filename].$LF$LF";
		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";

	my @ExpandedAtomSiteList = $Crystal->GetCExpandedAtomSiteListByOutputMode();

#Occupancyが1のサイト
	for(my $i = 0 ; $i < @ExpandedAtomSiteList ; $i++) {
		my $atom      = $ExpandedAtomSiteList[$i];
		my $atomname  = $atom->AtomNameOnly();
		my $charge    = $atom->Charge();
		my ($x,$y,$z) = $atom->Position(1);
		my $occupancy = $atom->Occupancy();
		my $mult      = $atom->Multiplicity();
		my $id        = $atom->Id();
		next if($occupancy < 0.9999);

		if($CoreShellType{$atomname} eq 'core-shell') {
			printf OUT "%2s core %12.8f  %12.8f  %12.8f  %9.6f %9.6f %7.4f  %d %d %d\n", 
					$atomname, $x, $y, $z, 0.0,     $occupancy, 0.0, 1,1,1;
			printf OUT "%2s shel %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;
		}
		else {
			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    = $atom->Charge();
		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);

		if($CoreShellType{$atomname} eq 'core-shell') {
			printf OUT "%2s core %12.8f  %12.8f  %12.8f  %9.6f %9.6f %7.4f  %d %d %d\n", 
					$atomname, $x, $y, $z, 0.0,     $occupancy, 0.0, 1,1,1;
			printf OUT "%2s shel %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;
		}
		else {
			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";

	my @Species;
	my @Buckingham;
	my @Spring;
	my %SpeciesCoreList;
	my %SpeciesShellList;
	my %SpringList;
	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 ($SpeciesCore1, $SpeciesShell1, 
			    $SpeciesCore2, $SpeciesShell2,
			    $Buckingham,
			    $Spring1, $Spring2)
			    	= $this->ReadBuckinghamPotential($name1, $name2, 
			    		$CoreShellType{$name1}, $CoreShellType{$name2});
			if($SpeciesCore1  ne '' and $SpeciesCoreList{$name1} == 0) {
				push(@Species, $SpeciesCore1);
				$SpeciesCoreList{$name1}++;
			}
			if($SpeciesShell1  ne '' and $SpeciesShellList{$name1} == 0) {
				push(@Species, $SpeciesShell1);
				$SpeciesShellList{$name1}++;
			}
			if($SpeciesCore2  ne '' and $SpeciesCoreList{$name2} == 0) {
				push(@Species, $SpeciesCore2);
				$SpeciesCoreList{$name2}++;
			}
			if($SpeciesShell2  ne '' and $SpeciesCoreList{$name2} == 0) {
				push(@Species, $SpeciesShell2);
				$SpeciesShellList{$name2}++;
			}
			push(@Buckingham, $Buckingham) if($Buckingham ne '');
			if($Spring1  ne '' and $SpringList{$name1} == 0) {
				push(@Spring, $Spring1);
				$SpringList{$name1}++;
			}
			if($Spring2  ne '' and $SpringList{$name2} == 0) {
				push(@Spring, $Spring2);
				$SpringList{$name2}++;
			}
		}
	}
	print OUT "species ", scalar @Species, "\n";
	for(my $i = 0 ; $i < @Species ; $i++) {
		print OUT "$Species[$i]\n";
	}
	print OUT "buck\n";
	for(my $i = 0 ; $i < @Buckingham ; $i++) {
		print OUT "$Buckingham[$i]\n";
	}
	if($nShellAtoms > 0) {
		print OUT "spring\n";
		for(my $i = 0 ; $i < @Spring ; $i++) {
			print OUT "$Spring[$i]\n";
		}
	}
	print OUT "\n";

	if($Function =~ /MD/i) {
		print OUT "#MD options\n";
		print OUT "supercell 1 1 1\n";
		print OUT "shellmass 0.1\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 "#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)> <step> <# of steps> <1st step of MD>\n";
		print OUT "temperature 300 0.1 1000\n";
		print OUT "#pressure P <GPa/kPa/MPa/Pa/atm/Nm-2/kbar (unit, Def=GPa)>\n";
		print OUT "pressure    0.1\n";
		print OUT "#Simulation time for equilibration prior to MD production phase\n";
		print OUT "equilibration 0.100 ps\n";
		print OUT "#Temperature scaling time: Def = equilibration\n";
		print OUT "tscale        0.100 ps\n";
		print OUT "#Simulation time to collect production data\n";
		print OUT "production    1.000 ps\n";
		print OUT "timestep      0.001 ps\n";
		print OUT "sample        0.001 ps\n";
		print OUT "write         0.005 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 "#Use 1 Angstrom taper on the short range potential for smoothness\n";
		print OUT "#cutp 12.0 1.0\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 300 0.1 1000\n";
		print OUT "#pressure P <GPa/kPa/MPa/Pa/atm/Nm-2/kbar (unit, Def=GPa)>\n";
		print OUT "#pressure      0.1\n";
		print OUT "#Simulation time for equilibration prior to MD production phase\n";
		print OUT "#equilibration 0.10 ps\n";
		print OUT "#Temperature scaling time: Def = equilibration\n";
		print OUT "#tscale        0.10 ps\n";
		print OUT "#Simulation time to collect production data\n";
		print OUT "#production    1.00 ps\n";
		print OUT "#timestep      0.001 ps\n";
		print OUT "#sample        0.001 ps\n";
		print OUT "#write         0.005 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;