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

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

use strict;
use File::Basename;

use ProgVars;
use Crystal::XCrySDen;
use Sci qw($torad $e0 $a0 $e $pi log10 erfc);
use Sci::Science;

#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 $GULPDir         = ProgVars::GULPDir();
my $GULPLibDir      = ProgVars::GULPLibDir();

#===============================================
# コンストラクタ、デストラクタ
#===============================================
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'} = LAMMPS::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 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 LAMMPS::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 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 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 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 ReadNextCrystalFromCustomDump
{
	my ($this, $idx, $in, $cry) = @_;

	my $k = 1.0;
	my @AtomTypeList= $cry->GetCAtomTypeList();

	my $Crystal = new Crystal();
#	$Crystal->SetCrystalName($Title);
#	$Crystal->SetSampleName($Title);

	for(my $i = 0 ; $i < @AtomTypeList ; $i++) {
		$Crystal->AddAtomType($AtomTypeList[$i]->AtomNameOnly());
	}

	my ($TimeStep, $nAtoms);
	my $IsRead = 0;
	while(1) {
		my $line = $in->ReadLine();
		if(!defined $line) {
			if(!defined $nAtoms or !$IsRead) {
				return (undef, 0);
			}
		}
#print "l:$line";

		if($line =~ /ITEM:\s*(\S.*)\s*$/) {
			my $key = $1;
print "key:$key\n";
			if($key eq 'TIMESTEP') {
				$TimeStep = $in->ReadLine() + 0.0;
print "  Time step: $TimeStep\n";
			}
			elsif($key eq 'NUMBER OF ATOMS') {
				$nAtoms = $in->ReadLine() + 0;
print "  nAtoms: $nAtoms\n";
			}
			elsif($key =~ /^BOX BOUNDS/) {
				my ($xlo, $xhi, $xy) = Utils::Split("\\s+", $in->ReadLine());
				my ($ylo, $yhi, $xz) = Utils::Split("\\s+", $in->ReadLine());
				my ($zlo, $zhi, $yz) = Utils::Split("\\s+", $in->ReadLine());
# Vectors ($xhi-$xlo, 0, 0)  ($xy, $yhi-$ylo, 0)  ($xz, $yz, $zhi-$zlo)
				$Crystal->SetLatticeVector($xhi-$xlo, 0.0, 0.0, $xy, $yhi-$ylo, 0.0, $xz, $yz, $zhi-$zlo);
			}
			elsif($key =~ /^ATOMS\s*(\S.*)\s*$/) {
				$IsRead = 1;
				my @list = Utils::Split("\\s+", $1);
print "  list: ", join(', ', @list), "\n";
				my %idx;
				for(my $i = 0 ; $i < @list ; $i++) {
					$idx{$list[$i]} = $i;
				}
				my %AtomCount;
				for(my $i = 0 ; $i < $nAtoms ; $i++) {
#					my ($id, $iType, $x, $y, $z, $ix, $iy, $iz) = Utils::Split("\\s+", $in->ReadLine());
#					my ($id, $iType, $xr, $yr, $zr, $ix, $iy, $iz) = Utils::Split("\\s+", $in->ReadLine());
					my @a = Utils::Split("\\s+", $in->ReadLine());
					my $id    = $a[$idx{id}];
					my $iType = $a[$idx{type}];
					my $x     = $a[$idx{xs}];
					my $y     = $a[$idx{ys}];
					my $z     = $a[$idx{zs}];

					my $atomname = $AtomTypeList[$iType-1]->AtomNameOnly();
					$AtomCount{$atomname}++;
					my $Label = $atomname . $AtomCount{$atomname};
#					my ($x, $y, $z) = $Crystal->CalFractionalCoordinate($k * $xr, $k * $yr, $k * $zr);
#printf " $i: $atomname: $Label (%8.3g, %8.3g, %8.3g) => (%8.3g, %8.3g, %8.3g)\n", $xr, $yr, $zr, $x, $y, $z;
printf " $i: $atomname: $Label (%12.4g, %12.4g, %12.4g)\n", $x, $y, $z;
					$Crystal->AddAtomSite($Label, $atomname, $x, $y, $z, 1.0);
				}
				last;
			}
			else {
				print "\nError: Invalid key in ITEM {$key]\n";
			}
		}
		else {
			;
		}
	}
	$Crystal->{TimeStep} = $TimeStep;
	$Crystal->{nAtom} = $nAtoms;
	
	$Crystal->ExpandCoordinates();
	my @AtomList = $Crystal->GetCExpandedAtomSiteList();
print "  nAtoms = ", scalar @AtomList, "\n";

	return ($Crystal, 1);
}

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

sub MakeXSFFromCustomDump
{
	my ($this, $pStep, $pTime, $pTemp, $DumpFile, $XSFFile, $IniCrystal, $HistoryFile) = @_;

	print("\n\n<b>Make XSF File from LAMMPS custom dump file:</b>\n");
	print("  Read from [$DumpFile].\n");
	print("  Save to   [$XSFFile].\n");

	my $nStep = $this->GetnStepInCustomDump($DumpFile);
print "nStep: $nStep\n";

	my $in = JFile->new($DumpFile, 'r') or die "Can not read [$DumpFile].\n";
	my $csv;
	if($HistoryFile) {
		$csv = JFile->new($HistoryFile, 'w') or die "Can not write to [$HistoryFile].\n";
		$csv->print("iStep,TimeStep,Temp,a,b,c,alpha,beta,gamma,Volume,Density\n");
#		$csv->print("iStep,TimeStep,Temp,a,b,c,alpha,beta,gamma,Volume,Density,msd\n");
	}

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

	my (@pos0, @prevpos, @pos, @offset);
	my ($Crystal, $FinalCrystal);
	my $ret;
	my $iStep = 0;
	my $msd = 0.0;
	while(1) {
#$site->SetVelocity($dx, $dy, $dz);
#$site->SetForce(0, 0, $tot);
		($Crystal, $ret) = $this->ReadNextCrystalFromCustomDump($iStep, $in, $IniCrystal);
		last if(!$ret or !$Crystal);

		$iStep++;
print "iStep=$iStep\n";
		$FinalCrystal = $Crystal->Duplicate();

		if($csv) {
			my @atoms = $Crystal->GetCExpandedAtomSiteList();
			if($iStep == 1) {
				for(my $i = 0 ; $i < @atoms ; $i++) {
					my ($x, $y, $z) = $atoms[$i]->Position();
					($prevpos[$i][0], $prevpos[$i][1], $prevpos[$i][2])
						= ($pos0[$i][0], $pos0[$i][1], $pos0[$i][2]) = ($x, $y, $z);
				}	
			}
			for(my $i = 0 ; $i < @atoms ; $i++) {
				my ($x, $y, $z) = $atoms[$i]->Position();
# wrapされた座標から、実際の変位量を計算する
# ただし、一単位格子以上の変位があると、変位量は不明になる。msdの計算はあきらめる
				if($prevpos[$i][0] - $x > 0.5) {
					$offset[$i][0] += 1.0;
				}
				elsif($x - $prevpos[$i][0] > 0.5) {
					$offset[$i][0] -= 1.0;
				}

				if($prevpos[$i][1] - $y > 0.5) {
					$offset[$i][1] += 1.0;
				}
				elsif($y - $prevpos[$i][1] > 0.5) {
					$offset[$i][1] -= 1.0;
				}

				if($prevpos[$i][2] - $z > 0.5) {
					$offset[$i][2] += 1.0;
				}
				elsif($z - $prevpos[$i][2] > 0.5) {
					$offset[$i][2] -= 1.0;
				}

				my $r = $Crystal->GetInterAtomicDistance(
						$pos0[$i][0], $pos0[$i][1], $pos0[$i][2], 
						$x+$offset[$i][0], $y+$offset[$i][1], $z+$offset[$i][2]);
				$msd += $r*$r;
#print "$i: r = $r ($pos0[$i][0], $pos0[$i][1], $pos0[$i][2]) - ($x+$offset[$i][0], $y+$offset[$i][1], $z+$offset[$i][2])\n";
				($prevpos[$i][0], $prevpos[$i][1], $prevpos[$i][2]) = ($x, $y, $z);
			}
#exit if($iStep >= 3);
			my ($a,$b,$c,$alpha,$beta,$gamma) = $Crystal->LatticeParameters();
			my $V = $Crystal->Volume();
			my $d = $Crystal->Density();
			my $temp = Algorism::Interpolate($pStep, $pTemp, $Crystal->{TimeStep});
			$csv->printf("$iStep,$Crystal->{TimeStep},$temp,$a,$b,$c,$alpha,$beta,$gamma,$V,$d\n");
#			$csv->printf("$iStep,$Crystal->{TimeStep},$temp,$a,$b,$c,$alpha,$beta,$gamma,$V,$d\n", 
#				$msd / $Crystal->{nAtom} / $iStep);
		}
		
		$xc->WriteXSFFileFromCrystal($XSFFile, $iStep, $Crystal);
	}

	$csv->Close() if($csv);
	
	print("\nMake XSF File from VASP *CAR Files: Finished\n");

	return $FinalCrystal;
}

sub SaveHistoryFromLog
{
	my ($this, $LogFile, $LogHistoryFile) = @_;

print "\n";
print "Save log hisotry to [$LogHistoryFile] from [$LogFile].\n";

	my @array;
	my $in  = JFile->new($LogFile, 'r') or return undef;
	my $out = JFile->new($LogHistoryFile, 'w') or return undef;
	
	my $iRun = 0;
	my $tc = 0;
	my $PrevLabelLine = '';
	while(1) {
		my $line = $in->SkipTo("^Step ");
		if($line) {
			$array[$iRun] = {};
		}
		else {
			last;
		}
		my @label = Utils::Split("\\s+", $line);
		last if($label[0] ne 'Step');
print "iRun: $iRun: Labels: ", join(', ', @label), "\n";

		if($PrevLabelLine ne $line) {
			$out->print(join(',', @label), "\n");
			$PrevLabelLine = $line;
		}
		$tc++;

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

			my @a = Utils::Split("\\s+", $line);
			last if($a[0] !~ /^[+\-\d\.eEdD]+$/);
			
			$array[$iRun]{iRun}[$c] = $iRun+1;
			for(my $i = 0 ; $i < @a ; $i++) {
				$array[$iRun]{$label[$i]}[$c] = $a[$i];
			}
			$out->print(join(',', @a), "\n");
			$c++;
			$tc++;
		}
		last if(!$line);
		$iRun++;
	}
	
	$in->Close();
	$out->Close();
print "\n";

	return \@array;
}

sub ReadLogFileToHashArray
{
	my ($this, $LogFile) = @_;

print "\n";
print "Read log file from [$LogFile].\n";

	my @array;
	my $in = JFile->new($LogFile, 'r') or return undef;
	
	my $iRun = 0;
	while(1) {
		my $line = $in->SkipTo("^Step ");
		if($line) {
			$array[$iRun] = {};
		}
		else {
			last;
		}
		my @label = Utils::Split("\\s+", $line);
		last if($label[0] ne 'Step');
print "iRun: $iRun: Labels: ", join(', ', @label), "\n";

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

			my @a = Utils::Split("\\s+", $line);
			last if($a[0] !~ /^[+\-\d\.eEdD]+$/);
			
			$array[$iRun]{iRun}[$c] = $iRun+1;
			for(my $i = 0 ; $i < @a ; $i++) {
				$array[$iRun]{$label[$i]}[$c] = $a[$i];
			}
			$c++;
		}
		last if(!$line);
		$iRun++;
	}
	
	$in->Close();
print "\n";

	return \@array;
}


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

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

	my $Crystal = new Crystal();
#	$Crystal->SetCrystalName($Title);
#	$Crystal->SetSampleName($Title);

	my $in = new JFile($InFile, "r");
	return 0 unless($in);
	
	my $line;
	while(1) {
		$line = $in->ReadLine();
		last if(!defined $line);
		last if($line =~ /^\s*pair_style\s/);
	}
	$line = $in->ReadLine();
	$line =~ s/^\s*#//;
	my @atoms = Utils::Split("\\s+", $line);	
print "atoms: ", join(', ', @atoms), "\n";
	for(my $i = 0 ; $i < @atoms ; $i++) {
		$Crystal->AddAtomType($atoms[$i]);
	}
	$in->Close();

	$in = new JFile($DataFile, "r");
	return 0 unless($in);
	
	my ($xhi, $yhi, $zhi, $xy, $xz, $yz);
	while(1) {
		$line = $in->ReadLine();
		last if(!defined $line);
		last if($line =~ /^\s*(Masses|Atoms)/i);

		my @a = Utils::Split("\\s+", $line);
		if($line =~ /xlo/) {
			$xhi = $k * ($a[1] - $a[0]);
		}
		elsif($line =~ /ylo/) {
			$yhi = $k * ($a[1] - $a[0]);
		}
		elsif($line =~ /zlo/) {
			$zhi = $k * ($a[1] - $a[0]);
		}
		elsif($line =~ /xy/) {
			$xy = $k * $a[0];
			$xz = $k * $a[1];
			$yz = $k * $a[1];
		}
	}
print "Lattice vector: ($xhi, 0.0, 0.0)  ($xy, $yhi, 0.0)  ($xz, $yz, $zhi)\n";
	$Crystal->SetLatticeVector($xhi, 0.0, 0.0, $xy, $yhi, 0.0, $xz, $yz, $zhi);

	$in->rewind();
	while(1) {
		$line = $in->ReadLine();
		last if(!defined $line);
		last if($line =~ /^\s*Atoms/);
	}

	my $c = 0;
	my %AtomCount;
	while(1) {
		$line = $in->ReadLine();
		last if(!defined $line);
		
		my ($iAtom, $iType, $Z, $xr, $yr, $zr) = Utils::Split("\\s+", $line);
		next if(!defined $zr);

		$c++;
		my $atomname = $atoms[$iType-1];
		$AtomCount{$atomname}++;
		my $Label = $atoms[$iType-1] . $AtomCount{$atomname};
		my ($x, $y, $z) = $Crystal->CalFractionalCoordinate($k * $xr, $k * $yr, $k * $zr);
printf " $c: $atomname: $Label (%8.3g, %8.3g, %8.3g) => (%8.3g, %8.3g, %8.3g)\n", $xr, $yr, $zr, $x, $y, $z;
		$Crystal->AddAtomSite($Label, $atomname, $x, $y, $z, 1.0);
	}
	$Crystal->ExpandCoordinates();
	my @AtomList = $Crystal->GetCExpandedAtomSiteList();
print "  nAtoms = ", scalar @AtomList, "\n";

	$in->Close();

	$Crystal->ExpandCoordinates();
print "\n";

	return $Crystal;
}

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

	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    = $pCharge->{$name1};
			for(my $j = $i ; $j < @$pAtomTypes ; $j++) {
				my $atom2 = $pAtomTypes->[$j];
				my $name2 = $atom2->AtomNameOnly();
				my $Z2    = $pCharge->{$name2};
				my $Aij   = $pBParams->{"$name1-$name2"}{Aij};
				my $bij   = $pBParams->{"$name1-$name2"}{bij};
				my $Cij   = $pBParams->{"$name1-$name2"}{Cij};
#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);
			}
		}
		$csv->print("\n");
	}
	$csv->Close();
#exit;

	return 1;
}

sub SaveInputFiles
{
	my ($this, $Crystal, $Function, $InFile, $DataFile, $CopyScript, $ScriptPath, $IsChooseRandomly) = @_;
	$CopyScript = 0 if(!defined $CopyScript);

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

	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 (%BParams, %Charge);
	my ($Aijmax, $bijmax, $Cijmax) = (0.0, 0.0, 0.0);
	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\n";
			my ($Zi, $Zj, $Aij, $bij, $Cij) = $this->ReadBuckinghamPotential(
				$name1, $name2, $CoreShellType{$name1}, $CoreShellType{$name2});
			$Charge{$name1} = $Zi;
			$Charge{$name2} = $Zj;
#print "i,j: $i,$j: $name1 $name2  Zi=$Zi  Zj=$Zj\n";
			$BParams{"$name1-$name2"} = $BParams{"$name2-$name1"} = {
				Aij => $Aij,
				bij => $bij,
				Cij => $Cij,
				};
print "i,j: $i,$j: $name1 $name2: $Charge{$name1}, $Charge{$name2}, $Aij, $bij, $Cij\n";
			$Aijmax = $Aij if($Aij > $Aijmax);
			$bijmax = $bij if($Aij != 0.0 and $bij > $bijmax);
			$Cijmax = $Cij if($Cij > $Cijmax);
		}
	}
#exit;
	my ($rmin, $rstep, $nr) = (0.01, 0.01, 400);
	$this->SavePotentials('Potentials.csv', $rmin, $rstep, $nr, \@AtomTypeList, \%BParams, \%Charge);

	print "\n";
	my $alpha = 0.4;
	my $EPS   = 1.0e-6;
	print "Estimate calculation range\n";
	print "  Aijmax = $Aijmax\n";
	print "  bijmax = $bijmax\n";
	print "  Cijmax = $Cijmax\n";
	print "  EWALD alpha is assumed to be $alpha\n";
	my $N = int(-log10($EPS) + 1.0e-6);
	print "  for EPS = $EPS = 10^$N \n";

	my $RDmax     = (2.26 + 0.26 * $N)/$alpha;
	my $erfcRDmax = erfc($alpha*$RDmax);
	my $RRmax    = 1.0 + $bijmax * $N * log(10.0);
	my $expRRmax = exp(-$RRmax / $bijmax);
	my $RPmax = $N * log(10.0);
	my $recRPmax6 = 1.0 / $RPmax**6;
	printf "    Direct term    : RDmax = %12.6g A  erfc(alpha*RDmax)  = %10.4g\n", $RDmax, $erfcRDmax;
	printf "    Repulsive term : RRmax = %12.6g A  exp(-RRmax/bijmax) = %10.4g\n", $RRmax, $expRRmax;
	printf "    Dispersion term: RPmax = %12.6g A  1.0/RPmax^6        = %10.4g\n", $RPmax, $recRPmax6;
	
	my $Rmax = $RDmax;
	$Rmax = $RRmax if($Rmax < $RRmax);
	$Rmax = $RDmax if($Rmax < $RDmax);
	my @NRmax;
	$NRmax[0] = int( $Rmax / sqrt($Crystal->{g11} * sin($torad*$beta ) * sin($torad*$gamma)) ) + 1;
	$NRmax[1] = int( $Rmax / sqrt($Crystal->{g22} * sin($torad*$alpha) * sin($torad*$gamma)) ) + 1;
	$NRmax[2] = int( $Rmax / sqrt($Crystal->{g33} * sin($torad*$alpha) * sin($torad*$beta )) ) + 1;
	print  "  Cal range estimated by Rmax = $Rmax A: ($NRmax[0], $NRmax[1], $NRmax[2])\n";
	my $pi2 = 2.0 * $pi;
	my $CalN = (int(2.0/3.0 * $pi2 * $Rmax**3.0 / $V) + 1.0) * $nExpandedAtomSiteList;
	print  "     Number of calculations in direct term: $CalN\n";
	
	my $max  = sqrt( $Crystal->{g11} * sin($torad*$beta ) * sin($torad*$gamma) ) * $NRmax[0];
	my $max2 = sqrt( $Crystal->{g22} * sin($torad*$alpha) * sin($torad*$gamma) ) * $NRmax[1];
	$max = $max2 if($max < $max2);
	   $max2 = sqrt( $Crystal->{g33} * sin($torad*$alpha) * sin($torad*$beta ) ) * $NRmax[2];
	$max = $max2 if($max < $max2);
	my $max3 = sqrt( $Crystal->{g11} * sin($torad*$beta ) * sin($torad*$gamma) ) * ($NRmax[1] - 1);
	   $max2 = sqrt( $Crystal->{g22} * sin($torad*$alpha) * sin($torad*$gamma) ) * ($NRmax[2] - 1);
	$max3 = $max2 if($max3 < $max2);
	   $max2 = sqrt( $Crystal->{g33} * sin($torad*$alpha) * sin($torad*$beta ) ) * ($NRmax[3] - 1);
		$max3 = $max2 if($max3 < $max2);
	print  "     Rmax(Calculation) = $max3 - $max\n";

	my $g2max = 4.0 * $alpha**2 * $N * log(10.0) + 1.0;
	my $expg2max = exp(-$g2max / 4.0 / $alpha**2);
	my ($Ra, $Rb, $Rc, @Rang) = $Crystal->ReciprocalLatticeParameters();
	my @lsin;
	for(my $i = 0 ; $i < 3 ; $i++) {
		$lsin[$i] = sin($torad * $Rang[$i]);
	}
	my @hgmax;
	$hgmax[0] = int( sqrt($g2max / ($Crystal->{Rg11} * $lsin[1] * $lsin[2]) ) / $pi2) + 1;
	$hgmax[1] = int( sqrt($g2max / ($Crystal->{Rg22} * $lsin[0] * $lsin[2]) ) / $pi2) + 1;
	$hgmax[2] = int( sqrt($g2max / ($Crystal->{Rg33} * $lsin[0] * $lsin[1]) ) / $pi2) + 1;
	print  "  Cal range estimated for reciprocal space: "
	      ."G2max = $g2max   exp(-G2max / 4.0 / alpha^2) = $expg2max: hmax = ($hgmax[0], $hgmax[1], $hgmax[2])\n";

	my $CalN = int( 2.0/3.0 * $pi2 * $g2max**1.5 / $RV / $pi2**3.0) + 1;
	print  "     Number of calculations in reciprocal term: $CalN\n";
	$max  = sqrt( $Crystal->{Rg11} * $lsin[1] * $lsin[2] ) * $hgmax[0];
	$max2 = sqrt( $Crystal->{Rg22} * $lsin[0] * $lsin[2] ) * $hgmax[1];
	$max = $max2 if($max < $max2);
	$max2 = sqrt( $Crystal->{Rg33} * $lsin[0] * $lsin[1] ) * $hgmax[3];
	$max = $max2 if($max < $max2);
	$max3 = sqrt ($Crystal->{Rg11} * $lsin[1] * $lsin[2] ) * ($hgmax[0] - 1);
	$max2 = sqrt ($Crystal->{Rg22} * $lsin[0] * $lsin[2] ) * ($hgmax[1] - 1);
	$max3 = $max2;
	$max2 = sqrt ($Crystal->{Rg33} * $lsin[0] * $lsin[1] ) * ($hgmax[2] - 1);
	$max3 = $max2 if($max3 < $max2);
	print  "     G^2max(Calculation) = ", ($max3*$pi2)**2.0, " - ", ($max*$pi2)**2.0, "\n";
#exit;

	print "\n";
	print("    Delete old input files [$InFile] and [$DataFile].\n");
	unlink($InFile);
	unlink($DataFile);
#print "C=$CopyScript $ScriptPath\n";
	if($CopyScript) {
		my $TargetScriptPath = 'GoLAMMPS.bat';
		print("    Delete old script file [$TargetScriptPath].\n");
		unlink($TargetScriptPath);
		if(!Utils::CopyFile($ScriptPath, $TargetScriptPath)) {
			print "Error: Can not copy [$ScriptPath] to [$TargetScriptPath].\n";
			exit;
		}
	}

	print "\n";
	print("    Save to [$DataFile].\n");
#print "f:$DataFile\n";
	my $out2 = JFile->new($DataFile, 'w') or die "Can not write to [$DataFile].\n";
	$out2->print("\n");
	$out2->printf("%10d atoms\n", $nExpandedAtomSiteList);
	$out2->print("\n");
	$out2->printf("%10d atom types\n", $nAtomType);
	$out2->print("\n");
# Vector: (xhi, 0, 0), (xy, yhi, 0), (xz, yz, zhi)
	$out2->printf("       %18.15g       %18.15g xlo xhi\n", 0.0, $a11);
	$out2->printf("       %18.15g       %18.15g ylo yhi\n", 0.0, $a22);
	$out2->printf("       %18.15g       %18.15g zlo zhi\n", 0.0, $a33);
	$out2->printf("       %18.15g       %18.15g       %18.15g xy xz yz\n", $a21, $31, $a32);
	$out2->print("\n");
	$out2->print("Masses\n");
	$out2->print("\n");
print "Atom types:\n";
	for(my $i = 0 ; $i < $nAtomType ; $i++) {
		my $p = $AtomTypeList[$i];
		my $atomname  = $p->AtomNameOnly();
		my $mass      = $p->AtomicMass();
		$out2->printf(" %4d %25.15f # %s\n", $i+1, $mass, $atomname);
printf(" %4d %25.15f # %s\n", $i+1, $mass, $atomname);
	}
	$out2->print("\n");
	$out2->print("Atoms\n");
	$out2->print("\n");
print "\n";
print "Atom sites:\n";
	for(my $i = 0 ; $i < $nExpandedAtomSiteList ; $i++) {
		my $p = $ExpandedAtomSiteList[$i];
		my $atomname  = $p->AtomNameOnly();
		my $occupancy = $p->Occupancy();
#		my $mult      = $p->Multiplicity();
#		my $id        = $p->Id();
		if($occupancy < 0.9999) {
			print "\nError: Partial occupancy (occ < 0.9999) is not allowed for atom #$i $atomname\n\n"; 
			exit;
		}
		my $charge = $Charge{$atomname};

		my $idAtom;
		my $pat;
		for(my $j = 0 ; $j < $nAtomType ; $j++) {
			if($atomname eq $AtomTypeList[$j]->AtomNameOnly()) {
				$idAtom = $j + 1;
				$pat = $AtomTypeList[$j];
				last;
			}
		}
		my ($x, $y, $z) = $p->Position();
		my $xr = $a11 * $x + $a21 * $y + $a31 * $z;
		my $yr = $a12 * $x + $a22 * $y + $a32 * $z;
		my $zr = $a13 * $x + $a23 * $y + $a33 * $z;
#                            "       1     1        -2.000000000000000         6.579577218202047         0.005786653198768         3.512676454470000
		$out2->printf("   %5d %5d       %18.15f       %18.15f       %18.15f       %18.15f  # %s\n", $i+1, $idAtom, $charge, $xr, $yr, $zr, $atomname);
print "  $atomname: ($x, $y, $z)  ($xr, $yr, $zr)\n";
	}
	$out2->print("\n");
	$out2->Close();

	print "\n";
	print("    Save to [$InFile].\n");
	unlink($InFile);
	my $out1 = JFile->new($InFile, 'w')   or die "Can not write to [$InFile].\n";

	$out1->print("#dimension 3\n");
	$out1->print("#newton on\n");
	$out1->print("units           metal\n");

#モデル構造
	$out1->print("\n");
	$out1->print("#atom_style atomic\n");
	$out1->print("atom_style      charge\n");
	$out1->print("#boundary x y z: [x,y,z]= p(eriodic)\n");
	$out1->print("boundary        p p p\n");
	$out1->print("box             tilt large\n");
	$out1->print("read_data       lmp_tmp.data\n");
#	$out1->print("read_data       $DataFile\n");

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

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

	$out1->print("pair_style      born/coul/long 10.0\n");
	for(my $i = 0 ; $i < $nAtomType ; $i++) {
#print "i=$i:$AtomTypeList[$i]\n";
		my $name1 = $AtomTypeList[$i]->AtomNameOnly();
		for(my $j = $i ; $j < $nAtomType ; $j++) {
			my $name2 = $AtomTypeList[$j]->AtomNameOnly();
# O  core O core    25.410 0.6937 32.32 0.0 12.0 0 0 0
#print "($i - $j) $name1 - $name2: ", $BParams{"$name1-$name2"}, "\n";
			my $Aij = $BParams{"$name1-$name2"}{Aij};
			my $bij = $BParams{"$name1-$name2"}{bij};
			my $Cij = $BParams{"$name1-$name2"}{Cij};
			$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";
			$out1->printf("pair_coeff   %3d %3d %12.4g %12.4f %12.4f %8.4f %8.4f\n", 
				$i+1, $j+1, $Aij, $bij, 0.0, $Cij, 0.0);
		}
	}

#計算条件
	$out1->print("\n");
	$out1->print("kspace_style    ewald 1.0e-4\n");
	$out1->print("#neighbor 2.0 bin\n");
	$out1->print("#neigh_modify    delay 10 every 1\n");
	$out1->print("#neigh_modify every 2 delay 10 check yes\n");
	$out1->print("neigh_modify    delay 0\n");
	$out1->print("#neighbor 0.3 bin\n");

#構造書き出し
	$out1->print("\n");
	$out1->print("#compute: com yes => subtract the drift of the center-of-mass\n");
	$out1->print("#         average yes => use the average as the reference position\n");
	$out1->print("#compute MSD1 all msd\n");
	$out1->print("#compute MSD1 all msd com no average no\n");
	$out1->print("dump            1 all custom 1000 lmp_tmp.dump id type xs ys zs ix iy iz\n");
	$out1->print("dump            2 all xtc 1000 lmp_tmp.xtc\n");
	$out1->print("#dump d0 all image 100 dump.*.jpg type type\n");
	$out1->print("#dump d1 mobile image 500 snap.*.png element element ssao yes 4539 0.6\n");
	$out1->print("#dump d2 all image 200 img-*.ppm type type zoom 2.5 adiam 1.5 size 1280 720\n");
	$out1->print("#dump m0 all movie 1000 movie.mpg type type size 640 480\n");
	$out1->print("#dump m1 all movie 1000 movie.avi type type size 640 480\n");
	$out1->print("#dump m2 all movie 100 movie.m4v type type zoom 1.8 adiam v_value s\n");

	$out1->print("#fix user-defined-ID group-ID style args\n");
	$out1->print("#  see https://lammps.sandia.gov/doc/fix.html\n");
	$out1->print("#fix 1 all nve\n");
	$out1->print("#fix 3 all nvt temp 300.0 300.0 0.01\n");
	$out1->print("#fix mine top setforce 0.0 NULL 0.0\n");
	$out1->print("#unfix 1\n");

#構造最適化
	my $pre;
	$pre = ($Function =~ /Relax/i)? '' : '#';
	$out1->print("\n");
	$out1->print("${pre}timestep        0.001\n");
	$out1->print("${pre}thermo_style    custom step etotal temp press lx vol\n");
	$out1->print("${pre}thermo 500\n");
	$out1->print("${pre}fix Relax1 all box/relax iso 0.0 vmax 0.01\n");
	$out1->print("#min_style: cg, sd, hftn, quickmin, fire\n");
	$out1->print("${pre}min_style cg\n");
	$out1->print("${pre}minimize 0.0 1.0e-20 1000 100000\n");
	if($Function =~ /Relax/i) {
		$out1->print("${pre}unfix Relax1\n");
	}
	else {
		$out1->print("#unfix Relax1\n");
	}
#MD
	$pre = ($Function =~ /MD/i)? '' : '#';
	$out1->print("\n");
	$out1->print("${pre}timestep        0.001\n");
	$out1->print("${pre}compute MSD1 all msd com no average no\n");
	$out1->print("${pre}thermo_style    custom step time temp epair emol etotal press vol density");
	for(my $i = 1 ; $i <= $nAtomType ; $i++) {
		$out1->print(" c_MSD1[$i]");
	}
	$out1->print("\n");
	$out1->print("${pre}thermo          10\n");
	$out1->print("${pre}velocity        all create 300.0 12345\n");
	$out1->print("${pre}fix             MD1 all npt temp  300.0 3000.0 0.2 iso 0.0 0.0 2.0\n");
	$out1->print("${pre}run             50000\n");
	$out1->print("${pre}unfix MD1\n");
	$out1->print("${pre}fix             MD1 all npt temp 3000.0 3000.0 0.2 iso 0.0 0.0 2.0\n");
	$out1->print("${pre}timestep        0.001\n");
	$out1->print("${pre}run             50000\n");
	$out1->print("${pre}unfix MD1\n");
	$out1->print("${pre}fix             MD1 all npt temp 3000.0  300.0 0.2 iso 0.0 0.0 2.0\n");
	$out1->print("${pre}timestep        0.001\n");
	$out1->print("${pre}run             50000\n");

	$out1->print("\n");
	$out1->print("write_restart   lmp_tmp.restart\n");
#	$out1->print("write_restart   $filenames[0].restart\n");

	$out1->print("\n");
	$out1->print("\n");
	$out1->print("# For reference, see blow, or search 'LAMMPS and keyword' \n");
	$out1->print("# URL   : https://lammps.sandia.gov/\n");
	$out1->print("# Manual: https://lammps.sandia.gov/doc/Manual.html\n");
	$out1->print("\n");
	$out1->print("#dump user-defined ID group-ID style nCycleStepForDump FileName args\n");
	$out1->print("# see https://lammps.sandia.gov/doc/dump.html\n");
	$out1->print("# style: atom atom/gz atom/mpiio cfg cfg/gz cfg/mpiio custom custom/gz custom/mpiio\n");
	$out1->print("#        dcd h5md image local molfile movie netcdf netcdf/mpiio vtk xtc xyz xyz/gz xyz/mpiio\n");
	$out1->print("# style=movie: mpg, avi, m4v\n");
	$out1->print("# style=custom: args=id(AtomId), (mol(MoleculeID), proc, procp1, type(AtomType), element, mass\n");
	$out1->print("#              x, y, z(unscaled atom coords), xs, ys, zs(scaled coords), xu, yu, zu(unwrapped coords)\n");
	$out1->print("#              xsu, ysu, zsu(scaled unwrapped), ix, iy, iz(box image)\n");
	$out1->print("#              vx, vy, vz(velocities), fx, fy, fz(forces on atom)\n");
	$out1->print("#              q(atom charge), mux, muy, muz, mu(orientation of dipoles)\n");
	$out1->print("#              radius, diameter, omegax, omegay, omegaz\n");
	$out1->print("#              angmomx, angmomy, angmomz, tqx, tqy, tqz\n");
	$out1->print("#              c_ID, c_ID[N], f_ID, f_ID[N], v_name\n");

# マニュアル表示
	$out1->print("\n");
	$out1->print("#units    : mass  len time energy   velocity force  Temp press  viscosity charge E-field dipole  density\n");
	$out1->print("# real    : g/mol A    fs  kcal/mol A/fs   kcal/mol.A K  atm      poise   e      V/A     e.A        g/cm3\n");
	$out1->print("# metal   : g/mol A    ps  eV       A/ps     eV/A     K  bars     poise   e      V/A     e.A        g/cm3\n");
	$out1->print("# si      : kg    m    s   J        m/s      N        K  Pa       Pa.s    C      V/m     C.m        kg/m3\n");
	$out1->print("# cgs     : g     cm   s   erg      cm/s     dynes    K  dyne/cm2 poise   esu  dyne/esu  10^18debye g/cm3\n");
	$out1->print("# electron: me    bohr fs  Hr       bohr/au  Hr/bohr  K  Pa               e      V/cm     debye\n");
	$out1->print("\n");
	$out1->print("#thermo_style style args (to write to log file)\n");
	$out1->print("# see https://lammps.sandia.gov/doc/thermo_style.html\n");
	$out1->print("#  style: step, elapsed, elaplong, dt, time\n");
	$out1->print("#         cpu, tpcpu, spcpu, cpuremain, part, timeremain\n");
	$out1->print("#         atoms, temp, press, pe, ke, etotal, enthalpy\n");
	$out1->print("#         evdwl, ecoul, epair, ebond, eangle, edihed, eimp\n");
	$out1->print("#         emol, elong, etail\n");
	$out1->print("#         vol, density, lx, ly, lz, xlo, xhi, ylo, yhi, zlo, zhi\n");
	$out1->print("#         xy, xz, yz, xlat, ylat, zlat\n");
	$out1->print("#         bonds, angles, dihedrals, impropers\n");
	$out1->print("#         pxx, pyy, pzz, pxy, pxz, pyz\n");
	$out1->print("#         fmax, fnorm, nbuild, ndanger\n");
	$out1->print("#         cella, cellb, cellc, cellalpha, cellbeta, cellgamma\n");
	$out1->print("#         c_ID, c_ID[I], c_ID[I][J], f_ID, f_ID[I], f_ID[I][J]\n");
	$out1->print("#         v_name, v_name[I]\n");

	$out1->Close();

	return 1;
}


1;