はじめに
VASPは密度汎関数理論(Density Functional Theory: DFT)に基づいて電子状態を計算するソフトウェアです。擬ポテンシャル法を採用しており、全電子計算を行うWIEN2kなどのソフトウェアに比べて高速に計算できます。また、精度に関しても良い成績を残しており、様々な分野で広く使われています。
しかし一方で、VASPの出力形式は状態密度(Density of state: DOS)・バンド分散・フェルミ面などの基本的な解析結果を可視化しづらく、計算実行後に出力ファイルを編集する手間が必要です。この点はQuantum ESPRESSOの方がユーザーフレンドリーな設計になっていると感じます。尤も、後処理を行う手段を確立しさえすれば問題ではなくなるのですが、初めてVASPを利用する方にとっては躓きやすい場所だと思います。
VASPを利用する際、事前・事後処理を行ってくれるソフトウェアを利用するか、自分で出力ファイルの編集を行う必要があります。現在、無償・有償含めてVASPの事前/事後処理計算を行うソフトウェアは幾つかあり、そのうち無償で提供されているものはpythonを利用したCLIツールです。
基本的にはこれらのツールの利用をお薦めしますが、はじめてVASPを利用される方にとって、出力ファイルの内容をある程度把握し、編集する方法を知ることは無駄にはならないと思います。また、で紹介されているp4vaspという無償のGUIツールがありますが、python2系で記述されており現在サポートされていないことや、フェルミ面の描画に対応していない点で、著者個人としては使用をお薦めできません。そのためチュートリアルに従ってVASPの使い方を勉強している際、p4vaspを利用せずにバンド分散・状態密度を出力する方法を紹介することは重要だと考えました。
そこでこの記事では、VASPの実行結果から状態密度・バンド分散・フェルミ面を可視化するには、どのように出力ファイルを編集すればよいかを紹介します。以下では、グラフ出力する方法としてgnuplotやexcelなどを利用することを想定していますが、他のグラフ出力ソフトにも適用できると考えています。
本記事の内容は基本的に非推奨の方法なので、難しいと感じた方は、で紹介しているソフトウェアの利用を参考にしてください。
VASPを実行後、DOSCARファイルに状態密度の情報が出力されます。このファイルを編集してgnuplotなどのグラフ出力ソフトで可視化できるようにすることがここでの目的となります。
DOSCARファイルの数値の詳しい説明はを確認してください。DOSCARの出力形式は、部分状態密度の解析を有効した場合、スピンを有効にした場合などで変化しますが、ここでは最も単純な場合についてのみ解説します。状態密度を可視化するうえで重要になる情報は、DOSCARファイルの6行目と7行目以降になります。
6行目は以下のような情報が出力されます。
E(max), E(min),NEDOS,E(fermi),1.0000
ここで重要になる量はNEDOSとE(fermi)です。NEDOSは状態密度のデータ点数を表しており、E(fermi)はフェルミエネルギーです。
7行目以降は、
energy DOS integrated-DOS
の順で出力されており、NEDOSの値だけこの行が続きます。 以上が、DOSCARファイルから状態密度を可視化するために必要な情報です。
では、DOSCARからグラフ出力するためのファイルを作成しましょう。 フェルミエネルギーを基準にしない場合、必要な作業は7行目以降のデータを別のファイルに貼り付けるだけで完了です。フェルミ準位を基準にしたい場合は、1列目のエネルギーの値をフェルミエネルギー分だけ引いたものを作る必要があります。編集の方法は様々ですが、ここでは一例としてFortranを使ったプログラムを作成しました。DOSCARファイルが存在するディレクトリでgfortranなどでコンパイルをし、実行すればdos.gnuファイルを出力するようにしています。このプログラムは可読性を重視するためにスピンを有効にした計算には対応していません。
PROGRAM MAIN implicit none integer file_DOS,file_OUT integer Ne double precision en,dos,sdos,Ef double precision dummy(2) integer i0 !open files open(newunit=file_DOS,file= "DOSCAR",status= "old") open(newunit=file_OUT,file="dos.gnu",status="replace") !skip header informations do i0 = 1,5 read(file_DOS,*) enddo read(file_DOS,*) dummy(1),dummy(2),Ne,Ef do i0 = 1,Ne read(file_DOS,*) en,dos,sdos write(file_out,*) en-Ef,dos,sdos enddo END PROGRAM
バンド分散
では次に、バンド分散の可視化について紹介します。バンド分散の可視化にはEIGENVALファイルを用います。まずは入力ファイルの一つであるKPOINTSで、k点の経路を指定してVASPを実行しましょう(方法は別の記事に譲ります)。EIGENVALファイルにはフェルミエネルギーの情報がないので、DOSCARファイルからフェルミエネルギーを読んでくる必要があります(この時のDOSCARファイルは、バンド分散の計算を実行した際に得られたものではなく、状態密度計算を実行した際に得られたDOSCARファイルを利用してください)。その他の値でバンド分散を可視化するために必要な情報は、EIGENVAL内の6行目と8行目以降です。
6行目は、
Ne,Nk,Nb
となっており、それぞれ、電子数・波数点の数・バンド数を表しています。このうち波数点の数とバンド数を可視化の際に利用します。
7行目以降は以下のような構造をしており、
k1x k1y k1z k-weight 1 eval(1,1) e-weight 2 eval(1,2) e-weight 3 eval(1,3) e-weight k2x k2y k2z k-weight 1 eval(2,1) e-weight 2 eval(2,2) e-weight 3 eval(2,3) e-weight ： ：
波数点(k1x,k1y,k1z)における固有エネルギーevalがNbだけ羅列されています(上記はNb=3を想定)。このブロックがNkだけ続く構造です(e-weight，k-weightは、単純な解析の場合バンド分散の可視化で使用しないためインデックスを省略しました)。スピンを有効にした場合は上記の構造が変わりますが、この記事ではそちらについては解説しません。
さて、グラフ出力ソフトを使って可視化するには、EIGENVALから以下のような形式のファイルを作成する必要があります。
1 eval( 1,1) 2 eval( 2,1) 3 eval( 3,1) : Nk eval(Nk,1) 1 eval( 1,2) 2 eval( 2,2) 3 eval( 3,2) : Nk eval(Nk,2)
k点の数はKPOINTSで指定した数だけ出るので、手作業で行うことは現実的ではありません。また、フェルミエネルギーを基準点に取りたいときは、DOSCARからフェルミエネルギーを持ってくる必要があります。編集の方法は様々ですが、ここでは一例としてFortranを使ったプログラムを作成しました。EIGENVAL，DOSCARファイルが存在するディレクトリでgfortranなどでコンパイルをし、実行すればband.gnuファイルを出力するようにしています(この時DOSCARファイルは、状態密度の計算で得られたものを利用してください)。このプログラムは可読性を重視するためにスピンを有効にした計算には対応していません。
PROGRAM MAIN IMPLICIT NONE integer file_EIGEN,file_DOS,file_OUT integer Nk,Nb double precision,allocatable:: eng(:,:) double precision Ef double precision dummy(4) integer i0,ik,ib !==================================== !get Fermi energy from DOSCAR file !==================================== !import DOSAR file open(newunit=file_DOS ,file= "DOSCAR",status= "old") !skip 5 lines do i0 = 1,5 read(file_DOS,*) enddo !get fermi energy (DOSCAR 6th line: Emax,Emin,NEDOS,Ef,1.0) read(file_DOS,*) dummy(1),dummy(2),dummy(3),Ef,dummy(4) close(file_DOS) !==================================== !get number of K-points and bands !==================================== !import EIGENVAL file open(newunit=file_EIGEN,file="EIGENVAL",status= "old") !skip 5 lines do i0 = 1,5 read(file_EIGEN,*) enddo !(EIGENVAL 6th line: Ne,Nk,Nb) read(file_EIGEN,*) dummy(1),Nk,Nb !==================================== !get energies !==================================== !construct energy array allocate(eng(Nk,Nb)) do ik = 1,Nk !skip blank line read(file_EIGEN,*) !read k coordinate (not used in this program) read(file_EIGEN,*) dummy(1),dummy(2),dummy(3) do ib = 1,Nb !read eigen energy read(file_EIGEN,*) dummy(1),eng(ik,ib),dummy(2) enddo enddo close(file_EIGEN) !==================================== !construct output file !==================================== !make output file open(newunit=file_OUT ,file="band.gnu",status="replace") do ib = 1,Nb do ik = 1,Nk write(file_OUT,*)ik,eng(ik,ib) - Ef enddo write(file_OUT,*) enddo close(file_OUT) !==================================== END PROGRAM
フェルミ面
フェルミ面の描画を自力で行うことは、初めて利用される方にはお薦めできません。また、事前・事後処理を行うソフトウェアでもフェルミ面の描画はサポートされていない場合があります。
フェルミ面を可視化するためのビューアとしてはXcrysdenやFermisurferがあり、これらを利用するためにはBXSF形式やFRMSF形式のファイルを作成する必要があります。EIGENVALファイルなどからBXSF形式などのファイルを作成することはできますが、基本的にはで紹介するソフトウェアを利用することをお薦めします。
おわりに
この記事では、出力ファイルを直接編集することで結果を可視化する方法について紹介しました。基本的にお薦めできない方法ですが、DOSCARやEIGENVALがどのような構造をしているかを知ることは重要だと思います。VASPから状態密度・バンド分散・フェルミ面を可視化する方法は上記のように煩雑なので、基本的にはで紹介するソフトウェアを利用することをお薦めします。 また、上述した事後処理に煩わされずに簡単に結果を可視化したい場合、Quantum ESPRESSOを利用することも一つの手段だと思います。
