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

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

use strict;
use File::Basename;

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 $LAMMPSDatabaseDir = Deps::MakePath($LAMMPSPerlDir, 'Potentials', 0);
my $GULPDir           = ProgVars::GULPDir();
my $GULPLibDir        = ProgVars::GULPLibDir();
my @PotentialLibs     = ForceField::PotentialLibs();

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

	return $this;
}

sub DESTROY
{
	my $this = shift;
}

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


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

if($pInitialCrystal) {
print "MSD initialized with $pInitialCrystal\n";
}
else {
print "MSD initialized without initial crystal structure\n";
}
	$this->{pInitialCrystal} = ($pInitialCrystal)? $pInitialCrystal->Duplicate() : undef;
	$this->{pPrevCrystal}    = ($pInitialCrystal)? $pInitialCrystal->Duplicate() : undef;
	$this->{nMSDCount} = 0;
	return if(!defined $pInitialCrystal);

	$this->{pAtomType} = $pInitialCrystal->GetpAtomType(0);
	$this->{pAtomSite} = $pInitialCrystal->GetpAtomSite(1);
if(0) {
	my $pPos           = $pInitialCrystal->GetpAtomPosition(1);
# $pPosは［ｘ，ｙ，ｚ］のポインタの配列のため、コピーして使う
	my @Pos1 = @$pPos;
	my @Pos2 = @$pPos;
	for(my $i = 0 ; $i < @$pPos ; $i++) {
		$Pos1[$i][0] = $Pos2[$i][0] = $pPos->[$i][0];
		$Pos1[$i][1] = $Pos2[$i][1] = $pPos->[$i][1];
		$Pos1[$i][2] = $Pos2[$i][2] = $pPos->[$i][2];
	}
	$this->{pInitialPos} = \@Pos1;
	$this->{pPrevPos}    = \@Pos2;
}
	my $pPos = $this->{pInitialPos} = $pInitialCrystal->GetpDuplicatedAtomPosition(1);
	$this->{pPrevPos}    = $pInitialCrystal->GetpDuplicatedAtomPosition(1);
	$this->{pMSD}          = [];
	$this->{pOffsetForMSD} = [];
	for(my $i = 0 ; $i < @$pPos ; $i++) {
		$this->{pMSD}[$i]          = 0.0;
		$this->{pOffsetForMSD}[$i] = [0.0, 0.0, 0.0];
	}
print "  nAtom0=", scalar @$pPos, "\n";

	return 1;
}

sub AddMSD
{
	my ($this, $pCrystal) = @_;
	
	if(!defined $this->{pInitialCrystal}) {
		$this->InitializeMSD($pCrystal);
	}

	my $pSite     = $pCrystal->GetpAtomSite(1);
#print "nAtom(site)=", scalar @$pSite, "\n";
	my $pPos     = $pCrystal->GetpAtomPosition(1);
#print "nAtom(pos)=", scalar @$pPos, "\n";
	my $pRefPos  = $this->{pInitialPos};
	my $pPrevPos = $this->{pPrevPos};
# $pPosは［ｘ，ｙ，ｚ］のポインタの配列のため、コピーして使う
	my @PrevPos;
	for(my $i = 0 ; $i < @$pPos ; $i++) {
		$PrevPos[$i][0] = $pPrevPos->[$i][0];
		$PrevPos[$i][1] = $pPrevPos->[$i][1];
		$PrevPos[$i][2] = $pPrevPos->[$i][2];
	}
	my $poff     = $this->{pOffsetForMSD};
#print "  pIni=$pRefPos  pPrev=$pPrevPos\n";
#print "  ini ref pos = ($pRefPos->[0][0], $pRefPos->[0][1], $pRefPos->[0][2])\n";

	for(my $i = 0 ; $i < @$pPos ; $i++) {
		my $p = $pPos->[$i];
		my ($x, $y, $z) = @$p;
		$x = Utils::Reduce01($x);
		$y = Utils::Reduce01($y);
		$z = Utils::Reduce01($z);

		my $px = $PrevPos[$i][0];
		my $py = $PrevPos[$i][1];
		my $pz = $PrevPos[$i][2];

# wrapされた座標から、実際の変位量を計算する
# ただし、1/2単位格子以上の変位があると、変位量は不定になる。msdの計算はあきらめる
		if(abs($px - $x) > 1.0) {
			print "Error: Displacement larger than 1.0: i=$i  dx = $px - $x\n";
			exit;
		}
		if(abs($py - $y) > 1.0) {
			print "Error: Displacement larger than 1.0: i=$i  dy = $py - $x\n";
			exit;
		}
		if(abs($pz - $z) > 1.0) {
			print "Error: Displacement larger than 1.0: i=$i  dz = $pz - $x\n";
			exit;
		}

		my @off = (0.0, 0.0, 0.0);
		if($px - $x > 0.5) {
			$off[0] = 1.0;
			$poff->[$i][0] += 1.0;
		}
		elsif($x - $px > 0.5) {
			$off[0] = -1.0;
			$poff->[$i][0] -= 1.0;
		}

		if($py - $y > 0.5) {
			$off[1] = 1.0;
			$poff->[$i][1] += 1.0;
		}
		elsif($y - $py > 0.5) {
			$off[1] = -1.0;
			$poff->[$i][1] -= 1.0;
		}

		if($pz - $z > 0.5) {
			$off[2] = 1.0;
			$poff->[$i][2] += 1.0;
		}
		elsif($z - $pz > 0.5) {
			$off[2] = -1.0;
			$poff->[$i][2] -= 1.0;
		}

		my $r = $pCrystal->GetInterAtomicDistance(
				$pRefPos->[$i][0], $pRefPos->[$i][1], $pRefPos->[$i][2], 
				$x+$off[0], $y+$off[1], $z+$off[2]);
#		my $r = $pCrystal->GetInterAtomicDistance(
#				$px, $py, $pz, $x+$poff->[$i][0], $y+$poff->[$i][1], $z+$poff->[$i][2]);
		$this->{pMSD}[$i] = $r*$r;
#print "i=$i: r=$r = ($pRefPos->[$i][0], $pRefPos->[$i][1], $pRefPos->[$i][2])"
#	." - ($x+$off[0], $y+$off[1], $z+$off[2])\n" if($i == 5);
#print "      msd=$this->{pMSD}[$i]\n" if($i == 0);

#print "  ini ref pos = ($pRefPos->[0][0], $pRefPos->[0][1], $pRefPos->[0][2]) $pRefPos\n" if($i ==0);
		($PrevPos[$i][0], $PrevPos[$i][1], $PrevPos[$i][2]) = ($x, $y, $z);
#print "  fin ref pos = ($pRefPos->[0][0], $pRefPos->[0][1], $pRefPos->[0][2]) \n" if($i ==0);
	}

	return (++$this->{nMSDCount}, $this->{pMSD});
}

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

	my $pAtomType = $this->{pAtomType};
	my $pAtomSite = $this->{pAtomSite};
	my $pPos      = $this->{pInitialPos};

	my $TotalMSD = 0.0;
	my @AtomSiteMSD;
	my @AtomTypeMSD;
	for(my $i = 0 ; $i < @$pAtomType ; $i++) {
		$AtomTypeMSD[$i] = 0.0;
	}
#print "nAtomType=", scalar @$pAtomType, "\n";
#print "nAtomSite=", scalar @$pAtomSite, ", ", scalar @$pPos, "\n";
	my @natom;
	for(my $i = 0 ; $i < @$pAtomSite ; $i++) {
		my $p = $pAtomSite->[$i];
		my $name = $p->AtomNameOnly();

		$AtomSiteMSD[$i] = $this->{pMSD}[$i];
#		$AtomSiteMSD[$i] = $this->{pMSD}[$i] / $this->{nMSDCount};
#print "i=$i: AtomSiteMSD[$i] = $this->{pMSD}[$i] / $this->{nMSDCount}\n";
		for(my $j = 0 ; $j < @$pAtomType ; $j++) {
			my $namej = $pAtomType->[$j]->AtomNameOnly();
#print "j=$j: $namej eq $name\n";
			if($namej eq $name) {
				$natom[$j]++;
				$AtomTypeMSD[$j] += $AtomSiteMSD[$i];
#print "i=$i,j=$j: AtomTypeMSD[$j] = $AtomTypeMSD[$j] += $AtomSiteMSD[$i]\n";
			}
		}
		$TotalMSD += $this->{pMSD}[$i];
	}
	for(my $j = 0 ; $j < @$pAtomType ; $j++) {
		$AtomTypeMSD[$j] /= $natom[$j] if($natom[$j] > 0);
#print "j=$j: $AtomTypeMSD[$j] /= $natom[$j]\n";
	}
	$TotalMSD /= scalar @$pPos;
#print "Total msd = $TotalMSD /= ", scalar @$pPos, "\n";

	return (\@AtomSiteMSD, \@AtomTypeMSD, $TotalMSD);
}


1;