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

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

use strict;
use Sci::Science;
use JFile;

use ProgVars;
#use Crystal::MyUtility;
use Crystal::Crystal;
use Crystal::AtomType;
use Crystal::AtomSite;
use Crystal::CIF;

my $WebRootDir     = ProgVars::WebRootDir();
my $ProgramDir     = ProgVars::ProgramDir();
my $KListDBDir     = ProgVars::KListDBDir();

#===============================================
# 変数取得関数
#===============================================
sub KListDBDir {
	my $this = shift;
	return $this->{'KListDBDir'} if($this->{'KListDBDir'});
	return $this->{'KListDBDir'} = $KListDBDir;
}
sub SetKListDBDir { 
	my($this, $d) = @_; 
	$this->{'KListDBDir'} = $d;
	return $d;
}

#===============================================
# コンストラクタ、デストラクタ
#===============================================
BEGIN {
}

sub new
{
	my ($module) = @_;
	my $this = {};
	bless $this;
	return $this;
}

sub DESTROY
{
	my $this = shift;
}

#===============================================
# 一般関数
#===============================================
sub AddXSFAnimationFileStepFromCIF
{
	my ($this, $outfile, $iStep, $CIF) = @_;

	my $Crystal = $CIF->GetCCrystal();
	$this->AddXSFAnimationFileStepFromCrystal($outfile, $iStep, $Crystal);
}

sub WriteXSFFileHeader
{
	my ($this, $outfile) = @_;
	return $this->WriteXSFAnimationFileHeader($outfile, '');
}

sub WriteXSFAnimationFileHeader
{
	my ($this, $outfile, $nStep) = @_;

	unless(open(OUT,">$outfile")) {
		return 0;
	}
	binmode(OUT);

	if($nStep ne '') {
		print OUT "ANIMSTEPS $nStep\n";
	}
	print OUT "CRYSTAL\n";
#	print OUT "\n";

	close(OUT);
	return 1;
}

sub WriteXSFAnimationFileHeaderFromCIF
{
	my ($this, $outfile, $nStep, $CIF) = @_;
	return $this->WriteXSFAnimationFileHeader($outfile, $nStep);
}

sub AddXSFAnimationFileStepFromCrystal
{
	my ($this, $outfile, $iStep, $Crystal) = @_;

	return 0 unless(open(OUT,">>$outfile"));
	binmode(OUT);

	my ($a11,$a12,$a13,$a21,$a22,$a23,$a31,$a32,$a33) = $Crystal->LatticeVectors();
	print  OUT "PRIMVEC $iStep\n";
	printf OUT " %16.10f %16.10f %16.10f\n", $a11,$a12,$a13;
	printf OUT " %16.10f %16.10f %16.10f\n", $a21,$a22,$a23;
	printf OUT " %16.10f %16.10f %16.10f\n", $a31,$a32,$a33;
	print  OUT "CONVVEC $iStep\n";
	printf OUT " %16.10f %16.10f %16.10f\n", $a11,$a12,$a13;
	printf OUT " %16.10f %16.10f %16.10f\n", $a21,$a22,$a23;
	printf OUT " %16.10f %16.10f %16.10f\n", $a31,$a32,$a33;

	my @AtomList = $Crystal->GetCExpandedAtomSiteList();
	my $nAtom = @AtomList;

	print  OUT "PRIMCOORD $iStep\n";
	printf OUT " %d  1\n", $nAtom;
	for(my $i = 0 ; $i < $nAtom ; $i++) {
		my $atom = $AtomList[$i];
		my $name = $atom->AtomNameOnly();
		my ($x,$y,$z) = $atom->Position(1);
#print "X: $i: ($x,$y,$z)<br>\n";
		my ($vx,$vy,$vz) = $atom->Velocity();
#print "V: $i: ($vx,$vy,$vz)<br>\n";
		my $AtomicNumber = AtomType::GetAtomicNumber($name);
		my ($xc,$yc,$zc) = $Crystal->FractionalToCartesian($x,$y,$z);

		printf OUT "%3d   %12.7f %12.7f %12.7f   %12.7f %12.7f %12.7f\n",
			$AtomicNumber, $xc, $yc, $zc, $vx, $vy, $vz;
	}

	close(OUT);
	return $iStep;
}

sub WriteXSFFileFromCrystal
{
	my ($this, $outfile, $iCycle, $Crystal, $Crystal2) = @_;

	return 0 unless(open(OUT,">>$outfile"));
	binmode(OUT);

#	my $Title = $Crystal->CrystalName();
#	$Title = 'undefined' if($Title eq '');
#	print OUT "# this is a specification\n";
#	print OUT "# of $Title crystal structure\n";
#	print OUT "\n";

#	print OUT "CRYSTAL\n";
#	print OUT "\n";

	my ($a11,$a12,$a13,$a21,$a22,$a23,$a31,$a32,$a33) = $Crystal->LatticeVectors();
	print OUT "PRIMVEC $iCycle\n";
	printf OUT " %16.10f %16.10f %16.10f\n", $a11,$a12,$a13;
	printf OUT " %16.10f %16.10f %16.10f\n", $a21,$a22,$a23;
	printf OUT " %16.10f %16.10f %16.10f\n", $a31,$a32,$a33;
#	print  OUT "\n";
	print  OUT "CONVVEC $iCycle\n";
	printf OUT " %16.10f %16.10f %16.10f\n", $a11,$a12,$a13;
	printf OUT " %16.10f %16.10f %16.10f\n", $a21,$a22,$a23;
	printf OUT " %16.10f %16.10f %16.10f\n", $a31,$a32,$a33;
#	print  OUT "\n";

	my @AtomList = $Crystal->GetCExpandedAtomSiteList();
	my @AtomList2;
	@AtomList2 = $Crystal2->GetCExpandedAtomSiteList() if($Crystal2 ne '');
	my $nAtom = @AtomList; #$Crystal->nTotalExpandedAtomSite();
	my $nAtom2 = @AtomList2;
	$nAtom = $nAtom2 if($nAtom2 > 0 and $nAtom > $nAtom2);
#print "nAtom: $nAtom, $nAtom2<br>\n";

	print  OUT "PRIMCOORD $iCycle\n";
	printf OUT " %d  1\n", $nAtom;
#$Crystal2 != ''の時、重心の変化を調べる
	my ($centx1, $centy1, $centz1) = (0.0, 0.0, 0.0);
	my ($centx2, $centy2, $centz2) = (0.0, 0.0, 0.0);
	if(defined $Crystal2) {
		for(my $i = 0 ; $i < $nAtom ; $i++) {
			my $atom  = $AtomList[$i];
			my $atom2 = $AtomList2[$i];
			my ($x,$y,$z) = $atom->Position(1);
			my ($x2,$y2,$z2) = $atom2->Position(1);
			$x2 = $Crystal->FindNearestEquivalentFractionCoordinate($x2, $x);
			$y2 = $Crystal->FindNearestEquivalentFractionCoordinate($y2, $y);
			$z2 = $Crystal->FindNearestEquivalentFractionCoordinate($z2, $z);
			$centx1 += $x;
			$centy1 += $y;
			$centz1 += $z;
			$centx2 += $x2;
			$centy2 += $y2;
			$centz2 += $z2;
		}
		$centx1 = ($centx2 - $centx1) / $nAtom;
		$centy1 = ($centy2 - $centy1) / $nAtom;
		$centz1 = ($centz2 - $centz1) / $nAtom;
printf "Center change: (%12.6f, %12.6f, %12.6f)<br>\n", $centx1, $centy1, $centz1;
	}

	for(my $i = 0 ; $i < $nAtom ; $i++) {
		my $atom = $AtomList[$i];
		my $name = $atom->AtomNameOnly();
		my ($x,$y,$z)    = $atom->Position(1);
		my ($fx,$fy,$fz) = $atom->Force();
#print "$i: force=($fx,$fy,$fz)\n";
		($fx,$fy,$fz) = $atom->Velocity() if(!defined $fx);
#print "$i: Velocity=($fx,$fy,$fz)\n";
		my $AtomicNumber = AtomType::GetAtomicNumber($name);
		my ($xc,$yc,$zc) = $Crystal->FractionalToCartesian($x,$y,$z);
#print "$i($nAtom): ($xc,$yc,$zc)<br>";
		if(!defined $Crystal2) {
			if(!defined $fz) {
				$fx = $fy = $fz = 0.0;
			}
			printf OUT "%3d   %12.7f %12.7f %12.7f   %12.7f %12.7f %12.7f\n",
				$AtomicNumber, $xc, $yc, $zc, $fx, $fy, $fz;
		}
		else {
			my $atom2 = $AtomList2[$i];
			my ($x2,$y2,$z2) = $atom2->Position(1);
#print "$i: x-x2: (", $y2-$y, ")<br>\n";
			$x2 = $Crystal->FindNearestEquivalentFractionCoordinate($x2, $x) - $centx1;
			$y2 = $Crystal->FindNearestEquivalentFractionCoordinate($y2, $y) - $centy1;
			$z2 = $Crystal->FindNearestEquivalentFractionCoordinate($z2, $z) - $centz1;
#print "$i: x-x2: (", $y2-$y, ")<br>\n";

			my ($xc2,$yc2,$zc2) = $Crystal2->FractionalToCartesian($x2,$y2,$z2);
			printf OUT "%3d   %12.7f %12.7f %12.7f   %12.7f %12.7f %12.7f\n",
				$AtomicNumber, $xc, $yc, $zc, $xc2-$xc, $yc2-$yc, $zc2-$zc;
		}
	}

	close(OUT);

	return 1;
}

sub WriteXSFFileFromCIF
{
	my ($this, $outfile, $CIF, $CIF2) = @_;

	my $Crystal = $CIF->GetCCrystal();
	my $Crystal2;
	$Crystal2 = $CIF2->GetCCrystal() if($CIF2 ne '');

	return $this->WriteXSFFileFromCrystal($outfile, '', $Crystal, $Crystal2) = @_;




	unless(open(OUT,">>$outfile")) {
		return 0;
	}
	binmode(OUT);

	my $Title = $Crystal->CrystalName();
	$Title = 'undefined' if($Title eq '');
	print OUT "# this is a specification\n";
	print OUT "# of $Title crystal structure\n";
	print OUT "\n";

	print OUT "CRYSTAL\n";
	print OUT "\n";

	my ($a11,$a12,$a13,$a21,$a22,$a23,$a31,$a32,$a33) = $Crystal->LatticeVectors();
	print OUT "PRIMVEC\n";
	printf OUT " %16.10f %16.10f %16.10f\n", $a11,$a12,$a13;
	printf OUT " %16.10f %16.10f %16.10f\n", $a21,$a22,$a23;
	printf OUT " %16.10f %16.10f %16.10f\n", $a31,$a32,$a33;
	print OUT "\n";
	print OUT "CONVVEC\n";
	printf OUT " %16.10f %16.10f %16.10f\n", $a11,$a12,$a13;
	printf OUT " %16.10f %16.10f %16.10f\n", $a21,$a22,$a23;
	printf OUT " %16.10f %16.10f %16.10f\n", $a31,$a32,$a33;
	print OUT "\n";

	print OUT "PRIMCOORD\n";
	my @AtomList = $Crystal->GetCExpandedAtomSiteList();
	my @AtomList2;
	@AtomList2 = $Crystal2->GetCExpandedAtomSiteList() if($CIF2 ne '');
	my $nAtom  = @AtomList; #$Crystal->nTotalExpandedAtomSite();
	my $nAtom2 = @AtomList2;
	$nAtom = $nAtom2 if($nAtom2 > 0 and $nAtom > $nAtom2);
#print "nAtom: $nAtom, $nAtom2<br>\n";
	printf OUT " %d  1\n", $nAtom;
#$CIF2!=''の時、重心の変化を調べる
	my ($centx1, $centy1, $centz1) = (0.0, 0.0, 0.0);
	my ($centx2, $centy2, $centz2) = (0.0, 0.0, 0.0);
	if($CIF2 ne '') {
		for(my $i = 0 ; $i < $nAtom ; $i++) {
			my $atom  = $AtomList[$i];
			my $atom2 = $AtomList2[$i];
			my ($x,$y,$z) = $atom->Position(1);
			my ($x2,$y2,$z2) = $atom2->Position(1);
			$x2 = $Crystal->FindNearestEquivalentFractionCoordinate($x2, $x);
			$y2 = $Crystal->FindNearestEquivalentFractionCoordinate($y2, $y);
			$z2 = $Crystal->FindNearestEquivalentFractionCoordinate($z2, $z);
			$centx1 += $x;
			$centy1 += $y;
			$centz1 += $z;
			$centx2 += $x2;
			$centy2 += $y2;
			$centz2 += $z2;
		}
		$centx1 = ($centx2 - $centx1) / $nAtom;
		$centy1 = ($centy2 - $centy1) / $nAtom;
		$centz1 = ($centz2 - $centz1) / $nAtom;
print "Center change: ($centx1,$centy1,$centz1)<br>\n";
	}

	for(my $i = 0 ; $i < $nAtom ; $i++) {
		my $atom = $AtomList[$i];
		my $name = $atom->AtomNameOnly();
		my ($x,$y,$z) = $atom->Position(1);
		my $AtomicNumber = AtomType::GetAtomicNumber($name);
		my ($xc,$yc,$zc) = $Crystal->FractionalToCartesian($x,$y,$z);
#print "$i($nAtom): ($xc,$yc,$zc)<br>";
		if($CIF2 eq '') {
			printf OUT "%3d   %12.7f %12.7f %12.7f   %12.7f %12.7f %12.7f\n",
				$AtomicNumber, $xc, $yc, $zc, 0.0, 0.0, 0.0;
		}
		else {
			my $atom2 = $AtomList2[$i];
			my ($x2,$y2,$z2) = $atom2->Position(1);
#print "$i: x-x2: (", $y2-$y, ")<br>\n";
			$x2 = $Crystal->FindNearestEquivalentFractionCoordinate($x2, $x) - $centx1;
			$y2 = $Crystal->FindNearestEquivalentFractionCoordinate($y2, $y) - $centy1;
			$z2 = $Crystal->FindNearestEquivalentFractionCoordinate($z2, $z) - $centz1;
#print "$i: x-x2: (", $y2-$y, ")<br>\n";
			
			my ($xc2,$yc2,$zc2) = $Crystal2->FractionalToCartesian($x2,$y2,$z2);
			printf OUT "%3d   %12.7f %12.7f %12.7f   %12.7f %12.7f %12.7f\n",
				$AtomicNumber, $xc, $yc, $zc, $xc2-$xc, $yc2-$yc, $zc2-$zc;
		}
	}

	close(OUT);

	return 1;
}

sub WriteBXSFFile
{
	my ($this, $OutFile, $PrimCrystal, $Title, $nBand, $nx, $ny, $nz, $EF, $BandIndex, $BandEnergy, $ProgName, 
			$OneByOne, $IsPrint) = @_;
	$IsPrint = 1 if(!defined $IsPrint);

print("Save to [$OutFile]\n") if($IsPrint);
	open(OUT, ">$OutFile") or return 0;

	my ($nx1, $ny1, $nz1) = ($nx+1, $ny+1, $nz+1);
	print OUT <<EOT;
 BEGIN_INFO
   #
   # This is a Band-XCRYSDEN-Structure-File
   # Produced by $ProgName
   #
   # Launch as: xcrysden --bxsf $OutFile
   #
   Fermi Energy: $EF
 END_INFO

 BEGIN_BLOCK_BANDGRID_3D
   $Title
   BEGIN_BANDGRID_3D_fermi
     $nBand
     $nx1 $ny1 $nz1
EOT
	printf(OUT "    %12.4g %12.4g %12.4g\n", 0.0, 0.0, 0.0);
	printf(OUT "    %12.4g %12.4g %12.4g\n", $PrimCrystal->{Ra11}, $PrimCrystal->{Ra12}, $PrimCrystal->{Ra13});
	printf(OUT "    %12.4g %12.4g %12.4g\n", $PrimCrystal->{Ra21}, $PrimCrystal->{Ra22}, $PrimCrystal->{Ra23});
	printf(OUT "    %12.4g %12.4g %12.4g\n", $PrimCrystal->{Ra31}, $PrimCrystal->{Ra32}, $PrimCrystal->{Ra33});

	for(my $ib = 0 ; $ib < $nBand ; $ib++) {
		printf(OUT "   BAND: %3d\n", $BandIndex->[$ib]);
		my $ie = 0;
		for(my $ix = 0 ; $ix <= $nx ; $ix++) {
			for(my $iy = 0 ; $iy <= $ny ; $iy++) {
				for(my $iz = 0 ; $iz <= $nz ; $iz++) {
					printf(OUT "    %8.5f", $BandEnergy->[$ib][$ix][$iy][$iz]);
					print(OUT "\n") if($OneByOne);
					$ie++;
					if($ie % 10 == 0) {
#					if($ie % 10 == 0 and $iz != $nz) {
						print(OUT "\n") if(!$OneByOne);
					}
				}
				print(OUT "\n");
			}
			if($OneByOne) {
			}
			else {
				print(OUT "\n\n");
			}
			$ie = 0;
		}
	}
	print OUT <<EOT;
   END_BANDGRID_3D
 END_BLOCK_BANDGRID_3D
EOT

	close(OUT);
	return 1;
}

sub ReadKList
{
	my ($this, $iSPG, $SPGName, $CellType) = @_;
	my $KListDBDir = $this->KListDBDir();
	my $fmask = Deps::MakePath($KListDBDir, "${iSPG}_*.klist", 0);
#print "KList DB fmask: $fmask<br>\n";
	my ($path) = glob($fmask);
	if(!defined $path) {
#		if($SPGName =~ /^\s*F/i) {
		if($CellType =~ /^\s*F-Prim/i) {
			$path = Deps::MakePath($KListDBDir, "FCC-PrimitiveCell.klist", 0);
		}
#		elsif($SPGName =~ /^\s*I/i) {
		elsif($SPGName =~ /^\s*I-Prim/i) {
			$path = Deps::MakePath($KListDBDir, "BCC-PrimitiveCell.klist", 0);
		}
		elsif($SPGName =~ /^\s*[HR]/i) {
			$path = Deps::MakePath($KListDBDir, "HCP.klist", 0);
		}
	}
	if(!defined $path) {
		$path = Deps::MakePath($KListDBDir, "SimpleCubic.klist", 0);
	}
#print "KList DB file: $path<br>\n";
	
	my $in = new JFile($path, "r");
	return ($path) unless($in);
	
	my @KList;
	my $nPoint = 0;
	while(!$in->eof()) {
		my $line = $in->ReadLine();
		my $label = substr($line, 0, 10);
		my $rest = substr($line, 10);
		Utils::DelSpace($label);
		$nPoint++;
		next if($label eq '');
		last if($label eq 'END');
#print "line: $line<br>\n";
#print "rest: $rest<br>\n";
		my ($kx, $ky, $kz, $div) = Utils::Split("\\s+", $rest);
#print "div=$div: $kx $ky $kz<br>\n";
		$kx /= $div;
		$ky /= $div;
		$kz /= $div;
#print "$label: $kx $ky $kz ($nPoint)<br>\n";
		my %k;
		$k{'Label'}  = $label;
		$k{'kx'}     = $kx;
		$k{'ky'}     = $ky;
		$k{'kz'}     = $kz;
		$k{'nPoint'} = $nPoint;
		push(@KList, \%k);

		$nPoint = 0;
	}

	$in->Close();

	return ($path, @KList);
}

1;