ソフトウェアニュース・レビュウ
結晶学を活用した結晶モデル作成支援ソフトウェア「壺車」の2024年更新版
2025 年 24 巻 1 号 p. A12-A17
詳細
2025 年 24 巻 1 号 p. A12-A17
The author's code for VASP POSCAR handling is updated with new methods developed in 2024. Size "1/N"supercells can be obtained to find identical or smaller size unit cells that may belong to a Bravais lattice with a lower degree of freedom if deviations from rigid requirements for each Bravais lattice are allowed. Sortable number "message digests" of the configuration of different species on a (sub)lattice can be derived, which is useful for eliminating identical configurations in exhaustive combinatorial searches of species distributions. A molecule containing a given atom can be extracted from, for example, a molecular crystal supercell with many molecules.
第一原理計算等用の結晶モデルの作成支援ソフトウェアとして，著者はツール「壺車(つぼぐるま)」を開発した． [1]結晶構造はVASP POSCAR形式である．マニュアルが日本語で，UNIXコマンドラインから実行する1ファイルで完結するためインストールと実行が容易な特徴がある．なお，コマンドラインからphonopy [2]等の別のプログラムの実行を要することもある．本論文では，著者が2024年に開発した手法を実装した更新版「壺車2.1」を紹介する．前版 [1]と更新版プログラムは両方とも，J-STAGE DATAにてCC-BYライセンスで公開されている．
結晶の基底ベクトルを与えた時，体積が1/N倍となる「1/Nスーパーセル」と，格子定数に関する理想値からの多少のずれを許した場合のブラべ格子を求める手法を開発した [3]．実験研究者は結晶構造解析で誤差を含む格子定数を得ることができる．ある程度のずれを許すことで，より自由度が低い，すなわち対称性の高いブラべ格子で結晶を表現できるのであれば，そのような基底ベクトルの取り方が知りたいこともある．電池材料を例に挙げると，層状正極材料LiCoO2のCo層をLi2/3Mn1/3層に置換したものが過剰系正極材料Li2MnO3である．一方，LiとCoの2つの原子を同一視すると概ね岩塩型構造になり，また全ての原子を同一視するとほぼ単純立方格子となる．ブラべ格子と，慣用単位胞(conventional cell)と基本単位胞(primitive cell)の原子数は以下のとおりである．単純立方格子:(単純立方格子, 1,1)，岩塩型構造:(面心立方格子, 8,2)，LiCoO2:(菱面体格子, 12,4)，Li2MnO3:(底心単斜格子, 24,12)．例えば12原子のLi2MnO3の基本単位胞の基底ベクトルをもとに，「1/3スーパーセル」として4原子のLiCoO2の基本単位胞に相当する単位胞があることを発見し，基底ベクトルの変換行列を求め，格子定数a, b, cがほぼ一定と定量的にみなすことができれば菱面体格子とみなせることを提示するのが，この手法の目的である．仕様上，原子位置の情報を無視して格子点のみに着目するため，ブラべ格子は判定できても空間群自体は考慮できない．また，単斜晶と直方晶では，結晶学的慣習に基づいた基底ベクトルの取り方は，原理的に必ずしも実現できない．まずは
tb2.pl --orthonormal_M_N2 N > orthonormal_M_N2.outを実行する．ここでは，「N2倍のスーパーセルで，さらにN倍すると元のセルのN×N×N倍になりうるもの」を全て求めた情報を出力する．任意の格子はアフィン変換により正規直交基底の単純立方格子に変換できる．丸め誤差が最大限に回避でき，また分率座標(fractional coordinate)系とデカルト(Cartesian)座標系の変換が自明なことから，N2スーパーセルの基底ベクトル探しは単純立方格子で探索する．最も直方体に近い(maximally orthogonalized)スーパーセルを探す手法は文献 [4]に記載されている．基底ベクトルa，b探しが概ねO((N2)6)= O(N12)となる上，並列計算が不可能であるため，無駄を減らす小技を駆使しても，なお時間がかかる．しかし，この情報はNのみに依存するため，プログラムとしてこの部分を切り分けている．J-STAGE DATAにて，N≤47の計算結果をファイル名orthonormal_M_N2_*.txtとして公開する．その後，
tb2.pl --1/N N [eps_a] [eps_d] [DOF1|DOF2|DOF3|DOF4] <POSCAR
tb2.pl --1/N orthonormal_M_N2.out [eps_a] [eps_d] [DOF1|DOF2|DOF3|DOF4] <POSCARを実行することで，実際のセルの情報を得ることができる．eps_a，eps_dは元論文 [3]中の，それぞれ主に格子定数の角度と長さに関する，理想値からのずれの指標εa，εdで，デフォルト値が0.1と0.01である．対称性が低すぎるブラべ格子の出力を防ぐため，自由度を制限するDOF1, DOF2, DOF3, DOF4のいずれかをeps_dの後に指定できる．自由度が1(立方晶)，2(正方晶，六方晶，三方晶)，3(直方晶)，4(単斜晶)以下の全てのブラべ格子の可能性が欲しい場合，それぞれDOF1，DOF2，DOF3，DOF4を指定すればよい．
2.2 元素種配列の同一性判定格子，もしくは一部の副格子に複数種の元素種(species)が入ることができる結晶を考慮する．ここでの「元素種」には空孔(vacancy)も含む．また，スピンの向きや価数を変えた原子を別元素種扱いすることもできる．指定した(副)格子上の安定な元素種の配列を得る場合に，組み合わせを総当たりしてエネルギーを計算することは有効である．しかし，対称性が高い場合は等価な配列が多く発生し，全て計算することは効率が悪い．体心立方格子の2×2×2スーパーセルでは原子が入るサイトが16個ある．ここでは原子サイトと格子点は合致する．しかし，格子点は基本単位胞につき1つであるが，サイトは基本単位胞に複数含まれることがあるので，必ずしも格子点とサイトは同一ではない．16のサイトに元素AとBを8個ずつ占有させる場合，素直に考えると合計16C8=12,870の組み合わせが存在するが(なお，一つのサイトをAに固定することで，容易に組み合わせを15C7=6,435に半減できる)，対称性が異なる組み合わせはわずか32種しかない [5]．この例では，対称性を考慮することで，計算時間を約400倍加速することができる．
情報通信(IT)分野では，長い文字列であり，多くの情報を含むファイルに対し，ハッシュ化の処理を行うことで，短い決まった文字数のメッセージダイジェストを決定する．ファイルが伝達したい情報であるメッセージであり，これを短い文字列に変換したものがメッセージダイジェスト(要約)である．ファイルの内容が少しでも違うと全く異なるメッセージダイジェストが生成されるため，メッセージダイジェストを比較することでファイルの内容の同一性をある程度確認することができる． [6] 元素種配列の同一性判定においては，結晶構造をIT分野のファイルと考え，結晶構造固有の短い情報をIT分野のメッセージダイジェストと捉える．格子定数も一種のメッセージダイジェストと捉えられるが，あまりにも情報が欠落しすぎて結晶の同一性判定の役に立たない．
結晶では，原子間にいかなる相互作用を仮定しても，同一の結晶に対しては必ず同じ物性値が得られる．原子i,jの二体間の原子間ポテンシャルは，例えばUij=U(qi, qj, rij)と表現することができる．ここで，qi, qjは元素種ごとに決まる値で，rijは原子間距離である．元素種ごとに適当なqの値を決定し，任意のでたらめな「ポテンシャル」Uを定義することで，浮動小数点数の「エネルギー」を求めることができる．結晶のメッセージダイジェストは，IT分野と異なり，固定長であることは求められない．むしろ，原子位置の僅かな数値誤差にほぼ影響されず，ソート可能な数字であることが求められる．
文献 [5]は2種類の「ポテンシャル」を用意した．すなわち，元素種ごとの値qとして事前に定められた0.5から1.5の値を使用し，結晶ごとに二体間「ポテンシャル」を計算する距離のカットオフを決定し，「ポテンシャル」としてFij= qiqj(rij2-rij+1)とGij= qiqj(0.5rij2-3.5rij+6.5)を定義している．二つの異なる結晶が偶然ほぼ同じメッセージダイジェストに至るハッシュ衝突現象が起きうるため，回避するために複数の「ポテンシャル」が必要となることがある．現実的な「ポテンシャル」は意外かもしれないが都合が悪い．通常のポテンシャルは原子間距離が大きくなると値が速く0に収束するが，メッセージダイジェストを求める際には，それなりに長い距離の原子配置の情報を拾いたいため，長い距離でもある程度の値を有す「ポテンシャル」が望ましい．
VASP5 POSCAR形式では6行目に各副格子に入る元素を指定するが，複数元素種が入る場合，A_Bのような形で指定する．まず
tb2.pl --digest_get_environment < POSCAR > env_fileを実行することで，相互作用をまとめたデータファイルenv_fileを作成する．
各元素種ごとに非負整数の元素種IDが定められ，結晶の各サイトの元素種IDを指定することで，メッセージダイジェストを求めることができる．例えば，4つのサイトの元素種IDが0,2,0,4の場合
tb2.pl --digest_use_environment env_file 0 2 0 4を用いてメッセージダイジェストを得る．この例では，サイトの元素種IDをもとに元素が正しく配置されたPOSCARを得るために，
tb2.pl --digest_recover_poscar env_file 0 2 0 4を実行する．メッセージダイジェスト導出の実装については，3節で詳述する．
2.3 分子の抽出特定の原子を含む分子に関する機能を追加した．原子間距離が指定した値distance以下の原子2つは，結合しているものとする．分子内の指定した原子と結合している全ての原子は，その分子に含まれているものとする．分子内の一つの原子を指定することで，その分子内の全ての原子を，各原子が結合している原子のリストを元に決定することができる．
グラフ理論の概念を借りると，原子はノード(頂点)であり，結合している原子ノードはエッジ(辺)で接合している．一つの原子ノードを指定することで，そのノードが含まれるグラフを特定できる．そのグラフに含まれるノードに対応する原子は，同じ分子に属する．
分子の特定の具体的なアルゴリズムは以下のとおりである．分子に含まれる原子の「調査済」配列と「未調査」配列を考える．
(1) 最初に指定した原子を「未調査」配列に入れる．
(2) 「未調査」配列から一つの原子Aを取り出し，「調査済」配列に移す．
(3) その原子Aと結合する全ての原子を求める．「調査済」配列に含まれていない原子を全て，「未調査」配列に入れる．
(4) 「未調査」配列に原子が残っている場合，(2)に戻る．残っていない場合，「調査済」配列の原子が，求める分子の構成原子である．
今回の実装では，原子の結合について，周期的境界条件を考えないものとする．これにより，無限に大きい分子の生成を避けることができる．分子が周期的境界条件の境界をまたがない限り，分子は周期的境界条件下で適切に特定される．
一般的な有機分子ではC，N，O，H間の共有結合の距離が1.6 Åを超えることがあまりなく，一方，イオン結晶では結合距離が1.6 Åを下回ることが稀である．このため，共有結合性の軽分子を扱う場合，1.6がdistanceの値の目安となる．
POSCAR内のatom_number番目の原子が含まれる分子のPOSCARを出力するには
tb2.pl --extract_molecule atom_number distance < POSCARを実行する．分子内の原子の番号(元素に関わらず1からの通し番号)が必要な場合は
tb2.pl --molecule_atoms_numbers atom_number distance < POSCARとなる．atom_number番目の原子を含む分子だけを取り除くには，
tb2.pl --remove_molecule atom_number distance < POSCARatom_number番目の原子を含む分子だけの座標の数字を(x,y,z)だけ動かすには，
tb2.pl --shift_molecule atom_number distance x y z < POSCARを実行すればよい．最後の機能は分子の脱離の計算に有用である．
関連する機能として，指定した原子だけの座標の数字を(x,y,z)だけ動かすコマンド，
tb2.pl --shift_these_atoms x y z A B C.. < POSCARが追加された．ここで，A B C… は動かす原子の番号であり，原子を好きな数だけ指定することができる．1..3のような連番指定も可能である．--shift_moleculeで吸着後の分子がうまく指定できない場合は，脱離状態のPOSCARで--molecule_atoms_number で分子内原子を特定し，--shift_these_atomsで動かせばよい．
2.4 その他の追加機能POSCAR内のN_atom番目の原子に対し，元素をN_species番目に変え，N_species番目の元素の最後に移動するには
tb2.pl --change_atom_species N_atom N_species < POSCARを実行する．例えば，元素A，B，Cの原子がそれぞれ2，4，6個あるPOSCARを考える．N_atomが9，N_speciesが1の場合，POSCARの9番目の原子，すなわち元素Cの3番目の原子の，元素名がAとなる．この原子は元素Aの末尾となるPOSCARの3番目の原子となり，元素A，B，Cの原子数がそれぞれ3，4，5となる．これはドーパント原子の設定と計算に便利である．
元素FROMの原子を全て既存の元素TOの原子にして(元素をマージし)，元素FROMをなくすには，
tb2.pl --merge_element FROM TO < POSCARを実行する．
POSCAR内の各元素内の原子をシャッフル(乱数を用いて順番を変える)するには
tb2.pl --shuffle <POSCARを実行する．ある元素に複数の原子があり，このうちN個ランダムで選んで削除したい場合，その元素の原子をシャッフルし，下からN個の原子を--rm_atoms等で削除すればよい．全ての元素の原子がシャッフルされる．シャッフルの結果，元素が変わることはない．
原子の動径分布関数(radial distribution function, pair correlation functionとも)を求めるには，
tb2.pl --rdf species1 species2 [length] [sigma] [binsize]を実行する．元素種species1とspecies2は番号で指定する．例えば，POSCARの6行目が Li Co O だとしたら，Li-Co結合の動径分布関数を見たかったらspecies1とspecies2に1と2，もしくは2と1を指定する．動径分布関数は距離lengthまで求める．セル内の原子を1つ以上含む結合を判定する．同じセル内の2つの原子間の結合は，ダブルカウントされる．結合距離はrに対しデルタ関数的であるため，スペクトルとして表示するにはピークを広げる必要がある．
動径分布関数は，
で定義される．ここで，σはピークの広がりを示す値sigmaであり，検討するrの最大長がlength，ri は結合iの長さである(最大考慮値がlength)．RDF(r)の値は，0からlengthの間，binsize間隔で求められる．このため，出力されるRDFは長さlength/binsizeの配列となる(実際にはrの値とペアで出力される)．exp項を全てのrに対し計算するのは非効率的であるため，ri-4σからri+4σの間でのみ考慮する．なお，exp{-(4σ)2/(2σ2)}≈0.0003である．RDFの大きさがある程度揃っていると，異なる元素ぺアのRDFのグラフ比較が容易になるため，RDFの最大値が1になるよう規格化係数cが決定される．パラメータのデフォルト値として，length=8，sigma=0.1，binsize=0.01が設定されている．
POSCARファイルの原子位置と元素の情報をGaussianのgjf形式に変換するコマンドとして，tb2.pl --gjf [valence] [spin_multiplicity] <POSCARが追加された．入力にPOSCAR計算条件の情報がないため，完全なgjfファイルは作成できないものの，計算条件#t opt=tight b3lyp/6-31++g(d,p) formcheckおよび価数(valence) 0，スピン多重度(spin_multiplicity) 1を仮定したgjfファイルが生成される．引数のvalenceとspin_multiplicity は任意入力でデフォルトはそれぞれ0と1，spin_multiplicityを入力する場合はvalenceの入力も必須である．このgjfファイルはそのまま実行できないものの，POSCARから原子位置と元素の情報を簡便にgjf形式に変換できる．
結晶のメッセージダイジェストの例として，閃亜鉛鉱(zincblende)型(InGa)Asの8原子ユニットセルを検討する．POSCARとして以下を用いる．As(Ga,In)15.8 0 00 5.8 00 0 5.8As Ga_In4 4D0 0 00.5 0.5 00 0.5 0.50.5 0 0.50.25 0.25 0.250.75 0.75 0.250.25 0.75 0.750.75 0.25 0.75tb2.pl --digest_get_environment <POSCAR を実行すると，以下が標準出力される．digest_get_environment site no kazu: 800001_21_21_21_2==POSCAR==As(Ga,In)1.05.79999999999999 0.00000000000000 0.000000000000000.00000000000000 5.79999999999999 0.000000000000000.00000000000000 0.00000000000000 5.79999999999999As Ga_In4 4Direct(8)0.00000000000000 0.00000000000000 0.00000000000000 As0.50000000000000 0.50000000000000 0.00000000000000 As0.00000000000000 0.50000000000000 0.50000000000000 As0.50000000000000 0.00000000000000 0.50000000000000 As0.25000000000000 0.25000000000000 0.25000000000000 Ga_In0.75000000000000 0.75000000000000 0.25000000000000 Ga_In0.25000000000000 0.75000000000000 0.75000000000000 Ga_In0.75000000000000 0.25000000000000 0.75000000000000 Ga_In==Interactions==0 4 17.8706889953715 13.79741148380010 5 17.8706889953715 13.79741148380010 6 17.8706889953715 13.79741148380010 7 17.8706889953715 13.79741148380011 4 17.8706889953715 13.79741148380011 5 17.8706889953714 13.79741148380011 6 17.8706889953715 13.79741148380011 7 17.8706889953715 13.79741148380012 4 17.8706889953715 13.79741148380012 5 17.8706889953714 13.79741148380012 6 17.8706889953714 13.79741148380012 7 17.8706889953714 13.79741148380013 4 17.8706889953715 13.79741148380013 5 17.8706889953715 13.79741148380013 6 17.8706889953714 13.79741148380013 7 17.8706889953714 13.79741148380014 4 18 94 5 12.6862915010152 20.40202025355334 6 12.6862915010152 20.40202025355334 7 12.6862915010152 20.40202025355335 5 18 95 6 12.6862915010152 20.40202025355335 7 12.6862915010152 20.40202025355336 6 18 96 7 12.6862915010152 20.40202025355337 7 18 9
2行目以降に0が4行，1_2が4行あるのは，最初の4つのサイトには元素種ID 0のAsが，次の4つのサイトには元素種ID 1と2のGaとInがそれぞれ入ることを示している．元素種IDはPOSCARの6行目の登場順で決まる．次にPOSCARが転記され，最後にGa_Inサイトが絡む2体間相互作用の情報が入っている．左2つの整数は0から始まるサイトIDである．複数元素種が入らないサイト同士の相互作用は無視されるため，As-As間の相互作用，すなわちサイトID 0-3同士の相互作用は入らない．複数元素種が入るサイトと入らないサイトの相互作用は，この例では全てのGa_Inサイトが等価なため同じであるが，一般的には等価でないため考慮する必要がある．二つの小数は，2つのメッセージダイジェストの相互作用にそれぞれ対応する．
tb2.pl --digest_get_environment < POSCAR > env.fileを実行して，上の結果をファイルに保存する．
As4Ga2In2の組成を考える場合，最初の4サイトには元素ID 0のAs，残り4サイトには元素種ID 0と1のGaとInが2つずつ入るため，以下の6つを実行し，考えられる配置のメッセージダイジェストを得る．
tb2.pl --digest_use_environment env_file 0 0 0 0 1 1 2 2
tb2.pl --digest_use_environment env_file 0 0 0 0 1 2 1 2
tb2.pl --digest_use_environment env_file 0 0 0 0 1 2 2 1
tb2.pl --digest_use_environment env_file 0 0 0 0 2 1 1 2
tb2.pl --digest_use_environment env_file 0 0 0 0 2 1 2 1
tb2.pl --digest_use_environment env_file 0 0 0 0 2 2 1 1後半の元素種IDの羅列は，原子数が多くなる場合などは特に，別途プログラムを作って生成するとよい．
実行すると2つのメッセージダイジェストが出力される．結果は全て53.3585182715025 45.0129670548652となり，全てのGa/In配置が等価なものであることが判明する．
元素種IDからPOSCARを書き起こす場合，例えばtb2.pl --digest_recover_poscar env_file 0 0 0 0 1 2 1 2を実行するとAs(Ga,In) 0 0 0 0 1 2 1 21.05.79999999999999 0.00000000000000 0.000000000000000.00000000000000 5.79999999999999 0.000000000000000.00000000000000 0.00000000000000 5.79999999999999As Ga In4 2 2Direct(8)0.00000000000000 0.00000000000000 0.00000000000000 As0.50000000000000 0.50000000000000 0.00000000000000 As0.00000000000000 0.50000000000000 0.50000000000000 As0.50000000000000 0.00000000000000 0.50000000000000 As0.25000000000000 0.25000000000000 0.25000000000000 Ga0.25000000000000 0.75000000000000 0.75000000000000 Ga0.75000000000000 0.75000000000000 0.25000000000000 In0.75000000000000 0.25000000000000 0.75000000000000 Inが得られる．外部プログラムで用意した元素種IDをPOSCARに変換するのは手間がかかるので，この機能を用意している．原子の位置は元素種ごとにまとめるためGa_Inサイトの原子の順番が元のPOSCARと変わっている．
J-STAGE DATAに，著者が2024年に開発した2つの手法，および分子の取扱に関する機能等を追加実装した「壺車」の更新版を公開した．今後も著者が手法を開発した場合，追加実装し，J-STAGE DATAで公開していく予定である．
プログラムはJ-STAGE DATAで公開される. DOI:10.60383/data.jccj.28069799
[Program Availability] Xeon CPU, Rocky Linux 8, Perl 5, phonopy 2.19.1で動作確認．プログラムはPerl 5で実行可能な単一ファイルで構成される．Perlの外部モジュールは使用しない．コマンドラインで引数をつけて実行するため，Unix/Linux環境を推奨する．一部機能はコマンドラインでphonopyを実行できる必要がある．テストに用いた入力ファイルと変数はドキュメントに記載されている．
本研究は科研費(24H00395) ，NEDO(JPNP21016) ，触媒科学計測共同研究拠点，共同利用・共同研究(the Joint Usage/Research Center for Catalysis)の助成金を受けて実施された．
