Go to Qiita Advent Calendar 2024 Top
LoginSignup
3
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

@PHVTuber

Pythonで結晶を描画するプログラムを作る~分極座標とアフィン変換#2~

Last updated at Posted at 2021-09-28

投稿日：2021/09/28
更新日：2021/10/04 (反転操作の改良)
更新予定内容：反転、らせん、映進操作のプログラムの改良($n$回、$m$回施すことで変化する点の軌跡の可視化)

はじめに

前回は対象操作と平行移動(並進操作)から構成されたアフィン変換によるアフィン空間(分極座標の3次元$+$ダミー次元)上の操作の可視化について学びました。本記事では前記事で取り扱わなかった反転、鏡映、回映、回反操作の可視化を目指します。これらの操作は対象操作$A$の部分のみを変更するだけで良いため、プログラムの変更がない限り前記事のものを流用します。また、ここらの話は群論のビジュアル的な理解が進むかと思うので色々パラメータをいじって遊んでいただけたらと思います。

注意
下記で使う記号における$C_{n},S_{n},C_{\bar{n}},i$や用語は一般的に使われます。しかし、任意回転軸、任意面、任意原点などの操作やその数式は群というものを直感的に理解してもらう、変数の扱いに専門知識が必要ではない等の理由で独自に構成したものです。細心の注意を払って数式などを作成しておりますが批判的に読んで頂けたらと思います。その関係上、記事はアップデートありきで作成しております。

各操作のプログラムの共通部分

基本的にはPythonで結晶を描画するプログラムを作る~分極座標とアフィン変換~で取り扱った関数を利用します。この節以降から対称操作とその行列表現について説明していきますが、前述したように操作$A$以外は同じ処理になりますのでseitzMtx関数の引数mtxに代入するための関数を作っていきます。また、様々な対称操作を扱うため、Mysimulation3という関数を新たに導入します。

def Mysimulation3(x_bef,seitz,size=40,fontsize=15,fig_size1=8,fig_size2=8,color='r',arrow_length_ratio=0.1,lim=8,fileName='test',save=False):
    """
    pos_bftからpos_aftに向かう矢印、x,y,z軸の描画
    ------------------------------------------------
    x_bef：list
         点の初期位置
    seitz：array
         ザイツの拡大行列
    """
    # ●計算
    x_aft = movePos(seitz,x_bef)[:-1]

    # ●描画
    fig = plt.figure(figsize=(fig_size1,fig_size2))
    ax = fig.add_subplot(111,projection='3d')
    # 移動前と移動後の点のプロットおよびその間に矢印を引く
    x,y,z = x_bef
    u,v,w = x_aft
    ax.quiver(x,y,z,u-x,v-y,w-z,color=color,arrow_length_ratio=arrow_length_ratio)
    ax.scatter(x,y,z,color='b',s=size)
    ax.text(x,y,z,s=r'$x_{bft}$',fontsize=fontsize)
    ax.scatter(u,v,w,color='r',s=size)
    ax.text(u,v,w,s=r'$x_{aft}$',fontsize=fontsize)
    # x,y,z軸の追加
    ax.quiver(-lim,0,0,2*lim,0,0,linestyle='-', arrow_length_ratio=0.1,color='black') # x軸
    ax.quiver(0,-lim,0,0,2*lim,0,linestyle='-', arrow_length_ratio=0.1,color='black') # y軸
    ax.quiver(0,0,-lim,0,0,2*lim,linestyle='-', arrow_length_ratio=0.1,color='black') # z軸
    # 出力範囲
    ax.set_xlim(-lim,lim)
    ax.set_ylim(-lim,lim)
    ax.set_zlim(-lim,lim)
    # 軸ラベル
    ax.set_xlabel('x',size=14)
    ax.set_ylabel('y',size=14)
    ax.set_zlabel('z',size=14)
    if save == True:
        plt.savefig(f'{fileName}.png',transparent=True,dpi=300)
    return x_aft

# 拡大行列の計算
x_bef = [-0.5,-0.5,-1]
vecB = [0,0,0] # b=(0,0,0)'の移動
mtx = # ここに作成した行列を入れ込みます
seitz = seitzMtx(mtx,vecB)

_ = Mysimulation3(x_bef,seitz)

1. 反転操作

反転操作とは点$\boldsymbol{x}=(x,y,z)'$を$\boldsymbol{x'}=-\boldsymbol{x}$に移す操作となります。すなわち、

\boldsymbol{x'}=i\boldsymbol{x}

となるような操作$i$が反転操作となります。従って、反転操作は以下の行列で記述されます。

i=\left(\begin{array}{ccc}
-1 & 0 & 0\\
0 &-1 & 0\\
0 & 0 & -1
\end{array}\right) = -I_{3}\tag{1}

この操作は原点$(0,0,0)'$を基準とした反転操作ですが、任意位置ベクトル$\boldsymbol{p}$を基準としたとき、

\boldsymbol{x'}=-\boldsymbol{x}+2\boldsymbol{p} \tag{1'}

と記述されます[1]。この形をよくみると、平行移動と線形変換の2つから構成されている事がわかります。これは前記事で扱ったようにザイツの記法で1つの行列にまとめられます。$(1')$式をザイツの記法に落とし込めば下記が得られます。

i_{p} = (-I_{3}|2\boldsymbol{p}) \tag{1''}

$\boldsymbol{p}=(0,0,0)'$のとき、$(1'')$は$(1)$に一致する事が確認できるかと思います。

余談ですが、2次の非線形光学効果はこの反転操作を元に持たない系でしか発現しません。そのため、非線形光学活性を有する有機化合物を反転対称心を持たない構造に結晶化させるという試みがあります。非線形光学効果についてはネットで公開されている『非線形光学デバイスの物理』を参考にして下さい。

def Invarse():
    """
    反転操作
    """
    return -np.eye(3)+0

mtx = Invarse()

上記の関数から生成した行列を用いて、Mysimulation3関数を実行すると以下が得られます。

test1.png

2. 鏡映操作

鏡映操作とは任意の平面に対して面対称な位置に移す操作です。例えば、2次元平面においては線対称と同義ですが、$x$軸に対しての鏡映操作とは点$\boldsymbol{x}=(x,y)'$を$\boldsymbol{x'}=(x,-y)'$とするような操作になります。この操作を担う行列の作成には以下の公式を利用します。

平面$\Pi$に対する単位法線ベクトルを$\boldsymbol{v}=(v_1,v_2,v_3)'$とおく。このとき、この面に対する鏡映操作$\Sigma_{v}$は

\Sigma_{v} = \left(
\begin{array}{ccc}
1-2v_1^2 & -2v_1 v_2 & -2v_1 v_3\\
-2v_1 v_2 & 1-2v_2^2 & -2v_2 v_3\\
-2v_1 v_3 & -2v_2 v_3 & 1-2v_3^2
\end{array}
\right) \tag{2}

である。

この行列の成分に注目すると$\Sigma_{ij}=\Sigma_{ji}$という性質を持っている事がわかります。これは専門には対称行列と呼ばれ、$N\times N$行列を考えた時独立した要素が$\frac{N(N+1)}{2}$個ある事を示します。つまり、今回に限れば6個の要素の計算のみで十分である事がわかります。

def Mirror(vecV,digit='0.01'):
    """
    鏡映操作
    ------------------------------
    vecV : array-like
        原点を通る平面の単位法線ベクトル
    """
    v1,v2,v3 = vecV
    
    # 回転行列の成分
    S11 = 1-2*v1**2
    S12 = -2*v1*v2
    S13 = -2*v1*v3
    S21 = S12
    S22 = 1-2*v2**2
    S23 = -2*v2*v3
    S31 = S13
    S32 = S23
    S33 = 1-2*v3**2

    return np.array([Round([S11,S12,S13],digit)
                     ,Round([S21,S22,S23],digit)
                     ,Round([S31,S32,S33],digit)])+0

vecV = [1,0,0] # yz平面の法線ベクトル
mtx = Mirror(vecV,digit='0.01')

上記の関数から生成した行列を用いて、Mysimulation3関数を実行すると以下が得られます。

test2.png

3. 回映操作

回映とは回転操作$C_{n}$の後にその軸に対する原点を通る垂直な平面に対する鏡映操作$\Sigma$を施す操作です。回転軸単位ベクトル$\boldsymbol{p}=(p_{1},p_{2},p_{3})'$に対する回映操作$S_{n}(\boldsymbol{p}) (n=\frac{2\pi}{\theta})$と置いたとき、

S_{n}(\boldsymbol{p})=\Sigma_{p}C_{n}(\boldsymbol{p})\tag{3}

と記述されます。ここで、$C_{n}(\boldsymbol{p})$は回転軸ベクトル$\boldsymbol{p}$に対する$n$回回転を表します。また、ここで採用している定義では$\Sigma$の添え字は鏡映操作を施す平面の法線ベクトルであるため、回転軸単位ベクトル$\boldsymbol{p}$になる事に留意する必要があります。対称性を議論するときに$n$が整数だと嬉しいのですが、遊び辛いので任意回転角用の回映操作に書き直します。

S(\theta,\boldsymbol{p})=\Sigma_{p}R_{p}(\theta)\tag{3'}

ここで、$S(\theta|\boldsymbol{p})$は回転軸単位ベクトル$\boldsymbol{p}$に対する$\theta [rad]$回転後に、その回転軸に対して垂直な平面に対する鏡映操作を施す操作、$R_{p}(\theta)$は前記事で取り扱ったロドリゲスの回転公式です。余談ですが、$S_{n}$軸を持たない分子はキラルである事と同義です。キラルなフェリ磁性体は時間反転対称性と空間反転対称性が破れており、面白い磁気光学効果を示す事が報告されています。

def impRot(vecP,theta,digit='0.01'):
    """
    回映操作. (3')式
    ------------------------------
    vecP : array-like
        回転軸単位ベクトル
    theta：float
        回転角度
    """
    impRot = np.dot(Mirror(vecP,digit='0.01'),
                    Rodrigues(vecP,theta,digit='0.01'))

    return impRot+0

x_bef = [-1,-1,0] # xy平面上の点
vecP = [0,1,0] # y軸回転
theta = 45.0 # 45度回転
mtx = impRot(vecP,theta)

上記の関数から生成した行列を用いて、Mysimulation3関数を実行すると以下が得られます。

test3.png

4. 回反操作

回反操作とは回転後に反転操作を施す操作です。回転軸ベクトル$\boldsymbol{p}$に対する$n$回回反操作を$C_{\bar{n}}$とおくと、

C_{\bar{n}}(\boldsymbol{p}) = iC_{n}(\boldsymbol{p})\tag{4}

と記述されます。$\boldsymbol{p}$に対して$\theta [rad]$回転させた後に反転させる操作を$C(\bar{\theta}|\boldsymbol{p})$とおいたとき、以下のように記述されます。

C(\bar{\theta},\boldsymbol{p}) = iR_{p}(\theta)\tag{4'}

def rotInv(vecP,theta,digit='0.01'):
    """
    回反操作. (4')式
    ------------------------------
    vecP : array-like
        回転軸単位ベクトル
    theta：float
        回転角度
    """
    rotInv = np.dot(invarse(),
                    Rodrigues(vecP,theta,digit='0.01'))

    return rotInv+0

x_bef = [-1,-1,0] # xy平面上の点
vecP = [0,1,0] # y軸回転
theta = 45.0 # 45度回転
mtx = rotInv(vecP,theta)

5. 映進操作

ある面に対して鏡映操作を施し、その後、その面に平行なある方向に$\frac{1}{2}$(軸映進面、二重映進面、対角映進面の場合)または$\frac{1}{4}$(ダイアモンド映進面の場合)並進する操作です。ここでいう映進面とはこの対称操作を定義する面、すなわち鏡映面になります。鏡映面に対して平行な軸というのは鏡映面に対する法線ベクトルと直交したベクトルの方向であるため、結晶の性質を無視すれば無数にとれます。結晶に限れば、$a$,$b$,$c$軸方向、面に平行な直交した2つのベクトルを合成した方向、単位格子の対角線方向などの限られた方向で意味を持ちます。映進面の分類は神戸大学の講義資料『対称軸、対称⾯の⽅向と操作』が詳しいです。これらを踏まえてこのアフィン変換を作っていきます。映進操作とは対称操作$A$と並進操作$\boldsymbol{b}$から構成されるため、ザイツの記法を使います。鏡映面の単位法線ベクトルを$\boldsymbol{v}$、鏡映面に平行なベクトルを$\boldsymbol{p}$と置き、$\frac{1}{m}$並進するとき、映進操作$g(\boldsymbol{v},\boldsymbol{p},m)$は以下で与えられます。

g(\boldsymbol{v},\boldsymbol{p},m) = \left(\Sigma_{v}|\frac{\boldsymbol{p}}{m}\right) \tag{5}

$\boldsymbol{p}=\boldsymbol{a}$,$\boldsymbol{b}$,$\boldsymbol{c}$,$\boldsymbol{a}+\boldsymbol{b}$,$\dots$などの単位格子の頂点位置ベクトルを考えたとき、$m$回同じ操作を施すと初期位置と等価な同一点に移動します。コードは既存のもので作成できるため目新しいものがありませんが、サンプルを下記です。vecPとvecVが直交でさえあれば映進操作なのでこの関係が保たれていれば問題はありません。

def orthoCheck(vecA,vecB,digit="0.01"):
    """
    ベクトルが直交しているかチェック. 内積0のとき直交
    --------------------------------
    vecA,vecB：array-like
        要素数が同じベクトル
    """
    inPro = Round(np.dot(vecA,vecB),digit)
    if inPro==0:
        print("Yes")
    else:
        print("No")
    
vecP = [1,0,0] # a軸
vecV = [0,1,0] # b軸
m = 2 # 並進
mtx = Mirror(vecV) # 鏡映操作
seitz = seitzMtx(mtx,np.array(vecP)/m)

6. らせん操作

らせん操作$n_{m}(\boldsymbol{p})$とは回転軸ベクトル$\boldsymbol{p}$に対する回転操作$C_{n}(\boldsymbol{p})$を施した後、その回転軸ベクトル$\boldsymbol{p}$の方向に$\frac{m}{n}$並進する操作です。イメージは螺旋階段です。結晶の周期性と整合がとれる回転は前述したように$n=2,3,4,6$のみで、そのうえで映進操作は$1\leq m<n$の整数値のみが許容となります。この操作はザイツの記法を用いると下記です。

n_{m}(\boldsymbol{p}) = \left(C_{n}(\boldsymbol{p})|\frac{m\boldsymbol{p}}{n}\right)

プログラムは既存のものを組み合わせるだけです。ここでは、$3_{1}(\boldsymbol{a})$を描いてみます。$n=3$のとき、$1\leq m <3$であるため、$3_{2}(\boldsymbol{a})$も許容となります。$n=3$とは、$\frac{2\pi}{3}=120[deg]$回転を表します。

# 3_1(a)の操作
vecP = [1,0,0] # a軸
n = 3
m = 1 
theta = 120 # 120度回転
mtx = Rodrigues(vecP,theta)
seitz = seitzMtx(mtx,np.array(vecP)*(m/n))

7. ここまでのお話の整理

結晶の単位格子の結晶パラメータ$a,b,c,\alpha,\beta,\gamma$に注目すると7つのタイプに分類できます。この区分をより詳細に行うため、実格子座標系(通常の直交座標系)から分極座標系を考えました。この座標系に変換する事で原子の相対位置がわかりやすくなり、また、ほかの結晶でも同様の変換が意味を持つ座標系という利点があります。さらにこの座標系にダミー次元を加えた空間をアフィン空間と呼び、この空間上の変換はアフィン変換と呼ばれます。このアフィン変換を考える事は結晶の分類を目的としたとき重要な役割を果たします。通常、アフィン変換とは線形変換と平行移動から構成されるため無数に考える事ができますが、アフィン変換後の移動点もまた結晶中の原子であるような変換に限ると1から十数個程度に限られます。つまり、ゼロから結晶を分類していくとなると、許容されるアフィン変換を探すという行為を最低230種類の結晶に対して施す必要があります。分極座標系上の点$x$に対してアフィン変換を施し移動した点$x'$は等価点と呼ばれ、この等価点はアフィン変換と1対1対応するため、結晶の空間群を規定するためには等価点のリストを列挙すれば良いことになります。

遊び方の指針

上記のコードを書いて、描くだけで満足してしまう方が大半かと思いますが大変遊び甲斐のあるプログラムです。アフィン変換の定義は分野によって若干の定義の違いはあるものの本質に違いはなく、群論、結晶学、コンピュータグラフィックスなどの多くの分野で取り扱われる基本的な操作です。他分野への応用に関しては他の記事を参考にして貰うとして結晶学と群論における話題を提示しておきます。これらを解明できるようパラメータをいじって遊んでみて下さい。

  • 結晶において回転は$n=2,3,4,6$のみ許容がされる
  • 結晶においてらせん操作は$1\leq m<n$のみが許容される
  • 正十二面体を持つ結晶はない

参考文献

[1] https://en.wikipedia.org/wiki/Point_reflection

3
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
3
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?