#イントロ
お疲れさまです。橋本です。
VASPはLDOSとPDOSが出力可能です。
LDOS(Local Density Of States)とは、ある原子の所属する電子の状態密度で、PDOSとは各軌道に所属する電子の状態密度の事です。すなわち、LDOSの和がPDOSとなります。VASPの結果出力ファイル「DOSCAR」には全原子のDOSと、PDOSが書かれています。
状態密度がわかれば、電子がエネルギーのどの区間をどのくらい占有しているかがわかるので、電気伝導に寄与する電子の数や比熱等の物理量を勘定したり、化学結合の様子等が知る事ができます。
このように、かなり重要な物理量なので、ファイルを読み解くついでに結果を整理するPythonスクリプトを自作しました。
#0. 前提
結果の中身について手動で整理する事に興味が無い方は、例えば、pymatgen「https://pymatgen.org/ 」やp4vasp「http://www.p4vasp.at/#/ 」を利用すれば中身をブラックボックスでDOSを出力する事ができます。
pymatgenについては公式のリファレンスやRKS WEBSITE様「http://ryokbys.web.nitech.ac.jp/pymatgen.html 」の記事が参考になります。
p4vaspはVASP公式サイト「https://www.vasp.at/wiki/index.php/CO_partial_DOS 」に方法が記載されております。
とは言っても、使い勝手が良いかと言われれば、正直微妙でして。本記事では「そんなに難しいフォーマットでは無いので、DOSCARをpythonで抽出して、わかりやすい形にする」を目標とします。
#1. DOSCARを得る
LDOSが含まれたDOSCARを得るには計算を2回(もしくは3回)走らせる必要があります。
##Step1.1回目の計算
計算に必要なファイルを準備します。
・ 「POSCAR」：
構造最適化後、又は計算したい構造を作成する等をし、構造が決まったファイルを用意します。
・「KPOINTS」：
計算の精度が収束する十分な数のk点を切ったファイルを用意します。
・「POTCAR」：
全電子を計算するため、PAW(Projector Augmented Wave)法のファイルを用意します。
・「INCAR」：
計算設定ファイルINCARに用いるそれぞれのタグの意味は「https://qiita.com/youkihashimoto3110/items/de92172e0b5e9f3872d3 」をご参考ください。
1回目の計算は電荷最適化(Electronic Relaxation)計算を新規で回します。構造最適化のパラメータは全てコメントアウトしておきましょう。
必須タグは
ISTART = 0
PREC = Accurate
EDIFF = 1E-5
LREAL = Auto
ALGO = VeryFast
です。なお、金属の場合は
ISMEAR = 1 #または2,-5
SIGMA = 0.5
絶縁体・半導体の場合は
ISMEAR = 0 #または-5
SIGMA = 0.01
を使用し、outcarのエントロピー項が原子当たり1meV未満になるようにSIGMA値を調整します。
絶縁体・半導体の場合にISMEAR>0を使用しないこと！
また、ISMEAR=-5(テトラヘドロン法)を使用すると精度が向上しますが、大抵の構造の場合適切ではありません。(???)
また、適宜「収束しない場合の有効措置」を使用します。
設定が完了したら1回目の計算を実行します。
##step2. 2回目の計算
計算が終了したら2回目の計算にはいります。2回目の計算では電荷最適化を継続ジョブで行うため、「INCAR」ファイルのみ書き換えます。
1回目の計算ファイルからの変更点は
ISTART = 1
ISYM = 0
LORBIT = 11
です。他は同じにして再計算しましょう。
また、細かくDOSのエネルギー刻み幅を細かくしたい場合は
EMIN = -10.0;EMAX = 17.0;NEDOS = 1001
を使用します。それぞれの意味はテンプレートに記載しています。
非自己無頓着な初期電荷配置を行うICHARG=11も有効です。
#2. DOSCARを読み解く
2回目の計算で出力された「DOSCAR」を任意のテキストエディタで開きます。
DOSCARの中身は下図ののようになっています。
見えづらくて申し訳ありません。（ルビ振りのため改行していますが、本物には改行はありません。）
VASP公式「https://www.vasp.at/wiki/index.php/DOSCAR 」や東北大学宮本先生「http://www.aki.che.tohoku.ac.jp/~taiko/vasp.html 」のHPにも情報がありますのでご参考ください。
また、上で示した例は計算でスピンを考慮していません。スピンを考慮した計算の場合、各スピンごとにも出力されます。
#3. pythonで結果を出力する
すなわち、
「固定文字列が6行あり、はじめに全DOSが。5行目と同じ文字列があり、１番目の原子のLPDOSが、5行目と同じ文字列があり、２番目の....」
の構造であるとわかります。また、空改行はなく、空白をスキップする必要があります。
これらを元にプログラミングを行うと以下になります。
import csv
atom_index = [1,2,3] #LPDOSを出力したい原子の番号(POSCARから何番目か)
doscar_name = "DOSCAR" #読み込むDOSCARの名前
save_csv_name = "doscar.csv" #DOSCARをcsvとして保存する名前
save_dos = "dos.csv" #全DOSの保存名
#DOSCARからリストに格納
with open(doscar_name, newline='') as doscar:
dos_read = csv.reader(doscar, delimiter=' ', skipinitialspace=True)
dos_list = [i for i in dos_read]
#DOSCARをCSVで出力
with open(save_csv_name, mode='w', newline='') as doscar_csv:
dos_write = csv.writer(doscar_csv)
dos_write.writerows(dos_list)
#DOSCARから全dosを抽出
i = 5 #読み込み開始インデックス
dos = []
count = 0
dos.append(dos_list[i])
while True:
i = i + 1
count = count + 1
dos.append(dos_list[i])
if dos[0] == dos_list[i+1]:
break
with open(save_dos, mode="w", newline='') as dos_csv:
dos_write = csv.writer(dos_csv)
dos_write.writerows(dos)
#DOSCARから指定したlpdosを抽出してそれぞれCSVに保存
lpdos_all = []
while True:
i = i + 1
pdos = []
count = 0
pdos.append(dos_list[i])
while True:
i = i + 1
count = count + 1
pdos.append(dos_list[i])
if len(dos_list)==i+1:
lpdos_all.append(pdos)
break
if len(pdos[count]) != len(dos_list[i+1]):
lpdos_all.append(pdos)
break
if len(dos_list)==i+1:
break
#lpdosを保存
for j in atom_index:
save_lpdos = "lpdos"+str(j)+".csv"
with open(save_lpdos, mode="w", newline='') as lpdos_csv:
lpdos_write = csv.writer(lpdos_csv)
lpdos_write.writerows(lpdos_all[j-1])
ひどいプログラムですが、読みにくくは無いと思います。
atom_index=[]にPOSCARの順番にある原子の番号を指定してあげると、その番号のLPDOSを"lpdos_番号.csv"として出力してくれます。あとはExcel等でプロットしてあげればそれぞれの原子、軌道におけるDOSをプロットすることができます。
また、空白を抜いて配列に格納する「#DOSCARからリストに格納」のメソッドはVASPの出力ファイルを整理する際に有効です。
上のDOSを用いた時、atom_index=[100,160,200]を指定した時の結果を示します。
これだけでどんな構造で、どんな特性をもつ物質なのか予想できますね。^^;