"""
結晶構造を管理するためのパッケージ。
このモジュールは、`tkCrystal` クラスを提供し、結晶の格子パラメータ、対称性、
原子サイト、物性値などのデータ管理と操作を可能にします。
CIF (Crystallographic Information File) などの結晶構造データを取り扱い、
情報を表示したり、3Dで可視化したりする機能を含みます。
関連リンク: :doc:`tkcrystal_usage`
"""
import numpy as np
from numpy import sin, cos, tan, arccos, arcsin, arctan, sqrt, exp, log, pi
from numpy import linalg as la
import csv
import sys
from pprint import pprint
import matplotlib.pyplot as plt
from tklib.tkobject import tkObject
from tklib.tkcrystal.tkcrystalobject import tkCrystalObject
from tklib.tkcrystal.tkspacegroup import tkSpaceGroup
from tklib.tkcrystal.tkatomtype import tkAtomType
from tklib.tkcrystal.tkatomsite import tkAtomSite
from tklib.tksci.tksci import a0, NA, pi, pi2, torad, todeg, round, Factorize, IsFactor, Factors
from tklib.tksci.tksci import degcos, degsin, degtan, acos, asin, atan, degacos, degasin, degatan
from tklib.tkutils import SplitFilePath, colors
import tklib.tkre
from tklib.tkre import DelQuote
[ドキュメント]
class tkCrystal(tkCrystalObject):
"""
結晶構造データを管理するクラス。
`tkCrystalObject` を継承し、結晶の格子パラメータ、対称性、原子サイト、
物性値などの詳細な情報を格納し、操作するための機能を提供します。
結晶構造の読み込み、情報の表示、構造の可視化などをサポートします。
"""
def __init__(self, **args):
"""
tkCrystalオブジェクトを初期化します。
親クラスの初期化を行い、独自の初期化メソッド `Initialize` を呼び出し、
引数で渡された情報をオブジェクトに設定します。
:param args: キーワード引数。オブジェクトの属性を初期化するために使用されます。
:returns: なし
"""
super(tkCrystal, self).__init__(**args)
self.Initialize()
self.update(**args)
def __del__(self):
"""
tkCrystalオブジェクトが破棄される際に呼び出されるデストラクタ。
親クラスのデストラクタを呼び出します。
:returns: なし
"""
super(tkCrystal, self).__del__()
def __str__(self):
"""
オブジェクトの文字列表現を返します。
クラスの完全なパスを文字列として返します。
:returns: クラスのパスを示す文字列 (str)。
"""
return self.ClassPath()
[ドキュメント]
def Initialize(self):
"""
オブジェクトの初期設定を行います。
現時点では特に処理はありませんが、将来的な拡張のためのプレースホルダーとして機能します。
:returns: なし
"""
pass
[ドキュメント]
def print_inf(self, key = None):
"""
結晶の各種情報を標準出力に表示します。
結晶名、格子系、空間群、格子パラメータ、体積、組成、密度、原子タイプ、
原子サイトなどの詳細情報をカテゴリ別に表示します。
`key` 引数で表示する情報のカテゴリを指定できます。
:param key: 表示する情報のカテゴリ (str または None) を指定します。
"header", "sample", "symmetry", "cell", "lattice",
"reciprocal", "properties", "atomtypes", "asymsites", "allsites"
などのカテゴリが利用可能です。Noneの場合、全てのカテゴリを表示します。
:returns: なし
"""
cry = self
if key is None or "header" in key:
print("")
print("====== tkcrystal.Print start ======")
try:
print("CIF path = ", cry.path)
except:
pass
if key is None or "sample" in key:
print("CrystalName = ", cry.CrystalName(Create = 1))
print("SampleName = ", cry.SampleName(Create = 1))
# print("CrystalName = ", cry.__CrystalName)
# print("SampleName = ", cry.__SampleName)
if key is None or "symmetry" in key:
print("")
print("Symmetry:")
ls = cry.LatticeSystem()
print("Lattice system: ", ls)
la = cry.LatticeAxis()
print("Lattice axis: ", la)
SPGName, iSPG, iSet = cry.GetSpaceGroupInf()
print("Space group: [{}] #{} set {}".format(SPGName, iSPG, iSet))
for i in range(self.nSymmetryOperation()):
sym = self.SymmetryOperation(i)
print(" {}: {}".format(i,sym))
if key is None or "cell" in key:
print("")
print("Real space:")
a, b, c, alpha, beta, gamma = cry.LatticeParameters()
print("cell: {:12.8f} {:12.8f} {:12.8f} A {:10.6f} {:10.6f} {:10.6f}".format(a, b, c, alpha, beta, gamma))
if key is None or "lattice" in key:
gij, Rgij = cry.Metrics()
print("Metrics gij:")
pprint(gij)
aij = cry.LatticeVectors()
print("Lattice vectors aij:")
pprint(aij)
# cry.CalculateVolume()
V = cry.Volume()
print("Volume: {} A^3".format(V))
if key is None or "reciprocal" in key:
print("")
print("Reciprocal space:")
ra, rb, rc, ralpha, rbeta, rgamm = cry.ReciprocalLatticeParameters()
print("cell: {:12.6f} {:12.6f} {:12.6f} A^-1 {} {} {}".format(ra, rb, rc, ralpha, rbeta, rgamm))
print("Metrics Rgij (A):")
pprint(Rgij)
Raij = cry.ReciprocalLatticeVectors()
print("Lattice vectors Raij (A^-1):")
pprint(Raij)
RV = cry.ReciprocalVolume()
print("Volume: {} A^-3".format(RV))
if key is None or "properties" in key:
print("")
print("Properties:")
cc = cry.ChemicalComposition()
sumcc = cry.SumChemicalComposition()
fu = cry.FormulaUnit()
MW = round(cry.UnitcellWeight(), 6)
d = cry.Density()
ad = cry.AtomDensity()
print(" Chemical composition: ", cc)
print(" Sum chemical composition: ", sumcc)
print(" Formula unit: ", fu)
print(" Unit cell weight: ", MW)
print(" Density: ", d, "g/cm3")
print(" Atom density: {:12.6e} /cm3".format(ad))
if key is None or "atomtypes" in key:
print("")
print("Atom types:")
AtomTypes = self.AtomTypeList()
for t in AtomTypes:
typea = t.AtomType()
typeo = t.AtomTypeOnly()
charge = t.Charge()
AN = t.AtomicNumber()
M = t.AtomicMass()
try:
print(" %4s (%2s) charge=%8.4f [Z=%3d M=%6.4f]" % (typea, typeo, charge, AN, M))
except:
pass
if key is None or "asymsites" in key:
print("")
print("Asymmetric atom sites:")
AtomSites = self.AtomSiteList()
for atom in AtomSites:
label = atom.Label()
atomname = atom.AtomName()
atomname0 = atom.AtomNameOnly()
charge = atom.Charge()
pos = atom.Position()
occ = atom.Occupancy()
id = atom.IdAsymmetricAtomSite()
m = atom.Multiplicity()
print(" %3d: %5s: %4s charge=%6.3f (%6.3f, %6.3f, %6.3f) occ=%8.4f m=%d"
% (id, label, atomname, charge, pos[0], pos[1], pos[2], occ, m))
if key is None or "allsites" in key:
print("")
print("Expanded atom sites:")
AtomSites = self.ExpandedAtomSiteList()
for atom in AtomSites:
label = atom.Label()
atomname = atom.AtomName()
atomname0 = atom.AtomNameOnly()
charge = atom.Charge()
pos = atom.Position()
occ = atom.Occupancy()
id = atom.IdAsymmetricAtomSite()
m = atom.Multiplicity()
print(" %3d: %5s: %4s charge=%6.3f (%6.3f, %6.3f, %6.3f) occ=%8.4f m=%d"
% (id, label, atomname, charge, pos[0], pos[1], pos[2], occ, m))
if key is None or "header" in key:
print("====== tkcrystal.Print end ======")
[ドキュメント]
def PrintInf(self, key = None):
"""
`print_inf` メソッドのエイリアス。
互換性のため、`print_inf` と同じ機能を提供します。
:param key: 表示する情報のカテゴリ (str または None) を指定します。
`print_inf` メソッドの`key`引数と同様です。
:returns: なし
"""
return self.print_inf(key)
[ドキュメント]
def to_dict(self):
"""
結晶構造データを辞書形式で返します。
Materials Project (MP) のStructureオブジェクトに似た形式で、
結晶の空間群、格子情報、原子サイト、組成などを辞書として構造化して返します。
:returns: 結晶構造データを表現する辞書 (dict)。
"""
crystal_name = self.CrystalName(Create = 1)
sample_name = self.SampleName(Create = 1)
SPGName, iSPG, iSet = self.GetSpaceGroupInf()
ls = self.LatticeSystem()
la = self.LatticeAxis()
spacegroup = { "source": "tklib.tkcrystal.tkspacegroup",
"symbol": SPGName,
"number": iSPG,
"set" : iSet,
"point_group": "null",
"crystal_system": ls,
"hall": "null"
}
a, b, c, alpha, beta, gamma = self.LatticeParameters()
a = round(a, 8)
b = round(a, 8)
c = round(a, 8)
alpha = round(a, 8)
beta = round(a, 8)
gamma = round(a, 8)
V = round(self.Volume(), 8)
gij, Rgij = self.Metrics()
aij = self.LatticeVectors()
cc = self.ChemicalComposition()
sumcc = self.SumChemicalComposition()
fu = self.FormulaUnit()
MW = round(self.UnitcellWeight(), 6)
d = round(self.Density(), 8)
ad = round(self.AtomDensity(), 8)
AtomTypes = self.AtomTypeList()
elements = []
for t in AtomTypes:
typea = t.AtomType()
typeo = t.AtomTypeOnly()
charge = t.Charge()
AN = t.AtomicNumber()
M = t.AtomicMass()
elements.append(typeo)
# AtomSites = self.AtomSiteList()
AtomSites = self.ExpandedAtomSiteList()
total_charge = 0.0
sites = []
for atom in AtomSites:
label = atom.Label()
atomname = atom.AtomName()
atomname0 = atom.AtomNameOnly()
charge = atom.Charge()
pos = atom.Position()
pos = [round(pos[i], 8) for i in range(3)]
xc, yc, zc = self.FractionalToCartesian(*pos)
f = atom.Force()
if f[0] is None:
f = ["null", "null", "null"]
v = atom.Velocity()
if v[0] is None:
v = ["null", "null", "null"]
occ = atom.Occupancy()
id = atom.IdAsymmetricAtomSite()
m = atom.Multiplicity()
total_charge += charge
species = { "element": atomname, "occu": occ, "charge": charge,
"xyz": [xc, yc, zc], "abc": pos, "force": f, "velocity": v,
"properties": {} }
site = { "label": label, "species": [species] }
sites.append(site)
lattice = {
"matrix": aij.tolist(),
"a": a,
"b": b,
"c": c,
"alpha": alpha,
"beta" : beta,
"gamma": gamma,
"volume": V,
}
d = {}
d["structure"] = {
"@module": "tklib.tkcrystal.tkcrystal",
"@class": "Structure",
"charge": total_charge,
"spacegroup": spacegroup,
"elements": elements,
"nelements": len(elements),
"nsites": len(AtomSites),
"lattice": lattice,
"sites": sites
}
return d
[ドキュメント]
def draw_box(self, ax, aij, nrange, linecolor = 'black'):
"""
3Dプロットに単位格子を描画します。
与えられた格子ベクトルに基づいて、matplotlibの3D軸上に単位格子の辺を描画します。
:param ax: matplotlibの3D軸オブジェクト (Axes3D)。
:param aij: 格子ベクトルを表す3x3行列 (numpy.ndarray)。
:param nrange: 描画範囲。現在は未使用ですが、引数として存在します (list)。
:param linecolor: 単位格子の辺の色 (str)。デフォルトは 'black' です。
:returns: なし
"""
# (0,0,0) -> ax
ax.plot([0.0, aij[0][0]],
[0.0, aij[0][1]],
[0.0, aij[0][2]], color = linecolor)
# (0,0,0) -> ay
ax.plot([0.0, aij[1][0]],
[0.0, aij[1][1]],
[0.0, aij[1][2]], color = linecolor)
# (0,0,0) -> az
ax.plot([0.0, aij[2][0]],
[0.0, aij[2][1]],
[0.0, aij[2][2]], color = linecolor)
# ax -> ax + ay
ax.plot([aij[0][0], aij[0][0] + aij[1][0]],
[aij[0][1], aij[0][1] + aij[1][1]],
[aij[0][2], aij[0][2] + aij[1][2]], color = linecolor)
# ax -> ax + az
ax.plot([aij[0][0], aij[0][0] + aij[2][0]],
[aij[0][1], aij[0][1] + aij[2][1]],
[aij[0][2], aij[0][2] + aij[2][2]], color = linecolor)
# ay -> ay + ax
ax.plot([aij[1][0], aij[1][0] + aij[0][0]],
[aij[1][1], aij[1][1] + aij[0][1]],
[aij[1][2], aij[1][2] + aij[0][2]], color = linecolor)
# ay -> ay + az
ax.plot([aij[1][0], aij[1][0] + aij[2][0]],
[aij[1][1], aij[1][1] + aij[2][1]],
[aij[1][2], aij[1][2] + aij[2][2]], color = linecolor)
# az -> az + ax
ax.plot([aij[2][0], aij[2][0] + aij[0][0]],
[aij[2][1], aij[2][1] + aij[0][1]],
[aij[2][2], aij[2][2] + aij[0][2]], color = linecolor)
# az -> ax + ay
ax.plot([aij[2][0], aij[2][0] + aij[1][0]],
[aij[2][1], aij[2][1] + aij[1][1]],
[aij[2][2], aij[2][2] + aij[1][2]], color = linecolor)
# ax + ay -> ax + ay + az
ax.plot([aij[0][0] + aij[1][0], aij[0][0] + aij[1][0] + aij[2][0]],
[aij[0][1] + aij[1][1], aij[0][1] + aij[1][1] + aij[2][1]],
[aij[0][2] + aij[1][2], aij[0][2] + aij[1][2] + aij[2][2]], color = linecolor)
# ax + az -> ax + ay + az
ax.plot([aij[0][0] + aij[2][0], aij[0][0] + aij[1][0] + aij[2][0]],
[aij[0][1] + aij[2][1], aij[0][1] + aij[1][1] + aij[2][1]],
[aij[0][2] + aij[2][2], aij[0][2] + aij[1][2] + aij[2][2]], color = linecolor)
# ay + az -> ax + ay + az
ax.plot([aij[1][0] + aij[2][0], aij[0][0] + aij[1][0] + aij[2][0]],
[aij[1][1] + aij[2][1], aij[0][1] + aij[1][1] + aij[2][1]],
[aij[1][2] + aij[2][2], aij[0][2] + aij[1][2] + aij[2][2]], color = linecolor)
[ドキュメント]
def draw_unitcell(self, ax, nrange, kr = 1000.0, linecolor = 'black', color = 'black', atomcolor = None):
"""
3Dプロットに単位格子と原子サイトを描画します。
`draw_box` メソッドで単位格子を描画し、結晶内の各原子サイトを球でプロットします。
原子タイプごとに異なる色を使用します。
:param ax: matplotlibの3D軸オブジェクト (Axes3D)。
:param nrange: 描画範囲。現在は未使用ですが、引数として存在します (list)。
:param kr: 原子サイズをプロットする際の係数 (float)。原子の半径にこの係数を乗じてプロットサイズを決定します。
:param linecolor: 単位格子の辺の色 (str)。デフォルトは 'black' です。
:param color: 原子にデフォルトの色を指定する場合の色 (str)。現在は未使用。
:param atomcolor: 原子タイプごとの色のリスト (list of str)。指定がない場合はデフォルトの色を使用します。
:returns: なし
"""
if atomcolor is None:
atomcolor = colors
ncolor = len(atomcolor)
aij = self.LatticeVectors()
self.draw_box(ax, aij, nrange, linecolor)
AtomTypes = self.AtomTypeList()
atom_types = {}
for t in AtomTypes:
# typea = t.AtomType()
typeo = t.AtomTypeOnly()
# charge = t.Charge()
# AN = t.AtomicNumber()
# M = t.AtomicMass()
atom_types[typeo] = 1
c = 0
for key in atom_types.keys():
atom_types[key] = c
c += 1
AtomSites = self.ExpandedAtomSiteList()
drawn_atoms = {}
for i in range(len(AtomSites)):
atom = AtomSites[i]
# label = atom.Label()
# atomname = atom.AtomName()
atomname0 = atom.AtomNameOnly()
# charge = atom.Charge()
pos = atom.Position()
xc, yc, zc = self.FractionalToCartesian(*pos)
# occ = atom.Occupancy()
# id = atom.IdAsymmetricAtomSite()
# m = atom.Multiplicity()
atom_types
ion_radius = 1.0
acolor = atomcolor[atom_types[atomname0] % ncolor]
if drawn_atoms.get(atomname0, None) is None:
ax.scatter([xc], [yc], [zc], label = atomname0, marker = 'o', c = acolor, linewidth = 0.5, edgecolors = 'white', s = kr * ion_radius)
drawn_atoms[atomname0] = 1
else:
ax.scatter([xc], [yc], [zc], marker = 'o', c = acolor, linewidth = 0.5, edgecolors = 'white', s = kr * ion_radius)
ax.legend()
[ドキュメント]
def draw(self, figsize = (10, 10), nrange = [[-0.1, 1.1], [-0.1, 1.1], [-0.1, 1.1]], rmin = 0.1, kr = 100.0):
"""
結晶構造を3Dプロットで可視化します。
matplotlibを使用して結晶の単位格子と原子サイトを3Dプロットで表示します。
描画範囲と原子サイズ調整係数を指定できます。
:param figsize: 描画する図のサイズ (幅, 高さ) (tuple)。デフォルトは (10, 10) です。
:param nrange: 描画する単位格子の範囲 (list)。各軸の開始と終了座標を指定します。
例: [[xmin, xmax], [ymin, ymax], [zmin, zmax]]。
:param rmin: 同一原子サイトと判断する距離の閾値 (オングストローム単位) (float)。
現在は未使用です。
:param kr: 原子サイズをプロットする際の係数 (float)。`draw_unitcell` メソッドに渡されます。
:returns: なし
"""
fig = plt.figure(figsize = figsize)
ax = fig.add_subplot(111, projection = '3d')
# Real space unit cell
self.draw_unitcell(ax, nrange, kr, linecolor = 'black')
# Note: set_aspect() is not implemented for 3D plots
# ax.set_aspect('equal','box')
xlim = ax.get_xlim()
ylim = ax.get_ylim()
zlim = ax.get_zlim()
lim = [min([xlim[0], ylim[0], zlim[0]]), max([xlim[1], ylim[1], zlim[1]])]
ax.set_xlim(lim)
ax.set_ylim(lim)
ax.set_zlim(lim)
plt.show()