VASPの出力ファイルDOSCARからLDOSとPDOSを整理するスクリプト

#イントロ
お疲れさまです。橋本です。

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の中身は下図ののようになっています。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/571701/28a1e49e-f8ea-3c11-bb29-9dd839cc8d52.png)
見えづらくて申し訳ありません。（ルビ振りのため改行していますが、本物には改行はありません。）
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行目と同じ文字列があり、２番目の....」
の構造であるとわかります。また、空改行はなく、空白をスキップする必要があります。

これらを元にプログラミングを行うと以下になります。

```python:Extract_dos.py
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]```を指定した時の結果を示します。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/571701/cc36d76f-b40a-fb7d-9143-a27299860d86.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/571701/87be13a3-9bd3-f4b1-47a7-c4d448ba2d43.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/571701/5c14e947-a8d3-7a3c-f6c7-4e8bd53f6a23.png)
これだけでどんな構造で、どんな特性をもつ物質なのか予想できますね。^^;