package CrystalObject;
use Exporter;
@ISA = qw(Exporter);
#公開したいサブルーチン
@EXPORT = qw();
use strict;
#use Math::Matrix;
use Math::MatrixReal;
use Math::Vector::Real;
use Sci qw($pi $a0 $torad $todeg);# acos asin);
use Sci::Algorism;
use Sci::MyMatrixReal;
use Math::Complex;
use Crystal::PointGroup;
use Crystal::SpaceGroup;
use Crystal::AtomType;
use Crystal::AtomSite;
#use Crystal::Rietan;
use Crystal::RietanFP;
#==================================================================
# 静的メンバー関数
#==================================================================
sub RoundSymmetricPosition
{
my ($x) = @_;
return 1.0 / 3.0 if($x == 0.3333 or $x == 0.333);
return -1.0 / 3.0 if($x == -0.3333 or $x == -0.333);
return 2.0 / 3.0 if($x == 0.6667 or $x == 0.667);
return -2.0 / 3.0 if($x == -0.6667 or $x == -0.667);
return 1.0 / 6.0 if($x == 0.1667 or $x == 0.167);
return -1.0 / 6.0 if($x == -0.1667 or $x == -0.167);
return 5.0 / 6.0 if($x == 0.8333 or $x == 0.833);
return -5.0 / 6.0 if($x == -0.8333 or $x == -0.833);
return $x;
}
sub PolarizationFactor
{
my ($Theta, $MonochromaterAlpha) = @_;
my $cos2a = 1.0;
if(defined $MonochromaterAlpha) {
$cos2a = cos(2.0 * $MonochromaterAlpha * $torad);
$cos2a = $cos2a * $cos2a;
}
my $cos2Q = cos($torad * 2.0 * $Theta);
return 1.0 + $cos2a * $cos2Q * $cos2Q;
}
sub LorentzFactor
{
my ($Theta) = @_;
my $cosQ = cos($torad * $Theta);
my $sinQ = sin($torad * $Theta);
# my $sinQ = sin($torad * 2.0 * $Theta);
return 1.0 / $cosQ / $sinQ / $sinQ;
# my $cosQ = cos($torad * $Theta);
my $sin2Q = sin($torad * 2.0 * $Theta);
return $cosQ / $sin2Q / $sin2Q;
}
sub TemperatureFactor
{
# Biso in Angstrom, $s in nm^-1
my ($Biso, $s) = @_;
$s = 0.1 * $s;
my $M = $Biso * $s * $s;
return exp(-$M);
}
#==================================================================
# 変数取得関数
#==================================================================
sub SampleName { return shift->{SampleName}; }
sub SetSampleName { my ($this,$s)=@_; return $this->{SampleName} = $s; }
sub SetCrystalName { my ($this,$name) = @_; return $this->{CrystalName} = $name; }
sub CrystalName { my ($this) = @_; return $this->{CrystalName}; }
sub SPGName { my ($this) = @_; my $SPG = $this->GetCSpaceGroup(); return $SPG->SPGName(); }
sub iSPG { my ($this) = @_; my $SPG = $this->GetCSpaceGroup(); return $SPG->iSPG(); }
sub iSet { my ($this) = @_; my $SPG = $this->GetCSpaceGroup(); return $SPG->iSet(); }
sub FormulaUnit { return shift->{'FormulaUnit'}; }
sub SetFormulaUnit { my($this,$Z)=@_; return $this->{'FormulaUnit'} = $Z; }
sub ChemicalComposition { return shift->{'ChemicalComposition'}; }
sub SetChemicalComposition { my($this,$c)=@_; return $this->{'ChemicalComposition'} = $c; }
sub SumChemicalComposition { return shift->{'SumChemicalComposition'}; }
sub SetSumChemicalComposition { my($this,$c)=@_; return $this->{'SumChemicalComposition'} = $c; }
sub SetnAtomType { my ($this,$n) = @_; return $this->{nAtomType} = $n; }
sub nAtomType { my ($this) = @_; return $this->{nAtomType}; }
sub SetnAsymmetricAtomSite { my ($this,$n) = @_; return $this->{nAsymmetricAtomSite} = $n; }
sub nAsymmetricAtomSite { my ($this) = @_; return $this->{nAsymmetricAtomSite}; }
sub SetnTotalExpandedAtomSite { my ($this,$n) = @_; return $this->{nTotalExpandedAtomSite} = $n; }
sub nTotalExpandedAtomSite { my ($this) = @_; return $this->{nTotalExpandedAtomSite}; }
sub LatticeSystem {
my ($this) = @_;
return $this->{LatticeSystem} if($this->{LatticeSystem});
my $SPG = $this->GetCSpaceGroup();
return $SPG->LatticeSystem();
}
sub SetLatticeSystem
{
my ($this, $latticesystem, $SearchByLatticeParameter, $tollat, $tolangle) = @_;
my $SPG = $this->GetCSpaceGroup();
my $ls;
# my $ls = $SPG->LatticeSystem();
# if($ls eq '') {
$SPG->SetLatticeSystem($latticesystem, $SearchByLatticeParameter, $tollat, $tolangle);
$ls = $SPG->LatticeSystem();
# }
$this->{LatticeSystem} = $ls;
return $ls;
}
sub SetTemperature
{
my ($this,$T) = @_;
return $this->{Temperature} = $T;
}
sub Temperature
{
my ($this) = @_;
return $this->{Temperature};
}
#sub SpaceGroupObject::SetP1
#{
# my ($this) = @_;
# $this->ClearSymmetryOperation();
# $this->SetSPGName("P 1");
# $this->SetiSPG(1);
# $this->SetiSet(1);
# $this->AnalyzeTranslation();
# $this->AddSymmetryOperation('x,y,z');
# $this->SetLatticeSystem('triclinic');
#}
sub BurstToP1
{
my ($this) = @_;
my @ExpandedAtomSiteList = $this->GetCExpandedAtomSiteList();
my $SPG = $this->GetCSpaceGroup();
$this->SetSpaceGroup("P 1", 1, 1);
$this->SetLatticeSystem('triclinic', 0);
$SPG->SetP1();
$this->SetCAsymmetricAtomSiteList(@ExpandedAtomSiteList);
$this->ExpandCoordinates();
#my $ls = $this->LatticeSystem();
#print "ls=$ls\n";
#my $SPG = $this->GetCSpaceGroup();
#my $bl = $SPG->GetBravaisLattice();
#print "bl=$bl\n";
#exit;
return $this;
}
sub SetSpaceGroup
{
my ($this, $SPGName, $iSPG, $iSet) = @_;
my $SPG = new SpaceGroup();
$SPG->SetSPGName($SPGName);
$SPG->SetiSPG($iSPG);
$SPG->SetiSet($iSet) if(!defined $iSet);
$SPG->AddSymmetryOperation("x,y,z");
my $ret = $this->SetCSpaceGroup($SPG);
return $ret;
}
sub SetCSpaceGroup
{
my ($this, $SPG) = @_;
$SPG->SetLatticeParameters($this->LatticeParameters()) if(!$SPG->{a});
$this->{SpaceGroup} = $SPG;
$this->SetLatticeSystem() if(!$SPG->LatticeSystem());
return $SPG;
}
sub GetCSpaceGroup
{
my ($this) = @_;
if(not defined $this->{"SpaceGroup"} or $this->{"SpaceGroup"} eq '') {
my $SPG = new SpaceGroup();
$SPG->SetP1();
$this->{"SpaceGroup"} = $SPG;
}
return $this->{"SpaceGroup"};
}
sub GetBravaisLattice { my ($this) = @_; my $SPG = $this->GetCSpaceGroup();
return $SPG->GetBravaisLattice();
}
sub SetCAtomTypeList
{
my ($this, @atomlist) = @_;
$this->{nAtomType} = @atomlist;
return $this->{RefAtomTypeList} = \@atomlist;
}
sub ClearAtomTypeList
{
my ($this) = @_;
$this->{RefAtomTypeList} = [];
# delete $this->{RefAtomTypeList};
#print "p=$this->{RefAtomTypeList}\n";
}
sub GetpAtomType
{
my ($this) = @_;
my $RefAtomList = $this->{RefAtomTypeList};
}
sub GetCAtomTypeList
{
my ($this) = @_;
#$this->{RefAtomTypeList} = [];
my $RefAtomList = $this->{RefAtomTypeList};
#print "Ref:$RefAtomList
\n";
return () if($RefAtomList eq '');
return @$RefAtomList;
}
sub GetCAtomType
{
my ($this,$i) = @_;
my $RefAtomList = $this->{RefAtomTypeList};
my $atomtype = $RefAtomList->[$i];
if(!$atomtype) {
return 0;
}
return $atomtype;
}
sub AddAtomType
{
my ($this, $atomname, $CheckRegistered, $occ) = @_;
$CheckRegistered = 1 if(!defined $CheckRegistered);
$occ = 1.0 if(!defined $occ);
#print "AddAtomType: occ=$occ ($atomname: $CheckRegistered)\n";
my $atomlistref = $this->{RefAtomTypeList};
if($atomlistref) {
if($CheckRegistered and $this->FindAtomTypeByName($atomname)) {
return;
}
#print "atname[$atomname]: passed\n";
my $at = new AtomType();
$at->SetAtomName($atomname);
$at->SetOccupancy($occ);
push(@$atomlistref, $at);
$this->{"nAtomType"} = scalar @$atomlistref;
}
else {
my @atomlist;
my $at = new AtomType();
$at->SetAtomName($atomname);
$at->SetOccupancy($occ);
push(@atomlist, $at);
$this->{"nAtomType"} = @atomlist;
$this->{"RefAtomTypeList"} = \@atomlist;
}
#my $ref = $this->{"RefAtomTypeList"} ;
#print "ref: ", $ref->[$this->{"nAtomType"}-1]->AtomName(), "
";
return $this->{"nAtomType"};
}
sub SetCAsymmetricAtomSiteList
{
my ($this, @atomlist) = @_;
$this->{nAsymmetricAtomSite} = @atomlist;
return $this->{RefAsymmetricAtomSiteList} = \@atomlist;
}
sub GetpAtomSite
{
my ($this, $mode) = @_;
return $this->GetpAtomSiteList($mode);
}
sub GetpAtomSiteList
{
my ($this, $mode) = @_;
if($mode == 0) {
return $this->{"RefAsymmetricAtomSiteList"};
}
return $this->{RefExpandedAtomSiteList};
}
sub GetpAsymmetricAtomSite
{
my ($this) = @_;
return $this->{"RefAsymmetricAtomSiteList"};
}
sub GetAtomSiteList
{
my ($this, $mode) = @_;
if(!defined $mode or $mode == 0) {
return $this->GetCAsymmetricAtomSiteList();
}
else {
return $this->GetCExpandedAtomSiteList();
}
}
sub GetCAsymmetricAtomSiteList
{
my ($this) = @_;
my $RefAtomList = $this->{"RefAsymmetricAtomSiteList"};
return () if($RefAtomList eq '');
return @$RefAtomList;
}
sub GetCAtomSite
{
my ($this,$i) = @_;
return $this->GetAtomSite($i);
}
sub GetAtomSite
{
my ($this,$i) = @_;
my $RefAtomList = $this->{RefAsymmetricAtomSiteList};
my $atomsite = $RefAtomList->[$i];
if(!$atomsite) {
return 0;
}
return $atomsite;
}
sub GetpExpandedAtomSite
{
my ($this) = @_;
return $this->{RefExpandedAtomSiteList};
}
sub GetCExpandedAtomSiteList
{
my ($this) = @_;
my $RefAtomList = $this->{RefExpandedAtomSiteList};
return () if($RefAtomList eq '');
return @$RefAtomList;
}
sub GetCExpandedAtomSite
{
my ($this,$i) = @_;
my $RefAtomList = $this->{RefExpandedAtomSiteList};
my $atomsite = $RefAtomList->[$i];
if(!$atomsite) {
return 0;
}
return $atomsite;
}
sub FindAtomSite
{
my ($this, $mode, $atomname) = @_;
my $p = $this->GetpAtomSiteList($mode);
for(my $i = 0 ; $i < @$p ; $i++) {
my $atom = $p->[$i];
my $name = $atom->AtomNameOnly();
if(uc $atomname eq uc $name) {
return $atom;
}
}
}
sub FindiAtom
{
my ($this, $mode, $atomname) = @_;
my $p = $this->GetpAtomSiteList($mode);
for(my $i = 0 ; $i < @$p ; $i++) {
my $atom = $p->[$i];
my $name = $atom->AtomNameOnly();
if(uc $atomname eq uc $name) {
return $i;
}
}
}
sub FindNearestSite
{
my ($this, $x0, $y0, $z0) = @_;
my $p = $this->GetpAtomSiteList(1);
my $minidx = -1;
my $mindis = 1.0e10;
for(my $i = 0 ; $i < @$p ; $i++) {
my $atom = $p->[$i];
my ($x1, $y1, $z1) = $atom->Position(1);
my $dis = $this->GetNearestInterAtomicDistance(
$x0, $y0, $z0, $x1, $y1, $z1);
if($dis < $mindis) {
$mindis = $dis;
$minidx = $i;
}
}
return $minidx;
}
sub GetpAtomPosition
{
my ($this, $mode) = @_;
my $pSite = $this->GetpAtomSite($mode);
my @Pos;
for(my $i = 0 ; $i < @$pSite ; $i++) {
my ($x, $y, $z) = $pSite->[$i]->Position(0);
$Pos[$i][0] = $x;
$Pos[$i][1] = $y;
$Pos[$i][2] = $z;
}
return \@Pos;
}
sub GetpDuplicatedAtomPosition
{
my ($this, $mode) = @_;
my @Pos = $this->DuplicateAtomPositions($mode);
return \@Pos;
}
sub DuplicateAtomPositions
{
my ($this, $mode) = @_;
my $pSite = $this->GetpAtomSite($mode);
my @Pos;
for(my $i = 0 ; $i < @$pSite ; $i++) {
my ($x, $y, $z) = $pSite->[$i]->Position(0);
$Pos[$i][0] = $x;
$Pos[$i][1] = $y;
$Pos[$i][2] = $z;
}
return @Pos;
}
sub Volume
{
my ($this,$UseAtomicUnit) = @_;
my $k = 1.0 / 0.529177;
$k = 1.0 unless($UseAtomicUnit);
$k = $k*$k*$k;
return $this->{"Volume"} * $k;
}
sub ReciprocalVolume
{
my ($this,$UseAtomicUnit) = @_;
my $k = 1.0 / 0.529177;
$k = 1.0 unless($UseAtomicUnit);
$k = $k*$k*$k;
$this->{RVolume} = $this->CalculateReciprocalVolume() if($this->{RVolume} == 0.0);
return $this->{RVolume} / $k;
}
sub Duplicate
{
my ($this) = @_;
my $Crystal = new Crystal;
$Crystal->SetSampleName($this->SampleName());
$Crystal->SetCrystalName($this->CrystalName());
$Crystal->SetLatticeSystem($this->LatticeSystem());
my $SPGName = $this->SPGName();
my $iSPG = $this->iSPG();
my $iSet = $this->iSet();
$Crystal->SetSpaceGroup($SPGName, $iSPG, $iSet);
$Crystal->SetCSpaceGroup($this->GetCSpaceGroup());
$Crystal->SetLatticeParameters($this->LatticeParameters());
$Crystal->SetCAtomTypeList($this->GetCAtomTypeList());
$Crystal->SetCAsymmetricAtomSiteList($this->GetCAsymmetricAtomSiteList());
$Crystal->SetTemperature($this->Temperature());
$Crystal->CalcMetrics();
my @v = $this->LatticeVectors();
if(@v > 0) {
$Crystal->SetLatticeVectors(@v);
}
$Crystal->ExpandCoordinates();
return $Crystal;
}
sub Symmetrize
{
my ($this, $tollatt, $tolangle, $tolpos) = @_;
$tollatt = 1.0e-5 if(!defined $tollatt);
$tolangle = $tollatt if(!defined $tolangle);
$tolpos = $tollatt if(!defined $tolpos);
my ($a,$b,$c,$alpha,$beta,$gamma) = $this->LatticeParameters();
$alpha = 90 if(abs($alpha-90.0) < $tolangle);
$beta = 90 if(abs($beta -90.0) < $tolangle);
$gamma = 90 if(abs($gamma-90.0) < $tolangle);
$this->LatticeParameters($a, $b, $c, $alpha, $beta, $gamma);
$this->SetLatticeSystem('', 1, $tollatt, $tolangle);
my $SPG = $this->GetCSpaceGroup();
my $ls = $SPG->LatticeSystem();
if($ls eq 'cubic') {
my $ar = ($a+$b+$c) / 3.0;
$this->SetLatticeParameters($ar,$ar,$ar,90,90,90);
}
elsif($ls eq 'rhombohedral' or $ls eq 'trigonal') {
my $ar = ($a+$b+$c) / 3.0;
my $al = ($alpha+$beta+$gamma) / 3.0;
$this->SetLatticeParameters($ar,$ar,$ar,$alpha,$alpha,$alpha);
}
elsif($ls eq 'tetragonal') {
my $ar = ($a+$b) / 2.0;
$this->SetLatticeParameters($ar,$ar,$c,90,90,90);
}
elsif($ls eq 'hexagonal') {
my $ar = ($a+$b) / 2.0;
$this->SetLatticeParameters($ar,$ar,$c,90,90,120);
}
elsif($ls eq 'orthorhombic') {
$this->SetLatticeParameters($a,$b,$c,90,90,90);
}
my @site = $this->GetCAsymmetricAtomSiteList();
for(my $i = 0 ; $i < @site ; $i++) {
my ($x, $y, $z) = $site[$i]->Position();
$x = $this->RoundParameter($x, $tolpos);
$y = $this->RoundParameter($y, $tolpos);
$z = $this->RoundParameter($z, $tolpos);
$site[$i]->SetPosition($x, $y, $z);
}
$this->ExpandCoordinates();
}
sub SymmetrizeParameter
{
my ($this, $ls, $x, $y, $z, $tol) = @_;
$ls = $this->LatticeSystem() if(!defined $ls or $ls eq '');
if($ls eq 'cubic' or $ls eq 'rhombohedral' or $ls eq 'trigonal') {
$x = $y = $z = ($x+$y+$z) / 3.0;
}
elsif($ls eq 'tetragonal' or $ls eq 'hexagonal') {
$x = $y = ($x+$y) / 2.0;
}
$x = $this->RoundParameter($x, $tol);
$y = $this->RoundParameter($y, $tol);
$z = $this->RoundParameter($z, $tol);
return ($x, $y, $z);
}
sub RoundParameter
{
my ($this, $x, $tol) = @_;
return $tol * int( ($x+0.1*$tol) / $tol );
}
sub GetCnMultiplicityExpandedAtomSiteList
{
my ($this) = @_;
my $RefnAtomList = $this->{"RefnExpandedAtomSite"};
return () if($RefnAtomList eq '');
return @$RefnAtomList;
}
#How to build Crystal (Ex. from CIF.pm)
# my $crystal = new Crystal();
# $crystal->SetLatticeParameters(
# $this->LatticeParameters() );
# my $SPG = $this->GetCSpaceGroup();
# $crystal->SetCSpaceGroup($SPG);
# $crystal->SetCAtomTypeList(
# $this->GetCAtomTypeList() );
# $crystal->SetCAsymmetricAtomSiteList(
# $this->GetCAsymmetricAtomSiteList() );
#==================================================================
# コンストラクタ、デストラクト
#==================================================================
sub new
{
my ($module) = @_;
my $this = {};
bless $this;
require Math::Matrix;
require Math::MatrixReal;
my $SPG = new SpaceGroup();
$SPG->SetP1();
$this->{"SpaceGroup"} = $SPG;
my @AtomTypeList;
# my $at1 = new AtomType();
# push(@AtomTypeList, $at1);
# pop(@AtomTypeList);
$this->{"RefAtomTypeList"} = \@AtomTypeList;;
$this->{"nAtomType"} = 0;
my @AsAtomSiteList;
my $at2 = new AtomSite();
push(@AsAtomSiteList, $at2);
pop(@AsAtomSiteList);
$this->{"RefAsymmetricAtomSiteList"} = \@AsAtomSiteList;
$this->{"nAsymmetricAtomSite"} = 0;
my @ExAtomSiteList;
my $at3 = new AtomSite();
push(@ExAtomSiteList, $at3);
pop(@ExAtomSiteList);
$this->{"RefExpandedAtomSiteList"} = \@ExAtomSiteList;
my @nAtomList;
$this->{"RefnExpandedAtomSite"} = \@nAtomList;
$this->{"nTotalExpandedAtomSite"} = 0;
return $this;
}
sub DESTROY
{
my $this = shift;
}
#==================================================================
# 一般メンバー関数
#==================================================================
sub PreferentialOrientation
{
my ($this, $IsSheet, $ho, $ko, $lo, $p1, $p2, $h, $k, $l) = @_;
return 1.0 if($p1 == 1.0 or $p2 == 0.0);
my @VO = $this->GetVectorFromHKL($ho, $ko, $lo);
my @Vhkl = $this->GetVectorFromHKL($h, $k, $l);
my $eVO = Sci::NormalizeVector(\@VO);
my $eVhkl = Sci::NormalizeVector(\@Vhkl);
my $c = Sci::InnerProduct($eVO, $eVhkl);
my $Q = acos(abs($c));
$Q = $pi / 2.0 - $Q if(!$IsSheet);
#print "PO [$IsSheet]: $ho,$ko,$lo - $h,$k,$l => $Q\n";
return $p1 + (1.0 - $p1) * exp(-$p2 * $Q * $Q);
}
sub IsVariablePosition
{
my ($this, $x, $y, $z, $pos)= @_;
my @x1 = ($x, $y, $z);
my $c0 = new Crystal();
$c0->SetLatticeParameters($this->LatticeParameters());
$c0->SetCSpaceGroup($this->GetCSpaceGroup());
$c0->AddAtomType('H');
$c0->AddAtomSite('H1', 'H', $x1[0], $x1[1], $x1[2]);
$c0->ExpandCoordinates();
my @sites = $c0->GetCExpandedAtomSiteList();
my $nSite = @sites;
my $VariablePos = 1.0;
$x1[$pos] += 0.14142;
my $c1 = new Crystal();
$c1->SetLatticeParameters($this->LatticeParameters());
$c1->SetCSpaceGroup($this->GetCSpaceGroup());
$c1->AddAtomType('H');
$c1->AddAtomSite('H1', 'H', $x1[0], $x1[1], $x1[2]);
$c1->ExpandCoordinates();
@sites = $c1->GetCExpandedAtomSiteList();
$VariablePos = 0.0 if($nSite < @sites);
return $VariablePos;
}
sub SetMinimalLatticeParameters
{
my ($this, @a) = @_;
my $ls = lc $this->LatticeSystem();
#print "ls: $ls\n";
if($ls eq 'cubic') {
$this->SetLatticeParameters($a[0], $a[0], $a[0], 90.0, 90.0, 90.0);
}
elsif($ls =~ /tetra/) {
$this->SetLatticeParameters($a[0], $a[0], $a[1], 90.0, 90.0, 90.0);
}
elsif($ls =~ /ortho/) {
$this->SetLatticeParameters($a[0], $a[1], $a[2], 90.0, 90.0, 90.0);
}
elsif($ls =~ /(rhomb|trig)/) {
$this->SetLatticeParameters($a[0], $a[0], $a[0], $a[1], $a[1], $a[1]);
}
elsif($ls =~ /hex/) {
$this->SetLatticeParameters($a[0], $a[0], $a[1], 90.0, 90.0, 120.0);
}
elsif($ls =~ /mono/) {
$this->SetLatticeParameters($a[0], $a[1], $a[2], 90.0, $a[3], 90.0);
}
else {
$this->SetLatticeParameters($a[0], $a[1], $a[2], $a[3], $a[4], $a[5]);
}
my ($a,$b,$c,$alpha,$beta,$gamma) = $this->LatticeParameters();
#print "latt: $a $b $c $alpha $beta $gamma\n";
return ($a,$b,$c,$alpha,$beta,$gamma);
}
sub SetLatticeParameters
{
my ($this, $a, $b, $c, $alpha, $beta, $gamma) = @_;
# return if(!defined $a or $a == 0);
$this->{a} = $a;
$this->{b} = $b;
$this->{c} = $c;
$this->{alpha} = $alpha;
$this->{beta} = $beta;
$this->{gamma} = $gamma;
#print "CrystalObject::SetLatticeParameters\n";
$this->CalcMetrics();
$this->GetCSpaceGroup()->SetLatticeParameters($a, $b, $c, $alpha, $beta, $gamma);
return ($a, $b, $c, $alpha, $beta, $gamma);
}
sub LatticeParameters
{
my ($this, $UseAtomicUnit) = @_;
my $a = $this->{'a'};
#$a = 1.0 if(not defined $a);
my $b = $this->{'b'};
#$b = $a if(not defined $b);
my $c = $this->{'c'};
#$c = $b if(not defined $c);
my $alpha = $this->{'alpha'};
#$alpha = 90.0 if(not defined $alpha);
my $beta = $this->{'beta'};
#$beta = 90.0 if(not defined $beta);
my $gamma = $this->{'gamma'};
#$gamma = 90.0 if(not defined $gamma);
my $k = ($UseAtomicUnit)? 1.0 / 0.529177 : 1.0;
return ($k*$a, $k*$b, $k*$c, $alpha, $beta, $gamma);
}
sub ReciprocalLatticeParameters
{
my ($this, $UseAtomicUnit) = @_;
my $Ra = $this->{Ra};
my $Rb = $this->{Rb};
my $Rc = $this->{Rc};
my $Ralpha = $this->{Ralpha};
my $Rbeta = $this->{Rbeta};
my $Rgamma = $this->{Rgamma};
my $k = ($UseAtomicUnit)? 0.529177 : 1.0;
return ($k*$Ra, $k*$Rb, $k*$Rc, $Ralpha, $Rbeta, $Rgamma);
}
sub GetReciprocalDistanceFromK
{
my ($this, $kx0, $ky0, $kz0, $kx1, $ky1, $kz1) = @_;
#printf "k: (%6.3g,%6.3g,%6.3g) - (%6.3g,%6.3g,%6.3g): ",
# $kx0,$ky0,$kz0,$kx1,$ky1,$kz1;
$kx0 = $kx1 - $kx0 if(defined $kx1);
$ky0 = $ky1 - $ky0 if(defined $ky1);
$kz0 = $kz1 - $kz0 if(defined $kz1);
#printf "(%6.3g,%6.3g,%6.3g)\n", $kx0, $ky0, $kz0;
if(0) {
my $Ra = $this->{'Ra'};
my $Rb = $this->{'Rb'};
my $Rc = $this->{'Rc'};
#print "Rec lat: $Ra $Rb $Rc\n";
$kx0 = $kx0 / $Ra;
$ky0 = $ky0 / $Rb;
$kz0 = $kz0 / $Rc;
}
my $Rg11 = $this->{"Rg11"};
my $Rg12 = $this->{"Rg12"};
my $Rg13 = $this->{"Rg13"};
my $Rg21 = $this->{"Rg21"};
my $Rg22 = $this->{"Rg22"};
my $Rg23 = $this->{"Rg23"};
my $Rg31 = $this->{"Rg31"};
my $Rg32 = $this->{"Rg32"};
my $Rg33 = $this->{"Rg33"};
#print "Rg: ($Rg11,$Rg22,$Rg33,$Rg12,$Rg23,$Rg13)
\n";
#print "($h,$k,$l)-($Rg11,$Rg22,$Rg33,$Rg12,$Rg23,$Rg13)
\n";
my $d = $Rg11*$kx0*$kx0 + $Rg22*$ky0*$ky0 + $Rg33*$kz0*$kz0
+ 2.0*$Rg12*$kx0*$ky0 + 2.0*$Rg23*$ky0*$kz0 + 2.0*$Rg13*$kx0*$kz0;
$d = 0.0 if($d < 0.0);
$d = sqrt($d);
#print "d=$d
\n";
#print "enter: $d ($h,$k,$l),($Rg11,$Rg22,$Rg33,$Rg12,$Rg23,$Rg13)
";
return $d;
}
sub CalculateHKLInterplanarSpacing
{
my ($this, $h, $k, $l) = @_;
my $Rg11 = $this->{"Rg11"};
my $Rg12 = $this->{"Rg12"};
my $Rg13 = $this->{"Rg13"};
my $Rg21 = $this->{"Rg21"};
my $Rg22 = $this->{"Rg22"};
my $Rg23 = $this->{"Rg23"};
my $Rg31 = $this->{"Rg31"};
my $Rg32 = $this->{"Rg32"};
my $Rg33 = $this->{"Rg33"};
#print "($h,$k,$l)-($Rg11,$Rg22,$Rg33,$Rg12,$Rg23,$Rg13)
\n";
my $d = $Rg11*$h*$h + $Rg22*$k*$k + $Rg33*$l*$l
+ 2.0*$Rg12*$h*$k + 2.0*$Rg23*$k*$l + 2.0*$Rg13*$l*$h;
#print "d=$d
";
#print "enter: $d ($h,$k,$l),($Rg11,$Rg22,$Rg33,$Rg12,$Rg23,$Rg13)
";
return 0.0 if($d == 0.0);
$d = sqrt(1.0 / $d);
#print "d=$d
";
return $d;
}
sub CalculateDiffractionAngleFromInterplanarSpacing
{
my ($this, $wl, $d) = @_;
return 0.0 if($d == 0.0);
my $sq = $wl * 0.5 / $d;
#print "sq: $sq
";
return 0.0 if(abs($sq) > 1.0);
my $Q = Sci::dasin($sq);
return $Q;
}
sub CalMetrics
{
my ($this) = @_;
return $this->CalcMetrics();
}
sub CalcMetrics
{
my ($this) = @_;
my $SPG = $this->GetCSpaceGroup();
my $SPGName = $SPG->SPGName();
my $ls = $this->LatticeSystem();
my $pT = $this->GetMatrixConventionalToPrimitiveCell(0);
my ($a,$b,$c,$alpha,$beta,$gamma) = $this->LatticeParameters();
if($a == 0.0) {
print "Error: CrystalObject.pm::CalcMetrics: Invalid Lattice Parameters.
\n";
print "Lattice: $a $b $c $alpha $beta $gamma
\n";
}
my $torad = Sci::torad(); #3.1415926535 / 180.0;
$this->{g11} = $a * $a;
$this->{g12} = $this->{g21} = $a * $b * cos($torad * $gamma);
$this->{g12} = $this->{g21} = 0.0 if(abs($this->{"g12"}) < 1.0e-7);
$this->{g13} = $this->{g31} = $a * $c * cos($torad * $beta);
$this->{g13} = $this->{g31} = 0.0 if(abs($this->{"g13"}) < 1.0e-7);
$this->{g22} = $b * $b;
$this->{g23} = $this->{g32} = $b * $c * cos($torad * $alpha);
$this->{g23} = $this->{g32} = 0.0 if(abs($this->{"g23"}) < 1.0e-7);
$this->{g33} = $c * $c;
my $V = $this->CalculateVolume();
#print "V=$V
";
#print "a=$a\n";
my ($ax, $ay, $az, $bx, $by, $bz, $cx, $cy, $cz);
#print "SPG:[$this->{ConventionalSPGName}][$SPGName] LS: $ls\n";
#exit;
if($ls eq 'trigonal' and abs($alpha - 60.0) < 1.0e-3) {
$this->{ConventionalSPGName} = 'F';
my $ac = $a * sqrt(2.0);
$this->{ConventionalLatticeParameters} = [$ac, $ac, $ac, 90.0, 90.0, 90.0];
}
if((defined $this->{ConventionalSPGName} and $this->{ConventionalSPGName} =~ /^\s*[FI]/i)) {
#print "CrystalObject::CalcMetrics: Lattice Vectors for primitive cell (SPG: $this->{ConventionalSPGName} / $SPGName)\n";
#print "********[$this->{ConventionalSPGName}][$SPGName]\n";
#exit;
my $pT = $this->GetMatrixConventionalToPrimitiveCell(1);
my $pl = $this->{ConventionalLatticeParameters};
my ($a, $b, $c, $alpha, $beta, $gamma) = @$pl;
#print "a=$a, $b, $c, $alpha\n";
#exit;
#Replace with those for [ABC] after check for [FI]
# $ax = $c * $pT->[0][0];
# $ay = $c * $pT->[0][1];
# $az = $c * $pT->[0][2];
# $bx = $a * $pT->[1][0];
# $by = $a * $pT->[1][1];
# $bz = $a * $pT->[1][2];
# $cx = $b * $pT->[2][0];
# $cy = $b * $pT->[2][1];
# $cz = $b * $pT->[2][2];
$ax = $a * $pT->[0][0];
$ay = $b * $pT->[0][1];
$az = $c * $pT->[0][2];
$bx = $a * $pT->[1][0];
$by = $b * $pT->[1][1];
$bz = $c * $pT->[1][2];
$cx = $a * $pT->[2][0];
$cy = $b * $pT->[2][1];
$cz = $c * $pT->[2][2];
#print "X: $pT->[0][0] $pT->[0][1] $pT->[0][2]\n";
#print "Y: $pT->[1][0] $pT->[1][1] $pT->[1][2]\n";
#print "Z: $pT->[2][0] $pT->[2][1] $pT->[2][2]\n";
#exit;
}
elsif(defined $this->{ConventionalSPGName} and $this->{ConventionalSPGName} =~ /^\s*[ABC]/i) {
#print "CrystalObject::CalcMetrics: Lattice Vectors for primitive cell (SPG: $this->{ConventionalSPGName} / $SPGName)\n";
#print "********[$this->{ConventionalSPGName}][$SPGName]\n";
my $pT = $this->GetMatrixConventionalToPrimitiveCell(1);
my $pl = $this->{ConventionalLatticeParameters};
my ($a, $b, $c, $alpha, $beta, $gamma) = @$pl;
$ax = $a * $pT->[0][0];
$ay = $b * $pT->[0][1];
$az = $c * $pT->[0][2];
$bx = $a * $pT->[1][0];
$by = $b * $pT->[1][1];
$bz = $c * $pT->[1][2];
$cx = $a * $pT->[2][0];
$cy = $b * $pT->[2][1];
$cz = $c * $pT->[2][2];
#print "latt: $a, $b, $c, $alpha, $beta, $gamma\n";
#print "pT ($pT->[0][0], $pT->[0][1], $pT->[0][2])\n";
#print "pT ($pT->[1][0], $pT->[1][1], $pT->[1][2])\n";
#print "pT ($pT->[2][0], $pT->[2][1], $pT->[2][2])\n";
#print "a=($ax, $ay, $az)\n";
#print "b=($bx, $by, $bz)\n";
#print "c=($cx, $cy, $cz)\n";
#print "X: $pT->[0][0] $pT->[0][1] $pT->[0][2]\n";
#print "Y: $pT->[1][0] $pT->[1][1] $pT->[1][2]\n";
#print "Z: $pT->[2][0] $pT->[2][1] $pT->[2][2]\n";
#exit;
}
elsif($ls =~ /hex/i) {
# elsif($SPGName =~ /^\s*P\s*6/) {
#print "CrystalObject::CalcMetrics: Lattice Vectors for hexagonal (SPG: $SPGName)\n";
# $ax = $a;
# $ay = 0.0;
# $az = 0.0;
# $bx = $ax / 2.0;
# $by = sqrt(3) / 2.0 * $ay;
# $bz = 0.0;
$ax = $a / 2.0;
$ay = -sqrt(3.0) / 2.0 * $a;
$az = 0.0;
$bx = $ax;
$by = -$ay;
$bz = 0.0;
$cx = 0.0;
$cy = 0.0;
$cz = $c;
}
elsif($ls =~ /^(trig|rhomb)/i) {
# elsif($SPGName =~ /^\s*R/ and abs($alpha - 90.0) > 1.0e-2 and abs($alpha - $beta) < 1.0e-2) {
#print "CrystalObject::CalcMetrics: Lattice Vectors for rhombohedral (SPG: $SPGName)\n";
my $cosa = cos($torad * $alpha);
my $axy = $a * sqrt( (1.0 - $cosa) / 6.0);
$ax = sqrt(3.0) * $axy;
$ay = $axy;
$az = $axy * sqrt(6.0 / (1.0 - $cosa) - 4.0);
$bx = -$ax;
$by = $ay;
$bz = $az;
$cx = 0.0;
$cy = -2.0 * $axy;
$cz = $az;
#print "alpha=$alpha, cosa=$cosa, axy=$axy, a=$a\n";
#print "a=sqrt(a2+b2+c2)=", sqrt($ax*$ax+$ay*$ay+$az*$az), "\n";
}
else {
#print "CrystalObject::CalcMetrics: Lattice Vectors for conventional cell (SPG: $SPGName)\n";
$ax = $a;
$ay = 0.0;
$az = 0.0;
$bx = $b * Sci::dcos($gamma);
$by = $b * Sci::dsin($gamma);
$bz = 0.0;
$cx = $c * Sci::dcos($beta);
$cy = ($c * Sci::dcos($alpha) - $cx * Sci::dcos($gamma)) / Sci::dsin($gamma);
$cz = sqrt($c*$c - $cx*$cx - $cy*$cy);
}
$this->{a11} = $ax;
$this->{a12} = $ay;
$this->{a13} = $az;
$this->{a21} = $bx;
$this->{a22} = $by;
$this->{a23} = $bz;
$this->{a31} = $cx;
$this->{a32} = $cy;
$this->{a33} = $cz;
my $sina = Sci::dsin($alpha);
my $sinb = Sci::dsin($beta);
my $sinc = Sci::dsin($gamma);
my $cosa = Sci::dcos($alpha);
my $cosb = Sci::dcos($beta);
my $cosc = Sci::dcos($gamma);
$this->{Rg11} = $b*$b*$c*$c*$sina*$sina / $V / $V;
$this->{Rg12} = $this->{Rg21} = $a*$b*$c*$c * ($cosa*$cosb - $cosc) / $V / $V;
$this->{Rg12} = $this->{Rg21} = 0.0 if(abs($this->{"g12"}) < 1.0e-7);
$this->{Rg13} = $this->{Rg31} = $a*$b*$b*$c * ($cosc*$cosa - $cosb) / $V / $V;
$this->{Rg13} = $this->{Rg31} = 0.0 if(abs($this->{"g13"}) < 1.0e-7);
$this->{Rg22} = $a*$a*$c*$c*$sinb*$sinb / $V / $V;
$this->{Rg23} = $this->{Rg32} = $a*$a*$b*$c * ($cosb*$cosc - $cosa) / $V / $V;
$this->{Rg23} = $this->{Rg32} = 0.0 if(abs($this->{"g23"}) < 1.0e-7);
$this->{Rg33} = $a*$a*$b*$b*$sinc*$sinc / $V / $V;
my $Ra = $this->{Ra} = sqrt($this->{Rg11});
my $Rb = $this->{Rb} = sqrt($this->{Rg22});
my $Rc = $this->{Rc} = sqrt($this->{Rg33});
my $Ralpha = $this->{Ralpha} = $todeg * acos($this->{Rg23} / $this->{Rb} / $this->{Rc});
my $Rbeta = $this->{Rbeta} = $todeg * acos($this->{Rg31} / $this->{Rc} / $this->{Ra});
my $Rgamma = $this->{Rgamma} = $todeg * acos($this->{Rg12} / $this->{Ra} / $this->{Rb});
#print "RLatt: $Ra, $Rb, $Rc, $Ralpha, $Rbeta, $Rgamma\n";
my $kV = 1.0 / $V;
my ($Rax, $Ray, $Raz) = Sci::VectorProduct([$bx, $by, $bz], [$cx, $cy, $cz], $kV);
my ($Rbx, $Rby, $Rbz) = Sci::VectorProduct([$cx, $cy, $cz], [$ax, $ay, $az], $kV);
my ($Rcx, $Rcy, $Rcz) = Sci::VectorProduct([$ax, $ay, $az], [$bx, $by, $bz], $kV);
#my $cc = Sci::ScalarProduct([$Rax, $Ray, $Raz], [$Rax, $Ray, $Raz]);
#print "cc=$cc, $this->{Rg11}\n";
#$cc = Sci::ScalarProduct([$Rax, $Ray, $Raz], [$Rbx, $Rby, $Rbz]);
#print "cc=$cc, $this->{Rg12}\n";
#printf "Rx*1=(%7.3f, %7.3f, %7.3f)\n", $Rax, $Ray, $Raz;
#printf "Ry*1=(%7.3f, %7.3f, %7.3f)\n", $Rbx, $Rby, $Rbz;
#printf "Rz*1=(%7.3f, %7.3f, %7.3f)\n", $Rcx, $Rcy, $Rcz;
if(0) {
my $Rax = $Ra;
my $Ray = 0.0;
my $Raz = 0.0;
my $Rbx = $Rb * Sci::dcos($Rgamma);
my $Rby = $Rb * Sci::dsin($Rgamma);
my $Rbz = 0.0;
my $Rcx = $Rc * Sci::dcos($Rbeta);
my $Rcy = ($Rc * Sci::dcos($Ralpha) - $Rcx * Sci::dcos($Rgamma)) / Sci::dsin($Rgamma);
my $Rcz = sqrt($Rc*$Rc - $Rcx*$Rcx - $Rcy*$Rcy);
}
$this->{Ra11} = $Rax;
$this->{Ra12} = $Ray;
$this->{Ra13} = $Raz;
$this->{Ra21} = $Rbx;
$this->{Ra22} = $Rby;
$this->{Ra23} = $Rbz;
$this->{Ra31} = $Rcx;
$this->{Ra32} = $Rcy;
$this->{Ra33} = $Rcz;
#printf "Rx*2=(%7.3f, %7.3f, %7.3f)\n", $Rax, $Ray, $Raz;
#printf "Ry*2=(%7.3f, %7.3f, %7.3f)\n", $Rbx, $Rby, $Rbz;
#printf "Rz*2=(%7.3f, %7.3f, %7.3f)\n", $Rcx, $Rcy, $Rcz;
return 1;
}
sub LatticeVectorsByMatrix
{
my ($this, $UseAtomicUnit) = @_;
my @v = $this->LatticeVectors($UseAtomicUnit);
my @a1 = ($v[0], $v[1], $v[2]);
my @a2 = ($v[3], $v[4], $v[5]);
my @a3 = ($v[6], $v[7], $v[8]);
return (\@a1, \@a2, \@a3);
}
# @v in Math::Vector::Real
sub LatticeVectorsByVectors
{
my ($this, $UseAtomicUnit) = @_;
my @v = $this->LatticeVectors($UseAtomicUnit);
my $a1 = V($v[0], $v[1], $v[2]);
my $a2 = V($v[3], $v[4], $v[5]);
my $a3 = V($v[6], $v[7], $v[8]);
return ($a1, $a2, $a3);
}
sub LatticeVectors
{
my ($this, $UseAtomicUnit) = @_;
#my $SPG = $this->GetCSpaceGroup();
#my $SPGName = $SPG->SPGName();
#print "SPG in CrystalObject.pm:[$SPGName]\n";
#exit;
my $k = ($UseAtomicUnit)? 1.0 / 0.529177 : 1.0;
return ($k*$this->{a11}, $k*$this->{a12}, $k*$this->{a13},
$k*$this->{a21}, $k*$this->{a22}, $k*$this->{a23},
$k*$this->{a31}, $k*$this->{a32}, $k*$this->{a33});
}
sub ReciprocalLatticeVectorsByMatrix
{
my ($this, $UseAtomicUnit) = @_;
my @v = $this->ReciprocalLatticeVectors($UseAtomicUnit);
my @a1 = $v[0,2];
my @a2 = $v[3,5];
my @a3 = $v[6,8];
return (\@a1, \@a2, \@a3);
}
sub ReciprocalLatticeVectorsByVectors
{
my ($this, $UseAtomicUnit) = @_;
my @v = $this->ReciprocalLatticeVectors($UseAtomicUnit);
my $a1 = V($v[0,2]);
my $a2 = V($v[3,5]);
my $a3 = V($v[6,8]);
return ($a1, $a2, $a3);
}
sub ReciprocalLatticeVectors
{
my ($this, $UseAtomicUnit) = @_;
my $k = ($UseAtomicUnit)? 1.0 / 0.529177 : 1.0;
return ($k*$this->{Ra11}, $k*$this->{Ra12}, $k*$this->{Ra13},
$k*$this->{Ra21}, $k*$this->{Ra22}, $k*$this->{Ra23},
$k*$this->{Ra31}, $k*$this->{Ra32}, $k*$this->{Ra33});
}
sub GetLatticeRange
{
my ($this, $MaximumDistance) = @_;
my ($a,$b,$c,$alpha,$beta,$gamma) = $this->LatticeParameters();
my $sinA = abs(sin($alpha * $torad));
my $sinB = abs(sin($beta * $torad));
my $sinG = abs(sin($gamma * $torad));
$sinA = $sinB if($sinB < $sinA);
$sinA = $sinG if($sinG < $sinA);
my $nx = int($MaximumDistance / $a / $sinA);
my $ny = int($MaximumDistance / $b / $sinA);
my $nz = int($MaximumDistance / $c / $sinA);
return ($nx+1, $ny+1, $nz+1);
}
sub GetYRangeForIteration
{
my ($this, $ix) = @_;
my $ls = lc $this->LatticeSystem();
#print "ls: $ls\n";
if($ls eq 'cubic' or $ls =~ /(rhomb|trig)/) {
return $ix;
}
elsif($ls =~ /(rhomb|trig)/) {
return $ix;
}
elsif($ls =~ /tetra/ or $ls =~ /hex/) {
return $ix;
}
elsif($ls =~ /ortho/) {
return 0;
}
elsif($ls =~ /mono/) {
return 0;
}
else {
return 0;
}
}
sub GetZRangeForIteration
{
my ($this, $ix, $iy) = @_;
my $ls = lc $this->LatticeSystem();
#print "ls: $ls\n";
if($ls eq 'cubic' or $ls =~ /(rhomb|trig)/) {
return $iy;
}
elsif($ls =~ /tetra/ or $ls =~ /hex/) {
return 0;
}
elsif($ls =~ /ortho/) {
return 0;
}
elsif($ls =~ /mono/) {
return 0;
}
else {
return 0;
}
}
sub ConvertLatticeIndexByLatticeSystem
{
my ($this, $ix, $iy, $iz, $SortDirection) = @_;
$SortDirection = 1 if(!defined $SortDirection);
if($SortDirection == 0) {
return (0, $ix, $iy, $iz);
}
my $ls = lc $this->LatticeSystem();
#print "ls: $ls\n";
my $IsConverted = 0;
if($ls eq 'cubic' or $ls =~ /(rhomb|trig)/) {
$IsConverted = 1 if($ix > $iy or $iy > $iz or $ix > $iz);
if($SortDirection > 0) {
return ($IsConverted, sort { return $a <=> $b; } ($ix, $iy, $iz) );
}
return ($IsConverted, sort { return $b <=> $a; } ($ix, $iy, $iz) );
}
elsif($ls =~ /tetra/ or $ls =~ /hex/) {
$IsConverted = 1 if($ix > $iy);
my @a;
if($SortDirection > 0) {
@a = sort { return $a <=> $b; } ($ix, $iy);
}
else {
@a = sort { return $b <=> $a; } ($ix, $iy);
}
return ($IsConverted, @a, $iz);
}
else {
# die "\nCrystaObject::ConvertLatticeIndexByLatticeSystem: Error: Lattice system [$ls] is not supported.\n";
# print "\nCrystaObject::ConvertLatticeIndexByLatticeSystem: Warning: Lattice system [$ls] is not supported.\n";
}
return ($IsConverted, $ix, $iy, $iz);
}
sub Metrics
{
my ($this) = @_;
my $g11 = $this->{g11};
my $g12 = $this->{g12};
my $g13 = $this->{g13};
my $g21 = $this->{g21};
my $g22 = $this->{g22};
my $g23 = $this->{g23};
my $g31 = $this->{g31};
my $g32 = $this->{g32};
my $g33 = $this->{g33};
return ($g11, $g12, $g13,
$g21, $g22, $g23,
$g31, $g32, $g33)
}
sub GetPerpendicularHKL
{
my ($this, $h1, $k1, $l1, $h2, $k2, $l2) = @_;
my @v1 = $this->GetVectorFromHKL($h1, $k1, $l1);
my @v2 = $this->GetVectorFromHKL($h2, $k2, $l2);
my @v3 = Sci::VectorProduct(\@v1, \@v2);
# my @v3 = ($v1[1]*$v2[2] - $v1[2]*$v2[1],
# $v1[2]*$v2[0] - $v1[0]*$v2[2],
# $v1[0]*$v2[1] - $v1[1]*$v2[0]);
#printf "v1=(%7.5f, %7.5f, %7.5f)\n", @v1;
#printf "v2=(%7.5f, %7.5f, %7.5f)\n", @v2;
#printf "v3=(%7.5f, %7.5f, %7.5f)\n", @v3;
my $R = new Math::MatrixReal(3, 3);
$R->assign(1, 1, $this->{Ra11});
$R->assign(1, 2, $this->{Ra21});
$R->assign(1, 3, $this->{Ra31});
$R->assign(2, 1, $this->{Ra12});
$R->assign(2, 2, $this->{Ra22});
$R->assign(2, 3, $this->{Ra32});
$R->assign(3, 1, $this->{Ra13});
$R->assign(3, 2, $this->{Ra23});
$R->assign(3, 3, $this->{Ra33});
#print "R=$R\n";
my $RR = $R->inverse();
#print "RR=$RR\n";
my @hkl = ($RR->element(1,1)*$v3[0] + $RR->element(1,2)*$v3[1] + $RR->element(1,3)*$v3[2],
$RR->element(2,1)*$v3[0] + $RR->element(2,2)*$v3[1] + $RR->element(2,3)*$v3[2],
$RR->element(3,1)*$v3[0] + $RR->element(3,2)*$v3[1] + $RR->element(3,3)*$v3[2]);
my ($min, $max) = Utils::CalMinMax(\@hkl, 1);
#printf "hkl=%6.3f %6.3f %6.3f\n", @hkl;
$min = 1.0e30;
for(my $i = 0 ; $i < 3 ; $i++) {
$hkl[$i] = 0.0 if(abs($hkl[$i]) < $max * 1.0e-15);
$min = $hkl[$i] if($hkl[$i] != 0 and $min > $hkl[$i]);
}
if($min < 1.0e30) {
for(my $i = 0 ; $i < 3 ; $i++) {
$hkl[$i] /= $min;
}
}
#printf "hkl=%6.3f %6.3f %6.3f\n", @hkl;
#my @v32 = $this->GetVectorFromHKL(@hkl);
#printf "v3=(%7.5f, %7.5f, %7.5f)\n", @v32;
return @hkl;
}
sub GetVectorFromHKL
{
my ($this, $h, $k, $l) = @_;
my @V = ($this->{Ra11} * $h + $this->{Ra21} * $k + $this->{Ra31} * $l,
$this->{Ra12} * $h + $this->{Ra22} * $k + $this->{Ra32} * $l,
$this->{Ra13} * $h + $this->{Ra23} * $k + $this->{Ra33} * $l);
#print "hkl=$h,$k,$l: ", join(',', @V), ")\n";
return @V;
}
sub RMetrics
{
my ($this) = @_;
my $Rg11 = $this->{Rg11};
my $Rg12 = $this->{Rg12};
my $Rg13 = $this->{Rg13};
my $Rg21 = $this->{Rg21};
my $Rg22 = $this->{Rg22};
my $Rg23 = $this->{Rg23};
my $Rg31 = $this->{Rg31};
my $Rg32 = $this->{Rg32};
my $Rg33 = $this->{Rg33};
return ($Rg11, $Rg12, $Rg13,
$Rg21, $Rg22, $Rg23,
$Rg31, $Rg32, $Rg33)
}
# @v in Math::Vector::Real
sub SetLatticeVectorsByVectors
{
my ($this, @v) = @_;
#print "v1=", $v[0], "\n";
#print "v2=", $v[1], "\n";
#print "v3=", $v[2], "\n";
$this->SetLatticeVector(
$v[0][0], $v[0][1], $v[0][2],
$v[1][0], $v[1][1], $v[1][2],
$v[2][0], $v[2][1], $v[2][2]);
}
sub SetLatticeVectors
{
my ($this, $a11, $a12, $a13,
$a21, $a22, $a23, $a31, $a32, $a33) = @_;
return $this->SetLatticeVector(
$a11, $a12, $a13,
$a21, $a22, $a23,
$a31, $a32, $a33);
}
sub SetLatticeVector
{
my ($this, $a11, $a12, $a13,
$a21, $a22, $a23, $a31, $a32, $a33) = @_;
#print "($this, $a11, $a12, $a13, $a21, $a22, $a23, $a31, $a32, $a33)\n";
$this->{a11} = $a11;
$this->{a12} = $a12;
$this->{a13} = $a13;
$this->{a21} = $a21;
$this->{a22} = $a22;
$this->{a23} = $a23;
$this->{a31} = $a31;
$this->{a32} = $a32;
$this->{a33} = $a33;
#print "SetLatticeVector: a1=$this->{a11},$this->{a12},$this->{a13}\n";
#print "SetLatticeVector: a2=$this->{a21},$this->{a22},$this->{a23}\n";
#print "SetLatticeVector: a3=$this->{a31},$this->{a32},$this->{a33}\n";
$this->CalMetricsFromLatticeVector();
$this->CalLatticeParametersFromMetrics();
return;
}
sub SetMetrics
{
my ($this,$g11,$g12,$g13,$g21,$g22,$g23,$g31,$g32,$g33)
= @_;
$this->{'g11'} = $g11;
$this->{'g12'} = $g12;
$this->{'g13'} = $g13;
$this->{'g21'} = $g21;
$this->{'g22'} = $g22;
$this->{'g23'} = $g23;
$this->{'g31'} = $g31;
$this->{'g32'} = $g32;
$this->{'g33'} = $g33;
return ($this,$g11,$g12,$g13,$g21,$g22,$g23,$g31,$g32,$g33);
}
sub CalLatticeParametersFromMetrics
{
my ($this) = @_;
my $a = sqrt($this->{'g11'});
my $b = sqrt($this->{'g22'});
my $c = sqrt($this->{'g33'});
my $alpha = $this->{'g23'} / $b / $c;
$alpha = Sci::dacos($alpha);
my $beta = $this->{'g13'} / $a / $c;
$beta = Sci::dacos($beta);
my $gamma = $this->{'g12'} / $a / $b;
$gamma = Sci::dacos($gamma);
$a = Sci::Round($a, 6);
$b = Sci::Round($b, 6);
$c = Sci::Round($c, 6);
#print "a,b,c=$a,$b,$c\n";
$alpha = Sci::Round($alpha, 4);
$beta = Sci::Round($beta , 4);
$gamma = Sci::Round($gamma, 4);
#print "alpha=$alpha,$beta,$gamma\n";
$this->SetLatticeParameters($a,$b,$c,$alpha,$beta,$gamma);
}
sub CalMetricsFromLatticeVector
{
my ($this) = @_;
my ($a11, $a12, $a13, $a21, $a22, $a23, $a31, $a32, $a33) = (
$this->{a11}, $this->{a12}, $this->{a13},
$this->{a21}, $this->{a22}, $this->{a23},
$this->{a31}, $this->{a32}, $this->{a33});
$this->{g11} = $a11*$a11 + $a12*$a12 + $a13*$a13;
$this->{g22} = $a21*$a21 + $a22*$a22 + $a23*$a23;
$this->{g33} = $a31*$a31 + $a32*$a32 + $a33*$a33;
$this->{g12} = $this->{g21} = $a11*$a21 + $a12*$a22 + $a13*$a23;
$this->{g23} = $this->{g32} = $a21*$a31 + $a22*$a32 + $a23*$a33;
$this->{g31} = $this->{g13} = $a31*$a11 + $a32*$a12 + $a33*$a13;
return ($this->{g11}, $this->{g12}, $this->{g13},
$this->{g21}, $this->{g22}, $this->{g23},
$this->{g31}, $this->{g32}, $this->{g33})
}
sub CalculateLatticeConstantFromVector
{
my ($this) = @_;
my $a11 = $this->{a11};
my $a12 = $this->{a12};
my $a13 = $this->{a13};
my $a21 = $this->{a21};
my $a22 = $this->{a22};
my $a23 = $this->{a23};
my $a31 = $this->{a31};
my $a32 = $this->{a32};
my $a33 = $this->{a33};
#print "a1 ($a11, $a12, $a13)\n";
#print "a2 ($a21, $a22, $a23)\n";
#print "a3 ($a31, $a32, $a33)\n";
my $a = sqrt($a11*$a11 + $a12*$a12 + $a13*$a13);
my $b = sqrt($a21*$a21 + $a22*$a22 + $a23*$a23);
my $c = sqrt($a31*$a31 + $a32*$a32 + $a33*$a33);
#print "a=$a, $b, $c\n";
my $ab = $a11*$a21 + $a12*$a22 + $a13*$a23;
my $bc = $a21*$a31 + $a22*$a32 + $a23*$a33;
my $ca = $a31*$a11 + $a32*$a12 + $a33*$a13;
#print "ab=$ab, $bc, $ca\n";
#exit;
my $cosg = $ab / $a / $b;
my $cosa = $bc / $b / $c;
my $cosb = $ca / $c / $a;
my $alpha = Sci::dacos($cosa);
my $beta = Sci::dacos($cosb);
my $gamma = Sci::dacos($cosg);
$a = Utils::Round($a+0.0, 6); #sprintf("%11.6f", $a) + 0.0;
$b = Utils::Round($b+0.0, 6); #sprintf("%11.6f", $b) + 0.0;
$c = Utils::Round($c+0.0, 6); #sprintf("%11.6f", $c) + 0.0;
$alpha = Utils::Round($alpha+0.0, 6); #sprintf("%12.6f", $alpha) + 0.0;
$beta = Utils::Round($beta+0.0, 6); #sprintf("%12.6f", $beta) + 0.0;
$gamma = Utils::Round($gamma+0.0, 6); #sprintf("%12.6f", $gamma) + 0.0;
#print "a=$a $b $c $alpha $beta $gamma\n";
return ($a,$b,$c,$alpha,$beta,$gamma);
}
sub FindNearestEquivalentFractionCoordinate
{
my ($this, $x2, $x) = @_;
my $xret = $x2;
my $mind = abs($x2 - $x);
if($mind > abs($x2+1.0 - $x)) {
$xret = $x2+1.0;
$mind = abs($x2+1.0 - $x);
}
if($mind > abs($x2-1.0 - $x)) {
$xret = $x2-1.0;
$mind = abs($x2-1.0 - $x);
}
return $xret;
}
sub FractionalToCartesianByOutputMode
{
my ($this, $x,$y,$z) = @_;
my ($a11, $a12, $a13, $a21, $a22, $a23, $a31, $a32, $a33);
if($this->{'OutputMode'} eq 'expanded' or
$this->{'OutputMode'} eq 'choose' ) {
my $SuperLattice = $this->{'CSuperLattice'};
$a11 = $SuperLattice->{"a11"};
$a12 = $SuperLattice->{"a12"};
$a13 = $SuperLattice->{"a13"};
$a21 = $SuperLattice->{"a21"};
$a22 = $SuperLattice->{"a22"};
$a23 = $SuperLattice->{"a23"};
$a31 = $SuperLattice->{"a31"};
$a32 = $SuperLattice->{"a32"};
$a33 = $SuperLattice->{"a33"};
}
else {
$a11 = $this->{"a11"};
$a12 = $this->{"a12"};
$a13 = $this->{"a13"};
$a21 = $this->{"a21"};
$a22 = $this->{"a22"};
$a23 = $this->{"a23"};
$a31 = $this->{"a31"};
$a32 = $this->{"a32"};
$a33 = $this->{"a33"};
}
my $xc = $a11*$x + $a21*$y + $a31*$z;
my $yc = $a12*$x + $a22*$y + $a32*$z;
my $zc = $a13*$x + $a23*$y + $a33*$z;
return ($xc,$yc,$zc);
}
sub FractionalToCartesian
{
my ($this, $x,$y,$z) = @_;
my $a11 = $this->{a11};
my $a12 = $this->{a12};
my $a13 = $this->{a13};
my $a21 = $this->{a21};
my $a22 = $this->{a22};
my $a23 = $this->{a23};
my $a31 = $this->{a31};
my $a32 = $this->{a32};
my $a33 = $this->{a33};
my $xc = $a11*$x + $a21*$y + $a31*$z;
my $yc = $a12*$x + $a22*$y + $a32*$z;
my $zc = $a13*$x + $a23*$y + $a33*$z;
return ($xc,$yc,$zc);
}
sub CartesianToFractional
{
my ($this, $xc,$yc,$zc) = @_;
my $A = new Math::MatrixReal(3, 3);
$A->assign(1, 1, $this->{"a11"});
$A->assign(1, 2, $this->{"a21"});
$A->assign(1, 3, $this->{"a31"});
# $A->assign(1, 2, $this->{"a12"});
# $A->assign(1, 3, $this->{"a13"});
$A->assign(2, 1, $this->{"a12"});
# $A->assign(2, 1, $this->{"a21"});
$A->assign(2, 2, $this->{"a22"});
$A->assign(2, 3, $this->{"a32"});
# $A->assign(2, 3, $this->{"a23"});
$A->assign(3, 1, $this->{"a13"});
$A->assign(3, 2, $this->{"a23"});
# $A->assign(3, 1, $this->{"a31"});
# $A->assign(3, 2, $this->{"a32"});
$A->assign(3, 3, $this->{"a33"});
# print "A:
\n";print $A;print "
\n";
my $RA = $A->inverse();
# print "RA:
\n";print $RA;print "
\n";
my $X = new Math::MatrixReal(3,1);
$X->assign(1, 1, $xc);
$X->assign(2, 1, $yc);
$X->assign(3, 1, $zc);
# print "X:
\n";print $X;print "
\n";
my $XX = $RA * $X;
# print "XX:
\n";print $XX;print "
\n";
return ($XX->element(1,1),
$XX->element(2,1),
$XX->element(3,1));
}
sub CalculateVolume
{
my ($this) = @_;
my $a = $this->{"a"};
my $b = $this->{"b"};
my $c = $this->{"c"};
my $alpha = $this->{"alpha"};
my $beta = $this->{"beta"};
my $gamma = $this->{"gamma"};
my $cosa = Sci::dcos($alpha);
my $cosb = Sci::dcos($beta);
my $cosg = Sci::dcos($gamma);
my $vol = 1.0 - $cosa*$cosa - $cosb*$cosb - $cosg*$cosg
+ 2.0 * $cosa*$cosb*$cosg;
$vol = $a*$b*$c*sqrt($vol);
$vol = sprintf("%14.6g", $vol) + 0.0;
$this->{"Volume"} = Sci::Round($vol, 4);
return $vol;
}
sub CalculateReciprocalVolume
{
my ($this) = @_;
my $a = $this->{Ra};
my $b = $this->{Rb};
my $c = $this->{Rc};
my $alpha = $this->{Ralpha};
my $beta = $this->{Rbeta};
my $gamma = $this->{Rgamma};
my $cosa = Sci::dcos($alpha);
my $cosb = Sci::dcos($beta);
my $cosg = Sci::dcos($gamma);
my $Rvol = 1.0 - $cosa*$cosa - $cosb*$cosb - $cosg*$cosg
+ 2.0 * $cosa*$cosb*$cosg;
$Rvol = $a*$b*$c*sqrt($Rvol);
$Rvol = sprintf("%14.6g", $Rvol) + 0.0;
$this->{RVolume} = Sci::Round($Rvol, 4);
return $Rvol;
}
sub FindAtomTypeByName
{
my ($this, $atomname) = @_;
my $atomlistref = $this->{"RefAtomTypeList"};
my $IsPartial = 0;
for(my $i = 0 ; $i < @$atomlistref ; $i++) {
my $type = $atomlistref->[$i];
my $name = $type->Label();
#print "i=$i: name=$name\n";
if($name =~ /\[.*?\]/) {
$IsPartial = 1;
last;
}
}
#print "IsP=$IsPartial\n";
if(!$IsPartial) {
$atomname =~ s/([\.\d+-]+)//;
}
for(my $i = 0 ; $i < @$atomlistref ; $i++) {
my $type = $atomlistref->[$i];
my $name;
if($IsPartial) {
$name = $type->Label();
}
else {
$name = $type->AtomNameOnly();
$name =~ s/([\d+-]+)//;
}
#print "name: $atomname <=> $name\n";
return $type if(uc $name eq uc $atomname);
}
return undef;
}
sub SplitAtomType
{
my ($this) = @_;
my @AtomTypeList = $this->GetCAtomTypeList();
my @AtomSiteList = $this->GetCAsymmetricAtomSiteList();
$this->ClearAtomTypeList();
$this->ClearAtomSiteList();
my %at;
#my $n = @AtomSiteList;
#print "nsite=$n\n";
for(my $i = 0 ; $i < @AtomSiteList ; $i++) {
my $site = $AtomSiteList[$i];
my $name = $site->AtomNameOnly();
my $label = $site->Label();
my $charge = $site->Charge();
my $occ = $site->Occupancy();
my ($x, $y, $z) = $site->Position();
my ($fx,$fy,$fz) = $site->Force();
$at{$name}++;
my $name2 = "$name\{" . Utils::IntToAlpha($at{$name} - 1) . "\}$charge";
#print "name2=$name2\n";
$this->AddAtomSite($label, $name2, $x, $y, $z, $occ, $fx, $fy, $fz);
# $this->AddAtomType($name2, 1);
#print "Site (indep type) $i: $name2 (Z=$charge) (occ=$occ): $key\n";
}
$this->ExpandCoordinates();
}
sub SplitPartialSites
{
my ($this) = @_;
my @AtomTypeList = $this->GetCAtomTypeList();
my @AtomSiteList = $this->GetCAsymmetricAtomSiteList();
# my @ExpandedAtomSiteList = $this->GetCExpandedAtomSiteList();
$this->ClearAtomTypeList();
$this->ClearAtomSiteList();
my %at;
for(my $i = 0 ; $i < @AtomSiteList ; $i++) {
my $site = $AtomSiteList[$i];
my $name = $site->AtomNameOnly();
my $label = $site->Label();
my $charge = $site->Charge();
my $occ = $site->Occupancy();
my ($x, $y, $z) = $site->Position();
my ($fx,$fy,$fz) = $site->Force();
my $key;
if(abs($occ - 1.0) < 1.0e-4) {
$key = "$name$charge";
}
else {
$key = sprintf("$name${charge}[%6.4f]", $occ);
}
$this->AddAtomSite($label, $key, $x, $y, $z, $occ, $fx, $fy, $fz);
next if($at{$key});
$at{$key}++;
$this->AddAtomType($key, 1);
print "Site (indep type) $i: $name (Z=$charge) (occ=$occ): $key\n";
}
my @ATL = $this->GetCAtomTypeList();
for(my $i = 0 ; $i < @ATL ; $i++) {
my $type = $ATL[$i];
my $nameraw = $type->Label();
my $name = $type->AtomNameOnly();
my $charge = $type->Charge();
my $occ = $type->Occupancy();
my $key = sprintf("$name${charge}[%6.4f]", $occ);
print "Type $i: {$nameraw} $name (Z=$charge) (occ=$occ): $key\n";
}
my @ASL = $this->GetCAsymmetricAtomSiteList();
for(my $i = 0 ; $i < @ASL ; $i++) {
my $site = $ASL[$i];
my $nameraw = $site->AtomName();
my $name = $site->AtomNameOnly();
my $charge = $site->Charge();
my $occ = $site->Occupancy();
my $key = sprintf("$name${charge}[%6.4f]", $occ);
print "c $i: {$nameraw}: $name (Z=$charge) (occ=$occ): $key\n";
}
$this->ExpandCoordinates();
#exit;
}
sub TotalCharge
{
my ($this) = @_;
$this->ExpandCoordinates();
my @ExpandedAtomSiteList = $this->GetCExpandedAtomSiteList();
my $nTotalExpandedAtomSite = $this->nTotalExpandedAtomSite();
my $Z = 0.0;
for(my $i = 0 ; $i < $nTotalExpandedAtomSite ; $i++) {
my $atom = $ExpandedAtomSiteList[$i];
my $charge = $atom->Charge();
my $occ = $atom->Occupancy();
$Z += $occ * $charge;
}
return $Z;
}
sub CheckTotalCharge
{
my ($this, $EPS) = @_;
$EPS = 1.0e-3 if($EPS eq '');
my $Z = $this->TotalCharge();
if(abs($Z) > $EPS) {
print "Error in CrystalObject::CheckTotalCharge: Total charge ($Z) is greater than EPS=$EPS.\n";
exit 0;
}
return 1;
}
sub GetCBravaisLatticeAsymmetricAtomSiteList
{
my ($this) = @_;
$this->ExpandCoordinates();
my $SPG = $this->GetCSpaceGroup();
my @ExpandedAtomSiteList = $this->GetCExpandedAtomSiteList();
my $nTotalExpandedAtomSite = $this->nTotalExpandedAtomSite();
my (@Tx, @Ty, @Tz);
my $nT = $SPG->nTranslation();
for(my $it = 1 ; $it <= $nT ; $it++) {
($Tx[$it-1], $Ty[$it-1], $Tz[$it-1]) = $SPG->TranslationVector($it);
}
my @AtomSite;
my $c = 0;
for(my $i = 0 ; $i < $nTotalExpandedAtomSite ; $i++) {
my $atom = $ExpandedAtomSiteList[$i];
my $label = $atom->Label();
my $atomname = $atom->AtomNameOnly();
my $charge = $atom->Charge();
my ($x,$y,$z) = $atom->Position();
my $occ = $atom->Occupancy();
#print "$i: $atomname $charge [$label]: ($x, $y, $z) [occ=$occ]\n";
my $IsExist = 0;
for(my $it = 1 ; $it < $nT ; $it++) {
my ($tx, $ty, $tz) = ($Tx[$it], $Ty[$it], $Tz[$it]);
my $x1 = Utils::Reduce01($x + $tx);
my $y1 = Utils::Reduce01($y + $ty);
my $z1 = Utils::Reduce01($z + $tz);
for(my $j = 0 ; $j < $c ; $j++) {
my $atom2 = $AtomSite[$j];
my ($x2, $y2, $z2) = $atom2->Position();
my $dis = $this->GetNearestInterAtomicDistance($x1, $y1, $z1, $x2, $y2, $z2);
#print "$i-$it-$j: ($x,$y,$z): ($x1,$y1,$z1)-($x2,$y2,$z2)=$dis\n";
if($dis < 1.0e-6) {
$IsExist = 1;
last;
}
}
last if($IsExist);
}
if(!$IsExist) {
$AtomSite[$c] = $atom;
$c++;
}
}
return @AtomSite;
}
sub ClearAtomSiteList
{
my ($this) = @_;
$this->{RefAsymmetricAtomSiteList} = [];
$this->{RefExpandedAtomSiteList} = [];
}
sub SetAtomSite
{
my ($this, $iatom, $label, $atomname, $x, $y, $z, $occ, $fx, $fy, $fz) = @_;
my $at = $this->{RefAsymmetricAtomSiteList}[$iatom];
return undef if(!$at);
$at->SetLabel($label) if(defined $label);
$at->SetAtomName($atomname) if(defined $atomname);
$at->SetPosition($x, $y, $z) if(defined $x);
$at->SetOccupancy($occ) if(defined $occ);
$at->SetForce($fx, $fy, $fz) if(defined $fz);
# $at->SetId(@$atomlistref + 1);
return $at;
}
sub AddAtomSite
{
my ($this, $label, $atomname, $x, $y, $z, $occ, $fx, $fy, $fz) = @_;
#print "AddAtomSite: F=$fx, $y, $fz\n";
# return $this->{nAsymmetricAtomSite} if($atomname eq '');
my $atomlistref = $this->{RefAsymmetricAtomSiteList};
if($atomlistref) {
my $at = new AtomSite();
$at->SetLabel($label);
$at->SetAtomName($atomname);
$at->SetPosition($x, $y, $z);
$at->SetForce($fx, $fy, $fz) if(defined $fz);
$at->SetOccupancy($occ);
$at->SetId(@$atomlistref + 1);
push(@$atomlistref, $at);
$this->{nAsymmetricAtomSite} = scalar @$atomlistref;
}
else {
my @atomlist;
my $at = new AtomSite();
$at->SetLabel($label);
$at->SetAtomName($atomname);
$at->SetPosition($x, $y, $z);
$at->SetForce($fx, $fy, $fz) if(defined $fz);
$at->SetOccupancy($occ);
$at->SetId(1);
push(@atomlist, $at);
$this->{nAsymmetricAtomSite} = @atomlist;
$this->{RefAsymmetricAtomSiteList} = \@atomlist;
}
$this->AddAtomType($atomname, undef, $occ) if($atomname ne '');
return $this->{nAsymmetricAtomSite};
}
sub GetimaxByRmax
{
my ($this, $Rmax) = @_;
my ($a,$b,$c,$alpha,$beta,$gamma) = $this->LatticeParameters();
return (int($Rmax/$a)+1, int($Rmax/$b)+1, int($Rmax/$c)+1);
}
sub FindCoordinationsForAllAtoms
{
# $iTargetSites: 0: Expanded 1: Asymmetric sites
my ($Crystal, $iTargetSites, $CoordRMax, $IsPrint) = @_;
$IsPrint = 0 if(!defined $IsPrint);
my (@SiteList, $nSite);
if($iTargetSites == 0) {
@SiteList = $Crystal->GetCExpandedAtomSiteList();
}
else {
@SiteList = $Crystal->GetCAsymmetricAtomSiteList();
}
$nSite = @SiteList;
for(my $i = 0 ; $i < $nSite ; $i++) {
my $pAtom0 = $SiteList[$i];
my $label0 = $pAtom0->Label();
my $type0 = $pAtom0->AtomName();
my ($x0,$y0,$z0) = $pAtom0->Position(1);
# my $occupancy0 = $pAtom0->Occupancy();
# my $iAtomType0 = $Crystal->FindiAtomType($type);
my $id0 = ($iTargetSites == 0)? $i : $pAtom0->{IdAsymmetricAtomSite} - 1;
print "Coordination arround $i-th atom [$label0: $type0 (id=$id0)] ($x0, $y0, $z0)\n" if($IsPrint);
$pAtom0->{idSite}= ($iTargetSites == 0)? $i : $pAtom0->{IdAsymmetricAtomSite} - 1;
$Crystal->FindCoordinationsAroundAtom($pAtom0, $iTargetSites, $CoordRMax);
my $idSite0 = $pAtom0->{idSite};
my $pCoord = $pAtom0->{pCoordAtoms};
for(my $i = 0 ; $i < @$pCoord ; $i++) {
my $pCA = $pCoord->[$i];
my $iSite1 = $pCA->{iSite};
my $idSite1 = $pCA->{idSite};
my $label1 = $pCA->Label();
my $type1 = $pCA->AtomName();
my ($x1, $y1, $z1) = $pCA->Position();
my ($X1, $Y1, $Z1) = @{$pCA->{pCoord}};
my $r = $pCA->{r};
printf " %03d: [%3s: %4s (id = %03d)] (%8.4f, %8.4f, %8.4f): r = %8.4f frac:(%6.3f, %6.3f, %6.3f)\n",
$iSite1, $label1, $type1, $idSite1, $X1, $Y1, $Z1, $r, $x1, $y1, $z1 if($IsPrint);
}
}
return \@SiteList;
}
sub FindCoordinationsAroundAtom
{
my ($Crystal, $pAtom0, $iTargetSites, $CoordRMax) = @_;
my ($x0,$y0,$z0) = $pAtom0->Position(1);
my ($ixmax, $iymax, $izmax) = $Crystal->GetimaxByRmax($CoordRMax);
my @ExpandedAtomSiteList = $Crystal->GetCExpandedAtomSiteList();
my $nTotalExpandedAtomSite = $Crystal->nTotalExpandedAtomSite();
my @CoordAtoms;
my $c = 0;
for(my $i = 0 ; $i < $nTotalExpandedAtomSite ; $i++) {
my $pAtom1 = $ExpandedAtomSiteList[$i];
my $label1 = $pAtom1->Label();
my $type1 = $pAtom1->AtomName();
my $id1 = ($iTargetSites == 0)? $i : $pAtom1->{IdAsymmetricAtomSite} - 1;
my ($x1,$y1,$z1) = $pAtom1->Position(1);
my @dislist = $Crystal->GetNearestInterAtomicDistances(
$x0, $y0, $z0, $x1, $y1, $z1, $CoordRMax, 0, $ixmax, $iymax, $izmax);
my $p = $ExpandedAtomSiteList[$id1];
my $n1 = $pAtom1->AtomName();
my $n = $p->AtomName();
my $l = $p->Label();
if($n1 ne $n) {
print "Error: Atom name mismatch ($n1 ne $n) for ($label1 vs $l) (i=$i vs id=$id1)\n";
exit;
}
for(my $j = 0 ; $j < @dislist ; $j++) {
my $piCell = $dislist[$j]->{piCell};
my $r = $dislist[$j]->{r};
$CoordAtoms[$c] = new AtomSite;
%{$CoordAtoms[$c]} = %$pAtom1;
$CoordAtoms[$c]->{pSite} = $pAtom1;
$CoordAtoms[$c]->{iSite} = $i;
$CoordAtoms[$c]->{idSite} = $id1;
$CoordAtoms[$c]->{piCell} = $piCell;
$CoordAtoms[$c]->{pCoord} = [$x1 + $piCell->[0], $y1 + $piCell->[1], $z1 + $piCell->[2]];
$CoordAtoms[$c]->{r} = $r;
my $l = $CoordAtoms[$c];
#print " $i: [$label: $type (id = $id1)] ($l->{pCoord}[0], $l->{pCoord}[1], $l->{pCoord}[2]): r = $l->{r}\n";
$c++;
}
}
$pAtom0->{pCoordAtoms} = \@CoordAtoms;
return $pAtom0;
}
sub GetNearestInterAtomicDistances
{
my ($this, $x0, $y0, $z0, $x1, $y1, $z1, $CoordRMax, $AllowZero, $ixmax, $iymax, $izmax) = @_;
$AllowZero = 1 if(!defined $AllowZero);
if(!defined $ixmax) {
($ixmax, $iymax, $izmax) = $this->GetimaxByRmax($CoordRMax);
}
my @list;
my $c = 0;
for(my $ix = -$ixmax ; $ix <= $ixmax ; $ix++) {
for(my $iy = -$iymax ; $iy <= $iymax ; $iy++) {
for(my $iz = -$izmax ; $iz <= $izmax ; $iz++) {
my $r = $this->GetInterAtomicDistance($x0, $y0, $z0, $x1+$ix, $y1+$iy, $z1+$iz);
#print "$ix,$iy,$iz: r=$r\n";
if($r < $CoordRMax) {
next if(!$AllowZero and $r < 1.0e-5);
$list[$c++] = {
piCell => [$ix, $iy, $iz],
r => $r,
};
}
}
}
}
return @list;
}
sub GetNearestInterAtomicDistance
{
my ($this, $x0, $y0, $z0, $x1, $y1, $z1, $AllowZero) = @_;
$AllowZero = 1 if(!defined $AllowZero);
$x0 = $this->FindNearestEquivalentFractionCoordinate($x0, $x1);
$y0 = $this->FindNearestEquivalentFractionCoordinate($y0, $y1);
$z0 = $this->FindNearestEquivalentFractionCoordinate($z0, $z1);
my $r = $this->GetInterAtomicDistance($x0, $y0, $z0, $x1, $y1, $z1);
#print "r=$r\n";
if(!$AllowZero and $r <= 1.0e-5) {
$r = 1.0e10;
for(my $ix = -1 ; $ix <= 1 ; $ix++) {
for(my $iy = -1 ; $iy <= 1 ; $iy++) {
for(my $iz = -1 ; $iz <= 1 ; $iz++) {
my $r2 = $this->GetInterAtomicDistance($x0, $y0, $z0, $x1+$ix, $y1+$iy, $z1+$iz);
#print "$ix,$iy,$iz: r2=$r2\n";
$r = $r2 if($r2 > 1.0e-5 and $r > $r2);
}
}
}
}
return $r;
}
sub GetInterAtomicDistance
{
my ($this, $x0, $y0, $z0, $x1, $y1, $z1) = @_;;
my $dx = $x1 - $x0;
my $dy = $y1 - $y0;
my $dz = $z1 - $z0;
my $g11 = $this->{"g11"};
my $g22 = $this->{"g22"};
my $g33 = $this->{"g33"};
my $g12 = 2.0* $this->{"g12"};
my $g13 = 2.0* $this->{"g13"};
my $g23 = 2.0* $this->{"g23"};
my $r2 = $g11*$dx*$dx + $g22*$dy*$dy + $g33*$dz*$dz
+ $g12*$dx*$dy + $g13*$dx*$dz + $g23*$dy*$dz;
return sqrt($r2);
}
sub GetInterAtomicAngle
{
my ($this, $x0, $y0, $z0, $x1, $y1, $z1,
$x2, $y2, $z2) = @_;;
my $dis01 = $this->GetInterAtomicDistance($x0,$y0,$z0,$x1,$y1,$z1);
return 0.0 if($dis01 == 0.0);
my $dis02 = $this->GetInterAtomicDistance($x0,$y0,$z0,$x2,$y2,$z2);
return 0.0 if($dis02 == 0.0);
#print "dis: $dis01, $dis02\n";
my $dx01 = $x1 - $x0;
my $dy01 = $y1 - $y0;
my $dz01 = $z1 - $z0;
my $dx02 = $x2 - $x0;
my $dy02 = $y2 - $y0;
my $dz02 = $z2 - $z0;
my $g11 = $this->{"g11"};
my $g22 = $this->{"g22"};
my $g33 = $this->{"g33"};
my $g12 = 2.0* $this->{"g12"};
my $g13 = 2.0* $this->{"g13"};
my $g23 = 2.0* $this->{"g23"};
my $ip = $g11*$dx01*$dx02 + $g22*$dy01*$dy02 + $g33*$dz01*$dz02
+ $g12*$dx01*$dy02 + $g13*$dx01*$dz02
+ $g12*$dy01*$dx02 + $g23*$dy01*$dz02
+ $g13*$dz01*$dx02 + $g23*$dz01*$dy02;
my $cosa = $ip / $dis01 / $dis02;
#print "cos: $cosa
";
my $angle = Sci::dacos($cosa);
#print "cos: $cosa $angle
";
$angle = 360.0 - $angle if($angle > 180.0001);
return $angle;
}
#iAtomTypeは1から始まる
sub FindiAtomType
{
my ($this, $atomname) = @_;
#print "F:$atomname\n";
$atomname =~ s/^(\D+)(\d.*)$/$1/;
$atomname = uc $atomname;
#print "F:$atomname\n";
my @AtomTypeList = $this->GetCAtomTypeList();
my $nAtomType = @AtomTypeList;
for(my $i = 0 ; $i < $nAtomType ; $i++) {
my $name = $AtomTypeList[$i]->AtomName();
$name =~ s/^(\D+)(\d.*)$/$1/;
#print "$i: $name\n";
return $i+1 if($atomname eq uc $name);
}
return -1;
}
sub CalCartesianCoordinate
{
my ($this, $x, $y, $z) = @_;
my ($ax, $ay, $az, $bx, $by, $bz, $cx, $cy, $cz) = $this->LatticeVectorsByOutputMode(0);
my $rx = $ax * $x + $bx * $y + $cx * $z;
my $ry = $ay * $x + $by * $y + $cy * $z;
my $rz = $az * $x + $bz * $y + $cz * $z;
return ($rx, $ry, $rz);
}
sub CalFractionalCoordinate
{
my ($this, $rx, $ry, $rz) = @_;
my ($ax, $ay, $az, $bx, $by, $bz, $cx, $cy, $cz) = $this->LatticeVectorsByOutputMode(0);
my $T = new Math::MatrixReal(3, 3);
$T->assign(1, 1, $ax);
$T->assign(1, 2, $bx);
$T->assign(1, 3, $cx);
$T->assign(2, 1, $ay);
$T->assign(2, 2, $by);
$T->assign(2, 3, $cy);
$T->assign(3, 1, $az);
$T->assign(3, 2, $bz);
$T->assign(3, 3, $cz);
my $RT = $T->inverse();
my $x = $RT->element(1,1) * $rx + $RT->element(1,2) * $ry + $RT->element(1,3) * $rz;
my $y = $RT->element(2,1) * $rx + $RT->element(2,2) * $ry + $RT->element(2,3) * $rz;
my $z = $RT->element(3,1) * $rx + $RT->element(3,2) * $ry + $RT->element(3,3) * $rz;
return ($x, $y, $z);
}
sub FindIdenticalSite
{
my ($this, $atomname, $x, $y, $z, $TolD, $iSiteStart) = @_;
$TolD = 0.01 if(!defined $TolD);
$iSiteStart = 0 if(!defined $iSiteStart);
$atomname =~ s/\d[\+\-]?//;
my @AtomSiteList = $this->GetCExpandedAtomSiteList();
my $nSite = @AtomSiteList;
for(my $ia = $iSiteStart ; $ia < $nSite ; $ia++) {
my $site = $AtomSiteList[$ia];
my $name = $site->AtomNameOnly();
next if(uc $name ne uc $atomname);
my ($x1, $y1, $z1) = $site->Position();
my $dis = $this->GetNearestInterAtomicDistance($x, $y, $z, $x1, $y1, $z1);
return $ia if($dis < $TolD);
}
return undef;
}
sub FindIdenticalAsymmetricSite
{
my ($this, $atomname, $x, $y, $z) = @_;
$atomname =~ s/\d[\+\-]?//;
my @AtomSiteList = $this->GetCAsymmetricAtomSiteList();
my $nAsymmetricAtomSite = @AtomSiteList;
for(my $ia = 0 ; $ia < $nAsymmetricAtomSite ; $ia++) {
my $site = $AtomSiteList[$ia];
my $name = $site->AtomNameOnly();
next if(uc $name ne uc $atomname);
my ($x1, $y1, $z1) = $site->Position();
my $dis = $this->GetNearestInterAtomicDistance($x, $y, $z, $x1, $y1, $z1);
return $ia if($dis < 0.01);
}
return undef;
}
sub ExpandCoordinates
{
my ($this, $DoTranslation) = @_;
$DoTranslation = 1 if(!defined $DoTranslation);
#初期化
$this->{nTotalExpandedAtomSite} = 0;
$this->{Density} = 0.0;
#同一原子の判断半径を決める
my ($a, $b, $c) = $this->LatticeParameters();
$a = $b if($a < $b);
$a = $c if($a < $c);
my $OverlapCriteriaRadius = 0.001 * $a;
#非対称構造中のパラメータ
my @AtomTypeList = $this->GetCAtomTypeList();
my $nAtomType = @AtomTypeList;
my @AtomSiteList = $this->GetCAsymmetricAtomSiteList();
my $nAsymmetricAtomSite = @AtomSiteList;
my $SPG = $this->GetCSpaceGroup();
my $nTranslation = $SPG->nTranslation();
$nTranslation = 1 if(!$DoTranslation);
my $nSymmetryOperation = $SPG->nSymmetryOperation();
my $nTotalSymmetry = $nTranslation * $nSymmetryOperation;
#print "nTotalSymmetry: $nTotalSymmetry\n";
#座標展開を開始
my @ExpandedAtomSiteList;
my @nExpandedAtomSite;
my $count = 0;
#print "nAAS=$nAsymmetricAtomSite\n";
for(my $ia = 0 ; $ia < $nAsymmetricAtomSite ; $ia++) {
$nExpandedAtomSite[$ia] = 0;
my $atom = $AtomSiteList[$ia];
my $label = $atom->Label();
my $AtomName = $atom->AtomName();
my $iAtomType = $this->FindiAtomType($AtomName);
my ($x,$y,$z) = $atom->Position(1);
my ($fx,$fy,$fz) = $atom->Force();
$x = CrystalObject::RoundSymmetricPosition($x);
$y = CrystalObject::RoundSymmetricPosition($y);
$z = CrystalObject::RoundSymmetricPosition($z);
my $Occupancy = $atom->Occupancy();
my ($vx,$vy,$vz) = $atom->Velocity();
$atom->SetId($ia+1);
$ExpandedAtomSiteList[$count] = new AtomSite();
$ExpandedAtomSiteList[$count]->SetIdAsymmetricAtomSite($ia+1);
$ExpandedAtomSiteList[$count]->SetiAtomType($iAtomType);
$ExpandedAtomSiteList[$count]->SetLabel($label);
$ExpandedAtomSiteList[$count]->SetAtomName($AtomName);
$ExpandedAtomSiteList[$count]->SetPosition($x,$y,$z);
$ExpandedAtomSiteList[$count]->SetForce($fx,$fy,$fz) if(defined $fx);
$ExpandedAtomSiteList[$count]->SetOccupancy($Occupancy);
$ExpandedAtomSiteList[$count]->SetVelocity($vx,$vy,$vz) if(defined $vx);
$nExpandedAtomSite[$ia] += 1;
$count++;
#print "($ia) x,y,z: $x,$y,$z
\n";
for(my $is = 0 ; $is < $nTotalSymmetry ; $is++) {
my ($xs, $ys, $zs) = $SPG->DoSymmetryOperation($is, $x, $y, $z, 1);
#print " ($is) xs,ys,zs: $xs,$ys,$zs
\n";
my ($fxs, $fys, $fzs) = $SPG->DoVelocitySymmetryOperation($is, $fx, $fy, $fz, 1);
my ($vxs, $vys, $vzs) = $SPG->DoVelocitySymmetryOperation($is, $vx, $vy, $vz, 1);
my $IsPassed = 1;
for(my $it = 0 ; $it < $count ; $it++) {
my $iat = $ExpandedAtomSiteList[$it]->iAtomType();
next if($iat != $iAtomType);
my ($x0,$y0,$z0) = $ExpandedAtomSiteList[$it]->Position(1);
my $dis = $this->GetNearestInterAtomicDistance(
$x0, $y0, $z0, $xs, $ys, $zs);
#print "$label(ia=$ia:is=$is:it=$it): ";
#printf "(%7.4f,%7.4f,%7.4f)-(%7.4f,%7.4f,%7.4f): dis: %f (%f)
\n",
# $x0, $y0, $z0, $xs, $ys, $zs, $dis, $OverlapCriteriaRadius;
if($dis < $OverlapCriteriaRadius) {
$IsPassed = 0;
last;
}
}
next unless($IsPassed);
$ExpandedAtomSiteList[$count] = new AtomSite();
$ExpandedAtomSiteList[$count]->SetIdAsymmetricAtomSite($ia+1);
$ExpandedAtomSiteList[$count]->SetiAtomType($iAtomType);
$ExpandedAtomSiteList[$count]->SetLabel($label);
$ExpandedAtomSiteList[$count]->SetAtomName($AtomName);
$ExpandedAtomSiteList[$count]->SetPosition($xs,$ys,$zs);
$ExpandedAtomSiteList[$count]->SetForce($fxs,$fys,$fzs);
$ExpandedAtomSiteList[$count]->SetOccupancy($Occupancy);
$ExpandedAtomSiteList[$count]->SetVelocity($vxs,$vys,$vzs);
#Idは1からはじまる
$ExpandedAtomSiteList[$count]->SetId($ia+1);
#print "v: ($vx,$vy,$vz) ($vxs, $vys, $vzs)
\n";
$nExpandedAtomSite[$ia] += 1;
$count++;
}
}
$this->{nTotalExpandedAtomSite} = $count;
$this->{RefExpandedAtomSiteList} = \@ExpandedAtomSiteList;
$this->{RefnExpandedAtomSite} = \@nExpandedAtomSite;
#多重度を、非対称原子に設定
for(my $ia = 0 ; $ia < @AtomSiteList ; $ia++) {
my $atom = $AtomSiteList[$ia];
#print "atom: ($ia): m=$nExpandedAtomSite[$ia]
\n";
$atom->SetMultiplicity($nExpandedAtomSite[$ia]);
}
#多重度を、全原子に設定
for(my $i = 0 ; $i < @ExpandedAtomSiteList ; $i++) {
my $atom = $ExpandedAtomSiteList[$i];
my $ia = $atom->Id() - 1;
#print "atom: ($i,$ia): m=$nExpandedAtomSite[$ia]
\n";
$atom->SetMultiplicity($nExpandedAtomSite[$ia]);
}
my $Density = $this->CalculateDensity();
$this->{Density} = $Density;
$this->AnalyzeChemicalComposition();
return;
}
sub AnalyzeChemicalComposition
{
my ($this) = @_;
my @ExpandedAtomSiteList = $this->GetCExpandedAtomSiteList();
my %AtomCount;
for(my $i = 0 ; $i < @ExpandedAtomSiteList ; $i++) {
my $site = $ExpandedAtomSiteList[$i];
my $name = $site->AtomNameOnly(1);
# $name =~ s/\{[^\{\}]*\}//g;
$AtomCount{$name} += $site->Occupancy();
#print "$i: $name [$AtomCount{$name}]\n";
}
my @Names = keys %AtomCount;
foreach my $key (@Names) {
$AtomCount{$key} = Utils::RoundParameter($AtomCount{$key}, 1.0e-4);
}
my $Composition = '';
my $MinN = 1000000;
for(my $i = 0 ; $i < @Names ; $i++) {
my $name = $Names[$i];
next if($name =~ /^(O|S|Se|Te|F|Cl|I|N|As)$/i);
my $n = $AtomCount{$name};
$MinN = $n if($MinN > $n);
$n = '' if($n == 1);
$Composition = "$Composition $name$n";
#print "i=$i: [$n] [$MinN] [$name]\n";
}
for(my $i = 0 ; $i < @Names ; $i++) {
my $name = $Names[$i];
next unless($name =~ /^(O|S|Se|Te|F|Cl|I|N|As)$/i);
my $n = $AtomCount{$name};
$MinN = $n if($MinN > $n);
$n = '' if($n == 1);
$Composition = "$Composition $name$n";
}
$Composition =~ s/^\s*//;
$this->SetSumChemicalComposition($Composition);
my @elements = Algorism::Factorize($MinN);
my $fac = 1;
#print "$MinN = ";
for(my $i = 0 ; $i < @elements ; $i++) {
my $num = $elements[$i];
#print " * " if($i > 0);
#print "$num";
my $IsCommonFactor = 1;
for(my $j = 0 ; $j < @Names ; $j++) {
my $name = $Names[$j];
my $n = $AtomCount{$name};
unless(Algorism::IsFactor($n, $num*$fac)) {
$IsCommonFactor = 0;
}
}
$fac *= $num if($IsCommonFactor);
#print "fac[$i]=$fac\n";
}
#print "\n";
$this->SetFormulaUnit($fac);
#print "Z = $fac\n";
$Composition = '';
for(my $i = 0 ; $i < @Names ; $i++) {
my $name = $Names[$i];
next if($name =~ /^(O|S|Se|Te|F|Cl|I|N|As|)$/i);
my $n = $AtomCount{$name} / $fac;
$n = '' if($n == 1);
$Composition = "$Composition $name$n";
}
for(my $i = 0 ; $i < @Names ; $i++) {
my $name = $Names[$i];
next unless($name =~ /^(O|S|Se|Te|F|Cl|I|N|As|)$/i);
my $n = $AtomCount{$name} / $fac;
$n = '' if($n == 1);
$Composition = "$Composition $name$n";
}
$Composition =~ s/^\s*//;
$this->SetChemicalComposition($Composition);
}
sub CalculateDensity
{
my ($this) = @_;
my $vol = $this->Volume();
# return $this->{Density} if($vol == 0.0);
my $Density = 0.0;
my @AtomTypeList = $this->GetCAtomTypeList();
my @AtomSiteList = $this->GetCAsymmetricAtomSiteList();
my $nAsymmetricAtomSite = $this->{'nAsymmetricAtomSite'};
#my $n = @AtomTypeList;
#print "n=$n, $nAsymmetricAtomSite\n";
my @nMultiplicityExpandedAtomSiteList
= $this->GetCnMultiplicityExpandedAtomSiteList();
for(my $i = 0 ; $i < $nAsymmetricAtomSite ; $i++) {
my $AtomSite = $AtomSiteList[$i];
my $AtomName = $AtomSite->AtomNameOnly(0);
# my $AtomName = $AtomSite->AtomNameOnly(1);
#print "an=$AtomName i=$i\n";
my $iAtomType = $this->FindiAtomType($AtomName);
if($iAtomType < 0) {
print "Error in CrystalObject::CalculateDensity: Invalid iAtomType=$iAtomType\n";
exit;
}
#print "iat=$iAtomType\n";
my $AtomType = $AtomTypeList[$iAtomType-1];
my $occupancy = $AtomSite->Occupancy();
#print "i=$i: $AtomSite, $AtomName, $iAtomType, $AtomType, $occupancy\n";
my $AtomicMass = $AtomType->AtomicMass();
my $mult = $nMultiplicityExpandedAtomSiteList[$i];
#print "$i: $AtomName ($iAtomType): $AtomicMass\n";
#print "$i: $AtomName ($AtomicMass) [$mult]\n";
$AtomicMass = 0.0 if(not defined $AtomicMass);
$occupancy = 1.0 if(not defined $occupancy);
$mult = 1.0 if(not defined $mult);
$Density += $AtomicMass*$mult*$occupancy;
}
my $NA = Sci::NA(); #6.0221367e23;
#print "vol: $vol NA: $NA\n";
$this->{"Density"} = ($vol == 0)? 0.0 : Sci::Round($Density / $NA / $vol * 1.0e24, 4);
return $this->{"Density"};
}
sub SetDensity
{
my ($this, $density) = @_;
return $this->{"Density"} = $density;
}
sub Density
{
my ($this) = @_;
return $this->{"Density"};
}
sub FindMass
{
my ($this, $atomname) = @_;
my @AtomTypeList = $this->GetCAtomTypeList();
for(my $i = 0 ; $i < @AtomTypeList ; $i++) {
my $atom = $AtomTypeList[$i];
my $name = $atom->AtomNameOnly();
next if(uc $name ne uc $atomname);
return $atom->AtomicMass();
}
return undef;
}
sub ConvertLattice
{
my ($this, $T, $LatticeParameterOnly, $CheckConsistency, $IsPrint) = @_;
my $nTotalExpandedAtomSite = $this->nTotalExpandedAtomSite();
my $SPG = $this->GetCSpaceGroup();
my $SPGName = $SPG->SPGName();
my ($a, $b, $c, $alpha, $beta, $gamma) = $this->LatticeParameters();
#実空間座標(x,y,z)の変換行列
my $tRT = new Math::MatrixReal(3, 3);
$tRT->transpose($T);
$tRT = $tRT->inverse();
my $g = new Math::MatrixReal(3, 3);
$g->assign(1, 1, $this->{'g11'});
$g->assign(1, 2, $this->{'g12'});
$g->assign(1, 3, $this->{'g13'});
$g->assign(2, 1, $this->{'g21'});
$g->assign(2, 2, $this->{'g22'});
$g->assign(2, 3, $this->{'g23'});
$g->assign(3, 1, $this->{'g31'});
$g->assign(3, 2, $this->{'g32'});
$g->assign(3, 3, $this->{'g33'});
my $g2 = new Math::MatrixReal(3, 3);
if(0) {
# if($SPGName =~ /^\s*[FIABC]/i) {
for(my $i = 1 ; $i <= 3 ; $i++) {
for(my $j = 1 ; $j <= 3 ; $j++) {
my $tij = $T->element($i,$j);
#print "$i,$j: $tij\n";
$g2->assign($i, $j, $tij);
}
}
}
else {
for(my $i1 = 1 ; $i1 <= 3 ; $i1++) {
for(my $i2 = 1 ; $i2 <= 3 ; $i2++) {
my $val = 0.0;
for(my $j1 = 1 ; $j1 <= 3 ; $j1++) {
for(my $j2 = 1 ; $j2 <= 3 ; $j2++) {
$val += $g->element($j1,$j2)
* $T->element($i1,$j1)
* $T->element($i2,$j2);
}
}
#print "val: $i1,$i2: $val\n";
$g2->assign($i1, $i2, $val);
}
}
}
#print "SPGName:$SPGName\n";
$this->{ConventionalSPGName} = $SPGName;
$this->{ConventionalLatticeParameters} = [$a, $b, $c, $alpha, $beta, $gamma];
print "***ConvL:$a, $b, $c, $alpha, $beta, $gamma\n";
my $NewCrystal = new Crystal();
$NewCrystal->{ConventionalSPGName} = $SPGName;
$NewCrystal->{ConventionalLatticeParameters} = [$a, $b, $c, $alpha, $beta, $gamma];
$NewCrystal->SetSpaceGroup("P 1", 1);
$NewCrystal->SetMetrics(
$g2->element(1,1), $g2->element(1,2), $g2->element(1,3),
$g2->element(2,1), $g2->element(2,2), $g2->element(2,3),
$g2->element(3,1), $g2->element(3,2), $g2->element(3,3));
$NewCrystal->CalLatticeParametersFromMetrics();
# Repeat to get proper lattice vectors for Primitive cell reduced by conventional cell
$NewCrystal->CalcMetrics();
#return $NewCrystal;
#my ($a2, $b2, $c2, $alpha2, $beta2, $gamma2) = $NewCrystal->LatticeParameters();
#die "**a'=$a2 b'=$b2 c'=$c2 alpha'=$alpha2 beta'=$beta2 gamma'=$gamma2\n";
return $NewCrystal if($LatticeParameterOnly);
$this->ExpandCoordinates();
my $amin = 0;
my $amax = 0;
for(my $i = 1 ; $i <= 3 ; $i++) {
my $e = $T->element($i,1);
$amin = int($e-0.999) if($e < $amin);
$amax = int($e-0.001) if($e > $amax);
}
my $bmin = 0;
my $bmax = 0;
for(my $i = 1 ; $i <= 3 ; $i++) {
my $e = $T->element($i,2);
$bmin = int($e-0.999) if($e < $bmin);
$bmax = int($e-0.001) if($e > $bmax);
}
my $cmin = 0;
my $cmax = 0;
for(my $i = 1 ; $i <= 3 ; $i++) {
my $e = $T->element($i,3);
$cmin = int($e-0.999) if($e < $cmin);
$cmax = int($e-0.001) if($e > $cmax);
}
if($CheckConsistency) {
$amax++; $bmax++; $cmax++;
}
#$bmin = 0;
print "Expansion range: ($amin,$bmin,$cmin) - ($amax,$bmax,$cmax)
\n" if($IsPrint);
my @ExpandedAtomSiteList = $this->GetCExpandedAtomSiteList();
$nTotalExpandedAtomSite = $this->nTotalExpandedAtomSite();
my $count = 0;
my %AtomTypeCount;
for(my $i = 0 ; $i < $nTotalExpandedAtomSite ; $i++) {
my $atom = $ExpandedAtomSiteList[$i];
my $atomnameraw = $atom->AtomName();
my $atomname = $atom->AtomNameOnly();
my $charge = $atom->Charge();
my ($x,$y,$z) = $atom->Position(1);
my $occ = $atom->Occupancy();
my $key;
if(abs($occ - 1.0) < 1.0e-4) {
$key = "$atomname$charge";
}
else {
$key = sprintf("$atomname${charge}[%6.4f]", $occ);
}
#print "i=$i atom=$atomname: key-name: $key : $atomnameraw\n";
for(my $ia = $amin ; $ia <= $amax ; $ia++) {
for(my $ib = $bmin ; $ib <= $bmax ; $ib++) {
for(my $ic = $cmin ; $ic <= $cmax ; $ic++) {
#print("ia,b,c=", $ia, $ib, $ic);
my $xyz = new Math::MatrixReal(3, 1);
$xyz->assign(1,1, $x+$ia);
$xyz->assign(2,1, $y+$ib);
$xyz->assign(3,1, $z+$ic);
$xyz = $tRT * $xyz;
my $x1 = Utils::Reduce01($xyz->element(1,1));
my $y1 = Utils::Reduce01($xyz->element(2,1));
my $z1 = Utils::Reduce01($xyz->element(3,1));
#print " (", $x+$ia, ", ", $y+$ib, ", ", $z+$ic, ") - ($x1,$y1,$z1)
\n";
my $WillAdd = 1;
my $MinimumDistance = 0.1;
if($NewCrystal->nAsymmetricAtomSite() == 0) {
$AtomTypeCount{$atomname}++;
my $label = $atomname . $AtomTypeCount{$atomname};
$NewCrystal->AddAtomSite($label, $key, $x1, $y1, $z1, $occ);
# $NewCrystal->AddAtomSite($label, $atomname, $x1, $y1, $z1, $occ);
#print " Passed0: $label: $key: ($x1,$y1,$z1) occ=$occ
\n";
next;
}
my @AtomList = $NewCrystal->GetCAsymmetricAtomSiteList();
my $nAtom = @AtomList;
#print(" *** nAtom=$nAtom\n");
for(my $j = 0 ; $j < @AtomList ; $j++) {
#print("j=$j\n");
my $atom0 = $AtomList[$j];
# my $atom0 = $NewCrystal->GetCAtomSite($j);
my $atomnameraw0 = $atom0->AtomName();
my $atomname0 = $atom0->AtomNameOnly();
my $charge0 = $atom0->Charge();
my ($x0,$y0,$z0) = $atom0->Position(1);
my $occ0 = $atom0->Occupancy();
#print("atom: $atomnameraw0 [Z=$charge0]\n");
#print " is=18 $atomnameraw0 $atomname0\n" if($i == 18);
my $key0;
if(abs($occ0 - 1.0) < 1.0e-4) {
$key0 = "$atomname0$charge0";
}
else {
$key0 = sprintf("$atomname0${charge0}[%6.4f]", $occ0);
}
#print("key: $key ne $key0?\n");
#print " keys(is=$i iatom=$j ($ia,$ib,$ic)): $key $key0\n" if($key =~ /Si/);
next if($key ne $key0);
# next if($atomname ne $atomname0);
my $dis = $NewCrystal->GetNearestInterAtomicDistance(
$x0, $y0, $z0, $x1, $y1, $z1);
#print " $nAtom: ($i,$j): ($ia,$ib,$ic): ($x0,$y0,$z0)-($x1,$y1,$z1) ";
#print " dis: $dis
minR=$MinimumDistance\n";
if($dis < $MinimumDistance) {
#print " Dropped: : $key: ($x1,$y1,$z1)-($x0,$y0,$z0) occ=$occ
\n";
$WillAdd = 0;
last;
}
}
if($WillAdd) {
$AtomTypeCount{$atomname}++;
my $label = $atomname . $AtomTypeCount{$atomname};
#print " Passed: $label: $key: ($x1,$y1,$z1) occ=$occ
\n";
$NewCrystal->AddAtomSite($label, $key, $x1, $y1, $z1, $occ);
# $NewCrystal->AddAtomSite($label, $atomname, $x1, $y1, $z1, $occ);
}
}
}
}
$count++;
}
$NewCrystal->ExpandCoordinates();
return $NewCrystal;
}
sub SetIsSelected
{
my($this, $index, $IsSelected) = @_;
my $RefAtomList = $this->{"RefExpandedAtomSiteList"};
my $atom = $RefAtomList->[$index];
if($atom ne '') {
return $atom->SetIsSelected($IsSelected);
}
return 0;
}
sub IsSelected
{
my($this, $index) = @_;
my $RefAtomList = $this->{"RefExpandedAtomSiteList"};
my $atom = $RefAtomList->[$index];
if($atom ne '') {
return $atom->IsSelected();
}
return 0;
}
sub SameIdSiteIndex
{
my ($this, $idm1, $sid) = @_;
my $SameIdSiteIndex = $this->{'SameIdSiteIndex'};
return $SameIdSiteIndex->{"[$idm1][$sid]"};
}
sub SortAtomTypeOrder
{
my ($this, $sorttype) = @_;
#print "SORT: $sorttype
";
if($sorttype =~ /ChargeAscend/i) {
my $ref = $this->{"RefAtomTypeList"};
@$ref = sort { $a->Charge() <=> $b->Charge() } @$ref;
}
elsif($sorttype =~ /ChargeDescend/i) {
my $ref = $this->{"RefAtomTypeList"};
@$ref = sort { -$a->Charge() <=> -$b->Charge() } @$ref;
}
return 1;
}
sub FillAtomTypeData
{
my ($this) = @_;
my $ref = $this->{"RefAtomTypeList"};
for(my $i = 0 ; $i < @$ref ; $i++) {
$ref->[$i]->FillAtomTypeData();
}
}
sub ReadASF
{
my ($this, $Source, $StopByError) = @_;
$StopByError = 1 if(!defined $StopByError);
my @AtomTypeList = $this->GetCAtomTypeList();
my $nAtomType = @AtomTypeList;
for(my $i = 0 ; $i < $nAtomType ; $i++) {
my $Rietan = new RietanFP;
my $AtomName;
my $AtomNameOnly = $AtomTypeList[$i]->AtomNameOnly(1);
my $Charge = AtomType::ChargeToVal($AtomTypeList[$i]->Charge());
#print "a1a: $AtomName (Charge=$Charge)\n";
if($Charge == 1) {
$AtomName = "$AtomNameOnly+";
}
elsif($Charge == -1) {
$AtomName = "$AtomNameOnly-";
}
elsif($Charge > 1) {
$AtomName = "$AtomNameOnly$Charge+";
}
elsif($Charge < -1) {
$Charge = -$Charge;
$AtomName = "$AtomNameOnly$Charge-";
}
else {
$AtomName = "$AtomNameOnly";
}
#print "a1: $AtomName (Charge=$Charge)\n";
#print "S=$StopByError\n";
my @a = $Rietan->ReadASFParameters($AtomName, undef, $Source, $StopByError);
if(@a < 11) {
$AtomName = $AtomTypeList[$i]->AtomNameOnly(1);
#print "a2: $AtomName\n";
@a = $Rietan->ReadASFParameters($AtomName, undef, $Source, $StopByError);
}
if(@a < 9) {
print "Error: Can not read ASFDC for [$AtomName].\n";
next;
}
#print "array: ", join(',', @a), "\n";
#print "i=$i: $AtomName\n";
$this->{"ASFDCRed$i"} = $Rietan;
$this->{"ASFDCRed$AtomName"} = $Rietan;
$this->{"ASFDCAtomName$i"} = $Rietan->{ASFDCAtomName};
$this->{"ASFDCAtomName$AtomName"} = $Rietan->{ASFDCAtomName};
my $AtomicNumber = $AtomTypeList[$i]->AtomicNumber();
$this->{"ASFDCAtomNameZ$i"} = $AtomicNumber;
$this->{"ASFDCAtomNameZ$AtomNameOnly"} = $AtomicNumber;
print "ASFDC AtomName[$i]: $Rietan->{ASFDCAtomName} [Zi=$AtomicNumber]\n";
#print " ao:$AtomNameOnly: $i\n";
}
$this->{ASFDCRed} = 1;
}
sub asfParameters
{
my ($this, $iAtomType, $Source,$StopByError) = @_;
$Source = 'Cu' if(!defined $Source);
$StopByError = 1 if(!defined $StopByError);
$this->ReadASF($Source,$StopByError) if(!defined $this->{ASFDCRed});
my $Rietan = $this->{"ASFDCRed$iAtomType"};
#print "R=$Rietan\n";
my $pP = $Rietan->{pASFDCParameters};
#print "pP=$pP\n";
# my ($A1,$B1,$A2,$B2,$A3,$B3,$A4,$B4,$C, $b, $df1,$df2) = @$pP;
return @$pP;
my ($A1, $B1, $A2, $B2, $A3, $B3, $A4, $B4, $A5, $B5, $C, $b,
$Crdf1, $Crdf2, $Fedf1, $Fedf2, $Codf1, $Codf2, $Cudf1, $Cudf2,
$Modf1, $Modf2, $Agdf1, $Agdf2, $Kbdf1, $Kbdf2) = @$pP;
#print "A=$A1 [$Source][$Crdf1][$Crdf2][b=$b]\n";
if($Source eq 'Cr') {# and defined $Crdf2) {
return ($A1, $B1, $A2, $B2, $A3, $B3, $A4, $B4, $A5, $B5, $C, $Crdf1, $Crdf2);
}
elsif($Source eq 'Fe') {# and defined $Fedf2) {
return ($A1, $B1, $A2, $B2, $A3, $B3, $A4, $B4, $A5, $B5, $C, $Fedf1, $Fedf2);
}
elsif($Source eq 'Co') {# and defined $Codf2) {
return ($A1, $B1, $A2, $B2, $A3, $B3, $A4, $B4, $A5, $B5, $C, $Codf1, $Codf2);
}
elsif($Source eq 'Cu') {# and defined $Cudf2) {
return ($A1, $B1, $A2, $B2, $A3, $B3, $A4, $B4, $A5, $B5, $C, $Cudf1, $Cudf2);
}
elsif($Source eq 'Mo') {# and defined $Modf2) {
return ($A1, $B1, $A2, $B2, $A3, $B3, $A4, $B4, $A5, $B5, $C, $Modf1, $Modf2);
}
elsif($Source eq 'Ag') {# and defined $Agdf2) {
return ($A1, $B1, $A2, $B2, $A3, $B3, $A4, $B4, $A5, $B5, $C, $Agdf1, $Agdf2);
}
elsif($Source eq 'Kb') {# and defined $Kbdf2) {
return ($A1, $B1, $A2, $B2, $A3, $B3, $A4, $B4, $A5, $B5, $C, $Kbdf1, $Kbdf2);
}
return ();
}
# $s = sin Q / lambda (in nm)
sub AtomicScatteringFactor
{
my ($this, $s) = @_;
return $this->asf($s);
}
# $s = sin Q / lambda (in nm)
sub asfElectron
{
my ($this, $iAtomType, $s,$StopByError) = @_;
$StopByError = 1 if(!defined $StopByError);
$s = 1.0e-5 if($s <= 0.0);
# the unit needs to be checked again
# my $s2 = 0.01 * $s * $s; # convert nm-1 to A-1
my $s2 = $s * $s;
my $k = 1.0e-2 * 8.0 * $pi * $pi / ($a0 * 1.0e9);
my $Z = $this->{"ASFDCAtomNameZ$iAtomType"};
my $asfXray = $this->asf($iAtomType, $s,$StopByError);
#print "s=$s, i=$iAtomType, Z=$Z, a=$asfXray\n";
#exit;
my $val = $k * ($Z - $asfXray) / $s2;
#print "va [$iAtomType]=$val = $k * ($Z - $asfXray) / $s2\n";
return $val;
}
# $s = sin Q / lambda (in nm)
sub asf
{
my ($this, $iAtomType, $s, $Source,$StopByError) = @_;
$Source = 'Cu' if(!defined $Source);
$StopByError = 1 if(!defined $StopByError);
$this->ReadASF($Source,$StopByError) if(!defined $this->{ASFDCRed});
my $Rietan = $this->{"ASFDCRed$iAtomType"};
my $pP = $Rietan->{pASFDCParameters};
my ($A1, $B1, $A2, $B2, $A3, $B3, $A4, $B4, $A5, $B5, $C, $b, $df1,$df2) = @$pP;
return 0.0 if(!defined $B5 or !defined $C);
$df1 = 0.0 if(!defined $df1);
$df2 = 0.0 if(!defined $df2);
my $s2 = 0.01 * $s * $s; # convert nm-1 to A-1
my $asf = $A1*exp(-$B1*$s2) + $A2*exp(-$B2*$s2) + $A3*exp(-$B3*$s2) + $A4*exp(-$B4*$s2) + $A5*exp(-$B5*$s2) + $C;
# my $asf = $A1*exp(-$B1*$s2) + $A2*exp(-$B2*$s2) + $A3*exp(-$B3*$s2) + $A4*exp(-$B4*$s2) + $C;
#print "i=$iAtomType, s=$s [$A1,$B1] [$A2,$B2] [$A3,$B3] [$A4,$B4] [$A5,$B5] [$C]\n";
#print "cal: [$A1*exp(-$B1*$s2) + $A2*exp(-$B2*$s2) + $A3*exp(-$B3*$s2) + $A4*exp(-$B4*$s2) + $A5*exp(-$B5*$s2) + $C]\n";
#print "i=$iAtomType, s=$s, asf=$asf [A1=$A1, C=$C]\n";
$asf = cplx($asf + $df1, $df2);
#print "i=$iAtomType, s=$s, asf=$asf [A1=$A1, C=$C] [df1=$df1, df2=$df2]\n";
return $asf;
my ($Crdf1,$Crdf2, $Fedf1,$Fedf2, $Codf1,$Codf2, $Cudf1,$Cudf2,
$Modf1, $Modf2, $Agdf1,$Agdf2, $Kbdf1, $Kbdf2);
($A1,$B1,$A2,$B2,$A3,$B3,$A4,$B4,$C, $b,
$Crdf1,$Crdf2, $Fedf1,$Fedf2, $Codf1,$Codf2, $Cudf1,$Cudf2,
$Modf1, $Modf2, $Agdf1,$Agdf2, $Kbdf1, $Kbdf2) = @$pP;
# if($Source eq 'n') {
# return $b;
# }
$s2 = 0.01 * $s * $s;
$asf = $A1*exp(-$B1*$s2)+$A2*exp(-$B2*$s2)+$A3*exp(-$B3*$s2)+$A4*exp(-$B4*$s2)+$C;
if($Source eq 'Cr' and defined $Crdf2) {
$asf = cplx($asf + $Crdf1, $Crdf2);
}
elsif($Source eq 'Fe' and defined $Fedf2) {
$asf = cplx($asf + $Fedf1, $Fedf2);
}
elsif($Source eq 'Co' and defined $Codf2) {
$asf = cplx($asf + $Codf1, $Codf2);
}
elsif($Source eq 'Cu' and defined $Cudf2) {
#print "d=$Cudf1,$Cudf2\n";
$asf = cplx($asf + $Cudf1, $Cudf2);
}
elsif($Source eq 'Mo' and defined $Modf2) {
$asf = cplx($asf + $Modf1, $Modf2);
}
elsif($Source eq 'Ag' and defined $Agdf2) {
$asf = cplx($asf + $Agdf1, $Agdf2);
}
elsif($Source eq 'Kb' and defined $Kbdf2) {
$asf = cplx($asf + $Kbdf1, $Kbdf2);
}
return $asf;
}
sub Fhkl
{
my ($this, $h, $k, $l, $IsElectron,$StopByError) = @_;
$IsElectron = 0 if(!defined $IsElectron);
$StopByError = 1 if(!defined $StopByError);
my @ExpandedAtomSiteList = $this->GetCExpandedAtomSiteList();
if(@ExpandedAtomSiteList == 0) {
$this->ExpandCoordinates();
@ExpandedAtomSiteList = $this->GetCExpandedAtomSiteList();
}
my $nTotalExpandedAtomSite = $this->nTotalExpandedAtomSite();
my $dhkl = $this->CalculateHKLInterplanarSpacing($h, $k, $l);
my $s = 1.0 / 2.0 / ($dhkl * 0.1); # in nm
#print "($h $k $l): d=$dhkl s=$s\n";
my $F = cplx(0, 0);
for(my $i = 0 ; $i < $nTotalExpandedAtomSite ; $i++) {
# my $AtomName = $ExpandedAtomSiteList[$i]->AtomName();
my $AtomName = $ExpandedAtomSiteList[$i]->AtomNameOnly(1);
my $iAtomType = $ExpandedAtomSiteList[$i]->iAtomType() - 1;
my ($x,$y,$z) = $ExpandedAtomSiteList[$i]->Position(1);
my $occupancy = $ExpandedAtomSiteList[$i]->Occupancy();
#print "$i: $iAtomType: $AtomName: ($x,$y,$z) occ=$occupancy [IsElectron=$IsElectron]\n";
my $asf;
if($IsElectron) {
$asf = $occupancy * $this->asfElectron($iAtomType, $s,$StopByError);
}
else {
$asf = $occupancy * $this->asf($iAtomType, $s,$StopByError);
}
#print "i=$i [$h,$k,$l]: pos($x,$y,$z) o=$occupancy asf=$asf ($AtomName:", $this->{"ASFDCAtomName$iAtomType"}, ")\n";
my $phi = 2.0 * $pi * ($h * $x + $k * $y + $l * $z);
$F += $asf * exp(i * $phi);
}
return (Re($F), Im($F));
# return ($Fr, $Fi);
}
sub OrthorhombicMultiplicity
{
my ($this, $m, $h, $k, $l) = @_;
$m *= 2 if($h != 0);
$m *= 2 if($k != 0);
$m *= 2 if($l != 0);
return $m;
}
sub TetragonalMultiplicity
{
my ($this, $m, $h, $k, $l) = @_;
$h = abs($h);
$k = abs($k);
$m *= 2 if($h != 0);
$m *= 2 if($k != 0);
$m *= 2 if($l != 0);
if($h != $k) {
$m *= 2;
}
return $m;
}
sub HexagonalMultiplicity
{
my ($this, $m, $h, $k, $l) = @_;
return 24 if($h != 0 and $k != 0 and $l != 0 and $h != $k);
return 12 if($h != 0 and $h == $k and $l != 0);
return 12 if($h == 0 and $k != 0 and $l != 0);
return 12 if($h != 0 and $k == 0 and $l != 0);
return 12 if($h != 0 and $k != 0 and $h != $k and $l == 0);
return 6 if($h != 0 and $h == $k and $l == 0);
return 6 if($h == 0 and $k != 0 and $l == 0);
return 6 if($h != 0 and $k == 0 and $l == 0);
return 2 if($h == 0 and $k == 0 and $l != 0);
print "Error in CrystalObject::HexagonalMultiplicity: Invalid hkl [$h,$k,$l]\n";
return undef;
}
sub TrigonalMultiplicity
{
my ($this, $m, $h, $k, $l) = @_;
return $this->HexagonalMultiplicity($m, $h, $k, $l);
}
sub RhombohedralMultiplicity
{
my ($this, $m, $h, $k, $l) = @_;
return $this->HexagonalMultiplicity($m, $h, $k, $l);
}
sub CubicMultiplicity
{
my ($this, $m, $h, $k, $l) = @_;
$h = abs($h);
$k = abs($k);
$l = abs($l);
$m *= 2 if($h != 0);
$m *= 2 if($k != 0);
$m *= 2 if($l != 0);
if($h != $k and $k != $l and $l != $h) {
$m *= 6;
}
elsif($h == $k and $k != $l) {
$m *= 3;
}
elsif($k == $l and $l != $h) {
$m *= 3;
}
elsif($l == $h and $h != $k) {
$m *= 3;
}
return $m;
}
sub Multiplicity
{
my ($this, $h, $k, $l) = @_;
my $ls = $this->LatticeSystem();
#print "ls: $ls\n";
my $m = 1;
if($ls eq 'cubic') {
return $this->CubicMultiplicity($m, $h, $k, $l);
}
elsif($ls eq 'tetragonal') {
return $this->TetragonalMultiplicity($m, $h, $k, $l);
}
elsif($ls eq 'trigonal') {
return $this->TrigonalMultiplicity($m, $h, $k, $l);
}
elsif($ls eq 'hexagonal') {
return $this->HexagonalMultiplicity($m, $h, $k, $l);
}
elsif($ls eq 'orthorhombic') {
return $this->OrthorhombicMultiplicity($m, $h, $k, $l);
}
elsif($ls eq 'triclinic') {
return 2;
}
else {
print "Error in CrystalObject::Multiplicity: Lattice system [$ls] is not supported.\n";
return undef;
}
}
sub ConvertByMatrix
{
my ($this, $x, $y, $z, $T) = @_;
my $x1 = $T->element(1,1)*$x + $T->element(1,2)*$y + $T->element(1,3)*$z;
my $y1 = $T->element(2,1)*$x + $T->element(2,2)*$y + $T->element(2,3)*$z;
my $z1 = $T->element(3,1)*$x + $T->element(3,2)*$y + $T->element(3,3)*$z;
return ($x1, $y1, $z1);
}
sub PrimitiveLatticeVectors
{
my ($this, $UseAtomicUnit) = @_;
my $k = ($UseAtomicUnit)? 1.0 / 0.529177 : 1.0;
my ($ax,$ay,$az,$bx,$by,$bz,$cx,$cy,$cz) = $this->LatticeVectors($UseAtomicUnit);
my ($sT, $T) = $this->GetBravaisToPrimitiveCellMatrix();
my ($ax1, $ay1, $az1) = $this->ConvertByMatrix($ax, $ay, $az, $sT);
my ($bx1, $by1, $bz1) = $this->ConvertByMatrix($bx, $by, $bz, $sT);
my ($cx1, $cy1, $cz1) = $this->ConvertByMatrix($cx, $cy, $cz, $sT);
#print "vector ($ax1, $ay1, $az1, $bx1, $by1, $bz1, $cx1, $cy1, $cz1)\n";
#exit;
return ($ax1, $ay1, $az1, $bx1, $by1, $bz1, $cx1, $cy1, $cz1);
}
sub GetCReducedAtomSiteList
{
my ($this, $pM, $pVT, $pAtomSiteList) = @_;
my @ExpandedAtomSiteList;
if(!defined $pAtomSiteList) {
@ExpandedAtomSiteList = $this->GetCExpandedAtomSiteListByOutputMode();
}
else {
@ExpandedAtomSiteList = @$pAtomSiteList;
}
if(@$pM == 1) {
return (1, @ExpandedAtomSiteList);
}
my $nTotalExpandedAtomSite = @ExpandedAtomSiteList;
my $IsMatched = 1;
my %iIdenticalSite;
for(my $i = 0 ; $i < $nTotalExpandedAtomSite ; $i++) {
next if(defined $iIdenticalSite{$i});
my $atom0 = $ExpandedAtomSiteList[$i];
my $atomname0 = $atom0->AtomNameOnly();
my $charge0 = $atom0->Charge();
my ($x0, $y0, $z0) = $atom0->Position(1);
# my $x0 = new MyMatrixReal([[$x0], [$y0], [$z0]]);
my $Vx0 = new MyMatrixReal([[$x0], [$y0], [$z0]]);
my $occupancy0 = $atom0->Occupancy();
#print "$i: $atomname0 ($x0, $y0, $z0) occ=$occupancy0\n";
my $IsFound = 0;
for(my $it = 1 ; $it < @$pM ; $it++) {
my $M = $pM->[$it];
my $V = $pVT->[$it];
my $x1 = $M * $Vx0 + $V;
# my $x1 = $M * $x0 + $V;
for(my $i1 = $i + 1 ; $i1 < $nTotalExpandedAtomSite ; $i1++) {
my $atom1 = $ExpandedAtomSiteList[$i1];
my $atomname1 = $atom1->AtomNameOnly();
my $charge1 = $atom1->Charge();
my ($x1, $y1, $z1) = $atom1->Position(1);
my $occupancy1 = $atom1->Occupancy();
#print " it=$it ($xt, $yt, $zt): $i1: $atomname1 ($x1, $y1, $z1): occ=$occupancy1\n";
next if($atomname0 ne $atomname1);
next if(abs($charge0 - $charge1 ) > 1.0e-3);
next if(abs($occupancy0 - $occupancy1) > 1.0e-3);
#print " candidate\n";
my $dis = $this->GetNearestInterAtomicDistance(
$x1->e(0), $x1->e(1), $x1->e(2), $x1, $y1, $z1);
if($dis < 1.0e-3) {
#print " identical: it=$it: $i1: $atomname1 ($x1, $y1, $z1)/($xt, $yt, $zt): $dis\n";
$iIdenticalSite{$i1} = $iIdenticalSite{$i1} . ", $i";
$IsFound = 1;
last;
}
else {
#print " not same : it=$it: $i1: $atomname1 ($x1, $y1, $z1)/($xt, $yt, $zt): $dis\n";
}
}
}
$IsMatched = 0 if(!$IsFound);
}
# foreach my $key (sort keys %iIdenticalSite) {
#print "Identical: $key - $iIdenticalSite{$key}\n";
# }
my @ReducedAtomSiteList;
if($IsMatched) {
for(my $i = 0 ; $i < $nTotalExpandedAtomSite ; $i++) {
next if($iIdenticalSite{$i} ne '');
push(@ReducedAtomSiteList, $ExpandedAtomSiteList[$i]);
}
}
else {
@ReducedAtomSiteList = @ExpandedAtomSiteList;
}
return ($IsMatched, @ReducedAtomSiteList);
}
sub GetCReducedAtomSiteListByTranslation
{
my ($this, $pTranslationalVectors, $pAtomSiteList) = @_;
if(!defined $pTranslationalVectors) {
$pTranslationalVectors = $this->GetTranslationalVectors();
}
if(0) {
my (@m, @v);
my $m = new MyMatrixReal('unit', 3, 3);
for(my $i = 0 ; $i < @$pTranslationalVectors ; $i++) {
my $p = $pTranslationalVectors->[$i];
my $v = new MyMatrixReal([[$p->[0]], [$p->[1]], [$p->[2]]]);
push(@m, $m);
push(@v, $v);
}
return $this->GetCReducedAtomSiteList(\@m, \@v, $pAtomSiteList);
}
my @TV = @$pTranslationalVectors;
my @ExpandedAtomSiteList;
if(!defined $pAtomSiteList) {
@ExpandedAtomSiteList = $this->GetCExpandedAtomSiteListByOutputMode();
}
else {
@ExpandedAtomSiteList = @$pAtomSiteList;
}
if(@TV == 1) {
return (1, @ExpandedAtomSiteList);
}
my $nTotalExpandedAtomSite = @ExpandedAtomSiteList;
my $IsMatched = 1;
my %iIdenticalSite;
for(my $i = 0 ; $i < $nTotalExpandedAtomSite ; $i++) {
next if(defined $iIdenticalSite{$i});
my $atom0 = $ExpandedAtomSiteList[$i];
my $atomname0 = $atom0->AtomNameOnly();
my $charge0 = $atom0->Charge();
my ($x0, $y0, $z0) = $atom0->Position(1);
my $occupancy0 = $atom0->Occupancy();
#print "$i: $atomname0 ($x0, $y0, $z0) occ=$occupancy0\n";
my $IsFound = 0;
for(my $it = 1 ; $it < @TV ; $it++) {
my $xt = $x0 + $TV[$it][0];
my $yt = $y0 + $TV[$it][1];
my $zt = $z0 + $TV[$it][2];
for(my $i1 = $i + 1 ; $i1 < $nTotalExpandedAtomSite ; $i1++) {
my $atom1 = $ExpandedAtomSiteList[$i1];
my $atomname1 = $atom1->AtomNameOnly();
my $charge1 = $atom1->Charge();
my ($x1, $y1, $z1) = $atom1->Position(1);
my $occupancy1 = $atom1->Occupancy();
#print " it=$it ($xt, $yt, $zt): $i1: $atomname1 ($x1, $y1, $z1): occ=$occupancy1\n";
next if($atomname0 ne $atomname1);
next if(abs($charge0 - $charge1 ) > 1.0e-3);
next if(abs($occupancy0 - $occupancy1) > 1.0e-3);
#print " candidate\n";
my $dis = $this->GetNearestInterAtomicDistance($xt, $yt, $zt, $x1, $y1, $z1);
if($dis < 1.0e-3) {
#print " identical: it=$it: $i1: $atomname1 ($x1, $y1, $z1)/($xt, $yt, $zt): $dis\n";
$iIdenticalSite{$i1} = $iIdenticalSite{$i1} . ", $i";
$IsFound = 1;
last;
}
else {
#print " not same : it=$it: $i1: $atomname1 ($x1, $y1, $z1)/($xt, $yt, $zt): $dis\n";
}
}
}
$IsMatched = 0 if(!$IsFound);
}
# foreach my $key (sort keys %iIdenticalSite) {
#print "Identical: $key - $iIdenticalSite{$key}\n";
# }
my @ReducedAtomSiteList;
if($IsMatched) {
for(my $i = 0 ; $i < $nTotalExpandedAtomSite ; $i++) {
next if($iIdenticalSite{$i} ne '');
push(@ReducedAtomSiteList, $ExpandedAtomSiteList[$i]);
}
}
else {
@ReducedAtomSiteList = @ExpandedAtomSiteList;
}
return ($IsMatched, @ReducedAtomSiteList);
}
sub GetCReducedAtomSiteListByInversion
{
my ($this, $pAtomSiteList) = @_;
if(0) {
my $m = new MyMatrixReal([[-1, 0, 0], [0, -1, 0], [0, 0, -1]]);
my $v = new MyMatrixReal([[0], [0], [0]]);
return $this->GetCReducedAtomSiteList([$m], [$v], $pAtomSiteList);
}
($this, $pAtomSiteList) = @_;
my @ExpandedAtomSiteList;
if(!defined $pAtomSiteList) {
@ExpandedAtomSiteList = $this->GetCExpandedAtomSiteListByOutputMode();
}
else {
@ExpandedAtomSiteList = @$pAtomSiteList;
}
my $nTotalExpandedAtomSite = @ExpandedAtomSiteList;
my $IsMatched = 1;
my %iIdenticalSite;
for(my $i = 0 ; $i < $nTotalExpandedAtomSite ; $i++) {
next if(defined $iIdenticalSite{$i});
my $atom0 = $ExpandedAtomSiteList[$i];
my $atomname0 = $atom0->AtomNameOnly();
my $charge0 = $atom0->Charge();
my ($x0, $y0, $z0) = $atom0->Position(1);
my $occupancy0 = $atom0->Occupancy();
#print "$i: $atomname0 ($x0, $y0, $z0) occ=$occupancy0\n";
my $IsFound = 0;
my ($xt, $yt, $zt) = (Utils::Reduce01(-$x0), Utils::Reduce01(-$y0), Utils::Reduce01(-$z0));
for(my $i1 = $i ; $i1 < $nTotalExpandedAtomSite ; $i1++) {
my $atom1 = $ExpandedAtomSiteList[$i1];
my $atomname1 = $atom1->AtomNameOnly();
my $charge1 = $atom1->Charge();
my ($x1, $y1, $z1) = $atom1->Position(1);
my $occupancy1 = $atom1->Occupancy();
#print " ($xt, $yt, $zt): $i1: $atomname1 ($x1, $y1, $z1): occ=$occupancy1\n";
next if($atomname0 ne $atomname1);
next if(abs($charge0 - $charge1 ) > 1.0e-3);
next if(abs($occupancy0 - $occupancy1) > 1.0e-3);
my $dis = $this->GetNearestInterAtomicDistance($xt, $yt, $zt, $x1, $y1, $z1);
if($dis < 1.0e-3) {
#print " identical: $i1: $atomname1 ($x1, $y1, $z1)/($xt, $yt, $zt): $dis\n";
$iIdenticalSite{$i1} = $iIdenticalSite{$i1} . ", $i";
$IsFound = 1;
last;
}
else {
#print " not same : $i1: $atomname1 ($x1, $y1, $z1)/($xt, $yt, $zt): $dis\n";
}
}
$IsMatched = 0 if(!$IsFound);
}
# foreach my $key (sort keys %iIdenticalSite) {
#print "Identical: $key - $iIdenticalSite{$key}\n";
# }
my @ReducedAtomSiteList;
if($IsMatched) {
for(my $i = 0 ; $i < $nTotalExpandedAtomSite ; $i++) {
next if($iIdenticalSite{$i} ne '');
push(@ReducedAtomSiteList, $ExpandedAtomSiteList[$i]);
}
}
else {
@ReducedAtomSiteList = @ExpandedAtomSiteList;
}
return ($IsMatched, @ReducedAtomSiteList);
}
sub AnalyzeaKProduct
{
my ($this, $aKProduct) = @_;
my ($a,$b,$c,$alpha,$beta,$gamma) = $this->LatticeParametersByOutputMode(0);
my ($nx, $ny, $nz);
my $an = $aKProduct * 10.0;
$nx = int($an / $a) + 1;
$ny = int($an / $b) + 1;
$nz = int($an / $c) + 1;
return ($nx, $ny, $nz);
}
sub GetTranslationalVectors
{
my ($this) = @_;
my $SPG = $this->GetCSpaceGroup();
my $SPGName = $SPG->SPGName();
my @VT = ([ 0, 0, 0]);
if($SPGName =~ /^\s*F/i) {
@VT = (
[ 0, 0, 0],
[ 0, 1/2, 1/2],
[1/2, 0, 1/2],
[1/2, 1/2, 0],
);
}
elsif($SPGName =~ /^\s*I/i) {
@VT = (
[ 0, 0, 0],
[1/2, 1/2, 1/2],
);
}
elsif($SPGName =~ /^\s*A/i) {
@VT = (
[ 0, 0, 0],
[ 0, 1/2, 1/2],
);
}
elsif($SPGName =~ /^\s*B/i) {
@VT = (
[ 0, 0, 0],
[1/2, 0, 1/2],
);
}
elsif($SPGName =~ /^\s*C/i) {
@VT = (
[ 0, 0, 0],
[1/2, 1/2, 0],
);
}
return \@VT;
}
sub IsFaceCenteredCell
{
my ($this, $pAS, $eps) = @_;
my @m = (
[0.0, 0.5, 0.5],
[0.5, 0.0, 0.5],
[0.5, 0.5, 0.0],
);
return $this->IsBravaisCellByMatrix($pAS, $eps, \@m);
}
sub IsBodyCenteredCell
{
my ($this, $pAS, $eps) = @_;
my @m = (
[0.5, 0.5, 0.5]
);
return $this->IsBravaisCellByMatrix($pAS, $eps, \@m);
}
sub IsACenteredCell
{
my ($this, $pAS, $eps) = @_;
my @m = (
[0.0, 0.5, 0.5]
);
return $this->IsBravaisCellByMatrix($pAS, $eps, \@m);
}
sub IsBCenteredCell
{
my ($this, $pAS, $eps) = @_;
my @m = (
[0.5, 0.0, 0.5]
);
return $this->IsBravaisCellByMatrix($pAS, $eps, \@m);
}
sub IsCCenteredCell
{
my ($this, $pAS, $eps) = @_;
my @m = (
[0.5, 0.5, 0.0]
);
return $this->IsBravaisCellByMatrix($pAS, $eps, \@m);
}
sub GuessBravaisLattice
{
my ($this, $tollatt, $tolangle, $EPSCoord) = @_;
my ($a, $b, $c, $alpha, $beta, $gamma) = $this->LatticeParameters();
my $SPG = $this->GetCSpaceGroup();
my $SPGName = $SPG->SPGName();
my $ls = lc $this->LatticeSystem();
#print "SPG=$SPGName LS=$ls\n";
my $LatticeSystem = ($ls eq '')? $SPG->SetLatticeSystem(undef, 1, $tollatt, $tolangle) : $ls;
my $BravaisLattice;
if($SPGName =~ /^\s*F/i) {
$BravaisLattice = 'F';
}
elsif($SPGName =~ /^\s*I/i) {
$BravaisLattice = 'I';
}
elsif($SPGName =~ /^\s*A/i) {
$BravaisLattice = 'A';
}
elsif($SPGName =~ /^\s*B/i) {
$BravaisLattice = 'B';
}
elsif($SPGName =~ /^\s*C/i) {
$BravaisLattice = 'C';
}
#print "SPG=$SPGName BL=$BravaisLattice\n";
if($BravaisLattice eq '') {
my @as = $this->GetCExpandedAtomSiteListByOutputMode();
if($this->IsFaceCenteredCell(\@as, $EPSCoord)) {
$BravaisLattice = 'F';
}
elsif($this->IsBodyCenteredCell(\@as, $EPSCoord)) {
$BravaisLattice = 'I';
}
elsif($this->IsACenteredCell(\@as, $EPSCoord)) {
$BravaisLattice = 'A';
}
elsif($this->IsBCenteredCell(\@as, $EPSCoord)) {
$BravaisLattice = 'B';
}
elsif($this->IsCCenteredCell(\@as, $EPSCoord)) {
$BravaisLattice = 'C';
}
else {
$BravaisLattice = 'P';
}
}
return ($LatticeSystem, $BravaisLattice);
}
sub GetLatticeConversionMatrix
{
my ($this, $PresetRule, $DirectionSelect, $ParametersSelect, $IsPrint) = @_;
$IsPrint = 1 if(!defined $IsPrint);
my $sT = new Math::MatrixReal(3, 3);
my $T = new Math::MatrixReal(3, 3);
my $tRT = new Math::MatrixReal(3, 3);
my $RT = new Math::MatrixReal(3, 3);
my $tR = new Math::MatrixReal(3, 3);
if($PresetRule eq 'No' or $PresetRule eq 'P1') {
$sT = Math::MatrixReal->new_from_cols(
[ [ 1, 0, 0],
[ 0, 1, 0],
[ 0, 0, 1] ] );
$sT->transpose($sT);
}
elsif($PresetRule eq 'RhombHex') {
$sT = Math::MatrixReal->new_from_cols(
[ [ 1, -1, 0],
[ 0, 1, -1],
[ 1, 1, 1] ] );
$sT->transpose($sT);
}
elsif($PresetRule eq 'HexRhomb') {
$sT = Math::MatrixReal->new_from_cols(
[ [ 2/3, 1/3, 1/3],
[-1/3, 1/3, 1/3],
[-1/3,-2/3, 1/3] ]);
$sT->transpose($sT);
}
elsif($PresetRule eq 'HexOrtho') {
$sT = Math::MatrixReal->new_from_cols(
[ [ 1, 0, 0],
[ 1, 2, 0],
[ 0, 0, 1] ] );
$sT->transpose($sT);
}
elsif($PresetRule eq 'FCCPrim' or $PresetRule eq 'PrimFCC') {
$sT = Math::MatrixReal->new_from_cols(
[ [ 0.5, 0.5, 0],
[ 0, 0.5, 0.5],
[ 0.5, 0, 0.5] ] );
if($PresetRule eq 'PrimFCC') {
$sT = $sT->inverse();
}
$sT->transpose($sT);
}
elsif($PresetRule eq 'BCCPrim' or $PresetRule eq 'PrimBCC') {
$sT = Math::MatrixReal->new_from_cols(
[ [-0.5, 0.5, 0.5],
[ 0.5,-0.5, 0.5],
[ 0.5, 0.5,-0.5] ] );
if($PresetRule eq 'PrimBCC') {
$sT = $sT->inverse();
}
$sT->transpose($sT);
}
elsif($PresetRule eq 'ACenterPrim') {
$sT = Math::MatrixReal->new_from_cols(
[ [ 1.0, 0.0, 0.0],
[ 0.0, 0.5, 0.5],
[ 0.0,-0.5, 0.5] ] );
$sT->transpose($sT);
}
elsif($PresetRule eq 'BCenterPrim') {
$sT = Math::MatrixReal->new_from_cols(
[ [ 0.5, 0.0, 0.5],
[ 0.0, 1.0, 0.0],
[-0.5, 0.0, 0.5] ] );
$sT->transpose($sT);
}
elsif($PresetRule eq 'CCenterPrim') {
$sT = Math::MatrixReal->new_from_cols(
[ [ 0.5, 0.5, 0.0],
[-0.5, 0.5, 0.0],
[ 0.0, 0.0, 1.0] ] );
$sT->transpose($sT);
}
else {
die "Error: Invalid conversion rule [$PresetRule].\n";
print " Rule is specified by Matrix:\n";
$sT->assign(1, 1, eval(&GetArg('t11')));
$sT->assign(1, 2, eval(&GetArg('t12')));
$sT->assign(1, 3, eval(&GetArg('t13')));
$sT->assign(2, 1, eval(&GetArg('t21')));
$sT->assign(2, 2, eval(&GetArg('t22')));
$sT->assign(2, 3, eval(&GetArg('t23')));
$sT->assign(3, 1, eval(&GetArg('t31')));
$sT->assign(3, 2, eval(&GetArg('t32')));
$sT->assign(3, 3, eval(&GetArg('t33')));
}
print "Conversion rule:\n" if($IsPrint);
my $IsChosen = 0;
if($ParametersSelect eq 'ai') {
if($DirectionSelect eq 'OriginalToConverted') {
print "Real space vectors: (a'i) = (tij)(ai)\n" if($IsPrint);
$T = $sT;
$IsChosen = 1;
}
elsif($DirectionSelect eq 'ConvertedlToOriginal') {
print "Real space vectors: (ai) = (tij)(a'i)\n" if($IsPrint);
$T = $sT->inverse();
$IsChosen = 1;
}
}
elsif($ParametersSelect eq 'xyz') {
if($DirectionSelect eq 'OriginalToConverted') {
print "Real space coordinates: (x'i) = (tij)(xi)\n" if($IsPrint);
$T->transpose($sT);
$T = $T->inverse();
$IsChosen = 1;
}
elsif($DirectionSelect eq 'ConvertedlToOriginal') {
print "Real space vectors: (xi) = (tij)(x'i)\n" if($IsPrint);
$T->transpose($sT);
$IsChosen = 1;
}
}
elsif($ParametersSelect eq 'a*i') {
if($DirectionSelect eq 'OriginalToConverted') {
print "Reciprocal space vectors: (a'*i) = (tij)(a*i)\n" if($IsPrint);
$T->transpose($sT);
$T = $T->inverse();
$IsChosen = 1;
}
elsif($DirectionSelect eq 'ConvertedlToOriginal') {
print "Reciprocal space vectors: (a*i) = (tij)(a'*i)\n" if($IsPrint);
$T->transpose($sT);
$IsChosen = 1;
}
}
elsif($ParametersSelect eq 'hkl') {
if($DirectionSelect eq 'OriginalToConverted') {
print "Reciprocal space coordinates: (h'i) = (tij)(hi)\n" if($IsPrint);
$T = $sT;
$IsChosen = 1;
}
elsif($DirectionSelect eq 'ConvertedlToOriginal') {
print "Reciprocal space coordinates: (hi) = (tij)(h'i)\n" if($IsPrint);
$T = $sT->inverse();
$IsChosen = 1;
}
}
unless($IsChosen) {
print "Error in CrystalObject::GetLatticeConversionMatrix: Selected option {Parameters: $ParametersSelect} in {$DirectionSelect} direction is not supported.\n";
return 1;
}
#実空間座標(x,y,z)の変換行列
$tRT->transpose($T);
$tRT = $tRT->inverse();
#逆空間ベクトル(a*i)の変換行列: $tRT
$RT->transpose($T);
$RT = $RT->inverse();
#逆空間座標(h k l)の変換行列: $T
$tR = $T;
return ($sT, $T, $tRT, $RT, $tR);
}
sub GetPrimitiveCrystal
{
my ($this, $LatticeParameterOnly, $CheckConsistency, $IsPrint) = @_;
my $T = new Math::MatrixReal(3, 3);
my $matrix = $this->GetMatrixConventionalToPrimitiveCell(0);
$T = Math::MatrixReal->new_from_cols($matrix);
$T->transpose($T);
return $this->ConvertLattice($T, $LatticeParameterOnly, $CheckConsistency, $IsPrint);
}
sub GetMatrixConventionalToPrimitiveCell
{
my ($this, $ReferToConventionalCell) = @_;
my $SPGName = ($ReferToConventionalCell and $this->{ConventionalSPGName})?
$this->{ConventionalSPGName} : $this->SPGName();
#print "2********[$this->{ConventionalSPGName}][$SPGName]\n";
if($SPGName =~ /^\s*F/i) {
return [ [ 0.5, 0.5, 0],
[ 0, 0.5, 0.5],
[ 0.5, 0, 0.5] ]
}
elsif($SPGName =~ /^\s*I/i) {
return [ [-0.5, 0.5, 0.5],
[ 0.5,-0.5, 0.5],
[ 0.5, 0.5,-0.5] ];
}
elsif($SPGName =~ /^\s*A/i) {
return [ [ 1.0, 0.0, 0.0],
[ 0.0, 0.5, 0.5],
[ 0.0,-0.5, 0.5] ];
}
elsif($SPGName =~ /^\s*B/i) {
return [ [ 0.5, 0.0, 0.5],
[ 0.0, 1.0, 0.0],
[-0.5, 0.0, 0.5] ];
}
elsif($SPGName =~ /^\s*C/i) {
return [ [ 0.5, 0.5, 0.0],
[-0.5, 0.5, 0.0],
[ 0.0, 0.0, 1.0] ];
}
return [ [ 1.0, 0.0, 0.0],
[ 0.0, 1.0, 0.0],
[ 0.0, 0.0, 1.0] ];
}
sub GetMatrixConventionalToPrimitiveReciprocalCell
{
my ($this) = @_;
my $SPGName = $this->SPGName();
if($SPGName =~ /^\s*F/i) {
return [ [-1.0, 1.0, 1.0],
[ 1.0,-1.0, 1.0],
[ 1.0, 1.0,-1.0] ];
}
elsif($SPGName =~ /^\s*I/i) {
return [ [ 1.0, 1.0, 0],
[ 0, 1.0, 1.0],
[ 1.0, 0, 1.0] ]
}
elsif($SPGName =~ /^\s*A/i) {
return [ [ 1.0, 0.0, 0.0],
[ 0.0, 1.0, 1.0],
[ 0.0,-1.0, 1.0] ];
}
elsif($SPGName =~ /^\s*B/i) {
return [ [ 1.0, 0.0, 1.0],
[ 0.0, 1.0, 0.0],
[-1.0, 0.0, 1.0] ];
}
elsif($SPGName =~ /^\s*C/i) {
return [ [ 1.0, 1.0, 0.0],
[-1.0, 1.0, 0.0],
[ 0.0, 0.0, 1.0] ];
}
return [ [ 1.0, 0.0, 0.0],
[ 0.0, 1.0, 0.0],
[ 0.0, 0.0, 1.0] ];
}
sub IsBravaisCellByMatrix
{
my ($this, $pAS, $eps, $pm) = @_;
$eps = 0.01 if(!defined $eps);
my @Checked;
for(my $i = 0 ; $i < @$pAS ; $i++) {
next if($Checked[$i]);
my $site1 = $pAS->[$i];
my $atom1 = $site1->AtomNameOnly();
my ($x1, $y1, $z1) = $site1->Position(1);
my $occ1 = $site1->Occupancy();
#Utils::Reduce01($x)
for(my $im = $0 ; $im < @$pm ; $im++) {
my $xm = $x1 + $pm->[$im][0];
my $ym = $y1 + $pm->[$im][1];
my $zm = $z1 + $pm->[$im][2];
my $IsFound = 0;
for(my $j = $i+1 ; $j < @$pAS ; $j++) {
my $site2 = $pAS->[$j];
my $atom2 = $site2->AtomNameOnly();
my ($x2, $y2, $z2) = $site2->Position(1);
my $occ2 = $site2->Occupancy();
next if($atom1 ne $atom2);
my $dis = $this->GetNearestInterAtomicDistance(
$xm, $ym, $zm, $x2, $y2, $z2);
if($dis < $eps) {
$IsFound = 1;
$Checked[$j] = 1;
last;
}
}
return 0 if(!$IsFound);
}
}
return 1;
}
sub GetBravaisToPrimitiveCellMatrix
{
my ($this) = @_;
my $SPG = $this->GetCSpaceGroup();
my $SPGName = $SPG->SPGName();
my $sT = new Math::MatrixReal(3, 3);
my $T;
if($SPGName =~ /^\s*F/i) {
$sT = Math::MatrixReal->new_from_cols(
[ [ 0.5, 0.5, 0],
[ 0, 0.5, 0.5],
[ 0.5, 0, 0.5] ] );
}
elsif($SPGName =~ /^\s*I/i) {
$sT = Math::MatrixReal->new_from_cols(
[ [-0.5, 0.5, 0.5],
[ 0.5,-0.5, 0.5],
[ 0.5, 0.5,-0.5] ] );
}
elsif($SPGName =~ /^\s*A/i) {
$sT = Math::MatrixReal->new_from_cols(
[ [ 1.0, 0.0, 0.0],
[ 0.0, 0.5, 0.5],
[ 0.0,-0.5, 0.5] ] );
}
elsif($SPGName =~ /^\s*B/i) {
$sT = Math::MatrixReal->new_from_cols(
[ [ 0.5, 0.0, 0.5],
[ 0.0, 1.0, 0.0],
[-0.5, 0.0, 0.5] ] );
}
elsif($SPGName =~ /^\s*C/i) {
$sT = Math::MatrixReal->new_from_cols(
[ [ 0.5, 0.5, 0.0],
[-0.5, 0.5, 0.0],
[ 0.0, 0.0, 1.0] ] );
}
$T = $sT;
$T->transpose($T);
return ($sT, $T);
}
sub GetPrimitiveCell
{
my ($this, $IsPrint) = @_;
#実空間べクトル(ai)の変換行列
my ($sT, $T) = $this->GetBravaisToPrimitiveCellMatrix();
if($IsPrint) {
print "Conversion matrix for real space vector: (a'i) = (Tij)(ai)\n";
printf " |%10.4f %10.4f %10.4f|\n", $T->element(1,1), $T->element(1,2), $T->element(1,3);
printf " |%10.4f %10.4f %10.4f|\n", $T->element(2,1), $T->element(2,2), $T->element(2,3);
printf " |%10.4f %10.4f %10.4f|\n", $T->element(3,1), $T->element(3,2), $T->element(3,3);
}
return $this->ConvertLattice($T, 0, 1);
}
sub FindTranslations
{
my ($this, $OverlapCriteriaRadius, $IsPrint) = @_;
my $SPG = $this->GetCSpaceGroup();
my @TranslationArrays = $SPG->GetTranslationArrays();
my (@hit, $pID, @ID);
for(my $i = 0 ; $i < @TranslationArrays ; $i++) {
($hit[$i], $pID) = $this->CheckATransration(
$TranslationArrays[$i], $OverlapCriteriaRadius, $IsPrint);
if($hit[$i]) {
# print "** Translation [TranslationArrays[$i]{Name}] 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 < @TranslationArrays ; $i++) {
if($hit[$i]) {
print "** Translation [$TranslationArrays[$i]{Name}] is found\n" if($IsPrint);
return ($TranslationArrays[$i]{Name}, \@ID);
}
}
return ('', \@ID);
}
sub CheckATransration
{
my ($this, $ppT, $OverlapCriteriaRadius, $IsPrint) = @_;
#$IsPrint=1;
my $pT = $ppT->{Matrix};
my @SiteList = $this->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 = $this->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 FindSiteSymmetries
{
my ($this, $xc, $yc, $zc, $OverlapCriteriaRadius, $IsPrint) = @_;
my @MArrays = PointGroup->new()->GetSymmetryOperationArraysInCartesian();
#print " Check site symmetry for $label [$type] at ($xc, $yc, $zc)...\n";
my ($rxc, $ryc, $rzc) = $this->CalCartesianCoordinate($xc, $yc, $zc);
my @AtomSiteList = $this->GetCAsymmetricAtomSiteList();
my $nAsymmetricAtomSite = @AtomSiteList;
my @ExpandedAtomSiteList = $this->GetCExpandedAtomSiteList();
my $nTotalExpandedAtomSite = $this->nTotalExpandedAtomSite();
# my @nMultiplicityExpandedAtomSiteList = $Crystal->GetCnMultiplicityExpandedAtomSiteList();
my (@hit);
for(my $is = 0 ; $is < @MArrays ; $is++) {
my $pM = $MArrays[$is]->{Matrix};
my $Name = $MArrays[$is]->{Name};
my $pF = $MArrays[$is]->{Freedom};
$hit[$is] = 1;
for(my $i = 0 ; $i < $nAsymmetricAtomSite ; $i++) {
my $p = $AtomSiteList[$i];
my $label = $p->Label();
my $type = $p->AtomName();
my ($x,$y,$z) = $p->Position(1);
my $occupancy = $p->Occupancy();
my $x0 = $x - $xc;
my $y0 = $y - $yc;
my $z0 = $z - $zc;
my ($rx0, $ry0, $rz0) = $this->CalCartesianCoordinate($x0, $y0, $z0);
# my $xs = $pM->[0][0] * $x0 + $pM->[0][1] * $y0 + $pM->[0][2] * $z0 + $pM->[0][3] + $xc;
# my $ys = $pM->[1][0] * $x0 + $pM->[1][1] * $y0 + $pM->[1][2] * $z0 + $pM->[1][3] + $yc;
# my $zs = $pM->[2][0] * $x0 + $pM->[2][1] * $y0 + $pM->[2][2] * $z0 + $pM->[2][3] + $zc;
my $rxs = $pM->[0][0] * $rx0 + $pM->[0][1] * $ry0 + $pM->[0][2] * $rz0 + $pM->[0][3] + $rxc;
my $rys = $pM->[1][0] * $rx0 + $pM->[1][1] * $ry0 + $pM->[1][2] * $rz0 + $pM->[1][3] + $ryc;
my $rzs = $pM->[2][0] * $rx0 + $pM->[2][1] * $ry0 + $pM->[2][2] * $rz0 + $pM->[2][3] + $rzc;
my ($xs, $ys, $zs) = $this->CalFractionalCoordinate($rxs, $rys, $rzs);
#if($Name =~ /3z/) {
#print "(x,y,z)=($x,$y,$z) 0=($x0,$y0,$z0)\n";
#print "zs = $xs = $pM->[0][0] * $x0 + $pM->[0][1] * $y0 + $pM->[0][2] * $z0 + $pM->[0][3] + $xc\n";
#print "ys = $ys = $pM->[1][0] * $x0 + $pM->[1][1] * $y0 + $pM->[1][2] * $z0 + $pM->[1][3] + $yc\n";
#print "zs = $zs = $pM->[2][0] * $x0 + $pM->[2][1] * $y0 + $pM->[2][2] * $z0 + $pM->[2][3] + $zc\n";
#}
my $id2 = $this->FindIdenticalSite($type, $xs, $ys, $zs, $OverlapCriteriaRadius);#, $i);
#$IsPrint=1;
if(defined $id2) {
printf " Sym Op[$is] ($Name): ($xs,$ys,$zs): Atom ($x,$y,$z) is equivalent to site $id2\n",
$xs, $ys, $zs, $x, $y, $z if($IsPrint);
next;
}
else {
printf " Sym Op[$is] ($Name): (%7.4f,%7.4f,%7.4f): Equivalent site of (%7.4f,%7.4f,%7.4f) was not found\n",
$xs, $ys, $zs, $x, $y, $z if($IsPrint);
$hit[$is] = 0;
last;
}
}
}
my ($PG) = PointGroup->new()->FindPointGroup(\@hit, \@MArrays);
print " Point Group: $PG\n";
print " Matched symmetries: ";
for(my $i = 0 ; $i < @MArrays ; $i++) {
next if(!$hit[$i]);
my $Name = $MArrays[$i]->{Name};
print "$Name ";
}
print "\n";
if(0) {
print " Incompatible symmetries: ";
for(my $i = 0 ; $i < @MArrays ; $i++) {
next if($hit[$i]);
my $Name = $MArrays[$i]->{Name};
print "$Name ";
}
print "\n";
}
return (1, \@hit);
}
sub CheckSPG
{
my ($this, $TSym, $pID, $iSPG, $iSet, $OverlapCriteriaRadius, $IsPrint) = @_;
my ($a,$b,$c,$alpha,$beta,$gamma) = $this->LatticeParameters();
# $Crystal->ExpandCoordinates();
my $SPG = $this->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 = $this->GetCExpandedAtomSiteList();
my $nTotalExpandedAtomSite = $this->nTotalExpandedAtomSite();
my @nMultiplicityExpandedAtomSiteList = $this->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(!$SPG2->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
\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 = $this->FindIdenticalSite($atomname, $xs, $ys, $zs, $OverlapCriteriaRadius);
my $id2 = $this->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 FindOrigin
{
my ($this, $xc, $yc, $zc, $OverlapCriteriaRadius,
$ShowAllCoordination, $RCoordMax, $IsPrint) = @_;
my $Rmax = $RCoordMax;
my $kRth = 0.1;
print "\n** Find origin from initial center ($xc, $yc, $zc)\n";
my @MArrays = PointGroup->new()->GetSymmetryOperationArraysInCartesian();
my ($la, $lb, $lc, $alpha, $beta, $gamma) = $this->LatticeParameters();
my @AtomSiteList = $this->GetCExpandedAtomSiteList();
my $nTotalSite = $this->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 = $this->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;
my @F = (1, 1, 1);
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;
if($n1stCoord > 1) {
$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;
print " Find site symmetry at ($xc1, $yc1, $zc1)\n";
my ($ret, $pMatchedSyms) = $this->FindSiteSymmetries($xc1, $yc1, $zc1, $OverlapCriteriaRadius);#, $IsPrint);
for(my $is = 0 ; $is < @MArrays ; $is++) {
next if(!$pMatchedSyms->[$is]);
# my $pM = $MArrays[$is]->{Matrix};
# my $Name = $MArrays[$is]->{Name};
my $pF = $MArrays[$is]->{Freedom};
$F[0] = 0 if(!$pF->[0]);
$F[1] = 0 if(!$pF->[1]);
$F[2] = 0 if(!$pF->[2]);
}
print " ** Total freedom: ($F[0], $F[1], $F[2])\n";
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, \@F);
}
sub FindPossibleSPGs
{
my ($this, $ShowCrystalInf, $ShowAllCoordination, $ShowFoundStructure, $ScanAll,
$FindSiteSymmetry, $xc0, $yc0, $zc0,
$OverlapCriteriaRadius, $ShowAllCoordination2, $RCoordMax, $IsPrint) = @_;
# $OverlapCriteriaRadius, $ShowAllCoordination, $RCoordMax, $IsPrint) = @_;
$this->ExpandCoordinates();
my $SPG = $this->GetCSpaceGroup();
my $SPGName = $SPG->SPGName();
my $iSPG = $SPG->iSPG();
my $LatticeSystem = $SPG->LatticeSystem();
my $nTranslation = $SPG->nTranslation();
my $nSymmetryOperation = $SPG->nSymmetryOperation();
my @AtomTypeList = $this->GetCAtomTypeList();
my $nAtomType = @AtomTypeList;
my @AtomSiteList = $this->GetCAsymmetricAtomSiteList();
my $nAsymmetricAtomSite = @AtomSiteList;
my @ExpandedAtomSiteList = $this->GetCExpandedAtomSiteList();
my $nTotalExpandedAtomSite = $this->nTotalExpandedAtomSite();
my ($TSym, $pID) = $this->FindTranslations($OverlapCriteriaRadius);
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 = $this->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, $pFreedom) = $this->FindOrigin($xc0, $yc0, $zc0, $OverlapCriteriaRadius,
$ShowAllCoordination, $RCoordMax);
$xc = 0 if(!$pFreedom->[0]);
$yc = 0 if(!$pFreedom->[1]);
$zc = 0 if(!$pFreedom->[2]);
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);
}
$this->SetCAsymmetricAtomSiteList(@ExpandedAtomSiteList);
# @ExpandedAtomSiteList = $this->GetCExpandedAtomSiteList();
my ($MostProbableSPG, $MostProbableiSPG, $MostProbableiSet, $MaxnSym) = (undef, -1, -1, 0);
my @PossibleSPGs;
my $c = 0;
for(my $iSPG = 230 ; $iSPG >= 1 ; $iSPG--) {
for(my $iSet = 1 ; ; $iSet++) {
my ($SPG2, $IsMatched, $mess, $pSiteId, $pAsymSites)
= $this->CheckSPG($TSym, $pID, $iSPG, $iSet, $OverlapCriteriaRadius, 0);
last if(!defined $SPG2);
if($IsMatched) {
my $nASSite = @$pAsymSites;
my $SPGName2 = $SPG2->SPGName();
my $LatticeSystem2 = $SPG2->LatticeSystem();
my $nTranslation = $SPG2->nTranslation();
my $nSymmetryOperation = $SPG2->nSymmetryOperation();
my $nTotalSymmetry = $nTranslation * $nSymmetryOperation;
if($MaxnSym < $nTotalSymmetry) {
$MaxnSym = $nTotalSymmetry;
$MostProbableSPG = $SPG2;
$MostProbableiSPG = $iSPG;
$MostProbableiSet = $iSet;
}
$PossibleSPGs[$c] = {
pSPG => $SPG2,
iSPG => $iSPG,
iSet => $iSet,
SPGName => $SPGName2,
LatticeSystem => $LatticeSystem2,
nTranslation => $nTranslation,
nSymmetryOperation => $nSymmetryOperation,
nTotalSymmetry => $nTotalSymmetry,
pAsymSites => $pAsymSites,
nASSite => $nASSite,
};
$c++;
}
else {
# print "SPG $SPGName2 (#$iSPG-$iSet LS=$LatticeSystem2) is not matched [$mess]\n";
}
}
}
my $SPGName2 = $MostProbableSPG->SPGName();
my $LatticeSystem2 = $MostProbableSPG->LatticeSystem();
$nTranslation = $MostProbableSPG->nTranslation();
$nSymmetryOperation = $MostProbableSPG->nSymmetryOperation();
my $nTotalSymmetry = $nTranslation * $nSymmetryOperation;
my %MostPossibleSPG = (
iSPG => $MostProbableiSPG,
iSet => $MostProbableiSet,
SPGName => $SPGName2,
LatticeSystem => $LatticeSystem2,
nTranslation => $nTranslation,
nSymmetryOperation => $nSymmetryOperation,
nTotalSymmetry => $nTotalSymmetry
);
return (\%MostPossibleSPG, \@PossibleSPGs);
}
1;