目次

前のトピックへ

2. インストール

次のトピックへ

4. 収束のチェック

このページ

3. グラフェンのバンド構造と状態密度

まず、簡単な計算例として、グラフェンのバンド構造と状態密度を求めてみましょう。 どちらも scf計算->non-scf計算->後処理という流れになります。

3.1. scf計算

何を計算するにしても、まず scf 計算により電荷密度(charge density)を求める必要があります。 scf 計算に必要なファイルはインプットファイル graphene.scf.in は以下の通りです。 とりあえずこのファイルをダウンロードして、適当な directory において下さい。

&control
 calculation = 'scf'
 prefix='graphene',
 tstress = .true.
 tprnfor = .true.
 pseudo_dir = './',
 outdir='./work/'
 disk_io='low'
 wf_collect=.true.
/
&system
 ibrav = 4,
 celldm(1) = 4.602,
 celldm(3) = 4,
 nat = 2,
 ntyp = 1,
 ecutwfc = 30.0,
 ecutrho = 150.0,
 occupations = 'smearing'
 smearing = 'm-p'
 degauss = 0.01
/
&electrons
 mixing_beta = 0.7
 conv_thr = 1.0d-8
/
ATOMIC_SPECIES
C  12.0107  C.pz-van_ak.UPF
ATOMIC_POSITIONS {alat}
C 0.00 0.00 0.00
C 0.00 0.57735026918962576451 0.00
K_POINTS {automatic}
12 12 1 0 0 0

特に重要なパラメータは以下の通りです。

  • 単位胞の構造

     ibrav = 4,
     celldm(1) = 4.602,
     celldm(3) = 4,
    

    hexagonal 構造(ibrav=4)をもち、

    a1 = a(1,0,0), a2 = a(-1/2,sqrt(3)/2,0), a3 = a(0,0,c/a)

    としたときの a=4.602 (単位はボーア半径:aB, 1A~0.529177aB), c/a=4 としています。

  • 元素の種類と原子数

     nat = 2,
     ntyp = 1,
    

    炭素1種類(ntyp=1)で、単位胞内の原子数は2(nat=2)

  • 原子位置

    ATOMIC_POSITIONS {alat}
    C 0.00 0.00 0.00
    C 0.00 0.57735026918962576451 0.00
    

    xyz座標でa(=4.602)を単位として(alat)原子位置を指定。

    (0,0,0), (0, 1/sqrt(3), 0)

    の2点におくと、graphene ができあがる。

  • カットオフエネルギー

     ecutwfc = 30.0,
     ecutrho = 150.0,
    

    ecutwfc は波動関数を平面波で展開するときのカットオフエネルギー。単位は Ry。ecutrho は電荷密度を展開するときのカットオフ。 ノルム保存の擬ポテンシャルを使っていれば ecutrho は自動的に ecutrho = 4*ecutwfc となるため設定しなくてよいが、ウルトラソフトの場合は 4*ecutwfc より 大きめにとったほうがよい。

  • k点の取り方

    K_POINTS {automatic}
    12 12 1 0 0 0
    

    12x12x1 のk点をとることを表してます。0 0 0 はメッシュを全体的にシフトさせるかどうかで、今の場合Γ点を含むメッシュになっています。 12 12 1 1 1 0 なら、メッシュの中央にΓ点がきます。

  • pseudopotential

    C  12.0107  C.pz-van_ak.UPF
    

    炭素に対して pseudopotential のファイル C.pz-van_ak.UPF を使うことを意味しています。 これは、quantum ESPRESSO の website からダウンロードできるファイルの一つです。 pseudopotential の内容については、このファイルの先頭に簡単な記述がありますが、C.pz-van_ak.UPF は Perdew-Zunger のパラメータを用いたLDAを使った ultrasoft型の pseudopotential になっています。 12.0107 は炭素の質量ですが、この計算では利用しません。

    注釈

    ここでは、炭素のpseudopotentialとしてこのファイルを用いますが、炭素に対してこの pseudopotential の使用を推奨しているわけではないのでご注意下さい。 pseudopotential の選択(作成)は計算しようとしている物理量や計算量にも関係してくるため、その都度吟味する必要があります。

以上の内容は第一原理計算を行っている論文なら大抵記載されており、これでほぼ計算の内容は決まります。

各パラメータの意味の詳細は INPUT_PW.html を参照して下さい。

注釈

ここで設定したパラメータ以外にも非常に多くのパラメータが設定できます。 この計算では default の値で問題ないので省略していますが、本格的に計算する場合は INPUT_PW.html を一読しておくことをお勧めします。

計算を流すためにはさらに pseudopotentialのファイル C.pz-van_ak.upf を pseudo_dirで設定しているディレクトリに置く必要があります。 今の場合、pseudo_dir = ‘./’ ですので current directory に C.pz-van_ak.UPF を置いて下さい。 この状態で以下のコマンドを実行すれば scf計算が行われます。

% ~/espresso-5.3.0/bin/pw.x < graphene.scf.in > graphene.scf.out

結果は、graphene.scf.out に出力されます。このファイルをみて最後の方に

     the Fermi energy is    -0.4661 ev

!    total energy              =     -22.84944326 Ry
     Harris-Foulkes estimate   =     -22.84944326 Ry
     estimated scf accuracy    <          7.1E-10 Ry

     The total energy is the sum of the following terms:

     one-electron contribution =     -89.64903894 Ry
     hartree contribution      =      46.53466683 Ry
     xc contribution           =      -6.97848167 Ry
     ewald contribution        =      27.24348888 Ry
     smearing contrib. (-TS)   =      -0.00007836 Ry

     convergence has been achieved in   6 iterations

とあればうまく計算ができています。

なお、この計算によって得られるグラフェンのエネルギーは

% grep -a "^\!" graphene.scf.out

とすればコマンドラインからも確認できます。

3.2. band計算

次に、scf で計算した charge density を用いて任意のk点に対するKohn-Sham軌道のエネルギーを計算します。 scf計算に利用した graphene.scf.in から以下のハイライトした行だけ変更した graphene.nscf.in というファイルを作成してください。

&control
 calculation = 'bands'
 prefix='graphene',
 tstress = .true.
 tprnfor = .true.
 pseudo_dir = './',
 outdir='./work/'
 disk_io='low'
 wf_collect=.true.
/
&system
 ibrav = 4,
 celldm(1) = 4.602,
 celldm(3) = 4,
 nat = 2,
 ntyp = 1,
 ecutwfc = 30.0,
 ecutrho = 150.0,
 occupations = 'smearing'
 smearing = 'm-p'
 degauss = 0.01
 nbnd = 16
/
&electrons
 mixing_beta = 0.7
 conv_thr = 1.0d-8
/
ATOMIC_SPECIES
C  12.0107  C.pz-van_ak.upf
ATOMIC_POSITIONS {alat}
C 0.00 0.00 0.00
C 0.00 0.57735026918962576451 0.00
K_POINTS {tpiba_b}
4
0.00000000      0.00000000      0.00000000   20
0.66666667      0.00000000      0.00000000   20
0.50000000      0.28867500      0.00000000   20
0.00000000      0.00000000      0.00000000   20

nbnd の行は計算するバンド数の指定で、指定しなくても構いません。(デフォルトの値については INPUT_PW.html を参照してください。)

最後の5行で、Γ(0,0,0)-K(2/3,0,0)-M(1/2,1/2sqrt(3),0)-Γのライン上に各20点ずつ計算するよう指定しています。 tpiba_b はこの座標指定に 2pi/a を単位として用いることを表しています。 代わりにcrystal_b を指定すれば以下のように逆格子ベクトル単位で指定することもできます。

K_POINTS {crystal_b}
4
0.00000000      0.00000000      0.00000000   20
0.66666667     -0.33333333      0.00000000   20
0.50000000      0.00000000      0.00000000   20
0.00000000      0.00000000      0.00000000   20

注釈

espresso のバージョン5.1以降では、この最後の5行によるk点の指定をより単純に

4
gG 20
K  20
M  20
gG 20

と書くこともできます。 ibrav = 4 を指定しているので、 gG はΓ点、KはK点と解釈されます。 この表記方法の詳細はespressoを展開したディレクトリの中にある、Doc/brillouin_zones.pdfを参照してください。

このファイルを scf計算した directory と同じところにおいて

% ~/espresso-5.3.0/bin/pw.x < graphene.nscf.in > graphene.nscf.out

とやればバンド計算終了です。

あとは、計算したKohn-Shamエネルギーをファイルにまとめます。それ専用のプログラム bands.x があるので、以下のようなインプットファイル graphene.band.in を用意します。

&bands
   outdir = './work/',
   prefix='graphene',
   filband='graphene.band',
   lsym=.true.
/

注釈

espresso の ver.4 以前では、最初の行は &bands の代わりに &inputpp になります。

このファイルを用いて以下のコマンドを実行します。

% ~/espresso-5.3.0/bin/bands.x < graphene.band.in > graphene.band.out

これで、graphene.band および graphene.band.gnu というファイルができるはずです。

注釈

Mac でgfortranを使ってコンパイルした場合、エラーで落ちてしまうかもしれません。 そのときはこちらの インプット を試してみて下さい。

graphene.band というファイルには、以下のように各k点ごとにKohn-Shamエネルギーのリストが出力されています。

 &plot nbnd=  16, nks=    61 /
            0.000000  0.000000  0.000000
 -20.113  -8.374  -3.513  -3.513   2.628   4.431   5.452   8.033   8.033  10.527
  11.345  12.549  13.122  19.981  22.156  23.138
            0.033333  0.000000  0.000000
 -20.093  -8.349  -3.594  -3.558   2.656   4.458   5.479   8.047   8.122  10.553
  11.273  12.576  13.153  19.996  22.223  23.163
            0.066667  0.000000  0.000000
 -20.033  -8.276  -3.827  -3.692   2.738   4.540   5.560   8.088   8.375  10.629
  11.060  12.658  13.246  20.040  22.417  23.240
            0.100000  0.000000  0.000000
 -19.933  -8.155  -4.190  -3.911   2.874   4.675   5.695   8.159   8.760  10.717
  10.756  12.794  13.411  20.102  22.714  23.377
            0.133333  0.000000  0.000000
 -19.794  -7.986  -4.652  -4.206   3.066   4.864   5.884   8.266   9.229  10.258
  10.932  12.983  13.662  20.161  23.057  23.598
            0.166667  0.000000  0.000000
 -19.615  -7.769  -5.186  -4.568   3.312   5.107   6.128   8.410   9.701   9.732
  11.156  13.223  14.017  20.163  23.362  23.863
            0.200000  0.000000  0.000000
 -19.396  -7.505  -5.765  -4.988   3.613   5.403   6.426   8.597   9.066  10.220
...

また、graphene.band.gnu というファイルには gnuplot でプロットしやすいよう加工されたデータが出力されます。 エネルギーの単位は graphene.band では eV で graphene.band.gnu では Ry ~ 13.6 eV であることに注意してください。

これを gnuplot で band.plt というファイルを使ってプロットします。(この辺は人によってそれぞれ流儀があるところだと思います。あくまで参考程度に。)

% gnuplot band.plt

これで、次のようなファイルが生成されます。

_images/band.png

注釈

このプロットでは、E=0がFermi energyになるようにエネルギーをシフトしています。 (もともとの出力におけるエネルギー原点には意味がないので、普通、Fermi energy あるいは valence top などをエネルギー原点に なるようシフトしてプロットします。)

Fermi energy は graphene.scf.out に出力されていますが、あくまで離散的なk点を使って推測したものになります。 今の場合、厳密に K点のエネルギーを用いればいいと分かっているので、そちらを優先しています。

3.3. DOS計算

次に状態密度を求めてみましょう。

バンド計算したディレクトリとは別のディレクトリで、scf 計算を行って下さい。

% ~/espresso-5.3.0/bin/pw.x < graphene.scf.in > graphene.scf.out

その後、scf計算に利用した graphene.scf.in から以下のハイライトした行だけ変更した graphene.nscf.in というファイルを作成してください。

&control
 calculation = 'nscf'
 prefix='graphene',
 tstress = .true.
 tprnfor = .true.
 pseudo_dir = './',
 outdir='./work/'
 disk_io='low'
 wf_collect=.true.
/
&system
 ibrav = 4,
 celldm(1) = 4.602,
 celldm(3) = 4,
 nat = 2,
 ntyp = 1,
 ecutwfc = 30.0,
 ecutrho = 150.0,
 occupations = 'tetrahedra'
/
&electrons
 mixing_beta = 0.7
 conv_thr = 1.0d-8
/
ATOMIC_SPECIES
C  12.0107  C.pz-van_ak.UPF
ATOMIC_POSITIONS {alat}
C 0.00 0.00 0.00
C 0.00 0.57735026918962576451 0.00
K_POINTS {automatic}
36 36 1 0 0 0

最後の2行で、36x36x1 のメッシュ状のk点を計算するよう指定しています。 実際には対称性を用いて計算するk点の数を減らしています。 また、occupationsとしてtetrahedron法を用いるように指定しています。 tetrahedron法の詳細については、 tetrahedron法について をご覧下さい。 第一原理計算に限らず、限られたk点メッシュの情報からk空間での積分の近似計算を行う上で非常に有用な手法です。

このファイルを用いて nscf 計算を行います。

% ~/espresso-5.3.0/bin/pw.x < graphene.nscf.in > graphene.nscf.out

あとは、計算したエネルギーをまとめてDOSの計算を行います。 専用のプログラム dos.x のインプットとして、以下のようなファイル graphene.dos.in を用意します。

&dos
   outdir = './work/',
   prefix='graphene',
   fildos='graphene.dos',
/

注釈

espresso の ver.4 以前では、最初の行は &dos の代わりに &inputpp になります。

このファイルを用いて dos.x を実行すると、

% ~/espresso-5.3.0/bin/dos.x < graphene.dos.in > graphene.dos.out

graphene.dos というファイルができるはずです。 これはそのままプロットすれば以下のような図が得られます。

_images/dos.png