Bader分析自動化 の履歴(No.1)
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
Quantum EspressoのDFT計算結果に対して、Baderの電荷分析を適用する際の手順を自動化したpythonコードを作りました。Linux上で動きます（Macはちょっとした変更で対応可能）。やっていることは簡単で、Quantum Espressoのpp.x(計算結果の処理プログラム)で電荷分布のcubeファイルを全電子と荷電子それぞれについて作成します（そのためのppの入力ファイルの作成も自動で行います）。出来たcubeファイルを使ってBader分析を行わせます。また、cubeファイルをVESTAでオープン、また個々の原子の電荷が出力されている*-ACF.datファイルをエディタで開くことも自動で行わせます（その結果、上のスクリーンショットのようになります）。Bader分析のやり方については、MateriAppが詳しいです。
Bader関係のところをコメントアウトすると、DFT計算後にVESTAで電子密度分布を見る際のpp.x計算の自動化にもなります…
Quantum Espresso（私が現在使っているのはVer. 7.0）と標準のPython環境が必要です。また、Bader分析のプログラムが必要です。ここではHenkelman groupの公開しているプログラムを使っています。ここからLinux用のbinaryをダウンロードして、解凍すれば問題なく使えました。私はQuantum Espressoのbinディレクトリーにコピーしました。
こちらからpythonコードをダウンロードしてください。１つのpythonのファイルだけです。コードの最初の方のbaderPath, ppPath, vestaPathなどをご自分のシステムに合わせて修正してください。ppPathはQuantum Espressoのbinディレクトリになります。
まずはQuantum Espressoでscf計算(relax/vc-relaxでもよい）を実行します。ただし、この計算にはPAWポテンシャルを使う必要があります。私はpslibraryのpbesolのPAWをよく使っています。PAWが必要な理由は、Bader分析では全電子を扱う必要があるからです。当然ながら通常の擬ポテンシャルには価電子の情報しか入ってません。なお、scf計算等では入力ファイルは*.inとなっていることを想定しています。
例としてMgO.inを使ってscf計算したとします。拡張子を除いたMgOをベースファイル名とします。実行する場合は、このベースファイル名を引数にします。（もし同じディレクトリーにauto-bader.pyがあるなら）、以下を実行します。auto-bader.pyが別のディレクトリにある場合はそこへのパスを加えてください。
./auto-bader.py MgO
これで自動的に全電子と荷電子のpp.xが実行され、その後、Bader解析が実施され、その結果（各原子の電荷）が印刷された*-ACF.datファイルがエディタ（このコードではgnome-text-editor）で表示されます。gnomeでない方は対応したエディタに変更してください。また、VESTAで電子分布が表示されます。
VESTAを起動したくない場合は、cmd = vestaPath...でVESTAを起動しているところ（２箇所あります）をコメントアウトしてください。その次の行もコメントアウトしてください。
エディタが起動しない場合は、使っているシステムのエディタに変更してください(editorPathを書き換える)。
なお、この計算でファイルが色々と出来るので、ディレクトリ(prefix-bader)を作成して、それらをその中に移動してます。もし同じ名前のディレクトリが存在する場合はエラーが出て正常に計算ができませんので、そのディレクトリをrenameするか、削除してから実行してください。
いくつか計算をしてみました。MgOでは電荷が1.64、ZnO(zincite)では1.19と、ZnOの方が小さい（共有性がより強い）となり、Kittelの固体物理に載っているイオン性もそれくらいなのでいいようですが、SiO2（石英）ではSiが3.19、Oが1.59となりました。最上部の図にはMg2SiO4 forsteriteの結果が見えていますが、Mgで1.66, Siで3.09、Oで-1.61になりました。これは私の認識しているSi-Oのイオン性よりもかなり高いので、何故なのか悩んでいるところです。
複雑な結晶構造で単位セル中にある同じサイトの原子で電荷が結構異なることがありました（Mg2SiO4 wadsleyite）。この場合は構造最適化計算を空間群を使って計算しました。体心格子だったので、自動的に単純格子に変換されて計算されます。どうもその辺りで問題があったらしく、空間群使わずに計算したら、電荷は正常になりました。どうもQuantum Espresso側の問題でした。