VASPで第一原理計算
1 VASPについて
この文章は研究室の学生がVASPを使って材料研究を行うためのメモです。VASPの使える計算機は、研究室にある並列計算機、千葉大の高速演算サーバ、研究室で借りている九大のITOを使うことができます。
VASPはPAW法で平面波第一原理計算プログラムです。3,000ユーロで大窪がライセンスを購入しており、研究室の学生は自由に使うことができます。千葉大に設置されている高速演算サーバ(IBM製SR24000)上でpatchを当 てたVASP5.4.1(VASP.5.4.1.05Feb16)を大窪がコンパイルして実行できること を確認しています。いろいろソースコードを修正せねばコンパイルが通らず数 日ハマりました。研究室で借りているITOではコンパイル済みの5.4.4使えます。
- SR24000上でVASP5.4.1をコンパイルするためのpatch patch.SR24000.5.4.1
- VASP5.4.1コンパイル用のmakeファイル makefile.include.5.4.1
このページを見たナノサイエンス学科の張潤南くんがVASP 5.4.4をコンパイルするためのパッチとmakefile.includeを送ってくれました。少し変更してvaspの配布しているpatchも含めてコンパイルするためのpatchを作りました。
(PATHは/usr/binを最初に参照するようしてください。またgmakeを使ってください)
export PATH=/usr/bin:$PATH gmake std
fftw3とlapackは自前でコンパイルしておく必要があります。
1.1 学習のコツ
すべてのパラメータの細かい背景は理解する必要はありませんが、計算条件が計算結果(電子密度、バンド、DOS…)に与える影響を理解しておく必要があります。VASPのマニュアルにある以下のURLがよくまとまっているので、計算したい物理量を精度よく計算するためにどのような注意点があるか理解しておくことが必要です。VASPのマニュアルはwikiで良くまとめられています。
http://cms.mpi.univie.ac.at/wiki/index.php/The_VASP_Manual
第一原理計算やPAW法の理論は計算しつつ学べば良いかと思います(本当はよくありませんが…応用化学科だし?)。まずは計算を流して出力結果を理解しながら眺めつつ、計算の素人に何か聞かれて答えられるレベルまで知識を身につけてください。大窪の大学院無機材料化学の講義も役に立ちます。
2 入力ファイル
ファイル名は固定で以下の4つのファイルが入力ファイルになります。計算毎にディレクトリを掘ってファイルを配置します。
- INCAR (計算条件)
- POSCAR (原子座標と結晶系)
- POTCAR (擬ポテンシャル)
- KPOINTS (Kサンプリング)
2.1 INCARの例
計算の種類や計算条件を指定するファイルです。自分が行う計算で使っている中身のtagを理解してください。1ノード(20コア=NCORE*NPAR)で計算を行ういくつかの例を示します。各tag(オプション)は、VASPのマニュアルやwikiを参照してきちんと意味を理解してください。
2.1.1 DOSの計算
# マニュアル: http://cms.mpi.univie.ac.at/vasp/vasp/vasp.html SYSTEM = HOGE NCORE = 4 # sqrt(COREs) NPAR = 5 # NCORE*NPAR=CORES NWRITE = 1 # OUTCARへの出力の量 [1] LREAL = Auto # # 精度 Low|Medium|High|Normal|Accurate| ENCUTを自動設定 PREC = Accurate # restart 0:from_scratch 1:restart from WAVECAR 3:MD ISTART = 0 ICHARG = 1 # Electronic Relaxation # ENCUT = 450 # cutoff: grep [EATOM POTCARした最大値] ALGO = Normal # Normal | VeryFast | Fast | EDIFF = 1E-4 # SCFの収束条件eV [1E-4] #バンド数 Non-spin:NELECT/2 + NIONS/2 Spin:0.6*NELECT+NMAG # NBANDS = 50 # 計算するバンドの数を設定(価電子の半分くらい) ISPIN = 2 # スピンを考慮(1:しない 2:する) [1] # PDOSを計算したいとき LORBIT=11
2.1.2 構造最適化の計算
SYSTEM = HOGE NCORE = 4 # sqrt(COREs) NPAR = 5 # NCORE*NPAR=CORES NWRITE = 1 # OUTCARへの出力の量 [1] LREAL = Auto # # 精度 Low|Medium|High|Normal|Accurate| ENCUTを自動設定 PREC = Accurate # restart 0:from_scratch 1:restart from WAVECAR 3:MD ISTART = 0 ICHARG = 1 # Electronic Relaxation # ENCUT = 450 # cutoff: grep [EATOM POTCARした最大値] ALGO = Normal # Normal | VeryFast | Fast | EDIFF = 1E-4 # SCFの収束条件eV [1E-4] #バンド数 Non-spin:NELECT/2 + NIONS/2 Spin:0.6*NELECT+NMAG # NBANDS = 16 ISPIN = 1 # スピンを考慮(1:しない 2:する) [1] # Ionic relaxation EDIFFG = 1E-3 # 構造緩和でのエネルギーの収束条件 (負で原子に働く最大の力) NSW = 200 # 構造最適化の場合はNSW=200にしてIBRION = 2にする IBRION = 1 # 0:MD 1:RMM-DIISで最適化 [-1](if:NSW>0 or -1) ISIF = 3 # 構造最適化の種類 # |------+-------+------------+--------+------------+-------------| # | ISIF | cal | cal | change | change | | # | | force | stress ten | ions | cell shape | cell volume | # |------+-------+------------+--------+------------+-------------| # | 0 | yes | no | yes | no | no | # | 1 | yes | trace only | yes | no | no | # | 2 | yes | yes | yes | no | no | # | 3 | yes | yes | yes | yes | yes | # | 4 | yes | yes | yes | yes | no | # | 5 | yes | yes | no | yes | no | # | 6 | yes | yes | no | yes | yes | # | 7 | yes | yes | no | no | yes | # |------+-------+------------+--------+------------+-------------|
2.2 POSCAR
単位胞、原子の種類、原子の座標を指定するファイルです。
VESTAでVASP用のPOSCARファイルを出力できます。スーパーセルを作成したい場合は、自作プログラムの cif2poscar.py を使えば良いでしょう。元素置換や欠陥、表面の計算を行いたい場合は、VESTAで確認しつつPOSCARをマニュアルで編集すれば良いでしょう。
2.3 POTCAR
擬ポテンシャルのデータが入ったファイルです。
1つの原子で複数の擬ポテンシャルファイルがあります。どの軌道の電子まで 陽に取り扱うか、電子相関の関数でファイルが異なります。
どの擬ポテンシャルを使ったが後で分からなくなるので、以下のような簡単なシェルスクリプトを作成してPOTCARを作成することをおすすめします。ここで、~/projects/vasp/ppに擬ポテンシャルファイル一式を置いていることとしています。POSCARの原子の並びとPOTCARのデータの並びば揃える必要があります。
#!/bin/sh
DIR=~/projects/vasp/pp/pbe
cat ${DIR}/O/POTCAR\
${DIR}/Sr_sv/POTCAR\
${DIR}/Ti/POTCAR \
> POTCAR
2.4 KPOINTS
K点を分割する方法と数を指定するファイルです。
分割サイズがa, b, cで同じようにします(例えばb=2a=4cなら、k=4, 2, 1)。
Automatic mesh # コメント行 0 # 0とすると自動でメッシュを切る Auto # 格子定数で均等に分割 20 # nx, ny, nz = 20/a+0.5, 20/b+0.5, 20/c+0.5に自動セット
3 出力ファイル
いろいろファイルが出力されるが主なデータが入っているファイルは以下のとおり。
- OSZCAR (SCFサイクルでのエネルギーと力)
- CHGCAR (電子密度、 CHGCAR.vaspにrenameすればvestaで開くことができる)
- XDATCAR (構造最適化の場合に出力された原子座標、POSCARと同じフォーマット)
- DOSCAR (DOSのデータ)
- CONTCAR (POSCARと同じフォーマット。構造最適化等で生成される最終ステップの構造)
- WAVECAR (リスタート用の波動関数、バイナリファイル)
4 SR24000上で実行
4.1 計算サーバ
巨大なメモリを必要とするので、千葉大の高速演算サーバ(SR24000)で計算を行います。千葉大生は、アカウント申請すればSR24000を無料で使えます。各自で申請書を作成して提出してください。数日後に学内便でアカウント情報が来ます。
https://www.imit.chiba-u.jp/new-system/sr/index.html
計算資源は18ノード(1ノード20コア)です。100個ぐらいの原子だったら1ノードの計 算資源で十分です。複数ノードの計算は、jobがなかなか流れないので、基本1 ノードの計算を行ってください。2017年8月現在、ジョブはスタックしておらず空いています。
silicon上でもVASPの実行は可能ですが、原子が50個を越えるとキツいです。
ディスククォータを越える場合、以下のコマンドで作業ディレクトリを作成できます。ファイルの最終アクセスから30日までファイルを保存されます。
/usr/local/bin/mkshortdir
/short以下にユーザ名のディレクトリが生成されます。
4.2 シェルの初期設定
デフォルトのログインシェルはcshです。cshの操作になれてないのでbashを使うのが良いかと思います。参考に大窪の.bashrcを晒します。
export TERM=vt100 export LANG=C export LC_ALL=C export PATH=/opt/freeware/bin/:$PATH export PATH=/opt/IBM/xlf/15.1.0/bin:$PATH export PATH=$HOME/bin:$PATH export PS1="\u@sr26k[\w] " alias ls='ls -F --color' alias less='less -fr'
ログインシェルを変更することは許されていないようなので、ログインしたらbashを立ち上げてください。
4.3 計算の実行
scp等でファイルをコピーした後、sshログインして以下のスクリプトでジョブを投げます。notify_userのメールアドレスは、自分のものに変えてください。~/bin/vasp_stdにvaspの実行ファイルを配置しています。cshなので、shと文法がちょっと違います。
1ノード(20コア)で計算する場合の例です。R2クラスは1ノード(20コア)を使う場合、最大 で24時間の計算時間を確保ができます。計算時間を計測したい場合は、VASPを実行するコ マンドの頭にtimeを追加してください。
#!/bin/csh -f #@class=R2 #@job_type=parallel #@output=$(jobid).out #@error =$(jobid).e #@node = 1 #@network.MPI=sn_all,,US,,instances=2 #@bulkxfer=yes #@tasks_per_node=20 #@resources=ConsumableCpus(1) #@notify_user=hoge@faculty.chiba-u.jp #@queue setenv MEMORY_AFFINITY MCM setenv MP_PROCS 20 time /usr/bin/poe ~/bin/vasp_std # POSCAR.basefileの2行目(セル定数の倍数)を変えて次々計算する場合 # エネルギーとセル定数の関係を調べる # set datfile = "relax.dat" # set base = "POSCAR.basefile" # echo "# " `date` >> relax.dat # foreach i (0.94 0.96 0.98 1.00 1.02 1.04 1.06) # sed -e "2,2s/.*/$i/" $base > POSCAR # /usr/bin/poe ~/bin/vasp_std # echo $i `tail -n1 OSZICAR` >> relax.dat # end
4.4 ジョブの操作
以下のコマンドを覚えておいてください。
- llsubmit hoge.sh (ジョブを投げる)
- llq (jobを確認)
- llcancel job-number (job-numberのjobを削除)
4.5 計算の実行例
40原子SrTiO3(2x2x2のスーパーセル)のインプットファイル一式は以下のとおりです。
SrTiO3_example.tgz
標準出力の計算結果です。
running on 20 total cores
distrk: each k-point on 20 cores, 1 groups
distr: one band on 4 cores, 5 groups
using from now: INCAR
vasp.5.4.1 05Feb16 (build Jul 15 2017 15:23:34) complex
POSCAR found type information on POSCAR O Sr Ti
POSCAR found : 3 types and 40 ions
LDA part: xc-table for Pade appr. of Perdew
POSCAR, INCAR and KPOINTS ok, starting setup
FFT: planning ...
WAVECAR not read
WARNING: chargedensity file is incomplete
entering main loop
N E dE d eps ncg rms rms(c)
DAV: 1 0.259834868679E+04 0.25983E+04 -0.97309E+04 14360 0.119E+03
DAV: 2 0.952461291265E+02 -0.25031E+04 -0.22724E+04 14475 0.353E+02
DAV: 3 -0.365028240917E+03 -0.46027E+03 -0.44447E+03 15375 0.157E+02
DAV: 4 -0.386812834690E+03 -0.21785E+02 -0.21541E+02 18130 0.370E+01
DAV: 5 -0.388212813649E+03 -0.14000E+01 -0.13966E+01 25180 0.788E+00 0.102E+02
DAV: 6 -0.311631880589E+03 0.76581E+02 -0.40559E+02 19075 0.456E+01 0.435E+01
DAV: 7 -0.318764358063E+03 -0.71325E+01 -0.77346E+01 19170 0.218E+01 0.392E+01
DAV: 8 -0.317555990587E+03 0.12084E+01 -0.11857E+01 18720 0.111E+01 0.617E+00
DAV: 9 -0.317597565974E+03 -0.41575E-01 -0.13406E+00 19765 0.255E+00 0.154E+00
DAV: 10 -0.317594272460E+03 0.32935E-02 -0.10754E-01 20260 0.699E-01 0.101E+00
DAV: 11 -0.317598447248E+03 -0.41748E-02 -0.25612E-02 18425 0.361E-01 0.592E-01
DAV: 12 -0.317600284385E+03 -0.18371E-02 -0.33105E-03 19480 0.119E-01 0.177E-01
DAV: 13 -0.317600525326E+03 -0.24094E-03 -0.15128E-03 16145 0.876E-02 0.936E-02
DAV: 14 -0.317600614784E+03 -0.89458E-04 -0.24963E-04 12100 0.619E-02
1 F= -.31760061E+03 E0= -.31761243E+03 d E =0.354493E-01 mag= 0.0046
writing wavefunctions
0.0u 0.0s 10:47 0% 0+48k 0+0io 1pf+0w
14回のSCF計算で収束し、エネルギーは-317.601 eV、計算時間は10分47秒でした。
5 自作プログラム
インプットファイルの作成やジョブの投入、出力結果の解析のために、pythonとCでツールを自作しています。プログラムは学内からのみアクセス可能です。-hでオプションのヘルプを参照できます。ほぼ研究室で行っている研究のためのツールです。汎用的なツールではありません。似たようなツールはweb上にたくさん転がっていますが、他人が作ったものの使い方を覚えるより自分で理解して作った方が早いので自作しています。
5.1 入力/実行支援
5.1.1 cif2poscar.py
対称性を消した(P1)cifファイルから、スーパーセルを作ってPOSCARを生成するプログラム。
5.1.2 make_vaspin.py
POSCARとPOTCARからvaspで計算を行うための様々なINPUTファイルのテンプレートを生成するプログラム。実際は、シェルで各種計算用のファイルをINPUTにコピーしつつ restart_vasp.py でコピーしながら計算を進めるようにして使います。また候補となる擬ポテンシャルファイルを標準出力に出力します。
5.1.3 restart_vasp.py
backupディレクトリを作成してINCAR, OUTCAR, XDATCAR, POSCARの後ろに3つの数字をつけて保存し、CONCARをPOSCARにコピーするプログラム。3つの数字は1つずつincrementします。submitするシェルの中で呼び出して、入力と出力を保存しながらINPUTを更新して連続してvaspを実行するためのものです。
5.1.4 incar2submit.py
ディレクトリにあるINCAR.*をINCARにコピーしながらITO上で連続してsubmitするためのシェルスクリプトを生成するプログラム。 restart_vasp.py で各ステップでの入力と出力を保存しながら計算するようなシェルスクリプトを生成します。
5.2 出力/解析支援
5.2.1 dosplot.py
DOSCARを読んでプロットするプログラム。オプションを与えれば選択した元素毎や1原子毎、または軌道毎にプロットできる。
5.2.2 relax2minimum.py
格子定数等を変えて評価したエネルギー変化を多項式フィッテングすることで、極値を求めるプログラム。
5.2.3 chgcar2line.py
chgcar(電子密度)ファイルから指定した分率座標1点(-c)に基づいて任意の方向(-m x, y, z)に走査した電子密度のラインデータを取得してプロットするプログラム。複数のファイルを渡せば、同じラインでデータをプロットします。任意の点の最近接の点でをサンプルする場合(-n 0)と任意の点の内挿(-n 500とか)によりラインデータを生成することが選択できます。
5.2.4 xdatcar2trans
XDATCARをlammpstrjに変換するプログラム。またステップ毎の格子定数もファイルに出力します。C言語で作成しているので、makeして実行ファイルを作成してください。
5.2.5 outcar2plot.py
relax計算等で出力した複数ステップのOUTCARからステップ毎のenergyや格子定数、全原子の中で最大の力(構造最適化の収束条件)等を拾ってプロットしてくれます。