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

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

use strict;
use File::Basename;

use ProgVars;
use Utils;
use GraphData;
use Sci::Science;

use Crystal::Crystal;
use Crystal::AtomType;
use Crystal::AtomSite;
use Crystal::XCrySDen;
use Crystal::CIF;

#===============================================
# スクリプト大域変数
#===============================================
my $LF = "<br>\n";
my $DirectorySeparator = Deps::DirectorySeparator();

#===============================================
# パス（読み込みDB）
# Web関係変数
# CGI の仮想パス
# プログラム名
#===============================================
my $Program      = ProgVars::Program();
my $ProgramDir   = ProgVars::ProgramDir();
my $CRYSTAL06Dir = ProgVars::CRYSTAL06Dir();
my $BasisSetDir  = Deps::MakePath($CRYSTAL06Dir, "Basis_Sets", 0);
my $WebRootDir   = ProgVars::WebRootDir();
#my $DVXaDir      = ProgVars::DVXaDir();
#my $PeriodicTable1DBPath = Deps::MakePath($DVXaDir, "PeriodicTable1.html", 0);

#============================================================
#　変数等取得、再定義関数
#============================================================
sub ClearAll      { my $this=shift; undef $this->{'FileType'}; undef $this->{'DataArray'}; }
sub FileType      { return shift->{'FileType'}; }
sub FileName      { return shift->{'FileName'}; }
sub SetFileName   { my ($this,$f)=@_; return $this->{'FileName'} = $f; }
#sub DataArray     { return shift->{'DataArray'}; }
#sub SetDataArray  { my ($this, $da) = @_; return $this->{'DataArray'} = $da; }
sub SetSampleName { my ($this,$n)=@_; return $this->{'SampleName'} = $n; }
sub SampleName    { return shift->{'SampleName'}; }
sub GetFileName {
	my($this, $path, $fname) = @_; 
	my ($drive, $dir, $filename, $ext, $lastdir, $filebody)
		= Deps::SplitFilePath($path);
	my $s = Deps::MakePath($drive, $dir);
	return Deps::MakePath($s, $fname);
}
sub GetD12FileName { my($this,$f)=@_; return Deps::ReplaceExtension($f, ".d12"); }
sub GetD3FileName  { my($this,$f)=@_; return Deps::ReplaceExtension($f, ".d12"); }
sub FindD12FileName
{
	my ($path) = @_;
	my @files = <*.d12>;
print "  *** d12 file: $files[0]\n";
	return $files[0];
}
sub ReadSampleNameFromD12File
{
	my ($path) = @_;
	my $in = new JFile($path, "rb");
	if(!$in) {
print("Error: Can not read [$path].\n");
		return undef;
	}
	my $line = $in->ReadLine();
	$in->Close();
	Utils::DelSpace($line);

	return $line;
}

#============================================================
#　継承クラスで定義しなおす関数
#============================================================
sub CheckFileType
{
	my ($path) = @_;
	my ($drive, $dir, $filename, $ext, $lastdir, $filebody)
		= Deps::SplitFilePath($path);
#print "f: $filename\n";
	if($filename =~ /\.out$/i) {
		return "CRYSTAL06 OUTPUT file";
	}
	elsif($filename =~ /\.outp$/i) {
		return "CRYSTAL06 PROPERTY OUTPUT file";
	}
	elsif($filename =~ /^DOSS\.DAT$/i) {
		return "CRYSTAL06 DOS file";
	}
	elsif($filename =~ /^BAND\.DAT$/i) {
		return "CRYSTAL06 BAND file";
	}
	return undef;
}

sub ReadStructureFromD12
{
	my ($this, $path) = @_;
	my $Crystal = new Crystal;
	my $SPG = $Crystal->GetCSpaceGroup();
	my $in = new JFile($path, "rb");
	if(!$in) {
print("Error in CRYSTAL06::ReadStructureFromD12: Can not read [$path].\n");
		return undef;
	}

	my $Sample = $in->ReadLine();
	Utils::DelSpace($Sample);
	my $line = $in->ReadLine();
	if($line =~ /CRYSTAL/i) {
		$in->ReadLine();
		my $iSPG = $in->ReadLine() + 0;
print "  iSPG: $iSPG\n";
		$SPG->ReadRietanSpaceGroupDB($iSPG);
#print "SPG:ls: ", $SPG->LatticeSystem(), "\n";
		$Crystal->SetLatticeSystem($SPG->LatticeSystem());

		$line = $in->ReadLine();
		my (@a) = Utils::Split("\\s+", $line);
#print "l: $a[0]\n";
		$Crystal->SetMinimalLatticeParameters(@a);
		my $nSites = $in->ReadLine() + 0;
print "  nSites: $nSites\n";
		my %AtomCount;
		for(my $is = 0 ; $is < $nSites ; $is++) {
			$line = $in->ReadLine();
			my ($AtomicNumber, $x, $y, $z) = Utils::Split("\\s+", $line);
			my $atomtype = $AtomicNumber;
			$AtomCount{$atomtype}++;
			my $Label = $atomtype . "-" . $AtomCount{$atomtype};
			$Crystal->AddAtomSite($Label, $atomtype, $x, $y, $z, 1.0);
print "    $Label: $atomtype: ($x, $y, $z)\n";
		}
	}
	$in->Close();
	$Crystal->ExpandCoordinates();

	return $Crystal;
}


sub ReadOUTPUT
{
	my ($this, $path, $target) = @_;
	my $HartreeToeV = Sci::HartreeToeV();

print "CRYSTAL06::ReadOUTPUT: Read [$path].\n";

	my $in = new JFile($path, "rb");
	return undef unless($in);

	$this->SetDataArray(new GraphDataArray);
	my $Data = new GraphData;

	my ($SampleName) = $in->ReadLine();
	Utils::DelSpace($SampleName);
#print "Sample: $SampleName\n";
	$Data->SetTitle($SampleName);

	if($target =~ /ForceConvergence/i) {
	}
	else { #'EnergyConvergence'
		my @Iteration;
		my @Energy;
		my $iter = 0;
		$in->rewind();
		while(!$in->eof()) {
			my $line = $in->SkipTo("^ CYC ");
#			my $line = $in->SkipTo("^\\s*CYC \\d+");
#print "line: $line\n";
			last if($in->eof());
			my ($DETot) = ($line =~ /DETOT\s+([\d\.+\-Ee]+)/i);
			$Energy[$iter] = $DETot * $HartreeToeV;
			$Iteration[$iter] = $iter+1;
print "i=$iter: dE(total)=$Energy[$iter]\n";
			$iter++;
		}
		$Data->{'x0'}      = \@Iteration;
		$Data->{'x0_Name'} = "Iteration";
		$Data->{'y0'}      = \@Energy;
		$Data->{'y0_Name'} = "Energy / eV";
	}

	$in->Close();

	$Data->CalMinMax();
	$this->{'DataArray'}->AddGraphData($Data);
	$this->{'DataArray'}->{'TargetData'} = $target;

print "CRYSTAL06::ReadOUTPUT: finished.\n";
	return $path;
}

sub ReadOUTPUTP
{
	my ($this, $filename, $TargetData) = @_;
	return undef;
}

sub ReadDOSS
{
	my ($this, $filename) = @_;
	my $HartreeToeV = Sci::HartreeToeV();

print "CRYSTAL06::ReadDOSS: Read [$filename].\n";
	my $d12File = FindD12FileName($filename);
	my $SampleName = ReadSampleNameFromD12File($d12File);

	my $EF = 0.0;
print "EF: $EF eV\n";

	my $in = new JFile($filename, "rb");
	if(!$in) {
print("Error: Can not read [$filename].\n");
		return undef;
	}

	$this->SetDataArray(new GraphDataArray);
## NEPTS   102 NPROJ     3 NSPIN     1
	my $line = $in->ReadLine();
#print "line: $line\n";
	my ($nData) = ($line =~ /NEPTS\s+(\d+)\s/i);
	my ($nDOS)  = ($line =~ /NPROJ\s+(\d+)\s/i);
	my ($nSpin) = ($line =~ /NSPIN\s+(\d+)\s/i);
print("  nData: $nData\n");
print("  nDOS : $nDOS\n");
print("  nSpin: $nSpin\n");

	my $Data0 = new GraphData;
	$Data0->SetTitle($SampleName);
	my $Data1;
	if($nDOS > 1) {
		$Data1 = new GraphData;
		$Data1->SetTitle($SampleName);
	}

##
#@ XAXIS LABEL "ENERGY (HARTREE)"
#@ YAXIS LABEL "DENSITY OF STATES (STATES/HARTREE/CELL)"
	$in->ReadLine();
	$in->ReadLine();
	$in->ReadLine();

	my @Energy;
	my @pDOS;
	my @Ne;
	for(my $j = 0 ; $j < $nDOS ; $j++) {
		my @DOS;
		$pDOS[$j] = \@DOS;
	}

	for(my $i = 0 ; $i < $nData ; $i++)  {
		$line = $in->ReadLine();
		last if(not defined $line);

		my @a = Utils::Split("\\s+", $line);
		last if(not defined $a[1]);

		$Energy[$i] = $HartreeToeV*($a[0] - $EF);
print "  $i: $Energy[$i]   ";
		for(my $j = 0 ; $j < $nDOS ; $j++) {
			$pDOS[$j]->[$i] = $a[$j+1];
print $pDOS[$j]->[$i], " ";
		}
print "\n";
		$Ne[$i]     = 2;
		last if($in->eof());
	}
	for(my $j = 0 ; $j < $nDOS ; $j++) {
		$Data0->{"x$j"}        = \@Energy;
		$Data0->{"x${j}_Name"} = "Energy / eV";
		$Data0->{"y$j"}        = $pDOS[$j];
		$Data0->{"y${j}_Name"} = "DOS";
	}
	$Data0->CalMinMax();
	$this->{'DataArray'}->AddGraphData($Data0);

print "CRYSTAL06::ReadDOSS: finished.\n";
	return $filename;
}

sub ReadBAND
{
	my ($this, $path) = @_;
	my $HartreeToeV = Sci::HartreeToeV();

print "CRYSTAL06::ReadBAND: Read [$path].\n";
	my $d12File = FindD12FileName($path);
	my $SampleName = ReadSampleNameFromD12File($d12File);
	my $Crystal = $this->ReadStructureFromD12($d12File);

	my $EF = 0.0;
print "EF: $EF eV\n";

	my $in = new JFile($path, "rb");
	return undef unless($in);

## NKPT    20 NBND     8 NSPIN     1
	my $line = $in->ReadLine();
#print "line: $line\n";
	my ($nKPoint) = ($line =~ /NKPT\s+(\d+)\s/i);
	my ($nBand)   = ($line =~ /NBND\s+(\d+)\s/i);
	my ($nSpin)   = ($line =~ /NSPIN\s+(\d+)\s/i);
print("  nKPoint: $nKPoint\n");
print("  nBand  : $nBand\n");
print("  nSpin  : $nSpin\n");

## NPANEL     3
	$line = $in->ReadLine();
	my ($nPanel) = ($line =~ /NPANEL\s+(\d+)\s/i);
print("  nPanel : $nPanel\n");

	my (@kbx, @kby, @kbz, @nKPoints);
	for(my $ip = 0 ; $ip < $nPanel+1 ; $ip++) {
##      1   (1,0,0)/2
		$line = $in->ReadLine();
#print "line: [$line]";
		my ($n, $x, $y, $z, $m) = ($line =~ /(\d+)\D+(\d+)\D+(\d+)\D+(\d+)\D+(\d+)/);
#print "$ip=$n,$x,$y,$z,$m\n";
		$kbx[$ip] = $x / $m;
		$kby[$ip] = $y / $m;
		$kbz[$ip] = $z / $m;
		$nKPoints[$ip] = $n;
print "Boundary $ip [$nKPoints[$ip]]: ($kbx[$ip],$kby[$ip],$kbz[$ip])\n";
	}

	$in->SkipTo("\@ VIEW ");
	my $EnergyBandArray = new EnergyBandArray;
	$EnergyBandArray->SetTitle($SampleName);
	$EnergyBandArray->SetCrystal($Crystal);
	my $Distance = 0.0;
	for(my $iK = 0 ; $iK < $nKPoint ; $iK++) {
		$line = $in->ReadLine();
		my ($k, @e) = Utils::Split("\\s+", $line);
		for(my $iL = 0 ; $iL < @e ; $iL++) {
			if($iK == 0) {
				$EnergyBandArray->AddEnergyBand();
			}
			$EnergyBandArray->AddByKE($iL, $iK, 0, 0, 
					$HartreeToeV*($e[$iL]-$EF));
		}
#print "k[$iK]: $kx $ky $kz\n";
#		for(my $iL = 0 ; $iL < $nLevels ; $iL++) {
#			$line = $in->ReadLine();
#			my ($ib,$e) = Utils::Split("\\s+", $line);
#			$EnergyBandArray->AddByKE($iL, $kx, $ky, $kz, $e-$EF);
#		}
	}
	$EnergyBandArray->CalMinMax();
	$EnergyBandArray->GetBandBoundary();

	my $DataArray = new GraphDataArray;
	$this->SetDataArray($DataArray);

	$EnergyBandArray->SetGraphDataArray($DataArray);

print "CRYSTAL06::ReadBAND: finished.\n";
	return $path;
}

sub ReadFiles
{
	my ($this, $filename, $TargetData) = @_;
	$TargetData = '' unless($TargetData);
	$this->ClearAll();

	my $FileType = $this->{'FileType'} = CheckFileType($filename);
	$this->SetFileName($filename);

	if($FileType eq "CRYSTAL06 OUTPUT file") {
		return $this->ReadOUTPUT($filename, $TargetData);
	}
	elsif($FileType eq "CRYSTAL06 PROPERTY OUTPUT file") {
		return $this->ReadOUTPUTP($filename, $TargetData);
	}
	elsif($FileType eq "CRYSTAL06 DOS file") {
		return $this->ReadDOSS($filename);
	}
	elsif($FileType eq "CRYSTAL06 BAND file") {
		return $this->ReadBAND($filename);
	}
	return undef;
}

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

sub DESTROY { my $this = shift; }

#============================================================
#　ファイル読み込み関数
#============================================================

#============================================================
#　一般関数
#============================================================

sub MakeD12File
{
	my ($this, $Crystal, $Function, $SpinPolarized, $fname) = @_;

#print "<H1>Make .d12</H1>\n";

	my $SPGName = $Crystal->SPGNameByOutputMode();
	my $iSPG    = $Crystal->iSPGByOutputMode();
	my ($a,$b,$c,$alpha,$beta,$gamma) 
		= $Crystal->LatticeParametersByOutputMode(0);
#結晶系を格子定数から探しなおす
	my $LatticeSystem = lc $Crystal->SetLatticeSystem(undef, 1);
print("LatticeSystem: $LatticeSystem<br>\n");

#ファイル作製開始
	my $out = new JFile($fname, "w");
	unless($out) {
		print "Can not write to $fname.$LF$LF";
		return 0;
	}

	my $Sample = $this->SampleName();
	$out->print("$Sample\n");
	$out->print("CRYSTAL\n");
	my $SPGMode = 0; # Specified by International Table Number
	my $IsRhombohedralAxis = 0;
	if($SPGName =~ /^R/i and $LatticeSystem eq 'trigonal') {
		$IsRhombohedralAxis = 1;
	}
	my $UseStandardOrigin = 0;
	$out->print("$SPGMode $IsRhombohedralAxis $UseStandardOrigin\n");

	$out->print(" $iSPG     ($SPGName)\n");

	if($LatticeSystem eq 'cubic') {
		$out->printf("   %12.6f\n", $a);
	}
	elsif($LatticeSystem eq 'tetragonal' or $LatticeSystem eq 'hexagonal') {
		$out->printf("   %12.6f %12.6f\n", $a, $c);
	}
	elsif($LatticeSystem eq 'orthorhombic') {
		$out->printf("   %12.6f %12.6f %12.6f\n", $a, $b, $c);
	}
	elsif($LatticeSystem eq 'trigonal') {
		$out->printf("   %12.6f   %12.6f\n", $a, $alpha);
	}
	elsif($LatticeSystem eq 'monoclinic') {
		if(abs($alpha-90.0)>1.0e-5) {
			$out->printf("   %12.6f %12.6f %12.6f   %12.6f\n", $a, $b, $c, $alpha);
		}
		elsif(abs($beta-90.0)>1.0e-5) {
			$out->printf("   %12.6f %12.6f %12.6f   %12.6f\n", $a, $b, $c, $beta);
		}
		elsif(abs($gamma-90.0)>1.0e-5) {
			$out->printf("   %12.6f %12.6f %12.6f   %12.6f\n", $a, $b, $c, $gamma);
		}
	}
	else {
		$out->printf("   %12.6f %12.6f %12.6f   %12.6f %12.6f %12.6f\n", 
					$a, $b, $c, $alpha, $beta, $gamma);
	}

	my @ExpandedAtomSiteList = $Crystal->GetCExpandedAtomSiteListByOutputMode();
	my $nAtomSite = @ExpandedAtomSiteList;
	$out->print("$nAtomSite\n");

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

		$out->printf(" %3d   %12.8f  %12.8f  %12.8f\n", 
				$AtomicNumber, $x, $y, $z);
	}

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

		$out->printf(" %3d   %12.8f  %12.8f  %12.8f  (Occ=%f)\n", 
				$AtomicNumber, $x, $y, $z, $occupancy);
	}
	$out->print("END\n");

print "Read Basis Sets from [$BasisSetDir].<br>\n";
#print "DVXaDir: [$DVXaDir].<br>\n";
#print "PeriodicTable1DBPath: [$PeriodicTable1DBPath].<br>\n";

	my @AtomTypeList = $Crystal->GetCAtomTypeList();
	for(my $i = 0 ; $i < @AtomTypeList ; $i++) {
		my $type   = $AtomTypeList[$i];
		my $name   = lc $type->AtomNameOnly();
		my $charge = $type->Charge();
		my $AtomicNumber = AtomType::GetAtomicNumber($name);
		$AtomicNumber = int($AtomicNumber + 0.001);

		my $EName = AtomType::GetEnglishName($name);
#		my ($AtmNum, $JName, $EName, $NameOrigin) = &GetIonInfFromPeriodicTable1($name);
		$name = lc $EName;
#		$name = 'silver' if($name eq 'ag');
#		$name = 'indium' if($name eq 'in');
#		$name = 'oxygen' if($name eq 'o');

		my $DBFile = Deps::MakePath($BasisSetDir, "${name}.html");
print "   Read [$name] Basis Sets from [$DBFile].<br>\n";
#$out->print("   Read [$name] Basis Sets from [$DBFile].\n");
		my $content;
		my $in = new JFile($DBFile, "r");
		unless($in) {
			print "<H2>Error: Can not read [$DBFile].</H2>\n";
			$out->print("***Insert [$name] Basis set here\n");
			next;
		}
		while(!$in->eof()) {
			$content .= $in->ReadLine();
		}
#print "content: $content<br>\n";
		my ($basis) = ($content =~ /\<pre\>(.*?)\<\/pre\>/si);
#print "bcontent: $basis<br>\n";
		my @lines = Utils::Split("\\n", $basis);
		my $Go = 0;
		my $Format = 0;
		for(my $i = 0 ; $i < @lines ; $i++) {
#print("a:$lines[$i] (Go=$Go)<br>\n");
#			if($lines[$i] =~ /^\s*\d+\s+\d+\s*$"/) {
			my @a = Utils::Split("\\s+", $lines[$i]);
			if($Go == 0 and (@a == 2 or @a == 5) 
			   and $lines[$i] =~ /^[\s\.\d]*$/) {
				if(@a == 5) {
					$Format = 1;
					$out->print("$AtomicNumber N\n");
				}
				$Go = 1;
			}
			next if($Go == 0);
			$lines[$i] =~ s/\s*(\r\n)$//;
			last if($lines[$i] =~ /^\s*$/);
			if($Go == 1) {
				if($Format == 0) {
					$out->print("$AtomicNumber $a[1]\n");
				}
				else {
					$out->print("$lines[$i]\n");
				}
#print("$lines[$i]<br>\n");
			}
			elsif($Go == 2) {
#				if($lines[$i] =~ /^\s*\d/) {
#					$out->print("$lines[$i]\n");
#				}
				$out->print("$lines[$i]\n");
			}
			else {
				last if($lines[$i] =~ /^\s*[a-zA-Z]/);
				$out->print("$lines[$i]\n");
			}
			$Go++;
		}
	}

	$out->print("99   0\n");
	$out->print("END\n");
	if($SpinPolarized eq 'Yes') {
		$out->print("UHF\n");
	}
	$out->print("SHRINK\n");
	$out->print("8 6 8\n");
	$out->print("END\n");
	$out->print("FMIXING\n");
	$out->print("30\n");
	$out->print("END\n");

	$out->Close();

	return 1;
}

sub MakeD3File
{
	my ($this, $Crystal, $Function, $D3) = @_;

	my $Sample = $this->SampleName();

	unless(open(OUTCRYSTAL06, ">$D3")) {
		print "<H2>Error!: Can not write to [$D3]</H2>\n";
		return -1;
	}
	print OUTCRYSTAL06 <<ENDofDATA;
EDFT
END
ANGSTROM
ECHG
0
65
ATOMS
1  -1 1 1
1  0 0 0
1  1 1 -1
MARGINS
1.9 1.9 1.9 1.9
END
POLI
4 0 4
BAND
$Sample
5 12 30 7 14 1 0
0 0 0  6 6 0
6 6 0  9 6 3
9 6 3  6 6 6
6 6 6  0 0 0
0 0 0  9 6 3
NEWK
6 6
1 0
DOSS
2 100 7 14 1   14  0
-1 1
-1 2
PROF
 2 0 1 0
 118 3
CP
 1  100.  100.  1
 1 0 0
 1.
BR
15. .3
END
PROF
 3 0 0 0
CP
 3  13.  13.  1
 1 0 0   1 1 0   1 1 1
 .05
DIFF
BR
15. .3
CONV
.41
END
PPAN
PATO
0 0
PPAN
ECHG
0
65
ATOMS
1  -1 1 1
1  0 0 0
1  1 1 -1
MARGINS
1.5 1.5 1.5 1.5
END
PROF
 2 0 1 0
 118 3
CP
 1  100.  100.  1
 1 0 0
 1.
BR
15. .3
END
PROF
 3 0 0 0
CP
 3  13.  13.  1
 1 0 0   1 1 0   1 1 1
 .05
DIFF
BR
15. .3
CONV
.41
END
END
ENDofDATA

	close(OUTCRYSTAL06);
	return 1;
}


sub SaveInputFiles
{
	my ($this, $Crystal, $Function, $SpinPolarized,
		$Dir, $IsChooseRandomly) = @_;
	my $LF = "<br>\n";
	
	my $SampleName = $this->SampleName();
	my $D12 = Deps::MakePath($Dir, "$SampleName.d12", 0);
	if(unlink($D12)) {
		print "$D12 was deleted.$LF";
	}
	my $D3  = Deps::MakePath($Dir, "$SampleName.d3", 0);
	if(unlink($D3)) {
		print "$D3 was deleted.$LF";
	}

	my $ret = $this->MakeD12File($Crystal, $Function, $SpinPolarized, $D12);
	return () if($ret <= 0);
	$ret = $this->MakeD3File($Crystal, $Function, $D3);
	return ($D12) if($ret <= 0);

	return ($D12, $D3);
}


1;