へっぽこ博士のなんちゃって研究室

PythonやR、分子軌道計算等に関する記事を書きます

格子定数と格子ベクトルの相互変換

- 事始め


Quantum ESPRESSO[*1]の入力ファイルに格子定数を入力するには、以下の2通りの宣言方法があります[*2]。

・SYSTEMカード内にibravとcelldmを記述
・CELL_PARAMETERSカード内に格子ベクトルを記述

ibravの説明では、どちらもBravais格子に合わせて計算する必要があるように見えますが、後者の方法では三斜晶系の方法を覚えるだけで良いので便利です(その他は三斜晶系の特殊例にすぎません)。
また、計算途中の格子定数を求めたい場合、出力ファイルにはCELL_PARAMETERSしか出力されないため、格子ベクトルから変換する必要があります。
そこで、ここではそれらの相互変換方法を紹介します。

- 格子定数と格子ベクトルの関係


Quantum ESPRESSOにおける格子ベクトルは、a軸を特殊軸とする下記の仕様です[*3]。

  
\begin{bmatrix}  
 a_x & 0 & 0 \\
 b_x & b_y & 0 \\
 c_x & c_y & c_z \\
\end{bmatrix}

格子定数をa, b, c, α, β, γとしたとき、格子ベクトルの各要素は下記式でそれぞれ求めることができます。

  
\begin{aligned}  
 a_x &= a \\
 b_x &= b\cos \gamma \\
 b_y &= b\sin \gamma \\
 c_x &= c\cos \beta \\
 c_y &= c\frac{(\cos \alpha - \cos \beta \cos \gamma)}{\sin \gamma} \\
 c_z &= \sqrt{c^2 - c_x^2 - c_y^2} \\
\end{aligned}


Pythonのコード例

格子定数は、カフェインとテオフィリンの共結晶(COD:7108131)[*4]から引用しました。

import numpy as np

# cell = [ a, b, c, α, β, γ ]
cell = np.array([7.0081, 8.7904, 13.0856, 95.648, 97.09, 91.667])
V = np.zeros((3,3))

cos_a = np.cos(np.deg2rad(cell[3]))
cos_b = np.cos(np.deg2rad(cell[4]))
cos_c = np.cos(np.deg2rad(cell[5]))
sin_c = np.sin(np.deg2rad(cell[5]))

V[0][0] = cell[0]
V[1][0] = cell[1] * cos_c
V[1][1] = cell[1] * sin_c
V[2][0] = cell[2] * cos_b
V[2][1] = cell[2] * (cos_a - cos_b * cos_c) / sin_c
V[2][2] = np.sqrt(cell[2] ** 2 - V[2][0] ** 2 - V[2][1] ** 2)


sin γで除算するので、本来ならZeroDivisionErrorをチェックしたほうが良いですが、そうなるのはγが0°、または180°の時だけで、そのような結晶は存在しえないため、わざわざ記述する必要はないように思います。


CELL_PARAMETERSを使用する場合は、SYSTEMカード内でibrav=0を指定します[*5]。

入力ファイルの一例

&SYSTEM
 ibrav = 0
/

CELL_PARAMETERS {angstrom}
  7.008100  0.000000  0.000000
 -0.255717  8.786680  0.000000
 -1.615133 -1.335391 12.916695


一方、格子ベクトルから格子定数へは下記式で変換することができます。


\begin{aligned}
a &= a_x \\
b &= \sqrt{b_x^2 + b_y^2} \\
c &= \sqrt{c_x^2 + c_y^2 + c_z^2} \\

\alpha &= \cos^{-1} \left( \frac{\overrightarrow{b} \cdot \overrightarrow{c}}{bc} \right) \\
\beta &= \cos^{-1} \left( \frac{\overrightarrow{c} \cdot \overrightarrow{a}}{ca} \right) \\
\gamma &= \cos^{-1} \left( \frac{\overrightarrow{a} \cdot \overrightarrow{b}}{ab} \right) \\
\end{aligned}


Pythonのコード例2

numpy:ndarrayでは、@で内積を計算できます[*6]。

import numpy as np

V = np.array([[1.0, 0.0, 0.0],[1.0, 1.0, 0.0],[1.0, 1.0, 1.0]])

a = V[0][0]
b = np.sqrt(V[1] @ V[1])
c = np.sqrt(V[2] @ V[2])

cos_a = (V[1] @ V[2]) / (b * c)
cos_b = (V[2] @ V[0]) / (c * a)
cos_c = (V[0] @ V[1]) / (a * b)

alpha = np.rad2deg(np.arccos(cos_a))
beta = np.rad2deg(np.arccos(cos_b))
gamma = np.rad2deg(np.arccos(cos_c))


因みに、セル体積は下記式で格子ベクトルから簡単に求めることができます。


V = a_x \times b_y \times c_z


*1:https://www.quantum-espresso.org/

*2:https://www.quantum-espresso.org/Doc/INPUT_PW.html#idm199

*3:https://www.quantum-espresso.org/Doc/INPUT_PW.html#idm1361
この確認を怠ったために、c軸を特殊軸として計算していまい、結果が一致しないと数日悩むことになったのは内緒です。

*4:Crystallography Open Database (COD)
http://www.crystallography.net/cod/index.php

*5:cif2cellの出力を見るとAも指定するのが作法のようですが、なくても問題ありません。

*6:Python 3.5、Numpy 1.10以降で採用された仕様です。