package ForceField;
@ISA = qw(Exporter);

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

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

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

use Deps;
use Utils;

use ProgVars;
use MyApplication;

use Crystal::Crystal;
use Crystal::SpaceGroup;
use Crystal::AtomType;
use Crystal::Rietan;
use Sci qw($torad $e0 $a0 $e $pi log10 erfc);
use Sci::Science;

#===============================================
# 大域変数
#===============================================
my $HOME = $ENV{'HOME'};

my $ProgramDir        = ProgVars::ProgramDir();
my $LAMMPSPerlDir     = ProgVars::LAMMPSPerlDir();
my $TemplateDir       = Deps::MakePath($LAMMPSPerlDir, 'Template', 0);
my $LAMMPSDatabaseDir = Deps::MakePath($LAMMPSPerlDir, 'Potentials', 0);
my $GULPDir           = ProgVars::GULPDir();
my $GULPLibDir        = ProgVars::GULPLibDir();
my @PotentialLibs   = (
	$GULPLibDir,
	Deps::MakePath($ProgramDir, ["MD", "gulp-OtherLibraries",
		"www.ucl.ac.uk", "klmc", "Potentials", "Library"], 0),
	Deps::MakePath($ProgramDir, ["MD", "gulp-OtherLibraries"], 0),
	'.', 
	);

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

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

	return $this;
}

sub DESTROY
{
	my $this = shift;
}

#============================================================
#　継承クラスで定義しなおす関数
#============================================================

#===============================================
# 変数取得関数
#===============================================
sub PotentialLibs
{
	return @PotentialLibs;
}

sub LibraryDir { 
	my $this = shift;
	return $this->{'LibraryDir'} if($this->{'LibraryDir'});
	return $this->{'LibraryDir'} = $GULPLibDir;
}
sub SetLibraryDir { my($this,$d)=@_; return $this->{'LibraryDir'} = $d; }

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

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

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

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

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

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

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

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

sub GetpPotential
{
	my ($this) = @_;
	return $this->{pPot};
}

sub GetpAtoms
{
	my ($this) = @_;
	my $pPot = $this->GetpPotential();
	return $pPot->{pAtoms};
}

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

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

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

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

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

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

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

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

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

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

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

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

	return $Spring;
}

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

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

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

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

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

	return $pHash;
}

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

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

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

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

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

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

sub 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 ReadPotentialLibrary
{
	my ($this, $LibraryPath) = @_;

	my %hash;
	
	my $in = JFile->new($LibraryPath, 'r') or return undef;

	my %atoms;
	if($LibraryPath =~ /\.table/i) {
		$hash{table} = 1;
		$hash{type}  = 'table';
		my ($name1, $name2);
		my $RcutMax = 0.0;
		my $RcutMin = 1.0e10;

#ヘッダ情報読み込み
		$in->SkipTo("#\\s*Potential type:");
		while(1) {
			my $line = $in->ReadLine();
			last if(!defined $line);
			last if($line =~ /^#\s*$/);

			if($line =~ /([A-Z][a-z]*)\s*-\s*([A-Z][a-z]*):(.*)$/) {
#			elsif($line =~ /([A-Z][a-z\-_\.\d]*)\s*-\s*([A-Z][a-z\-_\.\d]*):\s*(\S.*)\s*$/) {
				$this->AddPotentialLines($line);
				my $name1 = $1;
				my $name2 = $2;
				$atoms{$name1}++;
				$atoms{$name2}++;
				my $s = $3;
				if($s =~ /De=([\-+\d\.]+)/) {
					$hash{"$name1-core"}{De} = $1;
					$hash{"$name1-core-$name2-core"}{De} 
						= $hash{"$name2-core-$name1-core"}{De} 
						= $1;
print "  De($name1-$name2)=$1\n";
				}
				if($s =~ /a0=([\-+\d\.]+)/) {
					$hash{"$name1-core"}{a0}
						= $hash{"$name2-core-$name1-core"}{a0} 
						= $1;
print "  a0($name1-$name2)=$1\n";
				}
				if($s =~ /r0=([\-+\d\.]+)/) {
					$hash{"$name1-core"}{r0}
						= $hash{"$name2-core-$name1-core"}{r0} 
						= $1;
print "  r0($name1-$name2)=$1\n";
				}
			}
			elsif($line =~ /([A-Z][a-z\-_\.\d]*):\s*(\S.*)\s*$/) {
				$this->AddPotentialLines($line);
				my $name1 = $1;
				$atoms{$name1}++;
				my $s = $2;
				if($s =~ /Z=([\-+\d\.]+)/) {
					$hash{"$name1-core"}{charge} = $1;
print "  Z($name1)=$1\n";
				}
			}
			elsif($line =~ /fCoulomb\s*=\s*(\S+)/) {
				$this->AddPotentialLines($line);
				$hash{fsCoulomb} = $1;
print "  fsCoulomb=$1\n";
			}
			elsif($line =~ /SplitEwaldSum\s*=\s*(\S+)/) {
				$this->AddPotentialLines($line);
				$hash{SplitEwaldSum} = $1;
print "  SplitEwaldSum=$1\n";
			}
			elsif($line =~ /Ewald\s+parameter.*=\s*(\S+)/) {
				$this->AddPotentialLines($line);
				$hash{gewald} = $1;
print "  Ewald parameter (G vector)=$1\n";
			}

		}

#ポテンシャル読み込み
		while(1) {
			my $line = $in->ReadLine();
			last if(!defined $line);
			next if($line =~ /^\s*#/);
			Utils::DelSpace($line);
			next if($line eq '');

#まず、keyとなるatom1-atom2ラインを探し、{key}パラメータに代入
			if($line =~ /([A-Z][a-z_\-\d]?)-([A-Z][a-z_\-\d]?)/) {
				($name1, $name2) = ($1, $2);
				$atoms{$name1}++;
				$atoms{$name2}++;
				Utils::DelSpace($line);
				$hash{"$name1-core-$name2-core"}{key} 
					= $hash{"$name2-core-$name1-core"}{key} 
					= $line;
print "  key for $name1-$name2: $line\n";
#				$this->AddPotentialLines($line);
			}
			elsif($line =~ /^\s*N\s+(\d+)/) {
				my $key1 = "$name1-core-$name2-core";
				my $key2 = "$name2-core-$name1-core";
#print "l:$line\n";
# Nラインを探してrパラメータを設定
				my $N = $1 + 0;
				
				while(1) {
					$line = $in->ReadLine();
					last if(!defined $line);
					Utils::DelSpace($line);
					last if($line eq '');
				}
				
				my $pos = $in->tell();
#最初の２データを読み込む
				$line = $in->ReadLine();
				my ($idx1, $r1, $U1, $F1) = Utils::Split("\\s+", $line);
				$line = $in->ReadLine();
				my ($idx2, $r2, $U2, $F2) = Utils::Split("\\s+", $line);

				$hash{$key1}{nR}    = $hash{$key2}{nR}    = $N;
				$hash{$key1}{Rmin}  = $hash{$key2}{Rmin}  = $r1 + 0.0;
				$hash{$key1}{Rstep} = $hash{$key2}{Rstep} = $r2 - $r1;
				$hash{$key1}{Rmax}  = $hash{$key2}{Rmax}  = $r1 + ($N - 1) * ($r2 - $r1);
				$RcutMax = $hash{$key1}{Rmax} if($RcutMax < $hash{$key1}{Rmax});
				$RcutMin = $hash{$key1}{Rmax} if($RcutMin > $hash{$key1}{Rmax});
				$this->AddPotentialLines("$key1: "
					."R = $hash{$key1}{Rmin} - $hash{$key1}{Rmax}, $hash{$key1}{Rstep} step, "
					."$hash{$key1}{nR} mesh");
print "    R = $hash{$key1}{Rmin} - $hash{$key1}{Rmax}, $hash{$key1}{Rstep} step, "
	."$hash{$key1}{nR} mesh\n";

# ポテンシャルデータ読み込み
				$in->seek($pos);
				my (@r, @U, @F);
				for(my $i = 0 ; $i < $N ; $i++) {
					$line = $in->ReadLine();
					last if(!defined $line);
					Utils::DelSpace($line);
					next if($line eq '' or $line =~ /^#/);
					
					my ($idx, $r, $U, $F) = Utils::Split("\\s+", $line);
					$r[$i] = $r;
					$U[$i] = $U;
					$F[$i] = $F;
#print "$key1 [$idx][$i]: r=$r U=$U\n";
#exit;
				}
				$hash{$key1}{pR} = $hash{$key2}{pR} = \@r;
				$hash{$key1}{pU} = $hash{$key2}{pU} = \@U;
				$hash{$key1}{pF} = $hash{$key2}{pF} = \@F;
			}
		}
		$hash{RmaxMax} = $RcutMax;
		$hash{RmaxMin} = $RcutMin;
print " ** Rmax = $hash{RmaxMin} - $hash{RmaxMax}\n";

		$in->Close();

		my @atoms = sort keys %atoms;
		$hash{pAtoms} = \@atoms;
		$this->{pPot} = \%hash;

		return \%hash;
	}

	my $pos;
	while(1) {
		my $line = $in->ReadLine();
		last if(!defined $line);
		$line =~ s/#.*$//;
		Utils::DelSpace($line);
		next if($line eq '');
#print "l:$line\n";

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

	my @atoms = sort keys %atoms;
	$hash{pAtoms} = \@atoms;
	$this->{pPot} = \%hash;
print "atoms: ", join(', ', @atoms), "\n";

	return \%hash;
}

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

	my $alpha     = $args{alpha};
	$pPot   = $this->GetpPotential() if(!defined $pPot);

	my @Atoms;
	if($pAtomTypes) {
		for(my $i = 0 ; $i < @$pAtomTypes ; $i++) {
			my $atom1 = $pAtomTypes->[$i];
			$Atoms[$i] = $atom1->AtomNameOnly();
		}
	}
	else {
		my $pAtoms = $this->GetpAtoms();
		@Atoms = @$pAtoms;
	}

print "\n";
print "Save potentials for potential type [$pPot->{type}] (gewald=$alpha)\n";
print "  Atoms: ", join(', ', @Atoms), "\n";
#foreach my $key (sort keys %$pPot) {
#	my $p = $pPot->{$key};
#	foreach my $key2 (sort keys %$p) {
#		print "   $key:$key2:$p->{$key2} \n";
#	}
#}

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

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

	my ($pr, $pU, %Uij);
	for(my $i = 0 ; $i < @Atoms ; $i++) {
		my $name1 = $Atoms[$i];
		for(my $j = $i ; $j < @Atoms ; $j++) {
			my $name2 = $Atoms[$j];
			($pr, $pU) = $this->CalPotentialArray(
				$name1, $name2, $rmin, $rstep, $nr, %args);
			$Uij{$name1}{$name2} = $pU;
#print "$i-$j: $name1-$name2: pr=$pr  pU=$pU  U=$pU->[0]\n";
#exit;
		}
	}

	for(my $l = 0 ; $l < $nr ; $l++) {
		$csv->print("$pr->[$l]");

		for(my $i = 0 ; $i < @Atoms ; $i++) {
			my $name1 = $Atoms[$i];
			for(my $j = $i ; $j < @Atoms ; $j++) {
				my $name2 = $Atoms[$j];
#print "$l:r=$pr->[$l]  Uij[$name1][$name2]=$Uij{$name1}{$name2}[$l]\n";
				$csv->print(",", $Uij{$name1}{$name2}[$l]);
			}
		}
		$csv->print("\n");
	}
	$csv->Close();

	return 1;
}

sub U
{
	my ($this, $name1, $name2, $r, %args) = @_;

	my $pPot  = $this->GetpPotential();
	my $type  = $pPot->{type};
	my $alpha = (defined $args{alpha})? $args{alpha} : $pPot->{gewald};
	my $p1  = $pPot->{"$name1-core"};
	my $p2  = $pPot->{"$name2-core"};
	my $pij = $pPot->{"$name1-core-$name2-core"};
	my $Z1  = $p1->{charge};
	my $Z2  = $p2->{charge};

	if($type eq 'buckingham') {
		my $Aij = $pij->{A};
		my $bij = $pij->{rho};
		my $Cij = $pij->{C};
		my $UC  = $KC * $Z1 * $Z2 / $r * erfc($alpha * $r);
		my $UA  = ($bij == 0.0)? 0.0 : $Aij * exp(-$r / $bij);
		my $UD  = -$Cij / $r**6;
		return $UC + $UA + $UD;
	}
	elsif($type eq 'morse') {
		my $De = $pij->{De};
		my $a  = $pij->{a0};
		my $r0 = $pij->{r0};
		my $fC = $pij->{fsCoulomb};
		my $UC = (1.0 - $fC) * $KC * $Z1 * $Z2 / $r;
		my $v  = 1.0 - exp(-$a * ($r - $r0));
		my $UM = $De * ($v*$v - 1.0);
		return $UC + $UM;
	}
	elsif($type eq 'table') {
		my ($pr0, $pU0, $pF0);
		$pr0 = $pij->{pR};
		$pU0 = $pij->{pU};
		$pF0 = $pij->{pF};
		return Algorism::Interpolate($pr0, $pU0, $r);
	}
}

sub CalPotential
{
	my ($this, $name1, $name2, $rmin, $rstep, $nr, %args) = @_;
	return $this->CalPotential($name1, $name2, $rmin, $rstep, $nr, %args);
}

sub CalPotentialArray
{
	my ($this, $name1, $name2, $rmin, $rstep, $nr, %args) = @_;

print "Calculate potential for $name1-$name2  r=($rmin, $rstep, $nr)\n";
	my (@r, @U);
	for(my $l = 0 ; $l < $nr ; $l++) {
		$r[$l] = $rmin + $l * $rstep;
		$U[$l] = $this->U($name1, $name2, $r[$l], %args);
	}

	return (\@r, \@U);
}

sub SaveLibrary
{
	my ($this, $pPot, $fname, $comment) = @_;

	$pPot = $this->GetpPotential() if(!defined $pPot);
	my $pAtoms = $this->GetpAtoms();

	print "Save potential parameters to GULP Library file [$fname].\n";
	if($pPot->{type} eq 'table') {
		print "  Error: Potential type [$pPot->{type}] is not compatible with library.\n";
		print "\n";
		return;
	}

	my $out = JFile->new($fname, 'w') or die "Can not write to [$fname].\n";
	$out->print("$comment\n") if($comment ne '');
	$out->print("#\n");
	$out->print("#  $pPot->{type} potential for IGZO\n");
	$out->print("#\n");

	$out->print("species 4\n");
	for(my $i = 0  ; $i < @$pAtoms ; $i++) {
		my $name1 = $pAtoms->[$i];
		my $Z = $pPot->{"$name1-core"}{charge};
		$out->printf("%2s   core %12.6f\n", $name1, $Z);
	}
	$out->print("\n");
	$out->print("$pPot->{type}\n");
	for(my $i = 0  ; $i < @$pAtoms ; $i++) {
	for(my $j = $i ; $j < @$pAtoms ; $j++) {
		my $name1 = $pAtoms->[$i];
		my $name2 = $pAtoms->[$j];
		my $key1 = "$name1-core-$name2-core";
		my $key2 = "$name2-core-$name1-core";

		if($pPot->{type} eq 'buckingham') {
			$out->printf("%3s core %3s core  %10.6g  %10.6g %10.6g  %4.1f %4.1f  %d %d %d\n",
				$name1, $name2, 
				$pPot->{$key1}{A}, $pPot->{$key1}{rho}, $pPot->{$key1}{C},
				$pPot->{$key1}{rmin}, $pPot->{$key1}{rmax},
				$pPot->{$key1}{idA}, $pPot->{$key1}{rho}, $pPot->{$key1}{C});
		}
		elsif($pPot->{type} eq 'morse') {
			$out->printf("%3s core %3s core  %10.6g %10.6g %10.6g %8.4f %4.1f %4.1f  %d %d %d\n",
				$name1, $name2, 
				$pPot->{$key1}{De}, $pPot->{$key1}{a0}, $pPot->{$key1}{r0}, $pPot->{$key1}{fsCoulomb},
				$pPot->{$key1}{rmin}, $pPot->{$key1}{rmax},
				$pPot->{$key1}{idDe}, $pPot->{$key1}{ida0}, $pPot->{$key1}{idr0});
		}
	}
	}
	$out->Close();
}

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

	my $pPot   = $this->GetpPotential();
	my $pAtoms = $this->GetpAtoms();

	print "Potential parameters\n";
	print "  Atoms: ", join(', ', @$pAtoms, "\n");
	print "  type: $pPot->{type}\n";
	if($pPot->{table}) {
		print "  Tabulated potential\n";
	}
	else {
		print "  Parameter potential\n";
	}

	print "  Charges\n";
	for(my $i = 0  ; $i < @$pAtoms ; $i++) {
		my $name1 = $pAtoms->[$i];
		my $Z = $pPot->{"$name1-core"}{charge};
		printf "    %4s: Z=%12.6f\n", $name1, $Z;
	}

	print "  Other parameters\n";
	if($pPot->{type} eq 'table') {
		print "  fsCoulomb=$pPot->{fsCoulomb}\n";
		print "  SplitEwaldSum=$pPot->{SplitEwaldSum}\n";
		print "  gewald=$pPot->{gewald}\n";
		my $RcutMax = $pPot->{RmaxMax};
		my $RcutMin = $pPot->{RmaxMin};
			print "    Rmax=$RcutMin - $RcutMax\n";
	}
	for(my $i = 0  ; $i < @$pAtoms ; $i++) {
	for(my $j = $i ; $j < @$pAtoms ; $j++) {
		my $name1 = $pAtoms->[$i];
		my $name2 = $pAtoms->[$j];
		my $key1 = "$name1-core-$name2-core";
		my $key2 = "$name2-core-$name1-core";
		if($pPot->{type} eq 'buckingham') {
			my $Aij = $pPot->{"$name1-core-$name2-core"}{A};
			my $bij = $pPot->{"$name1-core-$name2-core"}{rho};
			my $Cij = $pPot->{"$name1-core-$name2-core"}{C};
			printf "    %8s: Aij=%12.6g  rho=%12.6f  C=%12.6f\n",
				"$name1-$name2", $Aij, $bij, $Cij;
		}
		elsif($pPot->{type} eq 'morse') {
			my $De = $pPot->{$key1}{De};
			my $a  = $pPot->{$key1}{a0};
			my $r0 = $pPot->{$key1}{r0};
#			my $fC = $pPot->{$key1}{fsCoulomb};
			printf "    %8s: De=%12.6g  a0=%12.6f  r0=%12.6f\n",
				"$name1-$name2", $De, $a, $r0;
		}
		elsif($pPot->{type} eq 'table') {
			my $De    = $pPot->{$key1}{De};
			my $a     = $pPot->{$key1}{a0};
			my $r0    = $pPot->{$key1}{r0};
			my $nR    = $pPot->{$key1}{nR};
			my $Rmin  = $pPot->{$key1}{Rmin};
			my $Rstep = $pPot->{$key1}{Rstep};
			my $Rmax  = $pPot->{$key1}{Rmax};
			my $pr    = $pPot->{$key1}{pR};
			my $pU    = $pPot->{$key1}{pU};
			my $pF    = $pPot->{$key1}{pF};
			print "    $name1-$name2: De=$De  a0=$a  r0=$r0\n";
			print "          r=($Rmin, $Rmax)  rstep=$Rstep  nR=$nR\n";
			print "          pointes: pr=$pr  pU=$pU  pF=$pF\n";
		}
	}
	}
	print "\n";

	print "Potential lines\n";
	my @l = $this->GetPotentialLines();
	for(my $i = 0 ; $i < @l ; $i++) {
		print "  $l[$i]\n";
	}
	print "\n";
}


1;