QE用pythonコード の変更点

#author("2024-08-15T08:23:08+09:00;2024-08-14T10:14:04+09:00","default:masami","masami")
#author("2024-08-19T15:54:16+09:00;2024-08-14T10:14:04+09:00","default:masami","masami")
*Quantum-Espresso用のpythonコード(2023/11/24)
----
（更新）
-20240814 dynG.pyを追加。
-20240515 pp-cutoff.pyを追加。
-20231203 pwout2in.inにバグありました。変換前の*.inファイルがibrav=0だと、cell matrixが以前のものと*.outから抽出したものの両方が*.new.inに出力されてました。直しました。
-20231126 vesta2pw.inにバグありました。space group名をコメントアウトして出力する時に#としてしまった。QEでは!がコメントアウトなので、pw.x実行時にエラーとなります。直して、ファイル更新しました。
----
　Quantum-Espresso(以下QE)をよく使うので、そのための入力ファイルを作ったり、結果を処理するプログラムをpythonでいくつか作ってます。それらは本wikiの所々で紹介してますが、ここでまとめて紹介します。
　これらは私の個人的利用で主に使っていて、エラー処理などもそれほどやっていません。ただ、コメントは結構入れてます（自分が忘れるので）。それから最初の行はpythonへのパスが入っています。これはWindowsでは普段使いませんが、Linux, Macではこの行が必要です。多くのコードではMac/Linuxどちらでも使えるように#!/usr/bin/python3にしてますが、numpy, matplotlib, pillowなどのライブラリーを使っている場合は私の使っているMacのanaconda版へのパスとなっています。この部分をそれぞれのシステムに合わせて変えてください。動かない等はご連絡ください。
**vesta2pw.py
[[download:https://mkanzaki.sakura.ne.jp/codes/vesta2pw.py]]
　このコードはVESTAの結晶構造を記述したファイル（*.vesta）から、pw.xのscf/relax/vc-relax計算用入力雛形ファイルを作るものです。結晶構造の記述に空間群を利用します。空間群が必要ない場合はxtl2pw.pyの方を使ってください。
　このコードを使う前に、結晶構造をCOD等からcifなどの形式で取得して、VESTA上で表示してみます。そしてEdit DataでUnit Cell...の部分を見ます。同じ空間群でも軸の取り方等で異なった表記が存在します。QEで空間群を使う場合は、標準の空間群を使う必要があります。Edit DataでUnit Cell...の部分で右のパネルで１番になっている表記の空間群が標準なので、もし１番が選択されていない場合は、１番に変更して、それで構造を一度変換して（Applyボタンをクリックする）、VESTA形式で保存します。たとえばオリビン（カンラン石）や直方晶系ペロブスカイトではPbnmが普段使われていますが、これは非標準なものですが、歴史的理由から使い続けられています。これらの場合には１番のPnmaに変更する必要があります。そしてApplyボタンを押して、VESTA形式で出力します。そして、以下を実行します(vc-relax計算する場合）、
 % ./vesta2pw.py test.vesta vc-relax
test部分は実際の名前に変えてください。なお「./vesta2pw.py」としているのは、vesta2pw.pyのあるディレクトリーでターミナルから実行する場合です。セキュリティの関係で./が必要です。私はpythonコードはまとめて１箇所に置いておいて、そこを.bashrcのPATHに追加しています。そうすると、どこからでも % vesta2pw.pyで使えます。vc-relax以外にrelax, scfが使えます。これでtest.inというQEのpw.xの入力ファイルが作成されます。場合によって止まる場合が出る場合がありますが、その指示に従ってやり直してください。なお、作成された*.inファイルは私のよく使うpseudopotentialや収束条件になっているので、それらを必要なものに変更してください。最後のK_Pointsも（セルが長いほど小さくなるように）自動で調整してますが、必要な精度が出るように変更してください。また、K_Pointsの最後、デフォルトでは0 0 0となってますが、ここも必要に応じて変えてください（1 1 1とか）。pseudopotentialはpythonコード内に書き込んであるので、それらを自分のよく使うものに変更すれば、毎回直す必要がなくなります。また時々使うコマンドはfile2.write("")文でコメントアウトとして出力しておくと変更する際に便利です（いくつか実際にコメントアウトしている文が入ってます）。普通はこれで出来た*.inファイルでpw.xを実行できるはずです。
　pw.xは空間群から自動的に複合セル(fcc, bcc,底面心など）は単純セルに直して計算されるので（扱う原子数が少なくなり）、短時間で計算されます。また、菱面体セルの場合に、六方晶セルとして取り扱うことが多いと思いますが、空間群を使った場合にちょっと問題があるようで、菱面体セルを使った方がいいようです（この場合だけ２番を使う）。この場合も原子数が1/3になるので、計算も早くなります。この変換でもVESTAが使えます。EditのUnit Cell...のところで２番のrhombohedral cellを選択します。
　空間群を使うと対称性（点群）を正しくプログラムに認識させるのに役立ちます。QEは空間群を使わなくても入力データから点群を求めるのですが、時に正しくない点群を使うことがあります。特にph.xの計算時には正しい点群で計算して欲しいので（各モードの表現がおかしくなるので）、空間群を使うと正しく認識されるようになります。ただ、QEのバグなのか、空間群を使った計算結果がたまにおかしくなることがあります。計算がおかしくなる場合は次のxtal2pw.pyを使ってください（複合セルでは計算時間が長くなりますが）。
**xtl2pw.py
[[download:https://mkanzaki.sakura.ne.jp/codes/xtl2pw.py]]
　これも上記と同じ目的ですが、空間群を使わない場合のコードです（QEは元々空間群を扱えなかったので、こちらの方を先に作ってました）。VESTAの出力ファイルの１つであるxtlファイル形式を使います。これには対称性は入ってなくて、単位格子内の全ての原子位置が収められています。そのため欠陥構造なども扱えます。xtlに記述している格子定数に基づいてibrav（ブラべセルの選択）を自動で設定しています。
　まずはVESTAで計算したいものの結晶構造を表示させてください。そこからExport Data...で、file typeとして、fractional coordinate (xtl)を選択して、出力します。そして以下を実行します(vc-relax計算する場合）、
 % xtl2pw.py test.xtl vc-relax 0
　２番目の引数にはvc-relax以外にrelax, scfが使えます。実行すると*.inというQEのpw.xの入力ファイルが作成されます。３番目の引数は0を指定すると全ての原子を最適化します、1の場合は原子の座標が特殊座標の場合（例えば0,0,0）はその原子の座標を固定します。rの場合は原子位置にランダムな揺らぎを加えます。rの場合の揺らぎの大きさはコード内のdeltaで変えられます。このオプションは構造が鞍点にあると疑われる場合に使えるように最近導入しました。
　原子位置が0に非常に近い場合に単位格子外の余分な原子を完全に排除できないことや逆に内部の原子を消してしまうことが稀にあります。実行した時に単位セル中の原子数がターミナルに出力されますので、それを見ておかしい場合は、VESTA上でxtlを一度読んで、そこで余分な単位格子外の原子を削除してください（このためにはまずはbondを全て解除します）。そしてboundaryで0から0.9999という感じで指定してやるとうまくいきます。そしてそれをまたxtl形式で出力します。どうしてもうまくいかない場合は、xtlを直接編集します（そういう例はすごく稀ですが）。
　vesta2pw.py同様に、作成された*.inファイルは私のよく使うpseudopotentialや収束条件になっているので、それらを必要なものに変更してください。K_pointsも（セルが長いほど小さくなるように）自動で調整しますが、必要な精度が出るように変更してください。また、デフォルトでは0 0 0としてますが、必要に応じて変えてください。デフォルトのpseudopotential等はpythonコード内に書き込んであるので、それらを自分のよく使うものに変更すれば、毎回直す必要がなくなります。質量は理科年表のものを使用しています。
**pwout2xtl.py
[[download:https://mkanzaki.sakura.ne.jp/codes/pwout2xtl.py]]
　QEのpw.xの出力結果には構造データも入ってますが、現在のところVESTAでは直接読めません。VESTAで読めるxtl形式に変換してやるのが、このコードの目的です。なので、まずはQEでpw.xを実行します。その出力ファイル*.outを変換します。なお、*.outがvc-relaxやrelaxで構造が最適化されていない場合（エラーが生じた場合など）には、一番最後の構造が使われます。そのため、あまり意味がありませんが、scf計算の出力でも実行可能です。以下のように実行します。
 % pwout2xtl.py test.out
　Macの場合はVESTAが自動的に起動されて、結果のxtlファイルを表示します。Linux, Windowsの場合はMac用のVESTAの起動部分（コードの最後のところ）をコメントアウトして、それぞれVESTAの起動コマンドに変更してください。コメントアウトしているところにそれぞれの例を示しています。インストール先によっては変更の必要があります。VESTAの起動が不要な場合は、全てコメントアウトしてください。xtlファイルはVESTAにdrag&dropしてもオープンできます。
　複合格子の場合は単純格子で計算して、その結果（構造）も単純格子になっています。このコードでは元の複合格子には戻しませんが、VESTA上で複合格子に戻すことが可能です。最後の格子の変換のところを見てください。
**pwout2in.py
(2023/12/3) update: ibrav=0の入力ファイルだと前のcell matrixが残る不具合を解消
[[download:https://mkanzaki.sakura.ne.jp/codes/pwout2in.py]]
　これはpw.xの出力結果の座標データから入力ファイルを作るコードです。元の入力ファイルの内容をそのまま使いますが、格子定数と原子座標部分を前の計算の出力結果で置き換えます。例えば構造最適化計算が途中エラーで止まった時に最後の構造を使って再計算する場合や、前の構造結果を使って、条件を変えて（圧力を変えるなど）新たな計算をする時に使います。その時に、*.inと*.outの両方を読み込みますので、両方が存在しないとエラーになります。両方読み込む関係で、ファイルの指定は.in, .outなど拡張子を取ったベース名を引数に使います。test.in, test.outがある場合は以下のように指定します。なお、２つのファイルでベース名が異なるとやはりエラーになります。
 % pwout2in.py test
これでtest.new.inが作成されます。必要に応じてrenameしてください。
　新しく作られる*.inは常にibarv=0となります（元の入力ファイル自体が空間群で指定していても）。ibarv=0の場合はセルのマトリックスが新たに*.inファイルのATOMIC_POSITIONSの前に挿入されます。構造最適化がうまく行ってない時は最後の構造を使います。scfの場合にも使えますが、構造は基本入力ファイルと同じままです（処理的には他と同様に*.outから構造を取ります）。
**pp-cutoff.py
[[download:https:/mkanzaki.sakura.ne.jp/codes/pp-cutoff.py]]
　pwscf計算(pw.x)の入力ファイルで設定するecutwfcとecutrhoは、使うpseudopotentialファイル内に記述されているsuggested minimumよりも大きい値を使う必要があります。このプログラムはそれをチェックします。なお、pseudopotentialファイルによってはsuggested minimumが提供されてない場合があり、その場合は有効な答えは返ってきません。私が普段使っているpslibraryのpseudopotentialにはsuggested minimumが記述されています。
 % pp-cutoff.py test.in
実行すると、各pseudopotentialファイルのsuggested minimumのリストとそれらの最大値と入力ファイルで使っている値を標準出力に出力します。また、それらを比較して、設定値が低い場合にはwarningを出します。意外に大きな値が必要な場合がありますので、初めて使うpseudopotentialの場合はこれを実行した方がいいかもしれません。もちろん、pseudopotentialファイル(UPF)はテキスト形式なのでヘッダー部分を読んだり、grepすればいいのですが、元素多く使っている場合は結構面倒です。
**auto-bader.py
[[download:https:/mkanzaki.sakura.ne.jp/codes/auto-bader.py]]
　これはBaderの電荷解析を行う時にプログラムです。Henkelman group（https://theory.cm.utexas.edu/henkelman/code/bader/）の作成したBader解析プログラムが必要です。コードの最初の方のVESTA、pp.x、bader.xのパスを正しく設定します。まずはscf計算を行いますが、全電子を扱う関係で、PAWポテンシャルを使う必要があります。scf計算の後で、以下を実行します。
 % ./auto-bader.py test
実行すると、pp.x用の全電子と価電子の入力ファイルを作り、それを使ってpp.xを実行して、それらの電子密度分布の*.cubeファイルができます。それらをVESTAでオープンします。また、それらのcubeファイルを使ってBader解析を行います。それらのプロセスを自動化してます。詳しくは[[Bader分析自動化]]をご覧ください。
**density.py
[[download:https://mkanzaki.sakura.ne.jp/codes/density.py]]
　これは先のauto-bader.pyから、Bader部分を除いて、単に電子密度分布をVESTAで表示するプログラムです。pw.xを一度実行しておく必要があります。コードの最初の方のVESTAとpp.xのパスを正しく設定します。まずはscf計算（relax/vc-relaxでもOK）を行います。scf計算の後で、以下を実行します。PAWでない場合には価電子密度のみを、PAWの場合はそれに加えて全電子密度も表示します。
 % ./density.py test
実行すると、pp.x用の入力ファイルを作り、それを使ってpp.xを実行して電子密度分布の*.cubeファイルができます。それらをVESTAでオープンします。
**phq0in.py
[[download:https://mkanzaki.sakura.ne.jp/codes/phq0in.py]]
　これはフォノン計算(ph.x)でq=0計算の入力雛形ファイルを作る簡単なコード。基準モード周波数の計算やIR、ラマンの強度計算時に使います。pw.x用の入力ファイルから必要な情報を読むので、そのファイルが必要。また、dynmat.xの入力ファイルも作ります。
 % pwout2xtl.py test.in
これでtest.ph.inが作成されます。これをph.xの入力ファイルにします。epsilとlramanはデフォルトでfalseになってますが、LO-TO分裂を扱いたい時はepsilはtrueにしてください。metalの場合はtrueは使えません。Raman強度の計算をする時はlramanをtrueにします。ただ、Raman強度の計算は現在のところ擬ポテンシャルにLDAかつノルム保存タイプを使った場合にのみ可能です。ちょっとトリックを使うとノルム保存タイプのGGAでも計算できます。それについては[[QEでの振動モード計算]]に書いてます。
　また、質量がコメントアウトされてますが、質量は指定しないとpw.xの計算の時に使った質量を使います。xtl2pw.pyとvesta2pw.pyを使った場合は、その中で定義されている質量が使われていますので、普段は指定する必要はありません。これらは理科年表から取った値です。同位体などで質量が異なるものの振動モードを計算したい場合は、#を取って質量を直接指定してください。例えば水素を重水素にして振動数を計算するなら値を2にします。
**ph2vesta.py
[[download:https://mkanzaki.sakura.ne.jp/codes/ph2vesta.zip]] これは例も含んだzipファイル
　これはフォノン計算(ph.x)でq=0計算の基準モード周波数の結果を使って、そのモードの変位を原子に矢印をつけて可視化するコードです。VESTAで見ることができます。使うためにはph.x計算でq=0の計算を行います(上のphq0in.pyで入力雛形ファイルが作れます）。そしてph.x計算の元となるpw.xの出力から、やはり上記のpwout2xtl.pyでxtlファイルを作成して、それをVESTAで読み込みます。方向、境界や大きさなどを好みの設定にして、VESTAファイルとして保存します。これで準備ができました。以下を実行します。
 % ph2vesta.py test.ph.out test.dyn test.out.vesta
１番目の引数はph.xの出力ファイル名で、２番目はdynファイル名、３番目はVESTAファイル名になります。実行するとモードのリストが出てきて、キー入力待ちになります。見たいモードの「番号」をキー入力します。そうするとそのモードの変位に対応した矢印のついたVESTAファイルが作成され、VESTAが起動して、今作ったファイルを表示します。矢印が短い、色変えたい等があれば「0」を入力すると変更可能です。全モードのVESTAファイルを一度に作成したい場合は「a」を入力してください（この場合、VESTAは起動しません）。終了には「q」をキー入力します。矢印の長さや色などはVESTAのVectors...からも変更可能です。
　なお、複合セルで計算して単純セルに変換された場合は、まずそのままで上記のように進めて、VESTAファイルを作成して、VESTA上で複合セルに戻します。最後の「格子の変換」を参照。
**ph2gif.py
[[download:https://mkanzaki.sakura.ne.jp/codes/ph2gif.zip]] これは例も含んだzipファイル
　上記と密接に関連しているのですが、こちらはモードの変位をgif動画にするコードです。基準モード周波数の結果を使って、そのモードの変位した構造を作って、それをVESTAで表示して、pngで出力します。これにはVESTAを外部ファイルで制御する機能を使っています。そしてそれらをまとめてgif動画を作ります。gif動画を作るところでpillowというライブラリーを使っています。これはpython標準のライブラリーではないので、pipやcondaでインストールする必要があります。その関係で最初の行がanacondaパスになってます。適宜変更してください。実行途中に沢山ファイル(*.vestaと*.png)が生じますが、最後にgifファイル以外は消去されます。それらが必要な場合は消去しているところをコメントアウトしてください（頭に#を追加）。変位した構造を使って構造最適化計算したいことが稀にありますが、その用途に使えます。
　使い方はph2vesta.pyと同じで、ph.x計算でq=0の計算を行います(上のphq0in.pyで入力ファイルが作れます）。そしてph.x計算の元となるpw.xの出力から、やはり上記のpwout2xtl.pyでxtlファイルを作成して、それをVESTAで読み込みます。方向、境界や大きさなどを好みの設定にして、VESTAファイルとして保存します。原子間距離が変化するために、ボンドを通常の原子間距離の設定よりも広めに設定しないと動画の途中でボンドが切れることがあります。これで準備ができました。以下を実行します。
 # ph2gif.py test.ph.out test.dyn test.out.vesta
１番目の引数はph.xの出力ファイルで、２番目はdynファイル、３番目はVESTAファイルになります。実行するとモードのリストが出てきますので、見たいモードの番号をキー入力します。そうするとgifファイルが作成されます（ちょっと時間がかかります）。終了にはqをキー入力します。出来たgifファイルをウェブブラウザーなどで見ると動画として見えます。PowerPointやKeynoteに貼り付けてプレゼンテーションにも使えます。ちょっとコードを変更すると（コメントアウトを外すと）、座標軸に対して回転させながら振動を見るようなことも可能です。
　もし、ファイルが見つからないというエラーが出る時は、VESTAからのpngファイル出力に時間がかかっているためだと思います。forループ中の'sleep 1.0'の時間を1.0から増加させてみてください。
　なお、こちらの場合は複合セルで計算して単純セルに変換された場合に、このコードでは複合セルに戻した状態でのgifファイルを作ることはできません。xtl2pw.pyを使って複合セルのまま計算させるのが１つの解法です。その場合は単純セルで最適化しておいて、それをVESTAで複合セルに戻してxtlで出力、それをxtl2pw.pyに使ってpw.xの入力ファイルを作って、それをscfするだけで、それをベースにph.x計算すればいいと思います。
　できたgif動画で余白が多い場合は、画像の編集プログラム（GraphicConverterとか）を使ってクロップ処理をすればいいです。
**ph2png.py
[[download:https://mkanzaki.sakura.ne.jp/codes/ph2png.py]] 
　これは上記のph2gif.pyからgifを作成するpillowに関連するところを取り除いたもので、それ以外は同じです。ph2gif.pyが使えている方は必要ないはずです。何らかの理由でpillowを使わない・使えない方の利用を想定しています。使い方自体も同じです。
 % ph2png.py test.ph.out test.dyn test.out.vesta
ph2gif.pyとの違いは、実行後にpngファイルが沢山生じていることです（デフォルトで37個）。これらからgif動画を作りたい場合は、フリーのgif動画作成のアプリケーションがあるので、それらを利用することで可能です。私はpillowを使う前は、PicGIFというアプリを使っていました。また、余白が多い場合は、画像の編集プログラム（GraphicConverterとか）でクロップ処理をすればいいです。
**dynmat2plot.py
[[download:https://mkanzaki.sakura.ne.jp/codes/dynmat2plot.zip]] これは例も含んだzipファイル
　これはdynmat.xの出力ファイルを読んで、ガウス関数を各モードにかけて、スペクトル的に変換して、csvに出力する簡単なコード。matplotlibによる簡易なプロットも行います（現状確認用）。利用にはdynmat.xの出力ファイルが必要です。また、matplotlib, numpyが必要です。その関係で最初の行がanacondaパスになってます。
 % dynmat2plot.py test.dm.out 1.05 1.0
第１引数はdynmat.xの出力ファイル名、第２引数は横軸（波数）への補正スケールで、一般に計算された波数（振動数）は実測よりも5~10%小さめになるために、それを補正するためです（まあ補正がリニアでいいかどうかはありますが）。1.0だと補正なしになります。第３引数はガウスピークの半値幅 (cm-1)になります。大きく設定するとピークが幅広になります。４程度が適当だと思います。
　実行するとスペクトルのmatplotlibによるプロットが出てきますが、赤外が赤で、もしあるならば非偏光ラマンが青、偏光ラマン（偏光解消度をかけたもの）が緑になります（通常重なってます）。これらは粉末パターンに相当します。単結晶の方位による計算も可能ですが、QEのプログラムをいじる必要があります。なお、プロットのwindowを閉じないとプログラムが終了しません。また、同じ内容がcsvファイルとしても作られてますので、他のプロットプログラムで見ることや実測値を加えて、合わせてプロットすることが可能です。
　なお、デフォルトでは0から4000 cm-1まで計算します（1 cm-1刻み）。私の場合、OH伸縮振動があるケースだと3700 cm-1くらいまでピークがあるためですが、普段そこまで必要ない場合はコード中のxwidthを変更してください。
　なお赤外吸収の強度などについてはPDielecという粉末試料の赤外線スペクトルをDFT計算の結果から計算してくれるソフトがあります。[[PDielecの使用法メモ]]をご覧ください。
**dynG.py
[[download:https://mkanzaki.sakura.ne.jp/codes/dynG.py]]
　これは上記PDielecでQuantum Espressoの計算結果を使う場合、pw計算をibrav = 0でやって、それからph計算を実施して、出来たdynファイルが必要になります。私はibrav = 0の計算を普段はやらないので、PDielecのためにibarv = 0での計算をやり直すのは面倒です。そこでibrav = 0でないQuantum Espressoの計算結果からPDielecで読めるような*.dynGファイルを作るpythonコードを作りました。一部情報をpwの構造最適化結果から取得するので、pw.xの出力ファイルとph.xの出力ファイルの１つの*.dynファイルが必要です。
　これは上記PDielecでQuantum Espressoの計算結果を使う場合、pw計算をibrav = 0でやって、それからph計算を実施して、出来たdynファイルが必要になります。私はibrav = 0を使っての計算を普段やらないので、PDielecのためにibarv = 0での計算をやり直すのは面倒です。そこでibrav = 0でないQuantum Espressoの計算結果からPDielecで読めるような*.dynGファイルを作るpythonコードを作りました。一部情報をpwの構造最適化結果から取得するので、pw.xの出力ファイルとph.xの出力ファイルの１つの*.dynファイルが必要です。
 % dynG.py test-pw.out test-ph.dyn
これでtest-ph.dynGが作成されます。これをPDielecから読みますが、PDielecは拡張子.dynGしか読まないので、拡張子は変更しないでください。なおph計算ではepsil = .true.にすることを忘れないように。
　まだ問題があって、面心格子などで格子定数がおかしくなる問題があります…１つの回避策は構造をVestaからxtlで出力して、それをxtl2pw.pyで出力ファイルに変換して計算を進めることです。
これでtest-ph.dynGが作成されます。これをPDielecから読みますが、PDielecは拡張子.dynGしか読まないので、拡張子は変更しないでください。なおph計算ではepsil = .true.にすることを忘れないように（デフォルトはfalseです）。
**格子の変換
　QEで空間群を使って複合セル（F,I,C）などを扱った場合に、単純格子(primitive cell)に自動的に変換して計算してくれるのですが（その方が計算時間は短くなります）、計算後に格子（座標）を元に戻してくれるほど親切ではありません…それで困らないことも多いのですが、構造の変化など見る時は分かりづらいので元のセルに戻すことが必要となります。これにはVESTAの機能を使えば元のセルに簡単に戻すことができます。まあ変換前後の座標系を考えて、変換（回転）マトリックスがどうなるかを考えたらいいだけなのですが、毎回考えるのも面倒なので、以下に回転マトリックスPをまとめています。これをVESTAのEdit DataでUnit Cell...、Transformボタンをクリックして、回転マトリックスのところを以下のように書き換えて、OKボタンをクリックします。これらの場合、体積が増加するとのwarningが出ますが、OKをクリック。原子を追加するをクリックします。そしてApplyをクリックします。これで変換されたはずです。右手系から左手系へ変換されるやmatrixがゼロになるwarningが出たら、それはどこかが間違っています。
***単純セルから面心セル(fcc)へ戻す回転マトリックス
,-1,-1,1
,1,-1,-1
,1,1,1
***菱面体から六方晶セルへ戻す回転マトリックス
,-1,1,1
,0,-1,1
,1,0,1
***単純セルからC底面心へ戻す回転マトリックス
,1,1,0
,-1,1,0
,0,0,1
***単純セルから体心セル(BCC)へ戻す回転マトリックス
,1,0,1
,-1,1,0
,0,-1,1
または
,1,1,0
,-1,0,1
,0,-1,1