package X3D; use Exporter; @ISA = qw(Exporter); #use Crystal::VRML; #@ISA = qw(VRML); #公開したいサブルーチン @EXPORT = qw(); use strict; use File::Basename; use Sci::Science; use Crystal::Crystal; use Crystal::AtomType; use Crystal::AtomSite; use Crystal::CIF; BEGIN { } sub new { my ($module) = @_; my $this = {}; bless $this; return $this; } sub DESTROY { my $this = shift; } sub SetSampleName { my ($this, $name) = @_; return $this->{'SampleName'} = $name; } sub SampleName { my ($this) = @_; return $this->{'SampleName'}; } sub PrintClosure { my ($this, $out) = @_; print $out "\n"; print $out " \n"; print $out "\n"; } sub PrintHeader { my ($this, $filename, $out, $color, $fogrange) = @_; $this->PrintX3DHeader($filename, $out, $color, $fogrange); } sub PrintX3DHeader { my ($this, $filename, $out, $color, $sky, $fogrange) = @_; $color = "1 1 1" if($color eq ''); $sky = 'white' if(!defined $sky); $fogrange = -1 if(!defined $fogrange); print $out < EOT if($color eq 'sky') { print $out < --> EOT } else { print $out < EOT } if($fogrange > 0) { print $out < EOT } else { print $out < --> EOT } print $out < EOT return 1; } sub CalTransform { my ($this, $x0, $y0, $z0, $x1, $y1, $z1, $scale, $Ver) = @_; $x0 *= $scale; $y0 *= $scale; $z0 *= $scale; $x1 *= $scale; $y1 *= $scale; $z1 *= $scale; my $centerx = ($x0 + $x1) / 2.0; my $centery = ($y0 + $y1) / 2.0; my $centerz = ($z0 + $z1) / 2.0; my $centerxm = -$centerx; my $centerym = -$centery; my $centerzm = -$centerz; my $dx = $x1 - $x0; my $dy = $y1 - $y0; my $dz = $z1 - $z0; my $r = $dx*$dx + $dy*$dy + $dz*$dz; $r = sqrt($r); my $cosQ = $dy / $r; my $Q = Sci::acos($cosQ); my $dxm = -$dx; my $dzm = -$dz; my $Qm = -$Q; my $centerstr = sprintf("%8.4f %8.4f %8.4f", $centerx, $centery, $centerz); # my $centermstr = sprintf("%8.4f %8.4f %8.4f", $centerxm, $centerym, $centerzm); my $rotationstr = sprintf("%8.4f %8.4f %8.4f %8.4f", $dz, 0, $dxm, $Q); # my $rotationmstr = sprintf("%8.4f %8.4f %8.4f %8.4f", $dz, 0, $dxm, $Qm); if(abs($dx) <= 1.0e-4 and abs($dz) <= 1.0e-4) { $rotationstr = sprintf("%8.4f %8.4f %8.4f %8.4f", 1.0, 0, $dxm, $Q); # $rotationmstr = sprintf("%8.4f %8.4f %8.4f %8.4f", 1.0, 0, $dxm, $Qm); } return ($r, $centerstr, $rotationstr); } sub PrintCylinder { my ($this, $out, $x0, $y0, $z0, $x1, $y1, $z1, $cellbarradius, $scale, $Ver) = @_; # $cellbarradius *= $scale; $cellbarradius = sprintf("%8.4f", $cellbarradius * $scale); my ($r, $centerstr, $rotationstr) = $this->CalTransform($x0, $y0, $z0, $x1, $y1, $z1, $scale, $Ver); $r = sprintf("%8.4f", $r); print $out < EOT1 return 1; } sub PrintSphere { my ($this, $out, $rx, $ry, $rz, $r, $color, $scale, $rscale, $Ver) = @_; $rscale = 1.0 if($rscale <= 0.0); $rx *= $scale; $ry *= $scale; $rz *= $scale; $r *= $scale*$rscale; my $rxm = -$rx; my $rym = -$ry; my $rzm = -$rz; $rx = sprintf("%8.4f", $rx); $ry = sprintf("%8.4f", $ry); $rz = sprintf("%8.4f", $rz); $rxm = sprintf("%8.4f", $rxm); $rym = sprintf("%8.4f", $rym); $rzm = sprintf("%8.4f", $rzm); $r = sprintf("%8.4f", $r); print $out < EOT1 return 1; } sub PrintText { my ($this, $out, $x, $y, $z, $text, $size, $scale, $Ver) = @_; my $FontFamily = "SERIF"; #"TYPEWRITER" $x *= $scale; $y *= $scale; $z *= $scale; my $xm = -$x; my $ym = -$y; my $zm = -$z; $size *= $scale; $y -= $size * 1.0; $x = sprintf("%8.4f", $x); $y = sprintf("%8.4f", $y); $z = sprintf("%8.4f", $z); $xm = sprintf("%8.4f", $xm); $ym = sprintf("%8.4f", $ym); $zm = sprintf("%8.4f", $zm); $size = sprintf("%9.4f", $size); print $out < EOT1 return 1; } sub PrintCellRod { my ($this, $out, $x0, $y0, $z0, $x1, $y1, $z1, $cellbarradius, $scale, $Ver) = @_; # $this->PrintCylinder($out, $x0, $y0, $z0, $x1, $y1, $z1, $cellbarradius, $scale, $Ver); my ($r, $centerstr, $rotationstr) = $this->CalTransform($x0, $y0, $z0, $x1, $y1, $z1, $scale, $Ver); print $out < EOT } sub PrintLattice { my ($this, $out, $Crystal, $cellbarradius, $scale, $Ver) = @_; my ($a11,$a13,$a12, $a31,$a33,$a32, $a21,$a23,$a22) = $Crystal->LatticeVectors(); #print "Latt: $a11,$a12,$a13
\n"; #print "Latt: $a21,$a22,$a23
\n"; #print "Latt: $a31,$a32,$a33
\n"; $this->PrintCellRod($out, 0, 0, 0, $a11, $a12, $a13, $cellbarradius, $scale, $Ver); $this->PrintCellRod($out, 0+$a21, 0+$a22, 0+$a23, $a11+$a21, $a12+$a22, $a13+$a23, $cellbarradius, $scale, $Ver); $this->PrintCellRod($out, 0+$a31, 0+$a32, 0+$a33, $a11+$a31, $a12+$a32, $a13+$a33, $cellbarradius, $scale, $Ver); $this->PrintCellRod($out, 0+$a21+$a31, 0+$a22+$a32, 0+$a23+$a33, $a11+$a21+$a31, $a12+$a22+$a32, $a13+$a23+$a33, $cellbarradius, $scale, $Ver); $this->PrintCellRod($out, 0, 0, 0, $a21, $a22, $a23, $cellbarradius, $scale, $Ver); $this->PrintCellRod($out, 0+$a11, 0+$a12, 0+$a13, $a21+$a11, $a22+$a12, $a23+$a13, $cellbarradius, $scale, $Ver); $this->PrintCellRod($out, 0+$a31, 0+$a32, 0+$a33, $a21+$a31, $a22+$a32, $a23+$a33, $cellbarradius, $scale, $Ver); $this->PrintCellRod($out, 0+$a11+$a31, 0+$a12+$a32, 0+$a13+$a33, $a21+$a11+$a31, $a22+$a12+$a32, $a23+$a13+$a33, $cellbarradius, $scale, $Ver); $this->PrintCellRod($out, 0, 0, 0, $a31, $a32, $a33, $cellbarradius, $scale, $Ver); $this->PrintCellRod($out, 0+$a11, 0+$a12, 0+$a13, $a31+$a11, $a32+$a12, $a33+$a13, $cellbarradius, $scale, $Ver); $this->PrintCellRod($out, 0+$a21, 0+$a22, 0+$a23, $a31+$a21, $a32+$a22, $a33+$a23, $cellbarradius, $scale, $Ver); $this->PrintCellRod($out, 0+$a11+$a21, 0+$a12+$a22, 0+$a13+$a23, $a31+$a11+$a21, $a32+$a12+$a22, $a33+$a13+$a23, $cellbarradius, $scale, $Ver); } sub PrintAtom { my ($this, $out, $label, $rx, $ry, $rz, $scale) = @_; # $this->PrintSphere($out, $rx, $ry, $rz, 0.2, '1 0 0', $scale, $scale, 3.0); $rx *= $scale; $ry *= $scale; $rz *= $scale; print $out < EOT } sub PrintChemicalBond { my ($this, $out, $plabel, $x0, $y0, $z0, $x1, $y1, $z1, $scale) = @_; # $this->PrintCylinder($out, $x0, $y0, $z0, $x1, $y1, $z1, 0.1, $scale, 3.0); my ($r, $centerstr, $rotationstr) = $this->CalTransform($x0, $y0, $z0, $x1, $y1, $z1, $scale, 3.0); print $out < EOT } sub SaveFile { my ($this, $Crystal, $filename, $nx, $ny, $nz, $MaxBondLength, $nMaxBonds, $Scale, $rScale, $VRMLVersion) = @_; my $Ver = 3.0; my $LF = "
\n"; my $PrintLabel = 0; #ファイル作製開始 my $out; unless(open($out,">$filename")) { print "Can not write to $filename.$LF$LF"; return; } # Header $this->PrintX3DHeader($filename, $out, "0.3 0.3 0.3"); my $SampleName = $this->SampleName(); my $SPGName = $Crystal->SPGNameByOutputMode(); my $iSPG = $Crystal->iSPGByOutputMode(); my @AtomTypeList = $Crystal->GetCAtomTypeList(); my $nAtomType = @AtomTypeList; my @AsymmetricAtomSiteList = $Crystal->GetCAsymmetricAtomSiteList(); my @ExpandedAtomSiteList = $Crystal->GetCExpandedAtomSiteListByOutputMode(); my $nExpandedAtomSite = @ExpandedAtomSiteList; my ($a,$b,$c,$alpha,$beta,$gamma) = $Crystal->LatticeParametersByOutputMode(1); my $ColorOffset = 0.3; my @AtomColor; my @AtomRadius; for(my $i = 0 ; $i < @AtomTypeList ; $i++) { my $AtomType = $AtomTypeList[$i]; my $atomname = $AtomType->AtomNameOnly(); $AtomRadius[$i] = $AtomType->IonRadius(); #AtomicRadius(); $AtomRadius[$i] = 0.6 if($AtomRadius[$i] <= 0); $AtomRadius[$i] *= $rScale; if($AtomColor[$i] eq '') { my $red = $ColorOffset + rand(1.0-$ColorOffset); my $green = $ColorOffset + rand(1.0-$ColorOffset); my $blue = $ColorOffset + rand(1.0-$ColorOffset); $red /= 3 if($i % 3 == 0); $green /= 3 if($i % 3 == 1); $blue /= 3 if($i % 3 == 2); $red = sprintf("%7.4f", $red); $green = sprintf("%7.4f", $green); $blue = sprintf("%7.4f", $blue); $AtomColor[$i] = "$red $green $blue"; } } my $Max = $a*$nx; $Max = $b*$ny if($Max < $b*$ny); $Max = $c*$nz if($Max < $c*$nz); $Max = 1.0 if($Max <= 0.0); $Scale = $Scale / $Max * 10.0; my $CellBarRadius = 0.05; my $BondRadius = 0.1; my $centerx = -(-$nx+0.98 + 1.02) * $a / 2.0 * $Scale; my $centery = -(-$ny+0.98 + 1.02) * $b / 2.0 * $Scale; my $centerz = -(-$nz+0.98 + 1.02) * $c / 2.0 * $Scale; #print "aa ($nx): $centerx, ", $a*$Scale, " (", -$nx+0.98+1.02, ")$LF"; $centerx = sprintf("%8.4f", $centerx); $centery = sprintf("%8.4f", $centery); $centerz = sprintf("%8.4f", $centerz); my $ScaledCellBarRadius = sprintf("%8.4f", $CellBarRadius * $Scale); my $CellBarColor = '0.5 0.5 0.5'; print $out < EOT print $out < EOT for(my $i = 0 ; $i < @AsymmetricAtomSiteList ; $i++) { my $site = $AsymmetricAtomSiteList[$i]; my $atomname = $site->AtomNameOnly(); my $label = $site->Label(); my $iAtom = $Crystal->FindiAtomType($atomname); my $color = $AtomColor[$iAtom-1]; my $radius = $AtomRadius[$iAtom-1] * $Scale; if($PrintLabel) { print $out < EOTATOM } else { print $out < EOTATOM } } print $out < EOT my $BondColor = "0.9 0.9 0.9"; my $ScaledBondRadius = sprintf("%8.4f", $BondRadius * $Scale); for(my $i = 0 ; $i < @AsymmetricAtomSiteList ; $i++) { my $site1 = $AsymmetricAtomSiteList[$i]; my $label1 = $site1->Label(); for(my $j = $i ; $j < @AsymmetricAtomSiteList ; $j++) { my $site2 = $AsymmetricAtomSiteList[$j]; my $label2 = $site2->Label(); my $plabel; if($label1 gt $label2) { $plabel = "$label2$label1"; } else { $plabel = "$label1$label2"; } print $out < PROTOBOND } } print $out < EOTLATT $this->PrintLattice($out, $Crystal, $CellBarRadius, $Scale, $Ver); print $out "\n"; my $ShowCrystal = new Crystal; $ShowCrystal = $Crystal; $ShowCrystal->SetnAsymmetricAtomSite(0); $ShowCrystal->ExpandCoordinates(); print $out < EOTATOMS my $inx = int($nx-0.01) + 1; my $iny = int($ny-0.01) + 1; my $inz = int($nz-0.01) + 1; my %nAtomTypes; for(my $ix = -$inx ; $ix <= $inx ; $ix++) { for(my $iy = -$iny ; $iy <= $iny ; $iy++) { for(my $iz = -$inz ; $iz <= $inz ; $iz++) { for(my $i = 0 ; $i < $nExpandedAtomSite ; $i++) { my $AtomSite = $ExpandedAtomSiteList[$i]; my $SiteLabel = $AtomSite->Label(); #print "Label: $SiteLabel
\n"; my $atomname = $AtomSite->AtomNameOnly(); my $charge = $AtomSite->Charge(); my ($x,$y,$z) = $AtomSite->Position(1); #print "$atomname: ($x,$y,$z)+($ix,$iy,$iz)$LF"; $x += $ix; $y += $iy; $z += $iz; next if($x < -$nx+0.98 or 1.02 < $x); next if($y < -$ny+0.98 or 1.02 < $y); next if($z < -$nz+0.98 or 1.02 < $z); #print "$atomname: ($x,$y,$z)$LF"; my $occupancy = $AtomSite->Occupancy(); my ($xc, $yc, $zc) = $Crystal->FractionalToCartesian($x, $y, $z); my $iAtom = $Crystal->FindiAtomType($atomname); my $AtomType = $AtomTypeList[$iAtom-1]; my $Radius = $AtomRadius[$iAtom-1]; my $Color = $AtomColor[$iAtom-1]; $nAtomTypes{$atomname}++; my $label = $SiteLabel; my $occscale = (0.3 + $occupancy) / 1.3; $this->PrintAtom($out, $SiteLabel, $xc, $zc, $yc, $Scale); $ShowCrystal->AddAtomSite($label, $atomname, $x, $y, $z, $occupancy); } } } } print $out "\n"; print $out < EOTATOMS my @alist = $ShowCrystal->GetCAsymmetricAtomSiteList(); my $count = 0; for(my $i = 0 ; $i < @alist ; $i++) { my $atom0 = $alist[$i]; my $SiteLabel0 = $atom0->Label(); #print "SiteLabel0: $SiteLabel0 - \n"; my $name0 = $atom0->AtomNameOnly(); my ($x0,$y0,$z0) = $atom0->Position(0); my ($xc0, $yc0, $zc0) = $ShowCrystal->FractionalToCartesian($x0, $y0, $z0); for(my $j = $i+1 ; $j < @alist ; $j++) { my $atom1 = $alist[$j]; my $SiteLabel1 = $atom1->Label(); #print "SiteLabel1: $SiteLabel1
\n"; my $name1 = $atom1->AtomNameOnly(); my ($x1,$y1,$z1) = $atom1->Position(0); my $dis = $ShowCrystal->GetInterAtomicDistance( $x0, $y0, $z0, $x1, $y1, $z1); #print "($i,$j) ($name0,$name1) ($x0,$y0,$z0)-($x1,$y1,$z1) $dis$LF"; next if($dis < 0.1 or $dis > $MaxBondLength); #print "...Show ($MaxBondLength)$LF"; my ($xc1, $yc1, $zc1) = $ShowCrystal->FractionalToCartesian($x1, $y1, $z1); my $plabel; if($SiteLabel0 gt $SiteLabel1) { $plabel = "$SiteLabel1$SiteLabel0"; } else { $plabel = "$SiteLabel0$SiteLabel1"; } $this->PrintChemicalBond($out, $plabel, $xc0, $zc0, $yc0, $xc1, $zc1, $yc1, $Scale); $count++; if($count > $nMaxBonds) { print "Number of bonds reaches maximum ($nMaxBonds).$LF"; last; } } last if($count > $nMaxBonds); } $this->PrintClosure($out); close($out); #原子名のX3Dファイルを作成 my $NameFile = $filename; $NameFile =~ s/\..*?$//; $NameFile = "$NameFile-name.x3d"; my $out; unless(open($out,">$NameFile")) { print "Can not write to $NameFile.$LF$LF"; return; } $this->PrintX3DHeader($NameFile, $out, "0.3 0.3 0.3"); my $TotalRad = 0.0; my $MaxRad = 0.0; for(my $i = 0 ; $i < $nAtomType ; $i++) { my $atom = $AtomTypeList[$i]; my $atomname = $atom->AtomNameOnly(); my $Radius = $atom->IonRadius(); #AtomicRadius(); $Radius = 0.6 if($Radius <= 0); $TotalRad += $Radius; $MaxRad = $Radius if($Radius > $MaxRad); } my $ShiftY = $MaxRad * $nAtomType / 2.0 * 0.7 ; $ShiftY = sprintf("%8.4f", $ShiftY); my $xs = 0.0; my $ys = 0.0; my $zs = 0.0; for(my $i = 0 ; $i < $nAtomType ; $i++) { my $atom = $AtomTypeList[$i]; my $atomname = $atom->AtomNameOnly(); my $Radius = $atom->IonRadius(); #AtomicRadius(); print "r($atomname)=$Radius
\n"; # $Radius = 1.4 if($Radius <= 0 and $atomname eq 'O'); $Radius = 0.6 if($Radius <= 0); my $Color = $AtomColor[$i]; $this->PrintSphere($out, $xs, $ys, $zs, $Radius, $Color, $Scale, $rScale, $Ver); $this->PrintText( $out, $xs+$MaxRad*1.5, $ys+2.0*$MaxRad, $zs, $atomname, 2.0, $Scale, $Ver); $ys -= $MaxRad * 2 * 0.7; } $this->PrintClosure($out); close($out); return ($filename, $NameFile); } 1;