2024 年 23 巻 1 号 p. A2-A8
Generation of crystal models is an unavoidable part of first-principles calculations and other computational simulations at the atomic level. The author developed a versatile single-file code that can make various models, such as crystallographic conventional cells, supercells, and slab-and-vacuum models, as a one-liner command. Computer-assisted generation of steps-and-terraces and grain boundary models is also possible.
結晶モデルの作成は，第一原理計算を始めとする原子レベルの計算シミュレーションでは避けて通れない．結晶モデル作成を支援するソフトウェアシステムとして，例えばPymatgen [1]やAtomic Simulation Environment (ASE) [2]がよく知られている．これらは多数の人が作り上げたシステムであり，多機能である．しかし，Pythonに慣れていないと利用がとても困難で，マニュアルが英語で書かれていることもありインストールにすら四苦八苦するケースも多い．他方，Dassault Systèmes 社のMaterials Studio，クロスアビリティ社のWinmostar，アドバンスソフト社のAdvance/Nanolaboのような，結晶構造を見ながら，選択肢を選びつつモデルを作成し，計算シミュレーションまで行うことができるGraphical user interface (GUI)システムも存在するが，システムの機能に実装されていない機能との統合が難しい．このため，スクリプトを用いた処理が手軽に行える，できれば国産の代替ソフトウェアが望まれる．
著者は結晶モデル作成支援ツール「壺車(つぼぐるま)」を作成した．プログラムはJ-STAGE DATAにてCC-BYライセンスにて公開されている．名前の由来は第一原理計算ソフトウェアVASP [3, 4]の主要ファイルである「POTCAR」を無理やり和訳したものである．プログラム本体はPerlで書かれた1つの実行ファイルのみでインストールは簡単，動作は原則としてUNIXコマンドラインでの1行入力，日本国民への研究結果の還元および国内産業振興の一助となるためマニュアルは日本語，と使いやすさを重視している．「壺車」は著者が単独で作成・管理しているため，問題が発覚した時の対処が比較的容易である．また，新機能を追加する際，著者はソフトウェア設計の思想，既存コードの書き方の癖，各サブルーチンの挙動を熟知しているため，迅速に実装できる．結晶の対称性に関する情報はphonopy [5]プログラムをコマンドラインで実行して結果を参照するため，別途phonopyのインストールが必要となる．入出力周りについてはユーザーにより好みが異なるであろうが，著者の研究スタイルにとって最適な方式を採用している．このため，結晶モデルの入力フォーマットはVASPのPOSCAR形式にのみ対応している．PymatgenやASE等の他のシステムとの直接の互換性はなく(コマンドラインから呼び出すことで実行することは可能である)，著者自身は今後も互換性に関する機能の追加および英語版マニュアルの作成を行う予定はない．なお，自分の手間をかけてこれらを行いたい者がいるのであれば，著者はあえて妨げるつもりはない．
本ソフトウェアに実装された機能のうち主要なものをTable 1に記す．表中に説明に〇がついているものはphonopyが必須となるコマンドである．また，この節では代表的な機能と実行方法を解説する．
| Classification | Command | Descriptions |
| 他形式出力 | --cif6 | cif形式に変換 |
| 格子情報 | --data | 格子定数等の情報表示 |
| ユニットセル作成 | --cc | 〇慣用単位胞の作成 (Make crystallographic conventional cell) |
| --cp | 〇結晶学的基本単位胞の作成 (Make crystallographic primitive cell) | |
| --niggli | 基底ベクトルをNiggli既約にする (Niggli-reduce basis vectors) | |
| --tidy | 基底ベクトルを三角行列にする | |
| 格子(スラブモデル向け) | --basearea | ab面の面積を表示 |
| --height | スラブ高さを表示 | |
| --maxortho | c軸をできるだけab面に垂直に取る | |
| 格子(発展) | --rotate | 指定軸まわりに指定角度回転する |
| 組成(情報) | --numatoms | 原子数表示.元素版もあり |
| --stoichiometry | 組成表示 | |
| --volperatom | 原子当たり体積表示 | |
| 組成(原子位置変更なし) | -d | Fractional (direct) 座標にする |
| -c | Cartesian座標にする | |
| --selective_true | Selective dynamicsで全部Trueにする.原子範囲指定に加え, Falseにする機能もあり | |
| 組成(原子位置変更あり) | --in_cell | 全座標の値をfractional (direct)で0と1の間にする |
| --shift | 原子位置を指定変位分シフト | |
| 原子追加削除 | --rmatom | 指定原子(複数可)削除.元素版もあり |
| --addatom | 指定原子追加 | |
| --swapatom | 指定2原子の位置を交換.元素版もあり | |
| 原子間距離 | --dist | 各原子に対し,一定距離以下の原子と距離を表示.原子範囲を指定する機能もあり |
| --dist_min | 最小結合距離を表示 | |
| --coord | 各原子に対し配位数を表示 | |
| スーパーセル | --supercell | スーパーセルを作成(3要素指定,9要素変換行列指定の両方可) |
| スーパーセル(発展) | --N_supercell | 〇できるだけ直方体のスーパーセルを探す(Find maximally orthogonalized supercells) |
| --3dalike | 〇2つの格子のスーパーセルを比較し,そっくりさを表示 | |
| --2D_rectangular | 〇長方格子に近いN倍の2次元スーパーセルを探す.正方・六方格子版もあり | |
| K点(一般) | --kpoints | 指定メッシュ間隔のKPOINTSファイルを作成 |
| --kpoints_slab | 指定メッシュ間隔のスラブ計算用KPOINTSファイルを作成 | |
| K点(バンド図) | --kpath | 〇バンドパス用k点の位置情報を表示 |
| --kpath_mass | 〇有効質量計算用k点の位置情報を表示 | |
| NEB | --nebimages | 指定数のNEBのイメージPOSCARを作成 |
| --superimpose | 複数のPOSCARの原子位置を1つのPOSCARに重ねる | |
| 対称性(基本) | --sgnumber | 〇空間群番号を表示.Hermann–Mauguin空間群記号,点群記号,ブラべ格子版もあり |
| 対称操作(isometry) | --isometry | 〇結晶を自分自身に写す対称操作(isometry)の表示 |
| --isometry_nonpolar | 〇無極性を担保するisometryの表示 | |
| --surface_isometry | 〇指定面指数のスラブが無極性になれるかを表示 | |
| --equivalent_site | 〇指定原子と等価な原子を表示 | |
| 表面モデル(基本) | --change_slab_vacuum_thickness_ortho | スラブの真空厚さを指定値に,c軸をab面に垂直に変更 |
| --make_slab | z座標が特定の範囲の原子のみ残す | |
| 表面モデル(発展) | --hkl_primitive | 〇慣用単位胞の(hkl)面がab面となるような基本単位胞を作成 (Find primitive cells where the (hkl) plane of the conventional cell is the ab-plane) |
| --slab_poscar | 〇指定スラブ厚さ以上,指定真空厚さ以上の(hkl)面の無極性スラブモデルを作成 | |
| --unique_nonpolar | 〇無極性になりうる面指数を表示(点群で判断,等価な面指数については代表1つのみ) | |
| --termination_polarity_list | 指定した面指数範囲で無極性なスラブが取れるかを表示 | |
| 表面モデル(ステップモデル) | (複数機能) | 〇Phonopy必須 |
| 表面モデル(ファセットモデル) | (複数機能) | 〇Phonopy必須 |
| 表面原子判定 | --solid_angle | 提供したスラブモデルの表面原子判定(qconvex使用) |
| --extract_surface | 提供したスラブモデルの表面原子のみ抽出(qconvex使用) | |
| 表面エネルギー記述子 | --proximity | 〇指定(hkl)面に対し,無極性終端を自動検出しAtom proximity functionを算出 |
| --roughness | 提供したスラブモデルのSurface roughness indexを算出(qconvex使用) | |
| 粒界モデル作成 | (複数機能) | 〇Phonopy必須 |
note: The JCCJ Office considered Table 1 to be equivalent to the text and so allowed Japanese descriptions.
結晶構造の表現の基本は結晶学的慣習に基づいた慣用単位胞(crystallographic conventional cell)であり，定義がInternal Tables of Crystallography Volume A 6th edition (ITA) [6]の3.1.1.4節に記載されている．結晶の対称性はHermann–Mauguin記号で表現できるが，一部の結晶については，同じHermann–Mauguin記号の対称性を複数の格子定数で表現できることがある．一例は空間群Pmmmで基底ベクトルの長さが3,4,5Åの結晶で，6通りの表現がある．(なお，原子位置については，使用するWyckoff 位置(position)の記号について優先順位を設けたとしても，空間群P1等，無限の取り方が存在するケースがあるのでここでは議論しない．)各空間群番号の代表的なHermann–Mauguin記号に対し，事前に決められた格子定数の大小関係を要求することで，一意に慣用単位胞の格子定数を選ぶことができる．制約が必要な単斜晶の空間群はP2, P21, Pm, P2/m,P21/mの5つであり [7]，これらはITAの表3.5.2.3において”a=c かつ 90°<β<120°”がEuclidean正規化群 (normalizer)の対称性をenhanceするものである．直方晶ではITAの表2.1.3.4のnumber of distinct projectionsの数に応じて制限が決まる．数が1,2,3,6の場合，それぞれa<b<c，aが最短，a<b，制限なし，となる．なお，ITAの表3.5.2.3のEuclidean 正規化群の，最も制限が厳しいcell metricと最も緩い(generalな)cell metricの「Index of G in NE(G) or NE+(G)」の比を使おうとすると，空間群Ibcaで例外処理が必要となる． [8]他の結晶系では制約は不要である．
この一通りの表現の慣用単位胞を
tb2.pl --cc <POSCARで生成することができる．なお，PhonopyではITAで定まる慣用単位胞のうち一つを出力するが，複数の表現がある場合はどれが出力されるかは決まっていない．また，文献 [8]では結晶学的基本単位胞(crystallographic primitive cell)が定義されており，
tb2.pl --cp <POSCARで得られる．この結晶学的基本単位胞は必ずしも既約(reduced)，すなわち基底ベクトルが最も短くなるよう取られていない．既約な基本単位胞が欲しい場合は任意の基本単位胞を
tb2.pl --niggli <POSCARを用いてさらにNiggli 既約 [9, 10]にする．この--niggliコマンドは，与えられたPOSCARが実際に基本単位胞かどうかは気にせず，単に与えられた基底ベクトルをNiggli 既約にする．
通常の計算のKPOINTSファイルは
tb2.pl --kpoints 0.3 <POSCAR>KPOINTSで生成できる．この場合，k-meshは逆格子の基底ベクトルに沿って0.3×2π Å-1間隔で求められ，これはバンドギャップの開いた物質の通常の計算で一般的な値である(バンドギャップのない物質では間隔を短くする)．
逆格子空間でのバンド分散を示すバンド図は，特に半導体分野で重要となる．バンド図を描く際のバンドパスは，一般的に逆格子空間の高対称点を繋いだものである．高対称点は結晶の形と対称性で決まる．高対称点の多くには，結晶学の慣習に基づいた記号が与えられている． [11, 12] 文献 [8]に任意の結晶における高対称点と推奨バンドパスが記載されている．これらはSeeK-pathプログラム [13]から入手できるが，推奨バンドパスは
tb2.pl --kpath <POSCARでも入手できる．逆格子空間の高対称点の座標は実空間のユニットセルの取り方に依存する．上述のコマンドでは標準出力に結晶学的基本単位胞が生成される．併せてKPOINTS.Ckpath_0.025 ファイルに推奨バンドパス上のk点と高対称点の記号が記述される．バンドパスに沿ってサンプルする点の間隔はデフォルトで約0.025×2π Å-1 であるが，この間隔および生成ファイル名は自由に変更できる．
表面特性を計算するときは，二次元方向に無限に広がるスラブを真空層で挟むことで，三次元周期境界条件を満たしたスラブモデルを使用することが一般的である．Tasker [14]によると，スラブの極性は，全ての原子の層の平均電荷が0であるType 1，Type 1でないがスラブが「厚さ方向に双極子モーメントがない繰り返しユニット(repeat unit)」で構成されているType 2，そしてこのような繰り返しユニットが存在しないType 3に分類される．Type 1と2が無極性(nonpolar) スラブで，Type 3が極性(polar)スラブである．無極性スラブではスラブの表面と裏面が等価となるため，表面エネルギーを求めることができる．また，真空層中に静電ポテンシャルが厚み方向に一定値となる領域が生じるため，イオン化ポテンシャルや電子親和力が適切に定義できる．一方，著者ら [15]は全ての層が化学量論的に正しい(stoichiometricな) nonpolar type A, 厚さ方向に双極子モーメントがない繰り返しユニットの境界が原子の層の間にあるnonpolar type Bと層上にあるnonpolar type C，そしてそれ以外の極性スラブのpolarに分類した．岩塩型構造の(111)スラブはTaskerの定義ではpolarだが，著者らの定義では無極性と判定され，実際に表面エネルギーを得ることが可能である．
著者らはまず，nonpolar AとBスラブについて，バルクをへき開(cleave)し化学量論的に正しいスラブモデルが作れるのであれば自動生成する手法を開発した [15]．無極性であることは，スラブが特定の形の「結晶を自分自身に写す対称要素」(isometry)を有せば保証できる．また，結晶の特定の面指数において無極性スラブが作れるか否かは，その結晶の点群のみで決まる．ITAの表3.2.3.2にcrystal
formの一覧がある．同じcrystal
formに含まれる面指数は全て等価であるため，(hkl)と(
へき開の後に再構成した表面モデルとしては，スラブ両側の最表面の原子を半分ずつ削って化学量論的正しさと無極性を両立させたnonpolar type Cスラブの自動生成が可能である [18]．この再構成が必ずしも安定であるとは限らないが，選択肢として提示される．他にも，ステップモデル(steps-and-terraces model)の計算機支援型生成手法が開発されている [18]．
面指数とスラブと真空層の厚さが決まっている場合，例えば
tb2.pl --slab_poscar 20 15 2 3 -1 <POSCARでは(23
無極性スラブモデルの大量生成ワークフロー例を以下に記す．
tb2.pl --unique_nonpolar 4 2 3 6 <POSCAR>tempfile1 tb2.pl --termination_polarity_list tempfile2 <POSCAR>tempfile3 egrep ‘A|B’ tempfile2 | head -3 >tempfile3 tb2.pl --slab_poscar_list tempfile3 20 15 <POSCAR1番目のコマンドで，|h|≤4, |k|≤2, |l|≤3を満たす面指数(hkl)に対し，nonpolarスラブを実現可能なuniqueな面指数を出力する．ここではスラブモデルのスラブ表面に平行な面内面積が，最小値の6倍以下の物のみ出力する設定にする．2番目のコマンドで，化学量論的に正しい，かつ無極性type A,B,Cのモデルが取れるかの判定が出力される．何も出力されない面指数があるかもしれない．3番目のコマンドで，得られるスラブモデルがtype Cのみとなる面指数を排除する．もちろん，type Cのみが出る面指数を残しても構わない(この場合はegrep ‘A|B|C’もしくはヘッダを削除するtail -n +2となる)．面指数は面積順に並んでいるが，ここではhead -3を入れることで上から3つのみ出力する．どれだけの面指数を評価するかは結晶の対称性と構造，そして計算資源の余裕で決まる．一般論として，対称性が低い結晶ほど同じ倍率であっても出力される面指数の数が多くなり，多くの面を調べることになる．最後に，評価対象の面指数について，スラブ厚さが20Å以上，真空層15Å以上のスラブモデルを出力する．
ある結晶について多くの面指数の表面エネルギーが得られている場合，ある面が他の面のファセットに分解するかどうかを判定することができる．たとえばβ-Ga2O3の(60
粒界モデルの作成は，二つの粒が格子点の一部を共有する，すなわちcoincidence site lattice (CSL)を有するように重ね合わせたものを起点とするのが一般的であり，2003年に小川浩が公開した，ウェブブラウザで機能する粒界工房GBstudioでも採用されていた [20]．傾角粒界では，厳密な3次元CSLは，結晶が立方晶か，正方晶か六方晶で回転軸(rotation axis)が<001>の場合のみ可能である．一方，粒界モデルを2つのスラブの重ね合わせと考えると，界面での2次元CSLが存在さえすればよい．スラブを重ねることで粒界モデルを作る詳細が文献 [21]に記されている．
スラブモデルに対し，どの原子が表面にあるかを示す定量的な指標があると便利である．ある中心原子から一定半径内の原子(周辺原子)を抽出し，周辺原子がない空間の立体角Sを指標とすることができる．平らな表面では半球，すなわち2π sr の領域に周辺原子が存在しない [22]．この立体角は
tb2.pl --solid_angle 5 <POSCARで各原子に対し求めることができる．調べる範囲はこの例では5Åだが，系により適切な値が変わってくる．なお，コマンドラインでqconvexプログラムを実行できる必要がある([23]より入手可能)．
ある結晶で考えられる様々な終端のうち，物理的・化学的な根拠に基づいて，スラブモデルを実際に計算せずに安定と思われる終端を予測する記述子が3つ提唱されている．表面で切れた結合の面密度であるunsaturated coordination index (σ) [24]を求めるためには各原子の結合の本数を調べることが必要であるが，各原子周辺にある原子と結合距離の情報は
tb2.pl --dist 2.7 <POSCARのように求められる．この2.7はカットオフ結合長で，単位はÅである．この値は系により大きく変わるが，典型的な場合では最近接結合(nearest neighbor bond)と第二近接結合(next nearest neighbor bond)の長さの間をとればいい．切れた結合の数は，スラブモデルの結合の数と，バルクにおける同じ数の原子の結合数の差をとれば算出できる．
仮想平面からの原子の近さの指標であるatom proximity function of a plane (p)は，原子層の中間にある仮想平面の片側にある各原子に対し，仮想平面からの距離に応じて減衰するガウシアン関数のスコアをつけ，これらを足しあわせて規格化する．値は
tb2.pl –proximity_scan 1 2 3 0.3 <POSCAR_bulkで求められる．このPOSCAR_bulkはバルクのものを使用する．ここでは面指数(hkl)の|h|,|k|,|l|の絶対値がそれぞれ1,2,3以下の全ての終端の値を出す．0.3は，pの導出に使うパラメーターλの推奨値である．この記述子は，バルクをへき開した終端のみにしか使用できない． [25]
表面の粗さを指標化したsurface roughness index (v)は，各原子に対し周辺原子のない空間の立体角Sを求め，閾値Scutを超えた原子におけるS-Scutの和の面積密度をとることで求める．この値は
tb2.pl --roughness 0.75 2 < POSCAR.slabのように求められる．POSCAR_slabはスラブモデルのPOSCARであり，0.75と2はそれぞれパラメーターScut (単位 π sr)とd (無次元)の推奨値である． [25]．
上述の3記述子は，元素の情報を考慮しない．また，導出に表面モデルの緩和を必要としない．pとvは結合の情報も顕に考慮しない上，推奨パラメーターを使うことで事実上パラメーターフリーとなる．
基底ベクトルは一般的に縦ベクトルでa,b,cのように表記されるが，VASPではスーパーセルは横ベクトルで表現され，POSCAR形式では3~5行目に記載される．本稿ではあえて
のように変換行列をかけることになる．単に基底ベクトルを整数倍しない(変換行列の非対角成分に0でない値が含まれる)スーパーセルは，
tb2.pl --supercell d e f g h i j k l < POSCARと入力するだけで求められる．なお，全ての非対角成分が0の場合は，
tb2.pl --supercell d h l < POSCARでよい．
できるだけ直方体に近い(maximally orthogonal) N倍のスーパーセルが欲しい場合，
tb2.pl --N_supercell 10 < POSCARを実行すると，この場合，N=10のスーパーセルの変換行列の要素の羅列が出力される．実際には詳細なオプションを指定することができる．基底ベクトルの各種内積も併せて出力されるので，これらをもとに望ましいスーパーセルを探すことになる．例えば，基底ベクトルが短くなるスーパーセルは
Nudged elastic band (NEB)計算を行う場合，始点と終点の結晶であるPOSCAR.startとPOSCAR.endを用意し，
tb2.pl --neb_images POSCAR.start 4 <POSCAR.endを実行すると，00から05までのディレクトリが作成され，NEBイメージのPOSCARが各ディレクトリに出力される．なお，4は途中のイメージの数で変更可能である．ここでは，POSCAR.startが00に，POSCAR.endが05に入る．ほぼ同じ位置の原子の座標が始点で0.01，終点で0.99と約1ずれることが多々あるが，このようなずれは検出して自動的に対応している．NEB計算後の結晶構造は
tb2.pl –superimpose 01/CONTCAR 02/CONTCAR 03/CONTCAR 04/CONTCAR 05/POSCAR < 00/POSCARのように実行すると，一つのPOSCARファイルに複数のPOSCARの原子位置を重ねることができる．ここでは，標準入力の00/POSCARに他のPOSCARファイルの原子の情報が付加される．
他の有用そうな機能として，格子定数などの結晶の形状に関する基本的な情報を出力する
tb2.pl --data <POSCARPOSCARの6行目の元素名を用いてcifファイルを作る
tb2.pl --cif6 <POSCAR全ての原子のdirect (fractional)座標を0と1の間に収める
tb2.pl --incell <POSCAR原子座標をdirect (fractional)にする
tb2.pl -d <POSCAR逆にcartesianにする
tb2.pl -c <POSCAR原子の座標をすべて (p,q,r) ずらす(数値はdirect とcartesianのうち使用されているもの)
tb2.pl --shift p q r <POSCARが挙げられる．
プログラムはJ-STAGE DATAで公開される.DOI:10.60383/data.jccj.25339036[Program Availability]
Xeon CPU, CentOS 7.9, Perl 5, phonopy 2.19.1で動作確認．プログラムはPerl 5で実行可能な単一ファイルで構成される．Perlの外部モジュールは使用しない．コマンドラインで引数をつけて実行するため，Unix/Linux環境を推奨する．一部機能はコマンドラインでphonopyを実行できる必要がある．テストに用いた入力ファイルと変数はドキュメントに記載されている．
本研究は新学術領域研究(25106005)，科研費(18K04692，21H05101),CREST(JPMJCR17J3),およびNEDO(JPNP21016)の助成金を受けて実施された．
