Pythonで結晶を描画するプログラムを作る~単位格子を描く~

投稿日：2021/09/14
更新日：2021/09/15 (グラフ理論の簡単な説明を追加)

はじめに

VESTAなどの高性能なツールにより、後述するCIFデータさえ手に入れば空間群などの理解なしに材料の結晶構造が描画できてしまいます。その便利さを一旦捨てて、ゼロから結晶を描画する事を目指します。本記事はその第一回目です。

1. CIFデータ

CIFとはCrystallographic Information Fileの略でWPCI(Working Party on Crystallographic Information)が主導となって開発された結晶構造データを共有するデータフォーマットです。このCIFは基本的に以下の6つから構成されています[1]。

1. データ公表者および出版情報
2. 基礎結晶学的情報
3. 実験・解析方法の詳細
4. 得られた結晶データおよび結晶化学的情報
5. 解析時に使用したパラメータファイル
6. 解析に用いられたデータ本体

"基本的に"という言葉を付けたのには訳があって、データの提供主を示す1や解析に必要な2や4,6以外は記述されていない事があります。今すぐCIFデータを入手したいときはAtomWorkか[Materials Project] (https://materialsproject.org/ )で探すと良いかと思います。本記事ではMaterials Projectで公開されているSrTiO$_3$ (https://materialsproject.org/materials/mp-5229/ )を使っていきます。このファイルは拡張子が".cif"になっていますが、メモ帳などで開くと中身が確認できます。このファイルを読み解いて行きます。

下記はここからSrTiO$_3$の結晶構造データですという宣言文になっています。これは複数のCIFデータを1つのファイルにまとめたときにその区別用として威力を発揮します。

data_SrTiO3

下記はSrTiO$_3$が所属する空間群("_symmetry_space_group_name_H-M")とその番号("_symmetry_Int_Tables_number")を表しています。"_symmetry_space_group_name_H-M"はこの結晶が属す空間群をHermann-Mauguin記号(以降, H-M記号)で表しており、この結晶は$Pm\bar{3}m$という空間群に属している事を意味します。結晶を扱う上で重要な230種の空間群はナンバリングされており、$Pm\bar{3}m$は$221$という数字が与えられています。空間群と番号の対応はWikipediaの[List of space groups] (https://en.wikipedia.org/wiki/List_of_space_groups )を確認してみて下さい。ここで、当然の疑問として浮かび上がるのが、「何故H-M記号と番号をどちらも載せるのか?」だと思います。明確な理由があるので、次回または次々回あたりで触れたいと思います。

_symmetry_space_group_name_H-M   'Pm-3m'
_symmetry_Int_Tables_number   221

下記の部分は結晶パラメータと呼ばれ、単位格子を構成する重要な6つのパラメータです。それぞれ、$a$,$b$,$c$軸の格子定数($Å$)、$b$-$c$,$c$-$a$,$a$-$b$軸の挟角(degree)を表します。

_cell_length_a   3.94513000
_cell_length_b   3.94513000
_cell_length_c   3.94513000
_cell_angle_alpha   90.00000000
_cell_angle_beta   90.00000000
_cell_angle_gamma   90.00000000

他の部分についても触れたい所ですが、本記事の目標である単位格子の描画からずれてしまうので割愛します。別の記事で適宜取り扱っていきます。

#2. 単位格子
結晶の単位となる平行六面体を単位格子 (unit cell) といい、結晶を分類する上でなるべく高い対称性をもつものの中で最小の体積を持つ平行六面体から選ばれます。CIFデータではこの単位格子の情報が結晶パラメータ$a,b,c,\alpha,\beta,\gamma$として保存されています。このパラメータを7つの分類に整理したものを晶系 (crystal system) と呼び、以下の表のように分類されます(Markdown形式の表作成になれたら直します)。

	$a\neq b\neq c$	$a= b,\neq c$	$a= b=c$
$\alpha\neq\beta\neq\gamma$	三斜晶 (triclinic;a)	三斜晶 (triclinic;a)	三斜晶 (triclinic;a)
$\alpha =\gamma =90^\circ \neq \beta$	単斜晶 (monoclinic;m)	単斜晶 (monoclinic;m)	単斜晶 (monoclinic;m)
$\alpha =\beta=90^\circ, \gamma=120^\circ $	単斜晶 (monoclinic;m)	六方晶 (hexagonal;h)、三方晶 (trigonal;r)	六方晶 (hexagonal;h)、三方晶 (Rhombohedral;r)
$\alpha =\beta =\gamma =90^\circ$	直方晶 (orthorhombic;o)	正方晶 (tetragonal;t)	立方晶 (cubic;c)

この表に従ってSrTiO$_3$の結晶パラメータを見ると、$a=b=c,\alpha =\beta =\gamma =90^\circ $であるため、立方晶に分類される事がわかります。この晶系はCIFデータにおける"_space_group_crystal_system"に記述されます。このCIFデータにはありませんが、"_space_group_crystal_system"の項目がある場合は下記のように記述されます。

_space_group_crystal_system    cubic

単位格子を描くためには結晶パラメータのみで良く、単位格子頂点の位置ベクトル$\boldsymbol{A},\boldsymbol{B},\boldsymbol{C}$と置いたとき以下の式で表されます。

\begin{eqnarray}
\boldsymbol{A}&=&\left(\begin{array}{c}
A_{x}\\
A_{y}\\
A_{z}
\end{array}\right)=\left(\begin{array}{c}
a\\
0\\
0
\end{array}\right), \\
\boldsymbol{B}&=&\left(\begin{array}{c}
B_{x}\\
B_{y}\\
B_{z}
\end{array}\right)=\left(\begin{array}{c}
b\cos{\gamma}\\
b\sin{\gamma}\\
0
\end{array}\right),\tag{1}\\
\boldsymbol{C}&=&\left(\begin{array}{c}
C_{x}\\
C_{y}\\
C_{z}
\end{array}\right)=\left(\begin{array}{c}
c\cos{\beta}\\
\frac{c(\cos{\alpha}-\cos\beta\cos\gamma)}{\sin{\gamma}}\\
c\sqrt{\sin^2\beta - \frac{(\cos{\alpha}-\cos\beta\cos\gamma)^2}{\sin^2\gamma}}
\end{array}\right)
\end{eqnarray}

ここで、定数$n_1,n_2,n_3$で結びつけたベクトル$\boldsymbol{R}$：

\boldsymbol{R}:=n_1 \boldsymbol{A}+ n_2 \boldsymbol{B}+ n_3\boldsymbol{C}=(n_1,n_2,n_3)\tag{2}

を考えてみます。$n_1,n_2,n_3$を$0$または$1$に限ってみると、$(n_1,n_2,n_3)$は$2^3=8$通りある事がわかります。$\boldsymbol{A},\boldsymbol{B},\boldsymbol{C}$は単位格子の位置を表している事を思い出すとこの8パターンは平行六面体である単位格子の各頂点に対応している事がわかります。

余談ですが単位格子の体積$V$は$\boldsymbol{A},\boldsymbol{B},\boldsymbol{C}$のスカラー三重積を利用すれば、

\begin{eqnarray}
V&=& \boldsymbol{A}\cdot(\boldsymbol{B}\times\boldsymbol{C})\\
&=& A_{x}B_{y}C_{z}\\
&=& abc\sqrt{\sin^2\beta\sin^2\gamma -(\cos\alpha-\cos\beta\cos\gamma)^2}
\end{eqnarray}

と求まります。ここで、右辺1式から2式目への変換は$\boldsymbol{A}\cdot(\boldsymbol{B}\times\boldsymbol{C})$が上三角行列の行列式になっている事から対角成分の積のみを考えれば良いという関係を使っています。CIFデータから記述子などを生成する際に使えそうな式だと考える人がいるかと思いますが、"_cell_volume"に書いてあるため使う機会はありません。

#3. Pythonで単位格子を作成する
今から作成するコードの流れは結晶パラメータを受け取った後、(i) 結晶パラメータの角度の単位を度数(degree)からラジアン(radian)に変換、(ii) (1)式に結晶パラメータを代入、(iii) (2)式から単位格子の頂点を導出、(iv)3次元で図示です。

import math
import numpy as np
import itertools
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 関数
def deg2rad(deg):
    """
    deg→rad
    """
    return deg*(math.pi/180)

def vecA(a):
    """
    位置ベクトルA
    """
    Ax = a
    Ay = 0
    Az = 0
    return np.array([Ax,Ay,Az])

def vecB(b,gamma):
    """
    位置ベクトルB
    """
    gamma = deg2rad(gamma)
    Bx = b*math.cos(gamma)
    By = b*math.sin(gamma)
    Bz = 0
    return np.array([Bx,By,Bz])

def vecC(c,alpha,beta,gamma):
    """
    位置ベクトルC
    """
    alpha,beta,gamma = deg2rad(alpha),deg2rad(beta),deg2rad(gamma)
    Cx = c*math.cos(beta)
    Cy = c*(math.cos(alpha)-math.cos(beta)*math.cos(gamma))
    Cz = c*math.sqrt((math.sin(beta))**2-(math.cos(alpha)-math.cos(beta)*math.cos(gamma))**2/(math.sin(gamma))**2)
    return np.array([Cx,Cy,Cz])

def unitCell(idx,a,b,c,alpha,beta,gamma):
    """
    idx=(n1,n2,n3)から頂点位置(x,y,z)を算出
    """
    Ax,Ay,Az = idx[0]*vecA(a)
    Bx,By,Bz = idx[1]*vecB(b,gamma)
    Cx,Cy,Cz = idx[2]*vecC(c,alpha,beta,gamma)
    x = Ax + Bx + Cx # x座標
    y = Ay + By + Cy # y座標
    z = Az + Bz + Cz # z座標
    return x,y,z

def vtxPos(a,b,c,alpha,beta,gamma):
    """
    単位格子の頂点の位置を計算
    """
    val = [0,1] # n1,n2,n3がとりうる値
    X = np.zeros(8)
    Y = np.zeros(8)
    Z = np.zeros(8)
    for i,idx in enumerate(itertools.product(val,repeat=3)):
        x,y,z = unitCell(idx,a,b,c,alpha,beta,gamma)
        X[i],Y[i],Z[i] = x,y,z
    return X,Y,Z

# SrTiO3の結晶パラメータ
a = 3.94513
b = 3.94513
c = 3.94513
alpha = 90
beta = 90
gamma = 90

# 頂点の計算
X,Y,Z = vtxPos(a,b,c,alpha,beta,gamma)

# 図示
fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111,projection='3d')
ax.scatter(X,Y,Z,s=40,c='b')
ax.set_xlabel('x',size=14)
ax.set_ylabel('y',size=14)
ax.set_zlabel('z',size=14)
plt.show()
# plt.savefig('SrTiO3_unitCell.png',transparent=True,dpi=300)

上記のプログラムを実行すると以下の図が得られます。

これを平行六面体としてプロットされるようにプログラムを書き直していきます。平行六面体を連結グラフ(頂点はノード、辺はエッジに対応)として見たとき、1つのノードに対して3つのエッジがあるグラフである事がわかります。つまり、オイラーグラフの定理より、オイラー・グラフ (Eulerian graph) ではない事が言えます(補足参照)。これより、"ax.plot"による"1筆書き"で平行六面体を描く事ができないため、泥臭く1つのグラフに辺を加えていくような処理で対処する必要があるかと思います(楽な方法をご存じの方はご教示ください)。つまり、直線の端点を2つ指定して、"ax.plot"で線を追加するという処理を12回繰り返して平行六面体を作成します。2頂点分の座標を保持するvtxMerge関数を定義し、これを"ax.plot"に入れていきます。

def vtxMerge(idx1,idx2):
    """
    idx1=(n11,n12,n13)とidx2=(n21,n22,n23)の頂点の座標を持つ配列を生成
    """
    X = np.zeros(2)
    Y = np.zeros(2)
    Z = np.zeros(2)
    X[0],Y[0],Z[0] = unitCell(idx1,a,b,c,alpha,beta,gamma)
    X[1],Y[1],Z[1] = unitCell(idx2,a,b,c,alpha,beta,gamma)
    return X,Y,Z

# 平行六面体の頂点
vtx1 = (0, 0, 0)
vtx2 = (0, 0, 1)
vtx3 = (0, 1, 0)
vtx4 = (0, 1, 1)
vtx5 = (1, 0, 0)
vtx6 = (1, 0, 1)
vtx7 = (1, 1, 0)
vtx8 = (1, 1, 1)

# 平行六面体の辺
X1,Y1,Z1 = vtxMerge(vtx1,vtx2)
X2,Y2,Z2 = vtxMerge(vtx1,vtx3)
X3,Y3,Z3 = vtxMerge(vtx1,vtx5)
X4,Y4,Z4 = vtxMerge(vtx4,vtx2)
X5,Y5,Z5 = vtxMerge(vtx4,vtx3)
X6,Y6,Z6 = vtxMerge(vtx4,vtx8)
X7,Y7,Z7 = vtxMerge(vtx6,vtx2)
X8,Y8,Z8 = vtxMerge(vtx6,vtx5)
X9,Y9,Z9 = vtxMerge(vtx6,vtx8)
X10,Y10,Z10 = vtxMerge(vtx7,vtx3)
X11,Y11,Z11 = vtxMerge(vtx7,vtx5)
X12,Y12,Z12 = vtxMerge(vtx7,vtx8)

# プロット
fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111,projection='3d')
ax.plot(X1,Y1,Z1,c='b')
ax.plot(X2,Y2,Z2,c='b')
ax.plot(X3,Y3,Z3,c='b')
ax.plot(X4,Y4,Z4,c='b')
ax.plot(X5,Y5,Z5,c='b')
ax.plot(X6,Y6,Z6,c='b')
ax.plot(X7,Y7,Z7,c='b')
ax.plot(X8,Y8,Z8,c='b')
ax.plot(X9,Y9,Z9,c='b')
ax.plot(X10,Y10,Z10,c='b')
ax.plot(X11,Y11,Z11,c='b')
ax.plot(X12,Y12,Z12,c='b')
ax.set_xlabel('x',size=14)
ax.set_ylabel('y',size=14)
ax.set_zlabel('z',size=14)
plt.show()
# plt.savefig('SrTiO3_unitCell.png',transparent=True,dpi=300)

上記のプログラムを実行すると以下の図が得られます。ここでは図示してもあまり面白くない立方晶を取り扱っていますが、色々パラメータをいじって単位格子を眺めてみてください。

補足

今回取り扱った平行六面体を多面体グラフとしてみて、一筆書きができないと判断したロジックについて補足します。8つの頂点に$a,b,c,d,e,f,g,h$とラベルを振り、この点(vertex)の集合を$V$と置きます。この点の集合$V$から2つの頂点$i,j$を組み合わせて$\bar{ij}$と記述したものを辺とします。このとき、頂点をどういった順番で配置したかで考え得る辺が変わりますが$\bar{ab},\bar{ac},\dots,\bar{hg}$の12個の辺がある事がわかります($\bar{ij}=\bar{ji}$で考えています。これは無向性グラフと呼ばれます)。この集合を辺(edge)の集合$E$とします。この$V$と$E$の2つの集合から構成されるものをグラフ$G$と呼び、$G=(V,E)$と記述されます。このグラフ$G$が以下の定理を満たすとき、連結グラフと呼ばれます。

グラフ$G$が連結グラフならば$|E|\geq |V|-1$である。

平行六面体は$|E|=12$、$|V|=8$ (集合を絶対値で囲む操作は集合の要素数を数える事を意味します)であるため、上記の定理を満たします。ある点を端点とする辺の本数を次数(degree)といい、$deg(\cdot)$という記号を使って表します。つまり、平行六面体の次数は以下のように記述されます。

deg(a)=3,\quad deg(b)=3\quad ,\dots,\quad deg(h)=3

この各点における次数が以下の定理を満たすとき、そのグラフ$G$はオイラー・グラフと呼ばれ、同じ辺を二度なぞらずに線図形を描く事ができる事を意味します。

オイラーグラフの定理：
連結グラフ$G$がオイラー・グラフとなる必要十分条件は$G$の点の次数が全て偶数である事

平行六面体の各頂点における次数は$3$であるため、オイラー・グラフではない事がわかります。"ax.plot"は点と点を結んで折れ線グラフなどを描く事ができるため、多角形なら各頂点の次数が$2$なので一筆書きできますが多面体となると限られたものにしか使えないようです。

引用文献

[1] 松下能孝, J. Surface Analysis, 21, 71 (2014).

Pythonで結晶を描画するプログラムを作るシリーズ

• [Pythonで結晶を描画するプログラムを作る~分極座標とアフィン変換~] (https://qiita.com/PHVTuber/items/40bb5b151c4404de7dab)