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

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

use strict;
use File::Basename;

use Deps;
use Utils;
use Sci::Science;

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

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


#===============================================
# パス（読み込みDB）
# Web関係変数
# CGI の仮想パス
# プログラム名
#===============================================
my $Program    = basename($0);
my $ProgramDir = ProgVars->ProgramDir();
my $LDDir      = ProgVars->LDDir();
my $MXDDir     = ProgVars->MXDDir();
my $MXDOrtoDir = ProgVars->MXDOrtoDir();
my $MXDExecDir = ProgVars->MXDExecDir();
my $SX1DBPath  = ProgVars->SX1DBPath();

BEGIN {
}

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

sub DESTROY
{
	my $this = shift;
}

sub SetSampleName
{
	my ($this, $name) = @_;
	return $this->{'SampleName'} = $name;
}

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

sub GetBusingParameter
{
	my ($this,$atomname) = @_;
	
	$atomname =~ s/\{.*\}//g;
	unless(open(IN, "<$SX1DBPath")) {
		print "Error in MXD.pm::GetBusingParameter: "
			."Can not read DB [$SX1DBPath].<br>\n";
		return ('');
	}
	<IN>;
	while(<IN>) {
		my ($name, $mass, $charge, $ai, $bi, $ci, $rad)
			= split(/\s+/, $_);
#print "name: $name, $atomname ($mass, $charge, $ai, $bi) <br>\n";
		if(uc $name eq uc $atomname) {
			return ($mass, $charge, $ai, $bi, $ci, $rad);
		}
	}
	
	close(IN);
	return ('');
}

sub MakeLDFiles
{
	my ($this, $Crystal, $Function, $fname, $pModifiedParams, $EWALD, $EPS) = @_;

	my $LDDATPath = "$fname.dat";
	my $LD01Path  = "$fname.f01";
	my $LD02Path  = "$fname.f02";
	my $LD07Path  = "$fname.f07";
	if($pModifiedParams) {
		$LDDATPath = "$fname-modified.dat";
		$LD02Path  = "$fname-modified.f02";
	}

#	print("($LDDATPath) is deleted.\n") if(unlink($LDDATPath));
#	print("($LD01Path) is deleted.\n")  if(unlink($LD01Path));
#	print("($LD02Path) is deleted.\n")  if(unlink($LD02Path));
#	print("($LD07Path) is deleted.\n")  if(unlink($LD07Path));

	if($this->MakeLD01File($Crystal, $Function, $LD01Path) <= 0) {
		print("Error in LD::MakeLDFiles: Can not write to [$LD01Path].\n");
		return 0;
	}
	if($this->MakeLD02File($Crystal, $Function, $LD02Path, $pModifiedParams) <= 0) {
		print("Error in LD::MakeLDFiles: Can not write to [$LD02Path].\n");
		return 0;
	}
	if($this->MakeLD07File($Crystal, $Function, $LD07Path, $EWALD, $EPS) <= 0) {
		print("Error in LD::MakeLDFiles: Can not write to [$LD07Path].\n");
		return 0;
	}

#ファイル名リストのファイル .DAT をつくる
	print("\n");
	print("Create [$LDDATPath]\n");
	unless(open(OUT,">$LDDATPath")) {
		print "Can't write to $LDDATPath.\n\n";
		return;
	}
	print OUT "[Files]\n";
	print OUT "$LD01Path\n";
	print OUT "$LD02Path\n";
	print OUT "$fname.f04\n";
	print OUT "$LD07Path\n";
	print OUT "$fname.f08\n";
	print OUT "$fname.f09\n";
	if($Function =~ /Phonon/i) {
		print OUT "\n";
		print OUT "[General]\n";
		print OUT "0.4    EwaldAlpha\n";
		print OUT "0.8    MinIonicRadius\n";
		print OUT "1.0e-6 ConvergenceEPS\n";
		print OUT "\n";
		print OUT "[CalRange for Phonon]\n";
		print OUT "Execute\n";
		print OUT "\n";
		print OUT "[CalDistance]\n";
		print OUT "6.0    MaxDistance\n";
		print OUT "Execute\n";
		print OUT "\n";
		print OUT "[SaveCrystal]\n";
		print OUT "Execute\n";
		print OUT "\n";
		print OUT "[EnergyTest]\n";
		print OUT "NotExecute\n";
		print OUT "\n";
		print OUT "[CalEnergies]\n";
		print OUT "NotExecute\n";
		print OUT "\n";
		print OUT "[Phonon]\n";
		print OUT "0       fPhononDebug\n";
		print OUT "1       fUsedu/dx\n";
		print OUT "0       iPrintMatrixLevel=0: none 1: Wij 2: Wij/Gij "
			 ."3: Wij/Gij/Normalized Wij\n";
		print OUT "0       iEigenSolver=0: ZHEGV(Hermitian) "
			 ."1: ZGGEV(Asymmetric) 2: EigAB(Real)\n";
		print OUT "2       iPrintEigenVectorLevel=0: none 1: Complex "
			 ."2: Norm 3: Complex/Norm\n";
		print OUT "10      nPrintEigenVectors\n";
		print OUT "3.0e-5  500 500 500      CalPhononPrecision  "
			 ."nULimitNRange(3)\n";
		print OUT "1.0e-14 EigCHEPS (for iEigenSolver=2)\n";
		print OUT "3.0 3.0 wWindowFunc, nSigmaUsed\n";
		print OUT "2\n";
		print OUT "0.0 1.0 0.0 0.0 0.0 0.0  21\n";
		print OUT "0.0 0.0 0.0 0.0 0.0 1.0  21\n";
		print OUT "Execute\n";
		print OUT "0.0 1.0 0.0 1.0 0.0 0.0  21\n";
		print OUT "0.0 1.0 0.0 1.0 0.0 1.0  21\n";
	}
	close(OUT);

	return 1;
}

sub MakeLD01File
{
	my ($this, $Crystal, $Function, $LD01Path) = @_;

#空間群ファイルを作る　FILE01.DAT
	return 0 unless(open(OUT,">$LD01Path"));

	my $SPG = $Crystal->GetCSpaceGroup();
	my $SPGName = $Crystal->SPGName();
	$SPGName =~ s/[\s']//g;
	print OUT "$SPGName\n";

	my @SymmetryOperation = $SPG->ExpandSymmetryOperation();
	my $nSymmetryOperation = @SymmetryOperation;
	print OUT "0 1 $nSymmetryOperation\n";
	print OUT " 0.0000  0.0000  0.0000\n";

	for(my $i = 0 ; $i < $nSymmetryOperation ; $i++) {
		my $SymOp = $SymmetryOperation[$i];
		Utils::DelSpace($SymOp);

		my @SymOpXYZ = split(/,\s*/, $SymOp);
		my @XMatrix = SpaceGroup::SymOpToMatrix($SymOpXYZ[0]);
		my @YMatrix = SpaceGroup::SymOpToMatrix($SymOpXYZ[1]);
		my @ZMatrix = SpaceGroup::SymOpToMatrix($SymOpXYZ[2]);
		printf OUT "%7.4f %7.4f %7.4f %7.4f %7.4f %7.4f "
			."%7.4f %7.4f %7.4f\n",
			$XMatrix[0], $XMatrix[1], $XMatrix[2],
			$YMatrix[0], $YMatrix[1], $YMatrix[2],
			$ZMatrix[0], $ZMatrix[1], $ZMatrix[2];
		printf OUT "%7.4f %7.4f %7.4f\n",
			eval($XMatrix[3]), eval($YMatrix[3]), 
			eval($ZMatrix[3]);
	}
	print OUT "\n";
	close(OUT);

	return 1;
}

sub MakeLD02File
{
	my ($this, $Crystal, $Function, $LD02Path, $pModifiedParams) = @_;

#原子ファイルを作る　FILE02.DAT
	return 0 unless(open(OUT,">$LD02Path"));

	my @AtomTypeList    = $Crystal->GetCAtomTypeList();
	my $nAtomType = @AtomTypeList;
	print OUT "$nAtomType\n";

	for(my $i = 0 ; $i < $nAtomType ; $i++) {
		my $atom      = $AtomTypeList[$i];
		my $atomname  = $atom->AtomNameOnly();
		my $charge0   = $atom->Charge();
		my ($zc, $sign) = ($charge0 =~ /([\d\.]+)([+\-])/);
		if($sign eq '-') {
			$charge0 = -$zc;
		}

		my ($Mass, $charge, $ai, $bi, $ci, $rad)
			= $this->GetBusingParameter($atomname);
		$Mass = AtomType::GetAtomInformationByKey($atomname, "原子量")
			if($Mass eq '');
		$Mass =~ s/\(.*?\)//g;
		
		$charge = $charge0 if($charge0 ne '');

		$ai = 0.9  if($ai == 0);
		$bi = 0.08 if($bi == 0);
		$ci = 0.0  if($ci eq '');

		my $iMP;
		if($pModifiedParams) {
			for(my $j = 0 ; $j < @$pModifiedParams ; $j++) {
				if($atomname eq $pModifiedParams->[$j]{Name}) {
					$iMP = $j;
					last;
				}
			}
#print "$atomname => index of Bader charge atom: $iMP\n";
		}

		if(defined $iMP) {
			my $p = $pModifiedParams->[$iMP];
			if(defined $p->{Charge}) {
				$charge = $p->{Charge};
			}
			if(defined $p->{ai}) {
				$ai = $p->{ai};
			}
			if(defined $p->{bi}) {
				$bi = $p->{bi};
			}
			if(defined $p->{ci}) {
				$ci = $p->{ci};
			}
			if(defined $p->{rad}) {
				$rad = $p->{rad};
			}
		}

		print      "$atomname\n";
		print  OUT "$atomname\n";
		printf     "%8.4f  %8.4f %9.6f %9.6f %9.6f\n", eval($Mass), $charge, $ai, $bi, $ci;
		printf OUT "%8.4f  %8.4f %9.6f %9.6f %9.6f\n", eval($Mass), $charge, $ai, $bi, $ci;
	}
#Aij, Bij, Cij and Morse, Shell parameters
	for(my $i = 1 ; $i <= $nAtomType ; $i++) {
		for(my $j = $i ; $j <= $nAtomType ; $j++) {
			print OUT " 0.0  0.0    0.000000E+00  0.000000E+00  0.020000E+00  0.000000E+00 0.0e0\n";
 		}
	}
#ai, bi IDs
	for(my $i = 0 ; $i < $nAtomType ; $i++) {
		print OUT "   1   0\n";
	}
#Spring constants
	print OUT "2\n";
	print OUT "  2.200  3.400  4.200  4.300\n";
	for(my $i = 1 ; $i <= $nAtomType ; $i++) {
		for(my $j = $i ; $j <= $nAtomType ; $j++) {
			print OUT "8.6e-19  1.6e-19  0.0e-19\n";
		}
	}

	print OUT "\n";
	close(OUT);

	return 1;
}

sub MakeLD07File
{
	my ($this, $Crystal, $Function, $LD07Path, $EWALD, $EPS) = @_;

#結晶構造ファイルを作る　FILE07.DAT
	unless(open(OUT,">$LD07Path")) {
		print "Can't write to $LD07Path.$LF$LF";
		return;
	}

	print OUT "1\n";
	my $SampleName = $this->SampleName();
	print OUT "$SampleName\n";

	my ($a,$b,$c,$alpha,$beta,$gamma) 
		= $Crystal->LatticeParametersByOutputMode(0);
	printf OUT " %9.4f %9.4f %9.4f  %9.4f %9.4f %9.4f\n",
		$a, $b, $c, $alpha, $beta, $gamma;

	my $SPG = $Crystal->GetCSpaceGroup();
	my $SPGName = $SPG->SPGName();
	$SPGName =~ s/\s//g;
	print OUT "$SPGName\n";

	my @AsymmetricAtomSiteList = $Crystal->GetCAsymmetricAtomSiteList();
	my $nAsymmetricAtomSite    = @AsymmetricAtomSiteList;
	print OUT "$nAsymmetricAtomSite\n";
#	my @ExpandedAtomSiteList = $Crystal->GetCExpandedAtomSiteListByOutputMode();
#	my $nExpandedAtomSite    = @ExpandedAtomSiteList;
#	print OUT "$nExpandedAtomSite\n";

	for(my $i = 0 ; $i < $nAsymmetricAtomSite ; $i++) {
#	for(my $i = 0 ; $i < $nExpandedAtomSite ; $i++) {
#		my $atom      = $ExpandedAtomSiteList[$i];
		my $atom      = $AsymmetricAtomSiteList[$i];
		my $atomname  = $atom->AtomNameOnly();
		my $iAtom     = $Crystal->FindiAtomType($atomname);
#		my $charge    = $atom->Charge();
		my ($x,$y,$z) = $atom->Position(1);
		my $occupancy = $atom->Occupancy();
#		my $mult      = $atom->Multiplicity();
#		$mult = 1 if($IsExpandCoordinate or $IsChooseRandomly);

		printf OUT "%5i  %7.4f  %7.4f  %7.4f  %7.4f\n",
			$iAtom, $occupancy, $x, $y, $z;
	}
# LIDs
	print OUT "  1  1  1  0  0  0\n";
#XIDs
	for(my $i = 0 ; $i < $nAsymmetricAtomSite ; $i++) {
#	for(my $i = 0 ; $i < $nExpandedAtomSite ; $i++) {
		print OUT "  1  1  1\n";
	}

	$EWALD = 0.4    if($EWALD eq '');
	$EPS   = 1.0e-5 if($EPS eq '');
	print OUT "$EWALD\n";
	print OUT "0.8     $EPS  1.0E-5\n";
	print OUT "1 3.0 1\n";
	print OUT "  10  0.100000E-05   1   0\n";

	print OUT "\n";
	close(OUT);

	return 1;
}

1;