#!/usr/bin/perl

use lib 'd:/Programs/Perl/lib';
#use lib 'd:/MyWebs/cgi-bin/lib';

use lib '/home/tkamiya/bin/lib';
use lib '/home/tkamiya/bin';

use MyApplication;
use Crystal::CIF;
use Crystal::Crystal;
use Crystal::RietanFP;

use strict;
#use warnings;

#===============================================
# Applicationクラス作製
#===============================================
my $App = new MyApplication;
exit if($App->Initialize() < 0);

#===============================================
# コマンドラインオプション読み込み
#===============================================
$App->AddArgument("--TolD",                "--TolD=val: Torelance for atomic position (A) [Def 0.01 A]", "0.01");
$App->AddArgument("--xc0",                 "--xc0=val [Def 0.0]",                 "0.0");
$App->AddArgument("--yc0",                 "--yc0=val [Def 0.0]",                 "0.0");
$App->AddArgument("--zc0",                 "--zc0=val [Def 0.0]",                 "0.0");
$App->AddArgument("--ScanAll",             "--ScanAll=[0|1] [Def 0]",             "0");
$App->AddArgument("--ShowCrystalInf",      "--ShowCrystalInf=[0|1] [Def 0]",      "0");
$App->AddArgument("--ShowAllCoordination", "--ShowAllCoordination=[0|1] [Def 0]", "0");
$App->AddArgument("--RCoordMax",           "--RCoordMax=val [Def 7.0]",           "7.0");
$App->AddArgument("--ShowFoundStructure",  "--ShowFoundStructure=[0|1] [Def 0]",  "0");
$App->AddArgument("--help",                "--help: Show this help");
exit 1 if($App->ReadArgs() != 1);

my $OverlapCriteriaRadius = $App->GetGetArg('TolD');
   $OverlapCriteriaRadius = 0.01 if(!defined $OverlapCriteriaRadius);
my $xc0                   = $App->GetGetArg('xc0');
   $xc0                   = 0.0 if(!defined $xc0);
my $yc0                   = $App->GetGetArg('zc0');
   $yc0                   = 0.0 if(!defined $yc0);
my $zc0                   = $App->GetGetArg('zc0');
   $zc0                   = 0.0 if(!defined $zc0);
my $ScanAll               = $App->GetGetArg('ScanAll');
   $ScanAll               = 0 if(!defined $ScanAll);
my $ShowCrystalInf        = $App->GetGetArg('ShowCrystalInf');
   $ShowCrystalInf        = 0 if(!defined $ShowCrystalInf);
my $ShowAllCoordination   = $App->GetGetArg('ShowAllCoordination');
   $ShowAllCoordination   = 0 if(!defined $ShowAllCoordination);
my $RCoordMax             = $App->GetGetArg('RCoordMax');
   $RCoordMax             = 7.0 if(!defined $RCoordMax);
my $ShowFoundStructure    = $App->GetGetArg('ShowFoundStructure');
   $ShowFoundStructure    = 0 if(!defined $ShowFoundStructure);
my $CIFFile               = $App->GetGetArg(0);
   $CIFFile               = "NaCl-converted.cif" if(!defined $CIFFile);

#CIFクラスを取得
my $cif = new CIF;
#CIFファイル名を引数に与えて読み込む
unless($cif->Read($CIFFile)) {
	print "Error: Can not read [$CIFFile]\n\n";
	exit 2;
}
print "CIF                  : ", $cif->FileName(), "\n";
print "Initial center       : ($xc0, $yc0, $zc0)\n";
print "ShowCrystalInf       : $ShowCrystalInf\n";
print "ShowAllCoordination  : $ShowAllCoordination\n";
print "ShowFoundStructure   : $ShowFoundStructure\n";
print "ScanAll              : $ScanAll\n";
print "OverlapCriteriaRadius: $OverlapCriteriaRadius A\n";
print "RCoordMax            : $RCoordMax A\n";

&FindPossibleSPGs($OverlapCriteriaRadius);

exit 0;

#===============================================
# スクリプト終了
#===============================================

#==========================================
# サブルーチン
#==========================================
sub FindPossibleSPGs
{
	my ($OverlapCriteriaRadius) = @_;

#CIFクラスの内容から、Crystalクラスを作成
#作成した時点で、Metricsの計算、イオン情報(原子量など)の読み込み、
#格子体積の計算は終わっている
	my $Crystal = $cif->GetCCrystal();

	print "\n";
	my $CrystalName = $Crystal->CrystalName();
	print "CrystalName: $CrystalName\n";
	my ($a, $b, $c, $alpha, $beta, $gamma) = $Crystal->LatticeParameters();
	print "Lattice paramters: $a $b $c  $alpha $beta $gamma\n";
#	my ($g11,$g12,$g13,$g21,$g22,$g23,$g31,$g32,$g33) = $Crystal->Metrics();
#	print "Metrics: $g11 $g22 $g33  $g12 $g23 $g31\n";
#	my $vol = $Crystal->Volume();
#	print "Volume: $vol A^3\n";
#対称要素を使い、非対称単位中の原子を展開する
#この際、Densityも計算される
	$Crystal->ExpandCoordinates();
#	my $Density = $Crystal->Density();
#	print "Density: $Density g/cm^3\n";

#Crystalクラスから、SpaceGroupクラスを取り出す
#この時点で、並進対称要素の抽出、LatticeSystemの抽出は完了している
	my $SPG = $Crystal->GetCSpaceGroup();
	my $SPGName = $SPG->SPGName();
	my $iSPG    = $SPG->iSPG();
	my $LatticeSystem = $SPG->LatticeSystem();
	print "Space Group: $SPGName ($iSPG)\n";
	print "Lattice System: $LatticeSystem\n";

	my $nTranslation = $SPG->nTranslation();
	print "nTranslation: $nTranslation\n";
	if($ShowCrystalInf) {
		for(my $i = 0 ; $i < $nTranslation ; $i++) {
			my ($x,$y,$z) = $SPG->TranslationVector($i+1);
			print "   ($x, $y, $z)\n";
		}
		print "\n";
	}

	my $nSymmetryOperation = $SPG->nSymmetryOperation();
	print "nSymmetryOperation: $nSymmetryOperation\n";
	if($ShowCrystalInf) {
		for(my $i = 0 ; $i < $nSymmetryOperation ; $i++) {
			my $symop = $SPG->SymmetryOperation($i+1);
			print "   $symop\n";
	
			my ($x1,$y1,$z1,$t1,
				$x2,$y2,$z2,$t2,
				$x3,$y3,$z3,$t3)
				= $SPG->SymmetryOperationByMatrix($i+1);
			print "      $x1 $y1 $z1  $t1\n";
			print "      $x2 $y2 $z2  $t2\n";
			print "      $x3 $y3 $z3  $t3\n";
		}
		print "\n";
	}

#Crystalクラス中の、原子の種類 AtomTypeクラスのリストをとる
	my @AtomTypeList = $Crystal->GetCAtomTypeList();
	my $nAtomType = @AtomTypeList;
	print "nAtomType: $nAtomType\n";
	if($ShowCrystalInf) {
		for(my $i = 0 ; $i < $nAtomType ; $i++) {
			my $label        = $AtomTypeList[$i]->Label();
			my $atomname     = $AtomTypeList[$i]->AtomName();
			my $charge       = $AtomTypeList[$i]->Charge();
			my $AtomicNumber = $AtomTypeList[$i]->AtomicNumber();
			my $AtomicMass   = $AtomTypeList[$i]->AtomicMass();
			my $i1 = $i+1;
			print "   #$i1: $label : $atomname [$AtomicNumber] ($charge) ($AtomicMass)\n";
		}
		print "\n";
	}
	
#Crystalクラス中の、非対称単位中の原子 AsymmetricAtomSiteクラスのリストをとる
	my @AtomSiteList = $Crystal->GetCAsymmetricAtomSiteList();
	my $nAsymmetricAtomSite = @AtomSiteList;
	print "nAsymmetricAtomSite: $nAsymmetricAtomSite\n";
	if($ShowCrystalInf) {
		for(my $i = 0 ; $i < $nAsymmetricAtomSite ; $i++) {
			my $label     = $AtomSiteList[$i]->Label();
			my $type      = $AtomSiteList[$i]->AtomName();
			my ($x,$y,$z) = $AtomSiteList[$i]->Position(1);
			my $occupancy = $AtomSiteList[$i]->Occupancy();
			my $iAtomType = $Crystal->FindiAtomType($type);

			print "   $label ($type)[#$iAtomType]: ($x, $y, $z) [$occupancy]\n";
		}
	}

	my @ExpandedAtomSiteList = $Crystal->GetCExpandedAtomSiteList();
	my $nTotalExpandedAtomSite = $Crystal->nTotalExpandedAtomSite();
	print "nExpanededAtomSite: $nTotalExpandedAtomSite\n";
	if($ShowCrystalInf) {
		for(my $i = 0 ; $i < $nTotalExpandedAtomSite ; $i++) {
			my $label     = $ExpandedAtomSiteList[$i]->Label();
			my $type      = $ExpandedAtomSiteList[$i]->AtomName();
			my ($x,$y,$z) = $ExpandedAtomSiteList[$i]->Position(1);
			my $occupancy = $ExpandedAtomSiteList[$i]->Occupancy();
			my $iAtomType = $Crystal->FindiAtomType($type);

			print "   $label ($type)[#$iAtomType]: ($x, $y, $z) [$occupancy]\n";
		}
		print "\n";
	}

	my ($TSym, $pID) = &FindTranslations($Crystal);
	if($TSym ne '') {
		print "** Translation [$TSym] is found\n";
		if($ShowCrystalInf) {
			print "  Asymmetric sites:\n";
			my $c = 0;
			for(my $i = 0 ; $i < $nTotalExpandedAtomSite ; $i++) {
				next if($pID->[$i]);
	
				my $p         = $ExpandedAtomSiteList[$i];
				my $label     = $p->Label();
				my $type      = $p->AtomName();
				my ($x,$y,$z) = $p->Position(1);
				my $occupancy = $p->Occupancy();
				my $iAtomType = $Crystal->FindiAtomType($type);
				printf "   %3d: %4s (%2s)[#%02d]: (%8.5f, %8.5f, %8.5f) [%f]\n",
							$c, $label, $type, $iAtomType, $x, $y, $z, $occupancy;
				$c++;
			}
		}
		print "\n";
	}

	$TSym = '' if($ScanAll);

	my ($xc, $yc, $zc) = &FindOrigin($Crystal, $xc0, $yc0, $zc0);
	print "   Translation vector: (", -$xc, ", ", -$yc, ", ", -$zc, ")\n";

	for(my $i = 0 ; $i < $nTotalExpandedAtomSite ; $i++) {
		my $p              = $ExpandedAtomSiteList[$i];
		my $atomname       = $p->AtomNameOnly();
		my ($x, $y, $z)    = $p->Position();
		my ($x1, $y1, $z1) = ($x-$xc, $y-$yc, $z-$zc);
		$p->SetPosition($x1, $y1, $z1);
		printf "  Atom[%3d]: %2s: (%8.5f, %8.5f, %8.5f) to (%8.5f,%8.5f,%8.5f)\n",
					$i, $atomname, $x, $y, $z,$x1,$y1,$z1;#  if($ShowCrystalInf);
	}
	$Crystal->SetCAsymmetricAtomSiteList(@ExpandedAtomSiteList);
#	@ExpandedAtomSiteList = $Crystal->GetCExpandedAtomSiteList();

	my ($drive, $dir, $filename, $ext, $lastdir, $filebody) = Deps::SplitFilePath($CIFFile);
	my $TranslatedCIFFile = "$filebody-Translated.cif";
	if(!CIF->new()->CreateCIFFileFromCCrystal($Crystal, $TranslatedCIFFile, 0, 'unix')) {
		$App->print("Error: Can not create [$CIFFile].\n");
		return 0;
	}

	for(my $iSPG = 230 ; $iSPG >= 1 ; $iSPG--) {
		for(my $iSet = 1 ; ; $iSet++) {
			my ($SPG2, $IsMatched, $mess, $pSiteId, $pAsymSites) 
					= &CheckSPG($Crystal, $TSym, $pID, $iSPG, $iSet, $OverlapCriteriaRadius, 0);
			last if(!defined $SPG2);

			if($IsMatched) {
				my $SPGName2           = $SPG2->SPGName();
				my $LatticeSystem2     = $SPG2->LatticeSystem();
				my $nTranslation       = $SPG2->nTranslation();
				my $nSymmetryOperation = $SPG2->nSymmetryOperation();
				my $nTotalSymmetry     = $nTranslation * $nSymmetryOperation;

				print "SPG $SPGName2 (#$iSPG-$iSet LS=$LatticeSystem2) is matched\n";
				print "   nSym=$nSymmetryOperation  nTrans=$nTranslation\n";
				my $nASSite = @$pAsymSites;
				print "   $nASSite asymmetric sites are found.\n";

				my $NewCrystal = new Crystal();
				$NewCrystal->SetCrystalName($CrystalName);
				$NewCrystal->SetLatticeParameters($a, $b, $c, $alpha, $beta, $gamma);
				$NewCrystal->SetCSpaceGroup($SPG2);

				for(my $i = 0 ; $i < @$pAsymSites ; $i++) {
#print "SiteId=$pSiteId->[$i]\n";
#					next if(defined $pSiteId->[$i]);

					my $site      = $pAsymSites->[$i];
					my $label     = $site->Label();
					my $atomname  = $site->AtomName();
					my ($x,$y,$z) = $site->Position(1);
					my $occupancy = $site->Occupancy();
					print "   $i: $atomname: ($x, $y, $z) occ=$occupancy\n" if($ShowFoundStructure);

					$NewCrystal->AddAtomSite($label, $atomname, $x, $y, $z, $occupancy);
				}

				my $NewCIFName = $CIFFile;
				$NewCIFName =~ s/\.cif//i;
				$SPGName2 =~ s/\s+//g;
				$NewCIFName .= "-$SPGName2.cif";
				print "** Save CIF file to [$NewCIFName].\n";
				CIF->new()->CreateCIFFileFromCCrystal($NewCrystal, $NewCIFName, undef, undef, $ShowFoundStructure);
			}
			else {
#				print "SPG $SPGName2 (#$iSPG-$iSet LS=$LatticeSystem2) is not matched [$mess]\n";
			}
		}
	}

	return;
}

sub FindOrigin
{
	my ($Crystal, $xc, $yc, $zc, $IsPrint) = @_;

	my $Rmax = $RCoordMax;
	my $kRth = 0.1;

	print "** Find origin from initial center ($xc, $yc, $zc)\n";

	my @MArrays = (
		[ [-1,  0,  0, 0],
		  [ 0, -1,  0, 0],
		  [ 0,  0, -1, 0] ],
		[ [-1,  0,  0, 0],
		  [ 0,  1,  0, 0],
		  [ 0,  0,  1, 0] ],
		[ [ 1,  0,  0, 0],
		  [ 0, -1,  0, 0],
		  [ 0,  0,  1, 0] ],
		[ [ 1,  0,  0, 0],
		  [ 0,  1,  0, 0],
		  [ 0,  0, -1, 0] ],
		[ [-1,  0,  0, 0],
		  [ 0, -1,  0, 0],
		  [ 0,  0,  1, 0] ],
		[ [-1,  0,  0, 0],
		  [ 0,  1,  0, 0],
		  [ 0,  0, -1, 0] ],
		[ [ 1,  0,  0, 0],
		  [ 0, -1,  0, 0],
		  [ 0,  0, -1, 0] ],
		[ [ 0,  1,  0, 0],
		  [-1,  0,  0, 0],
		  [ 0,  0,  1, 0] ],
		[ [ 0,  0, -1, 0],
		  [ 0,  1,  0, 0],
		  [ 1,  0,  0, 0] ],
		[ [ 1,  0,  0, 0],
		  [ 0,  0,  1, 0],
		  [-1,  0,  0, 0] ],
		[ [ 1,  0,  0, 0],
		  [ 0, -1,  0, 0],
		  [ 0,  0, -1, 0] ],
		[ [ 0, -1,  0, 0],
		  [ 1,  0,  0, 0],
		  [ 0,  0,  1, 0] ],
		[ [ 0,  0,  1, 0],
		  [ 0,  1,  0, 0],
		  [-1,  0,  0, 0] ],
		[ [ 1,  0,  0, 0],
		  [ 0,  0, -1, 0],
		  [ 1,  0,  0, 0] ],
		);
	my @MLabel = (
		"i",
		"mx",
		"my",
		"mz",
		"2z",
		"2y",
		"2x",
		"4z",
		"4y",
		"4x",
		"4*3z",
		"4*3y",
		"4*3x",
		);

	my ($la, $lb, $lc, $alpha, $beta, $gamma) = $Crystal->LatticeParameters();
	my @AtomSiteList = $Crystal->GetCExpandedAtomSiteList();
	my $nTotalSite   = $Crystal->nTotalExpandedAtomSite();

if(0) {
	my (@site, %pSiteArray);
	for(my $i = 0 ; $i < $nTotalSite ; $i++) {
		my $site      = $AtomSiteList[$i];
		my $atomname  = $site->AtomNameOnly();
		my ($x,$y,$z) = $site->Position(1);
		my $occupancy = $site->Occupancy();
		my $ID        = sprintf("%s%04d", $atomname, int($occupancy*1000.0+1e-3));

		my %h = (
			ID          => $ID,
			'name'      => $atomname,
			'x'         => $x,
			'y'         => $y,
			'z'         => $z,
			'occupancy' => $occupancy,
			);

		$pSiteArray{$ID} = [] if(!$pSiteArray{$ID});
		my $p = $pSiteArray{$ID};
		push(@$p, \%h);
#		print "   $i: $ID ($x, $y, $z)\n";
	}

	foreach my $key (sort keys %pSiteArray) {
		my $p = $pSiteArray{$key};
		print "$key: \n";
		for(my $i = 0 ; $i < @$p ; $i++) {
			print "   ($p->[$i]{x}, $p->[$i]{y}, $p->[$i]{z}) occ=$p->[$i]{occupancy}\n";
		}
	}
}

	my %pRArray;

	my $nxmax = int($Rmax / $la) + 1;
	my $nymax = int($Rmax / $lb) + 1;
	my $nzmax = int($Rmax / $lc) + 1;
	for(my $nz = -$nzmax ; $nz <= $nzmax ; $nz++) {
	for(my $ny = -$nymax ; $ny <= $nymax ; $ny++) {
	for(my $nx = -$nxmax ; $nx <= $nxmax ; $nx++) {
		for(my $i = 0 ; $i < $nTotalSite ; $i++) {
			my $site      = $AtomSiteList[$i];
			my $atomname  = $site->AtomNameOnly();
			my ($x,$y,$z) = $site->Position(1);
			my $occupancy = $site->Occupancy();
			my $ID        = sprintf("%s%04d", $atomname, int($occupancy*1000.0+1e-3));

			my $r = $Crystal->GetInterAtomicDistance($xc, $yc, $zc, $x+$nx, $y+$ny, $z+$nz);
			next if($r > $Rmax);

			my %h = (
				ID          => $ID,
				'name'      => $atomname,
				'x'         => $x+$nx,
				'y'         => $y+$ny,
				'z'         => $z+$nz,
				'occupancy' => $occupancy,
				'r'         => $r,
				);

			$pRArray{$ID} = [] if(!$pRArray{$ID});
			my $p = $pRArray{$ID};
			push(@$p, \%h);
		}
	}
	}
	}

	my ($dx, $dy, $dz) = (0, 0, 0);
	my $c = 0;
	foreach my $key (sort keys %pRArray) {
		my $p = $pRArray{$key};
		@$p = sort { $a->{r} <=> $b->{r} } @$p;

		my ($i0, $Rmin, $Rth, $n1stCoord);
		if($p->[0]{r} < 3.0*$kRth) {
			$Rmin = $p->[1]{r}-$p->[0]{r};
			$i0 = 2;
		}
		else {
			$Rmin = $p->[0]{r};
			$i0 = 1;
		}
		$Rth  = $Rmin * $kRth;
		my ($xc1, $yc1, $zc1) = (0.0, 0.0, 0.0);
		for(my $i = 0 ; $i < @$p ; $i++) {
			if($i >= $i0 and $p->[$i]{r} - $p->[$i-1]{r} > $Rth) {
				$n1stCoord = $i;
				last;
			}
#			($xc1, $yc1, $zc1) += ($p->[$i-1]{x}, $p->[$i-1]{y}, $p->[$i-1]{z});
			$xc1 += $p->[$i]{x};
			$yc1 += $p->[$i]{y};
			$zc1 += $p->[$i]{z};
		}

#		($xc1, $yc1, $zc1) /= $n1stCoord;
		$xc1 /= $n1stCoord;
		$yc1 /= $n1stCoord;
		$zc1 /= $n1stCoord;

		$dx += $xc1;
		$dy += $yc1;
		$dz += $zc1;
		$c++;
		
		printf "%7s: # of atoms in 1st coordination = %2d: Center (%8.5f, %8.5f, %8.5f)\n", 
			$key, $n1stCoord, $xc1, $yc1, $zc1;
		if($ShowAllCoordination) {
			for(my $i = 0 ; $i < @$p ; $i++) {
				printf "   $i: (%8.5f, %8.5f, %8.5f)) r=%8.5f  occ=%f\n",
							$p->[$i]{x}, $p->[$i]{y}, $p->[$i]{z}, $p->[$i]{r}, $p->[$i]{occupancy};
			}
		}
	}
	print "\n";

	$dx /= $c;
	$dy /= $c;
	$dz /= $c;
#exit;

	return ($dx, $dy, $dz);
}


sub FindTranslations
{
	my ($Crystal, $IsPrint) = @_;
	
	my @TArrays = (
		[0.5, 0.5, 0.0],
		[0.5, 0.0, 0.5],
		[0.0, 0.5, 0.5],
		[0.5, 0.5, 0.5],
		);
	my @TLabel = (
		"C",
		"B",
		"A",
		"I",
		);

	my (@hit, $pID, @ID);
	for(my $i = 0 ; $i < @TArrays ; $i++) {
		($hit[$i], $pID) = &CheckTransration($Crystal, $TArrays[$i], $IsPrint);
		if($hit[$i]) {
#			print "** Translation [$TLabel[$i]] is found\n";
			for(my $j = 0 ; $j < @$pID ; $j++) {
				$ID[$j] = 1 if($pID->[$j]);
			}
		}
	}

	if($hit[0] and $hit[1] and $hit[2]) {
		print "** Translation [F] is found\n" if($IsPrint);
		return ("F", \@ID);
	}
	for(my $i = 0 ; $i < @TArrays ; $i++) {
		if($hit[$i]) {
			print "** Translation [$TLabel[$i]] is found\n" if($IsPrint);
			return ($TLabel[$i], \@ID);
		}
	}
	return ('', \@ID);
}

sub CheckTransration
{
	my ($Crystal, $pT, $IsPrint) = @_;
#$IsPrint=1;

	my @SiteList = $Crystal->GetCExpandedAtomSiteList();
	my $nSite = @SiteList;

	my @SiteId;
	for(my $i = 0 ; $i < $nSite ; $i++) {
		next if($SiteId[$i]);

		my $p         = $SiteList[$i];
		my $label     = $p->Label();
		my $atomname  = $p->AtomName();
		my ($x,$y,$z) = $p->Position(1);
		my $occupancy = $p->Occupancy();
print "   $i: $label ($atomname): ($x, $y, $z) [$occupancy]\n" if($IsPrint);

		my $xs = $x + $pT->[0];
		my $ys = $y + $pT->[1];
		my $zs = $z + $pT->[2];
		my $id2 = $Crystal->FindIdenticalSite($atomname, $xs, $ys, $zs, $OverlapCriteriaRadius, $i);

		if(defined $id2) {
print " Tr Op ($pT->[0], $pT->[1], $pT->[2]): ($xs,$ys,$zs): atom site $i is equivalent to site $id2\n" if($IsPrint);
			$SiteId[$id2] = 1 if($i < $id2);
			next;
		}
		else {
print " Tr Op ($pT->[0], $pT->[1], $pT->[2]): ($xs,$ys,$zs): Equivalent site was not found for atom site $i\n" if($IsPrint);
			return (0);
		}
	}
	
	return (1, \@SiteId);
}

sub CheckSPG
{
	my ($Crystal, $TSym, $pID, $iSPG, $iSet, $OverlapCriteriaRadius, $IsPrint) = @_;

	my ($a,$b,$c,$alpha,$beta,$gamma) = $Crystal->LatticeParameters();
#	$Crystal->ExpandCoordinates();
	my $SPG = $Crystal->GetCSpaceGroup();
#	my $SPGName = $SPG->SPGName();
#	my $iSPG    = $SPG->iSPG();
	my $LatticeSystem = $SPG->LatticeSystem();
#	print "Space Group: $SPGName ($iSPG)\n";
#	print "Lattice System: $LatticeSystem\n";
	my @ExpandedAtomSiteList = $Crystal->GetCExpandedAtomSiteList();
	my $nTotalExpandedAtomSite = $Crystal->nTotalExpandedAtomSite();
	my @nMultiplicityExpandedAtomSiteList = $Crystal->GetCnMultiplicityExpandedAtomSiteList();
#	print "nTotalExpandedAtomSite: $nTotalExpandedAtomSite\n";

	my $r              = new RietanFP;
	my $SPG2           = $r->ReadSpaceGroup($iSPG, $iSet);
#print "i=$iSPG, $iSet SPG2=$SPG2\n";
	return (undef) if(!$SPG2);
	my $SPGName        = $SPG2->SPGName();
	if($TSym ne '' and $SPGName !~ /^$TSym/i) {
		return ($SPG2, 0, "[$SPGName] in SPG is not compatible with translation symmetry [$TSym]");
	}
	my $LatticeSystem2 = $SPG2->LatticeSystem();

#print "LS: [$LatticeSystem][$LatticeSystem2]\n";
	if(!&IsCompatibleLatticeSystem($LatticeSystem, $LatticeSystem2)) {
		return ($SPG2, 0, "[$LatticeSystem2] in SPG is not compatible with [$LatticeSystem]");
	}

#非対称構造中のパラメータ
	my $nTranslation       = $SPG2->nTranslation();
	my $nSymmetryOperation = $SPG2->nSymmetryOperation();
	my $nTotalSymmetry     = $nTranslation * $nSymmetryOperation;
print "nTotalSymmetry: $nTotalSymmetry = $nSymmetryOperation * $nTranslation\n" if($IsPrint);

	my @SiteId;
	for(my $i = 0 ; $i < $nTotalExpandedAtomSite ; $i++) {
		next if($pID->[$i]);
		next if($SiteId[$i]);

		my $label     = $ExpandedAtomSiteList[$i]->Label();
		my $atomname  = $ExpandedAtomSiteList[$i]->AtomName();
		my ($x,$y,$z) = $ExpandedAtomSiteList[$i]->Position(1);
		my $occupancy = $ExpandedAtomSiteList[$i]->Occupancy();
		my $id        = $ExpandedAtomSiteList[$i]->IdAsymmetricAtomSite();
		my $mult      = $nMultiplicityExpandedAtomSiteList[$id-1];
print "   $i: [$id]$label ($atomname): ($x, $y, $z) [$occupancy][$mult]\n" if($IsPrint);

		my $Found = 0;
		for(my $is = 0 ; $is < $nTotalSymmetry ; $is++) {
			my ($xs, $ys, $zs) = $SPG2->DoSymmetryOperation($is, $x, $y, $z, 1);
#print "        ($is) xs,ys,zs: $xs,$ys,$zs<br>\n";
#			my ($fxs, $fys, $fzs) = $SPG->DoVelocitySymmetryOperation($is, $fx, $fy, $fz, 1);
#			my ($vxs, $vys, $vzs) = $SPG->DoVelocitySymmetryOperation($is, $vx, $vy, $vz, 1);

#			my $id2 = $Crystal->FindIdenticalSite($atomname, $xs, $ys, $zs, $OverlapCriteriaRadius);
			my $id2 = $Crystal->FindIdenticalSite($atomname, $xs, $ys, $zs, $OverlapCriteriaRadius, $i);
			if(defined $id2) {
print " Sym Op $is ($xs,$ys,$zs): atom site $i is equivalent to site $id2\n" if($IsPrint);
				$SiteId[$id2] = 1 if($i < $id2);
				$Found = 1;
				next;
			}
			else {
print " Sym Op $is ($xs,$ys,$zs): Equivalent site was not found for atom site $i\n" if($IsPrint);
				return ($SPG2, 0, "No equivalent site found for site $i");
			}
		}
		if(!$Found) {
print "   .Equivalent site was not found for atom site $i\n" if($IsPrint);
			return ($SPG2, 0, "No equivalent site found for site $i");
		}
	}

	my @AsymSites;
	for(my $i = 0 ; $i < $nTotalExpandedAtomSite ; $i++) {
		next if(defined $SiteId[$i]);
		push(@AsymSites, $ExpandedAtomSiteList[$i]);
	}

	return ($SPG2, 1, "OK", \@SiteId, \@AsymSites);
}

sub IsCompatibleLatticeSystem
{
	my ($CrystalLS, $SPGLS) = @_;

	if($CrystalLS =~ /cub/i) {
		return 1;
	}
	if($CrystalLS =~ /tet/i) {
		return 0 if($SPGLS =~ /cub/i);
		return 1;
	}
	if($CrystalLS =~ /ort/i) {
		return 0 if($SPGLS =~ /cub/i);
		return 0 if($SPGLS =~ /tet/i);
		return 1;
	}
	if($CrystalLS =~ /hex/i) {
		return 0 if($SPGLS =~ /cub/i);
		return 0 if($SPGLS =~ /tet/i);
		return 0 if($SPGLS =~ /ort/i);
		return 1;
	}
	if($CrystalLS =~ /trig/i) {
		return 0 if($SPGLS =~ /cub/i);
		return 0 if($SPGLS =~ /tet/i);
		return 0 if($SPGLS =~ /ort/i);
		return 1;
	}
	if($CrystalLS =~ /rho/i) {
		return 0 if($SPGLS =~ /cub/i);
		return 0 if($SPGLS =~ /tet/i);
		return 0 if($SPGLS =~ /ort/i);
		return 1;
	}
	if($CrystalLS =~ /mon/i) {
		return 0 if($SPGLS =~ /cub/i);
		return 0 if($SPGLS =~ /tet/i);
		return 0 if($SPGLS =~ /ort/i);
		return 0 if($SPGLS =~ /hex/i);
		return 0 if($SPGLS =~ /trig/i);
		return 0 if($SPGLS =~ /rho/i);
		return 1;
	}
	return 1 if($SPGLS =~ /trig/i or $SPGLS eq '');
	return 0;
}