sub LatticeConversion
{
my ($Action) = @_;
print "Lattice conversion
\n";
my @ParametersSelect = &GetArg('Parameters');
my $ParametersSelect = $ParametersSelect[0];
my @DirectionSelect = &GetArg('Direction');
my $DirectionSelect = $DirectionSelect[0];
my $CheckConsistency = &GetArg('CheckConsistency');
my $a = &GetArg('a'); $a = eval($a);
my $b = &GetArg('b'); $b = eval($b);
my $c = &GetArg('c'); $c = eval($c);
my $alpha = &GetArg('alpha'); $alpha = eval($alpha);
my $beta = &GetArg('beta'); $beta = eval($beta);
my $gamma = &GetArg('gamma'); $gamma = eval($gamma);
my $sxyz = new Math::MatrixReal(3, 1);
my $v = &GetArg('sx'); $v = eval($v);
$sxyz->assign(1, 1, $v);
$v = &GetArg('sy'); $v = eval($v);
$sxyz->assign(2, 1, $v);
$v = &GetArg('sz'); $v = eval($v);
$sxyz->assign(3, 1, $v);
my $shkl = new Math::MatrixReal(3, 1);
$v = &GetArg('sh'); $v = eval($v);
$shkl->assign(1, 1, $v);
$v = &GetArg('sk'); $v = eval($v);
$shkl->assign(2, 1, $v);
$v = &GetArg('sl'); $v = eval($v);
$shkl->assign(3, 1, $v);
my $sT = new Math::MatrixReal(3, 3);
my @PresetRule = &GetArg('PresetRule');
my $PresetRule = $PresetRule[0];
print "PresetRule: $PresetRule$LF";
if($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') {
$sT = Math::MatrixReal->new_from_cols(
[ [ 0.5, 0.5, 0],
[ 0, 0.5, 0.5],
[ 0.5, 0, 0.5] ] );
$sT->transpose($sT);
}
elsif($PresetRule eq 'BCCPrim') {
$sT = Math::MatrixReal->new_from_cols(
[ [-0.5, 0.5, 0.5],
[ 0.5,-0.5, 0.5],
[ 0.5, 0.5,-0.5] ] );
$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 {
print " Rule is specified by Matrix:$LF";
$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 "Direction: $DirectionSelect$LF";
print "Check consistency:$CheckConsistency$LF";
#基本ベクトルの変換行列
my $T = new Math::MatrixReal(3, 3);
print "Conversion rule:
\n";
my $IsChosen = 0;
if($ParametersSelect eq 'ai') {
if($DirectionSelect eq 'OriginalToConverted') {
print "Real space vectors: (a'i) = "
."(tij)(ai)$LF";
$T = $sT;
$IsChosen = 1;
}
elsif($DirectionSelect eq 'ConvertedlToOriginal') {
print "Real space vectors: (ai) = "
."(tij)(a'i)$LF";
$T = $sT->inverse();
$IsChosen = 1;
}
}
elsif($ParametersSelect eq 'xyz') {
if($DirectionSelect eq 'OriginalToConverted') {
print "Real space coordinates: (x'i) = "
."(tij)(xi)$LF";
$T->transpose($sT);
$T = $T->inverse();
$IsChosen = 1;
}
elsif($DirectionSelect eq 'ConvertedlToOriginal') {
print "Real space vectors: (xi) = "
."(tij)(x'i)$LF";
$T->transpose($sT);
$IsChosen = 1;
}
}
elsif($ParametersSelect eq 'a*i') {
if($DirectionSelect eq 'OriginalToConverted') {
print "Reciprocal space vectors: (a'*i) = "
."(tij)(a*i)$LF";
$T->transpose($sT);
$T = $T->inverse();
$IsChosen = 1;
}
elsif($DirectionSelect eq 'ConvertedlToOriginal') {
print "Reciprocal space vectors: (a*i) = "
."(tij)(a'*i)$LF";
$T->transpose($sT);
$IsChosen = 1;
}
}
elsif($ParametersSelect eq 'hkl') {
if($DirectionSelect eq 'OriginalToConverted') {
print "Reciprocal space coordinates: (h'i) = "
."(tij)(hi)$LF";
$T = $sT;
$IsChosen = 1;
}
elsif($DirectionSelect eq 'ConvertedlToOriginal') {
print "Reciprocal space coordinates: (hi) = "
."(tij)(h'i)$LF";
$T = $sT->inverse();
$IsChosen = 1;
}
}
unless($IsChosen) {
print "Selected option {Parameters: $ParametersSelect} in {$DirectionSelect}"
." direction is not supported.$LF";
return 1;
}
print "Input matrix: (tij)$LF";
print "\n";
printf " |%10.4f %10.4f %10.4f|\n",
$sT->element(1,1), $sT->element(1,2), $sT->element(1,3);
printf " |%10.4f %10.4f %10.4f|\n",
$sT->element(2,1), $sT->element(2,2), $sT->element(2,3);
printf " |%10.4f %10.4f %10.4f|\n",
$sT->element(3,1), $sT->element(3,2), $sT->element(3,3);
print "
\n";
print "
\n";
#実空間べクトル(ai)の変換行列
print "Real space vector: (a'i) = "
."(Tij)(ai)$LF";
print "\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);
print "
\n";
#実空間座標(x,y,z)の変換行列
my $tRT = new Math::MatrixReal(3, 3);
$tRT->transpose($T);
$tRT = $tRT->inverse();
print "Real space coordinate conversion matrix (x') = "
."((Tij)t)-1(x)$LF";
print "\n";
printf " |%10.4f %10.4f %10.4f|\n",
$tRT->element(1,1), $tRT->element(1,2), $tRT->element(1,3);
printf " |%10.4f %10.4f %10.4f|\n",
$tRT->element(2,1), $tRT->element(2,2), $tRT->element(2,3);
printf " |%10.4f %10.4f %10.4f|\n",
$tRT->element(3,1), $tRT->element(3,2), $tRT->element(3,3);
print "
\n";
#逆空間ベクトル(a*i)の変換行列: $tRT
my $RT = new Math::MatrixReal(3, 3);
$RT->transpose($T);
$RT = $RT->inverse();
print "Reciprocal space vector conversion matrix (a*') = "
."((Tij)t)-1(a*)$LF";
print "\n";
printf " |%10.4f %10.4f %10.4f|\n",
$RT->element(1,1), $RT->element(1,2), $RT->element(1,3);
printf " |%10.4f %10.4f %10.4f|\n",
$RT->element(2,1), $RT->element(2,2), $RT->element(2,3);
printf " |%10.4f %10.4f %10.4f|\n",
$RT->element(3,1), $RT->element(3,2), $RT->element(3,3);
print "
\n";
print "$LF";
#逆空間座標(h k l)の変換行列: $T
my $tR = new Math::MatrixReal(3, 3);
$tR = $T;
print "Reciprocal space coordinate conversion matrix (h') = (Tij)(h)$LF";
print "\n";
printf " |%10.4f %10.4f %10.4f|\n",
$tR->element(1,1), $tR->element(1,2), $tR->element(1,3);
printf " |%10.4f %10.4f %10.4f|\n",
$tR->element(2,1), $tR->element(2,2), $tR->element(2,3);
printf " |%10.4f %10.4f %10.4f|\n",
$tR->element(3,1), $tR->element(3,2), $tR->element(3,3);
print "
\n";
print "$LF";
print "
\n";
#座標(sxyz)の変換: tRT
if($sxyz->element(1,1) ne '') {
my $a1 =
my $txyz = $tRT * $sxyz;
print "Real space coordinate conversion:$LF";
printf "(%10.4f, %10.4f, %10.4f) => (%10.4f, %10.4f, %10.4f)$LF",
$sxyz->element(1,1), $sxyz->element(2,1), $sxyz->element(3,1),
$txyz->element(1,1), $txyz->element(2,1), $txyz->element(3,1);
print "$LF";
}
#逆空間座標(shkl)の変換: tR
if($shkl->element(1,1) ne '') {
my $thkl = $tR * $shkl;
print "Reciprocal space coordinate conversion:$LF";
printf "(%10.4f, %10.4f, %10.4f) => (%10.4f, %10.4f, %10.4f)$LF",
$shkl->element(1,1), $shkl->element(2,1), $shkl->element(3,1),
$thkl->element(1,1), $thkl->element(2,1), $thkl->element(3,1);
print "$LF";
}
#格子定数の変換
if($a ne '') {
print "Original lattice parameters:$LF";
printf "a=%12.6f b=%12.6f c=%12.6f alpha=%8.4f beta=%8.4f gamma=%8.4f$LF",
$a, $b, $c, $alpha, $beta, $gamma;
my $Crystal = new Crystal;
$Crystal->SetLatticeParameters($a,$b,$c,$alpha,$beta,$gamma);
my $NewCrystal = $Crystal->ConvertLattice($T, 1);
my ($a2,$b2,$c2,$alpha2,$beta2,$gamma2) = $NewCrystal->LatticeParameters();
print "Converted lattice parameters:$LF";
printf "a'=%12.6f b'=%12.6f c'=%12.6f alpha'=%8.4f beta'=%8.4f gamma'=%8.4f$LF",
$a2, $b2, $c2, $alpha2, $beta2, $gamma2;
print "$LF";
my ($g11, $g12, $g13, $g21, $g22, $g23, $g31, $g32, $g33)
= $Crystal->Metrics();
print "Real space metrics conversion:$LF";
print "\n";
printf " |%10.4f %10.4f %10.4f|\n", $g11, $g12, $g13;
printf " |%10.4f %10.4f %10.4f|\n", $g21, $g22, $g23;
printf " |%10.4f %10.4f %10.4f|\n", $g31, $g32, $g33;
print "
\n";
print "=>$LF";
($g11, $g12, $g13, $g21, $g22, $g23, $g31, $g32, $g33)
= $NewCrystal->Metrics();
print "\n";
printf " |%10.4f %10.4f %10.4f|\n", $g11, $g12, $g13;
printf " |%10.4f %10.4f %10.4f|\n", $g21, $g22, $g23;
printf " |%10.4f %10.4f %10.4f|\n", $g31, $g32, $g33;
print "
\n";
}
my $ret = &UploadFile('UploadCIFFilePath', $WorkingDir);
if($ret < 0) {
print "Error in UploadFile: Upload failed due to Code $ret.
\n";
return;
}
else {
my $InputFile = $ret;
print "
\n";
print "Convert CIF File ($InputFile).$LF";
print "Links:$LF";
&ShowFileLink($WebRootDir, $InputFile, "Input CIF File");
if($InputFile !~ /\.cif$/i) {
print "Error: ($InputFile) is not a CIF file.
\n";
return;
}
my $CIF = new CIF();
unless($CIF->Read($InputFile)) {
print "Error: Can not read $InputFile.
\n";
return -2;
}
my $SimplifiedCIFFile = $InputFile;
$SimplifiedCIFFile =~ s/\.cif/-simple\.cif/i;
if($CIF->WriteSimpleCIFFile($SimplifiedCIFFile)) {
&ShowFileLink($WebRootDir, $SimplifiedCIFFile, "Simplified CIF File");
}
my $Crystal = $CIF->GetCCrystal();
$Crystal->ExpandCoordinates();
my ($a,$b,$c,$alpha,$beta,$gamma) = $Crystal->LatticeParameters();
print "Lattice parameters:$LF";
print "a=$a b=$b c=$c alpha=$alpha beta=$beta gamma=$gamma$LF";
print "Volume of original lattice: ", $Crystal->Volume(), " A3$LF";
print "Density of original lattice: ", $Crystal->Density(),
" g/cm3$LF";
print "$LF";
my $NewCrystal = $Crystal->ConvertLattice($T, 0, $CheckConsistency);
my ($a2,$b2,$c2,$alpha2,$beta2,$gamma2) = $NewCrystal->LatticeParameters();
my ($g11, $g12, $g13, $g21, $g22, $g23, $g31, $g32, $g33)
= $Crystal->Metrics();
print "Real space metrics conversion:$LF";
print "\n";
printf " |%10.4f %10.4f %10.4f|\n", $g11, $g12, $g13;
printf " |%10.4f %10.4f %10.4f|\n", $g21, $g22, $g23;
printf " |%10.4f %10.4f %10.4f|\n", $g31, $g32, $g33;
print "
\n";
print "=>$LF";
($g11, $g12, $g13, $g21, $g22, $g23, $g31, $g32, $g33)
= $NewCrystal->Metrics();
print "\n";
printf " |%10.4f %10.4f %10.4f|\n", $g11, $g12, $g13;
printf " |%10.4f %10.4f %10.4f|\n", $g21, $g22, $g23;
printf " |%10.4f %10.4f %10.4f|\n", $g31, $g32, $g33;
print "
\n";
my $OutputFile = $InputFile;
$OutputFile =~ s/\.cif/-converted\.cif/i;
print "Converted Lattice parameters:$LF";
print "a'=$a2 b'=$b2 c'=$c2 alpha'=$alpha2 beta'=$beta2 gamma'=$gamma2$LF";
print "Volume of converted lattice: ", $NewCrystal->Volume(), " A3$LF";
print "Density of converted lattice: ", $NewCrystal->Density(),
" g/cm3$LF";
my $NewCIF = new CIF();
$NewCIF->SetFileName($OutputFile);
$NewCIF->SetCrystalName($Crystal->CrystalName());
$NewCIF->SetLatticeParameters($a2, $b2, $c2, $alpha2, $beta2, $gamma2);
$NewCIF->SetSpaceGroup("P 1", 1);
$NewCIF->SetVolume($NewCrystal->Volume());
$NewCIF->AddSymmetryOperation("x,y,z");
my @ExpandedAtomSiteList = $NewCrystal->GetCExpandedAtomSiteList();
my $nTotalExpandedAtomSite = $NewCrystal->nTotalExpandedAtomSite();
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();
$NewCIF->AddAsymmetricAtomSite($label, $atomname, $x, $y, $z,$occ);
}
$NewCIF->FillCIFData();
unless($NewCIF->WriteSimpleCIFFile($OutputFile)) {
print "Error: Could not write to ($OutputFile).
\n";
return -1;
}
&ShowFileLink($WebRootDir, $OutputFile, "Converted CIF File");
if(abs($Crystal->Density() - $NewCrystal->Density()) / $Crystal->Density()
> 0.05) {
print "Causion!!.
\n";
print "Densities of the original and converted cells are "
."different.
\n";
print "The conversion rule employed would not be compatible with "
."the original cell.
\n";
}
if($PresetRule eq 'HexOrtho' or
$PresetRule eq 'FCCPrim' or
$PresetRule eq 'BCCPrim') {
print "Find symmetry for the converted lattice
";
my ($drive, $dir, $filename, $ext, $lastdir, $filebody)
= Deps::SplitFilePath($InputFile);
my $SampleName = $filebody;
#my $OutputFile = $InputFile;
# $OutputFile =~ s/\.cif/-converted\.cif/i;
my $StructPath = $InputFile;
$StructPath =~ s/\.cif/\.struct/i;
my $SymStructPath = $InputFile;
$SymStructPath =~ s/\.cif/\-Symmetrized.struct/i;
my $SymCIFPath = $InputFile;
$SymCIFPath =~ s/\.cif/\-Symmetrized.cif/i;
my ($a0,$b0,$c0,$alpha0,$beta0,$gamma0) = $NewCrystal->LatticeParameters();
if($PresetRule eq 'HexOrtho') {
$NewCrystal->SetLatticeParameters($a0, $a0*1.3, $c0,$alpha0,$beta0,$gamma0);
}
elsif($PresetRule eq 'FCCPrim' or $PresetRule eq 'BCCPrim') {
$NewCrystal->SetLatticeParameters($a0,$b0, $c0, 65,65,65);
}
print "MakeStruct: Save to $StructPath
";
my $WIEN = new WIEN2k;
$WIEN->SetSampleName($SampleName);
$WIEN->SaveStructFile($NewCrystal, $StructPath, 0, 0, 0, 0);
&ShowFileLink($WebRootDir, $StructPath, "Struct File");
print "Execute sgroup: Save to $SymStructPath
";
print("$SGROUPPath -wi $StructPath$LF");
open(PIPEIN, "$SGROUPPath -wi $StructPath |");
open(OUT, ">$StructPath.out");
while() {
print OUT $_;
}
close(OUT);
close(PIPEIN);
print("$SGROUPPath -wi -wo $StructPath $SymStructPath$LF");
system("$SGROUPPath -wi -wo $StructPath $SymStructPath");
&ShowFileLink($WebRootDir, "$StructPath.out", "sgroup output");
&ShowFileLink($WebRootDir, $SymStructPath, "Symmetrized Struct File");
my $WIEN2k2 = new WIEN2k;
my $SymCrystal = $WIEN2k2->ReadStructFile($SymStructPath, 0);
unless($SymCrystal) {
print("Error: Can not read [$SymStructPath].$LF");
return 0;
}
my $SPG = $SymCrystal->GetCSpaceGroup();
my $SPGName = $SPG->SPGName();
if($SPGName =~ /^R/i or $SPGName =~ /^P\s*\-?3/i) {
$SPG->SetSPGName("$SPGName R");
my ($a,$b,$c,$alpha,$beta,$gamma) = $SymCrystal->LatticeParameters();
print(" Hexagonal axis: a=$a c=$c$LF");
my $ar = sqrt(3.0*$a*$a + $c*$c) / 3.0;
my $alphah = 1.5 / sqrt(3.0 + ($c/$a)*($c/$a));
$alphah = 2.0 * Sci::asin($alphah) * Sci::todeg();
$alphah = $alpha0 if(abs($alphah -65.0) < 1.0e-3);
print(" Rhombohedral axis: a=$ar alpha=$alphah$LF");
$SymCrystal->SetLatticeParameters($ar,$ar,$ar,$alphah,$alphah,$alphah);
}
my $CIF2 = new CIF;
if(!$CIF2->CreateCIFFileFromCCrystal($SymCrystal, $SymCIFPath, 0, 'unix')) {
print("Error: Can not create [$SymCIFPath].$LF");
return 0;
}
&ShowFileLink($WebRootDir, $SymCIFPath, "Symmetrized CIF File");
#ここでは使わないが、念のため、後始末でもとにもどす
$NewCrystal->SetLatticeParameters($a0,$b0,$c0,$alpha0,$beta0,$gamma0);
}
}
return 1;
}