Siは半導体であることはよく知られています。 このSiを使って、エネルギーバンドを計算方法を習得し、PWscfの基本的な使い方を身に付けます。 また、対称点の指定を行うことで、バンド分散をどのように描くのかについても勉強します。 最後に、実空間における電荷密度の分布の計算も行います。

計算を実行するにあたり、事前に集めておく情報は、

• 結晶のブラべー格子の種類と格子定数
• 結晶の単位胞に含まれる原子の種類・数・原子量・位置

です。 原子の種類に応じて、擬ポテンシャルのデータを含んだファイルを別途用意する必要があります。

単体のSiは最もよく知られている半導体で、その結晶構造はダイヤモンド構造です。 Siの化合物は岩石の成分として含まれているので、身の回りに最もありふれた元素の一つとも言えます。 現在のエレクトロニクスの成功は、 Si単体の結晶が高純度（99.999999999%）で作成できることもさるものながら、ハートリー・フォック近似を用いる第一原理計算によってその物性が計算できるところも大きいと考えられます。

以上のことから、Siは第一原理計算の入門の題材として適した素材であるといえます。

まず作業ディレクトリSiを作成し、そこで作業することにします。

$ mkdir Si
$ cd Si

Quantum ESPRESSOは擬ポテンシャル法に基づく第一原理計算ソフトなので、電子状態の計算には擬ポテンシャルのデータ・ファイルが必要です。 擬ポテンシャルは、その作成の方法によっていくつかの種類が用意されていますが、今回はSi.pz-vbc.UPFをダウンロードします。

• 今回は単体なので擬ポテンシャルのファイルの種類はあまり問題になりませんが、複数の原子を使う化合物だと、いろいろな原子の擬ポテンシャルのファイルを組み合わせる必要が出てきます。
• その際に、WIEN2kのような全電子近似を用いた第一原理計算の結果を使って、本当に正しいバンドになっているかを確認します
• それならQuantumESPRESSOいらんやんけ…という話になりますが、バンド計算をしたあとの処理に関しては擬ポテンシャル法の方が高速に計算できるので向いています。
• バンド分散、状態密度、フェルミ面、電荷密度分布などの一粒子状態に基づく量を計算したい場合はWIEN2kの方が便利です。
• また、角度分解光電子分光やde Haas-van Alphen効果、Shubnikov-de Haas効果などの実験から推定されるバンド構造やフェルミ面の情報を元に、最適な擬ポテンシャルの組み合わせを探すこともあります。

ディレクトリSiの中に次のファイルを作成します（scfはself-consistent fieldの略）。

Si.scf.in

&control
  prefix      = 'Si',
  calculation = 'scf',
  pseudo_dir  = './pseudo/',
  outdir      = './tmp/',
/
&system
  ibrav = 2, celldm(1) = 10.26, nat = 2, ntyp = 1,
  ecutwfc = 12.0,
/
&electrons
  conv_thr = 1.0d-8
/
ATOMIC_SPECIES
 Si 28.086 Si.pz-vbc.UPF
ATOMIC_POSITIONS
 Si 0.00 0.00 0.00
 Si 0.25 0.25 0.25
K_POINTS automatic
 6 6 6 1 1 1

入力ファイルの各パラメータは、Quantum ESPRESSOのインストール・ディレクトリの中にある PW/Doc/INPUT PW.html というファイルを参照すれば書いてあります。 もしパラメータを指定しなければ、デフォルトの値が使われます。 今回は入門編ということで、最低限、指定しなければいけないもののみ指定します。 指定する順番は任意です。

計算データの入出力に関する情報を指定します。

結晶構造に関する情報を指定します。 以下のパラメータはすべて初期値が設定されていませんので、ユーザーが必ず入力する必要があります。

実際の計算にあたってのオプションです。

使用する原子の種類、原子量、擬ポテンシャルのファイル名を指定します。

alat (celldim(1)あるいはAのこと)を1としたときの原子の単位胞内における位置をデカルト座標で指定します。

計算に使うk点を指定します。

• automaticオプションにより自動的に逆格子を6×6×6に分割したものを使います (Monkhorst-Pack grid)。分割の数が多ければ多きほど計算が正確になりますが、それだけ計算に時間がかかるようになります。
• 後ろの3つの数字は、分割した点のk点を使うなら0、分割した点と点の中間点を使うなら1を指定します（つまり0を指定したらガンマ点が含まれますが、1を指定すると含まれません）。
  • 半分だけずらした方が、少ないk点数でより精度が高くなることが知られています

端末からpw.xを実行します。

$ pw.x < Si.scf.in > Si.scf.out

実行したら、計算の途中経過が出力された Si.scf.out というファイルと、電荷密度などの計算結果のファイルが入った tmp/Si.save/ というディレクトリができます。 計算に失敗したら crash という名前のファイルもできるので、それを開いて原因を確認します。 計算が収束する様子は Si.scf.out の estimated scf accuracy を見ることで確認できます。

$ grep -e 'estimated scf' Si.scf.out
estimated scf accuracy < 0.04383890 Ry
estimated scf accuracy < 0.00208552 Ry
estimated scf accuracy < 0.00003904 Ry
estimated scf accuracy < 0.00000224 Ry
estimated scf accuracy < 0.00000008 Ry

Siは半導体ですので Si.scf.out にフェルミ準位が表示されません。 その代わりに、価電子帯のトップのエネルギーが表示されます。 あとのバンド計算で必要ですので確認しておきます。

$ grep -e 'highest' Si.scf.out
highest occupied level (ev): 4.4741

以下の手順で入力ファイルを作成します。

1. Si.scf.in をコピーして Si.bands.in というファイル名にします。
2. &control の calculation を scf から bands にします。
3. &systemにnbnd=8という記述を加えます。nbndは計算するバンドの数です。
   • Siの電子配置は[Ne]3s23p2なので、価電子は4つあります。
   • 単位胞に2つのSi原子が含まれるので、価電子帯のバンドの数は8本になります。
   • ただし、非磁性の場合の計算はエネルギーはスピンについて縮退しているので、バンドの数は結局4本になります。
   • nbndは指定しなければ自動的に価電子帯のバンドの数（今の場合は4）になります。
   • バンド計算の際は、伝導帯のバンドも見えないと不便なので、多めにバンドの数をとってやるわけです。
   • 今回は、3sと3pのすべてのバンドが見えるように8としました。
4. K_POINTS にエネルギーを計算するk点を指定します。

最終的にファイルは次のようになります。

Si.bands.in

&control
  prefix = 'Si',
  calculation = 'bands',
  pseudo_dir = './pseudo/',
  outdir = './tmp/',
/
&system
  ibrav = 2, celldm(1) = 10.26, nat = 2, ntyp = 1,
  ecutwfc = 12.0, nbnd = 8
/
&electrons
  conv_thr = 1.0d-8
/
ATOMIC_SPECIES
 Si 28.086 Si.pbe-rrkj.UPF
ATOMIC_POSITIONS
 Si 0.00 0.00 0.00
 Si 0.25 0.25 0.25
K_POINTS
 3
 0.0 0.0 0.0 1
 1.0 0.0 0.0 2
 0.5 0.5 0.5 3

K_POINTSの下にある数字の3は、3つのk点でエネルギーを計算することを指定します。 その下に各k点の座標と、1から3の番号付けを書きます。 座標はデカルト座標で、2π/aを単位とします。

再び端末からPWscfを実行します。

$ pw.x < Si.bands.in > Si.bands.out

プログラムの実行により作成された Si.bands.out から指定したk点におけるエネルギーがわかります。

Si.bands.out

（中略）
          k = 0.0000 0.0000 0.0000 (   229 PWs)   bands (ev):

    -6.3225   4.5875   4.5875   4.5875   6.3731   7.0827   7.0827   7.0827

          k = 1.0000 0.0000 0.0000 (   222 PWs)   bands (ev):

    -2.6903  -2.6903   2.0923   2.0923   5.5112   5.5112  13.3516  13.3516

          k = 0.5000 0.5000 0.5000 (   228 PWs)   bands (ev):

    -4.3138  -1.7254   3.5130   3.5130   5.5922   7.9810   7.9810  11.8807
（中略）

計算がうまく実行できることがわかったら、今度は逆格子空間で対称点を結んだ線上のバンド構造を計算しましょう。 Quantum ESPRESSOにはバンド構造を描くプログラムが含まれているので、それを利用します。

Si.bands.inを開いて、バンドを描く経路を指定します。 例としてK_POINTSで

• L点2π/a(1/2, 1/2, 1/2) → Γ点2π/a(0, 0, 0) → X点2π/a(0, 0, 1) → K点2π/a(0, 1, 1) → Γ点

という経路のk点の座標を入力します。 ここでaは格子定数です。

Si.bands.in

&control
  prefix = 'Si',
  calculation = 'bands',
  pseudo_dir = './pseudo/',
  outdir = './tmp/',
/
&system
  ibrav = 2, celldm(1) = 10.26, nat = 2, ntyp = 1,
  ecutwfc = 12.0, nbnd = 8
/
&electrons
  conv_thr = 1.0d-8
/
ATOMIC_SPECIES
 Si 28.086 Si.pbe-rrkj.UPF
ATOMIC_POSITIONS
 Si 0.00 0.00 0.00
 Si 0.25 0.25 0.25
K_POINTS tpiba_b
 5
 0.5  0.5  0.5   5
 0.0  0.0  0.0  10
 0.0  0.0  1.0  10
 0.0  1.0  1.0  10
 0.0  0.0  0.0   0

• tpiba_bはバンドをプロットするためのk点を自動生成するオプションです。
• 2π/aを単位としたデカルト座標で指定します。
• 各行の4つ目の数字は、次のk点までの間に何個点を取るかを指定します（最終行の4つ目の数字は無視されます）。

作成できたら、PWscfを使ってエネルギー分散を計算します。

$ pw.x < Si.bands.in > Si.bands.out

次に、バンド構造をプロットするために必要なファイルを作成します。

bands.in

&bands
  prefix = 'Si',
  outdir = './tmp/',
  filband = 'bands.dat',
/

prefix や outdir はさきほど指定したものと同じものにします。 これを bands.x というプログラムで実行します。

$ bands.x < bands.in > bands.out

bands.outとbands.datが生成されたら成功です。 これらのファイルを開いて中身にエラーが出ているようなら、入力ファイルに間違いがないか確認してやり直しましょう。

さて、以上で得られたデータを元に、バンド分散をplotband.xというプログラムを使って描いてみます。

指示に従って入力します。

$ plotband.x
    Input file > bands.dat
Reading 8 bands at 36 k-points
Range: -5.6940 16.4680eV Emin, Emax > -7.0 14.0
high-symmetry point: 0.5000 0.5000 0.5000 x coordinate 0.0000
high-symmetry point: 0.0000 0.0000 0.0000 x coordinate 0.8660
high-symmetry point: 0.0000 0.0000 1.0000 x coordinate 1.8660
high-symmetry point: 0.0000 1.0000 1.0000 x coordinate 2.8660
high-symmetry point: 0.0000 0.0000 0.0000 x coordinate 4.2802
output file (xmgr) > Si.bands.xmgr
bands in xmgr format written to file si.bands.xmgr
output file (ps) > Si.bands.ps
Efermi > 4.4741
deltaE, reference E (for tics) 1.0 4.4741
bands in PostScript? format written to file si.bands.ps

• Range … どこまでのエネルギーの範囲をプロットするか、エネルギーの最高値と最低値の表示を参考に決めます。
• output file (xmgr) … Graceというグラフ描画ソフト用の出力ファイルの名前を指定します。 また、output file (ps)は、Postscript用のファイルです。
• Efermi … フェルミ・エネルギーを指定しますが、今の場合は半導体なので、さっき調べた価電子帯のトップのエネルギーを入力します。
• deltaE … グラフのエネルギーの刻み幅
• reference E … 刻みの基準となるエネルギーの値

実は、バンドの経路はXCrysDenを使って指定すると非常に楽です。

1. XCrysDenを起動し、[File]→ [Open PEscf …] → [Open PWscf Input File] から Si.scf.in を開きます。
2. [Tools] → [k-path Selection]を選ぶことで、ブリルアンゾーンから経路を選ぶことができます。
3. 選び終えたら[OK]を押し、“Total number of k-points along the path”でk点数を指定した後、“File of type”でpwscfを選んで適当な名前をつけて保存します。
4. 得られたk点の経路を Si.bands.in の K_POINTS に貼り付けます。

ミラー面[1, −1, 0]で結晶を切り、そこでの電荷の分布を調べます。

電荷密度の計算には pp.x というプログラムを使用します。 まず、以下のようなファイルを作成します。

Si.pp_rho.in

&inputpp
  prefix = 'silicon',
  outdir = './tmp/',
  filplot = 'sicharge',
  plot_num= 0,
/
&plot
  nfile = 1
  filepp(1) = 'sicharge'
  weight(1) = 1.0
  iflag = 2
  output_format = 2
  fileout = 'Si.rho.dat'
  e1(1) = 1.0, e1(2) = 1.0, e1(3) = 0.0,
  e2(1) = 0.0, e2(2) = 0.0, e2(3) = 1.0,
  nx=56, ny=40
/

pp.xはまず&inputppを読み込んでプログラムを実行して必要な情報を計算した後、&plotを読み込んでグラフにプロットできる形式にします。 どちらか片方の計算を行いたい場合は、そのブロックのみを残して他は消してください。

入力ファイルと出力ファイルの設定を行います。

電荷密度をプロットするのに必要な部分です。

最終的な結果は \begin{equation} \sum^{\mathrm{nfile}}_{n=1} \mathrm{weight}(n) \times \mathrm{rho}(n) \end{equation} で表されます。 詳しくは PP/Doc/INPUT_PP.html を参照してください。

プログラムの実行は次のように行います。

$ pp.x < Si.pp_rho.in

すると sicharge と Si.rho.dat とというファイルができます。 得られた電荷密度は plotrho.x というプログラムを使って次のようにプロットします。

$ plotrho.x
    Input file > Si.rho.dat
r0 : 0.0000 0.0000 0.0000
tau1 : 1.0000 1.0000 0.0000
tau2 : 0.0000 0.0000 1.0000
read 2 atomic positions
output file > Si.rho.ps
Read 56 * 40 grid
Logarithmic scale (y/n)? > n
Bounds: 0.002452 0.075120
min, max, # of levels > 0 0.08 6

するとSi.rho.psというファイルができるので、電荷密度を確認します。

Si原子同士の間で電荷密度が大きくなっていて、Si原子間の共有結合の様子がわかります。