Abinitビルドの備忘録の続編。Abinitの使い方を書く。いつも大体書きかけ。
入力ファイル共通事項
Abinitの基本的な使い方は
abinit < files.files > log.log
mpirun -np 2 abinit < files.files > log.log #MPIを使う場合
という感じで、必要なファイル群が書かれたfiles.filesを標準入力として与えて実行する。もちろん入力するファイル名は他のものでもよい。
ただし、この与え方はabinit10で廃止予定の事項で、abinit9では事実上取り除かれているらしいが、一応中で色々やってくれているらしく、動く。一先ず、abnit10を使い始めるまでは、これを使う。
標準入力の中身
input.in #計算の指示を与える入力ファイルの名前
output.out #計算の出力が主に書き込まれる出力ファイルの名前
Hogei #計算にあたって、必要なファイル群(密度とか波動関数)のprefix
Hogeo #計算にあたって、出力されるファイル群(密度とか波動関数)のprefix
Hoge #何か必要。役目はあんまりよくわかってない。
/path/4/pseudopotential1 #一つ目の使う擬ポテンシャル/PAWのファイルのPATH
/path/4/pseudopotential2 #二つ目の使う擬ポテンシャル/PAWのファイルのPATH
/path/4/pseudopotential3 #三つ目の使う擬ポテンシャル/PAWのファイルのPATH
入力ファイル
上で書いたところのinput.inのこと。
# Crystalline Cu-FCC
#
# SCF calculation
kptopt 1 # Option for the automatic generation of k points,
# taking into account the symmetry
ngkpt 8 8 8 # 8x8x8 Brillouin zone sampling
prtden 1 # Print the density, for latter use
prtdos 1 # DoS with Gaussian smearing
prtwf 0 # No WFK print out
toldfe 1.0d-12 #
ecut 30.0 # Maximal plane-wave kinetic energy cut-off, in Hartree
occopt 7 # Specification of metallic occupation distribution, Gaussian
tsmear 0.0001 #Energy width for the metallic occupation disctribution
nstep 200 # Maximal number of SCF cycles
nband 8 # number of band to be prepared
#************************************************************************************
#* Generated by cif2cell 1.2.10 2020-11-14 1:31 *
#* T. Bjorkman, Comp. Phys. Commun. 182, 1183-1186 (2011). Please cite generously. *
#* *
#* () *
#* Van Ingen R.P. et al., J. Appl. Phys. 76, 1871-1883 (1994) *
#************************************************************************************
# Structural parameters
acell 3*6.8313598515
rprim 0.000000000000000 0.500000000000000 0.500000000000000
0.500000000000000 0.000000000000000 0.500000000000000
0.500000000000000 0.500000000000000 0.000000000000000
natom 1
ntypat 1
typat 1
znucl 29
xred 0.000000000000000 0.000000000000000 0.000000000000000
こんな感じで、[input variable] [value]の形式で書く。input variableについて、調べるときはAll Variables - abinitから調べる。
ここにはないけれど、Datasetというキーワードもあって、これを使うとSCF, DoS, Band, DFPTのような計算を一連のプロセスとして計算することができる。ただ、これは計算コストがものによって、違って、要求する計算資源も全然違うので、個人的には別々のプロセスで処理する方が、私は好み。
エネルギーのカットオフとか、Brillouin域のサンプリングはどうせ収束を見るので気にしなくてよいのだけれど、原子の配置とかは、間違えているのでどうしても期待した応えはでない。複雑な結晶の場合、原子配置は決して手で入力してはいけない。cif2cellことはじめに書いたように、cifから斜方格子での原子配置を書き出すアプリを使う。cif2cellはひとまずabinit形式の出力に対応しているので、それを使った。(cif2cellあんまり開発は進んでるような気がしないので、そのうち他のものに乗り換えないといけない気もしている)
SCF(Cu)
input.in #計算の指示を与える入力ファイルの名前
output.out #計算の出力が主に書き込まれる出力ファイルの名前
SCFi #計算にあたって、必要なファイル群(密度とか波動関数)のprefix
SCFo #計算にあたって、出力されるファイル群(密度とか波動関数)のprefix
SCF #何か必要。役目はあんまりよくわかってない。
./29-Cu.GGA.fhi #銅の擬ポテンシャル
SCFについていは、特に事前に必要なファイルはない。
# Crystalline Cu-FCC
#
# SCF calculation
ngkpt 8 8 8 # 8x8x8 Brillouin zone sampling
prtden 1 # Print the density, for latter use
prtdos 1 # DoS with Gaussian smearing
prtwf 0 # No WFK print out
toldfe 1.0d-12 # Energy tolerance for SCF
nstep 200 # Maximal number of SCF cycles
ecut 30.0 # Maximal plane-wave kinetic energy cut-off, in Hartree
occopt 7 # Specification of metallic occupation distribution, Gaussian
tsmear 0.0001 # Energy width for the metallic occupation disctribution
nband 8 # number of band to be prepared
#************************************************************************************
#* Generated by cif2cell 1.2.10 2020-11-14 1:31 *
#* T. Bjorkman, Comp. Phys. Commun. 182, 1183-1186 (2011). Please cite generously. *
#* *
#* () *
#* Van Ingen R.P. et al., J. Appl. Phys. 76, 1871-1883 (1994) *
#************************************************************************************
# Structural parameters
acell 3*6.8313598515
rprim 0.000000000000000 0.500000000000000 0.500000000000000
0.500000000000000 0.000000000000000 0.500000000000000
0.500000000000000 0.500000000000000 0.000000000000000
natom 1
ntypat 1
typat 1
znucl 29
xred 0.000000000000000 0.000000000000000 0.000000000000000
金属のSCFを計算するときに必要そうな入力をもれなく書いたもの。
Brillouin域サンプリングの収束性を確認する場合は、ngkptを大きな数にしていく。
DoSとバンド計算に使いたいので、prtdenを1にして、電子密度を書き出すようにしてある。特に意味はないけれど、計算も時間がかからないしファイル容量も小さいのでprtdosは1にして書き出すようにしてある。prtwfkは書き出されるファイルのサイズが大きいので、DFPTを後でやるときしか基本的には1にしないことにしている。
SCFの収束条件を決めるtoldfeは一応SCFの再現性を取るために小さめの値にしてある。単位胞のなかにたくさん原子、電子がある系とかだと1.0e-13とかにすると、精度の限界で最大反復回数nstepまで回ってしまうことがあるので、大体いつも1.0e-12くらいにしている。必要に応じて小さくする。
名前からわかる通りecutがエネルギーカットオフ。出来るだけ大きな値を取るのが良いけれど、小さい値でも大体それらしい答えが得られる(実空間基底を使う計算だとグリッドが粗いと全く話にならない固有値が出たりする)ので、初期は粗目のエネルギーカットオフを使うことが多い。Valence top/Fermi面がs,pなら10 Hartreeとか。d軌道の場合は30 Hartreeにすることが多い。もちろん最後に大きな値を入力して収束性を確認する。
金属の場合はoccoptを4,5,6,7のいずれかにする。絶縁体と違って、軌道エネルギーを各 $k$ 毎に電子数(の半分)だけ詰めるというアプローチができないので、それを指定するため。こうした占有数分布を導入する際は、付随してエネルギーの幅を決める必要があるので、その幅をtsmearで決める。occoptが4でFermi-Dirac分布の時は、これは温度に対応する。他の分布の場合は、一応エネルギーの幅くらいの意味しかない。
非整数値の占有数分布を導入したので、計算として事前に占有の可能性のある軌道の数を決めなければならない。それがnband。プリミティブセルで、スピン縮重でエネルギー幅が十分小さければ(0.01以下くらい)、セル内電子数の半分+1か+2で取っておけば十分のはず。エネルギー幅が0.1くらいになると、結構な数のバンドを用意しないといけない。あと、スーパーセルを使っている場合とかは、Fermi面近傍に縮退があったりするので、多めの数にしておいた方が良いはず。
acellは$\bf{a}_1, \bf{a}_2, \bf{a}_3$で作られる斜方格子の、それぞれのベクトルの大きさを指定する。三辺の長さがすべて同じときはサンプルのように3*10.0のように書くと、すべてに10が代入される。
rprimは$\bf{a}_1, \bf{a}_2, \bf{a}_3$を単位ベクトル(実は単位ベクトルにする必要はないけれど)にした際に、デカルト成分で表示したときの成分表示。かっこの中に書いたように、実際には単位ベクトルにする必要はなくて、とにかく以下のようになる:
$$
\begin{align}
\mathbf{a}_i = \mathrm{acell}[i-1] \times \left(\mathrm{rprim[i-1,0]}\mathbf{x} + \mathrm{rprim[i-1,1]}\mathbf{y} + \mathrm{rprim[i-1,2]}\mathbf{z}\right), i = 1,2,3.
\end{align}
$$
natomは原子数。ntypatは扱う原子の種類を指定する。typatは一次元配列で、各原子がどの原子種か指定する。znuclは原子番号。この計算ではCuのプリミティブセルの計算で、原子は一種類で一つしかセルの中に無いので、それぞれ、1, 1, 1, 39となっている。xredはプリミティブベクトル$\mathbf{a}_i, i=1,2,3$を単位とする相対座標を指定する。このケースではどこにおいても良いので、[0,0,0]が指定されている。
output.outに書かれた色々な量
力とか、バンドギャップとか、FFTに使ったグリッドの数とかその辺をどうやってみつけるか。
密度を可視化したい
abinitをビルドした時に同梱されているcut3Dをつかう。PATHはabinitの場所と同じ。
cut3Dを起動(コンソールで)と、入力するファイルを聞かれるので、例えばSCFo_DENのように密度のファイルを渡す。ここで、Tab補完は使えないので、注意。
この操作で、どの形式で書き出すか聞かれるので、ここではxcrysden用に*.xsfを指定する。
書き出すファイルの名前を聞かれて、座標を平行移動するかどうか聞かれるので、適当に回答して操作をすすめると、晴れてファイルが出来る。
この*.xsfはVESTAでも読み込めるらしい。
参考ABINIT 7 - RKS website
DoS
DoSの計算はSCFo_DENが必要。
バンド計算
DoSの計算はSCFo_DENが必要。
DFPT
DoSの計算はSCFo_WFKが必要。
Input variableのリスト(自分がよくつかうものを中心)
kptopt: Default: 4 if nspden == 4, 1 otherwise
chkprim: Default: 1
prtwf: Usually: 1
prtden: Usually: 1
prtdos: Default: 0
ecut: Default: None
toldfe: Default: 0.0
nstep: Default: 30
occopt: Default: 1
tsmear: Default: 0.01
nband: Default: None
ixc: Default: 1, PW-LDA
xc_tb09_c: Default: 99.99
chksymbreak: Default: 1
Libxcを使ったときの汎関数の指定の仕方
ixc = -XXXCCCのように、負の整数で与える。XXX, CCCはLibxc - a library of exchange-correlation functionals for density-functional theoryから選んだ交換汎関数と相関汎関数のid。二桁のものは0を先頭につけて、例えば012のようにする。
実は、交換汎関数の相関汎関数のidはかぶっていないので、ixc = -CCCXXXのように書いても正しいものが読み込まれる。また、交換相関汎関数がセットになっているものもあって、そうしたケースではその三桁のidをixc = -XCXのように書けばよい。
バンド計算のデータの取得
*_EIGから$k,\epsilon(b,k)$のデータに変換するためのpythonスクリプト(固有値の数が8の倍数以外の時のテストはまだしていない)
import numpy as np
import math
#Mathematical constants
pi = np.pi
tpi = 2.0*pi
fpi = 4.0*pi
zI = 1.0j
def abinitband(dirname):
fname = dirname + '/Bando_EIG'
f = open(fname,'r')
lines = f.readlines()
f.close()
Nlen = len(lines)
print(Nlen, ': No. of lines of the file, ',fname)
temp = lines[0].split()
Nk = np.int(temp[4])
print(Nk)
temp = lines[1].split(',')
temp = temp[1]
temp = temp.split()
Nb = np.int(temp[1])
Nbrow = Nb//8
if (Nb%8 != 0):
Nbrow += 1
print(Nb)
eigval = np.zeros([Nb,Nk],dtype = 'float64')
kred = np.zeros([Nk,3],dtype = 'float64')
for ik in range(Nk):
j = (Nbrow + 1)*ik + 1
temp = lines[j].split()
kred[ik,0] = np.float(temp[7])
kred[ik,1] = np.float(temp[8])
kred[ik,2] = np.float(temp[9])
ib = 0
for i in range(Nbrow - 1):
j = (Nbrow + 1)*ik + 2 + i
temp = lines[j].split()
eigval[i*8 + 0,ik] = np.float(temp[0])
eigval[i*8 + 1,ik] = np.float(temp[1])
eigval[i*8 + 2,ik] = np.float(temp[2])
eigval[i*8 + 3,ik] = np.float(temp[3])
eigval[i*8 + 4,ik] = np.float(temp[4])
eigval[i*8 + 5,ik] = np.float(temp[5])
eigval[i*8 + 6,ik] = np.float(temp[6])
eigval[i*8 + 7,ik] = np.float(temp[7])
j = (Nbrow + 1)*ik + 1 + (Nbrow - 1)
temp = lines[j].split()
for i in range(Nb-(Nbrow-1)*8):
eigval[(Nbrow-1)*8 + i,ik] = np.float(temp[i])
#
fname = dirname + '/output.out'
f = open(fname,'r')
lines = f.readlines()
f.close()
Nlen = len(lines)
print(Nlen, ': No. of lines of the file, ',fname)
for i in range(Nlen):
temp = lines[i]
fd1 = temp.find('R(1)')
fd2 = temp.find('R(2)')
fd3 = temp.find('R(3)')
if (fd1 > 0):
R1 = np.zeros(3)
temp = temp.split()
R1[0] = np.float(temp[1])
R1[1] = np.float(temp[2])
R1[2] = np.float(temp[3])
if (fd2 > 0):
R2 = np.zeros(3)
temp = temp.split()
R2[0] = np.float(temp[1])
R2[1] = np.float(temp[2])
R2[2] = np.float(temp[3])
if (fd3 > 0):
R3 = np.zeros(3)
temp = temp.split()
R3[0] = np.float(temp[1])
R3[1] = np.float(temp[2])
R3[2] = np.float(temp[3])
print(R1)
print(R2)
print(R3)
vol = np.dot(R1,np.cross(R2,R3))
G1 = np.cross(R2,R3)/vol*tpi
G2 = np.cross(R3,R1)/vol*tpi
G3 = np.cross(R1,R2)/vol*tpi
print(G1)
print(G2)
print(G3)
k = np.zeros([Nk],dtype = 'float64')
k[0] = 0.0
for ik in range(1,Nk):
dkred = kred[ik,:] - kred[ik-1,:]
dk = dkred[0]*G1 + dkred[1]*G2 + dkred[2]*G3
k[ik] = k[ik-1] + np.linalg.norm(dk)
return k, eigval
dirnameに該当ファイルのあるpathを文字列で代入して、k, eigval = abinitband(dirname)でデータが得られる。
トラブルシューティング
abinit 8.10.1でmetaGGAを汎関数として呼んだときにこける
こけた際に、以下のメッセージが出ている場合
--- !ERROR
src_file: m_vtowfk.F90
src_line: 1138
mpi_rank: 4
message: |
The eigenvector with band 1 has zero norm.
This usually happens when the number of bands (nband) is comparable to the number of planewaves (mpw)
Action: Check the parameters of the calculation. If nband ~ mpw, then decrease nband or, alternatively, increase ecut
...
なるエラーで止まることがある。一見指示通りにパラメーターを設定すれば解決できるような気がするが、これはそもそもabinitで運動エネルギー密度を評価してないため。以下のキーワードを入力ファイル(input.inとする)に足せば、期待通りに動く。一応公式マニュアルixcにも最後の方にきちんと書いてある。
usekden = 1 # For mGGA