グラフェンのバンド構造と状態密度 ================================ まず、簡単な計算例として、グラフェンのバンド構造と状態密度を求めてみましょう。 どちらも scf計算->non-scf計算->後処理という流れになります。 scf計算 ------- 何を計算するにしても、まず scf 計算により電荷密度(charge density)を求める必要があります。 scf 計算に必要なファイルはインプットファイル `graphene.scf.in `_ は以下の通りです。 とりあえずこのファイルをダウンロードして、適当な directory において下さい。 .. literalinclude:: files/graphene_band/graphene.scf.in 特に重要なパラメータは以下の通りです。 * 単位胞の構造 .. literalinclude:: files/graphene_band/graphene.scf.in :lines: 12-14 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 としています。 * 元素の種類と原子数 .. literalinclude:: files/graphene_band/graphene.scf.in :lines: 15-16 炭素１種類(ntyp=1)で、単位胞内の原子数は２(nat=2) * 原子位置 .. literalinclude:: files/graphene_band/graphene.scf.in :lines: 29-31 格子ベクトルa1, a2, a3を単位として(crystal)原子位置を指定。 0*a1 + 0*a2 + 0*a3, 1/3*a1 + 2/3*a2 + 0*a3 の２点におくと、graphene ができあがる。 * カットオフエネルギー .. literalinclude:: files/graphene_band/graphene.scf.in :lines: 17-18 ecutwfc は波動関数を平面波で展開するときのカットオフエネルギー。単位は Ry。ecutrho は電荷密度を展開するときのカットオフ。 ノルム保存の擬ポテンシャルを使っていれば ecutrho は自動的に ecutrho = 4*ecutwfc となるため設定しなくてよいが、ウルトラソフトの場合は 4*ecutwfc より 大きめにとったほうがよい。 * k点の取り方 .. literalinclude:: files/graphene_band/graphene.scf.in :lines: 32-33 12x12x1 のk点をとることを表してます。0 0 0 はメッシュを全体的にシフトさせるかどうかで、今の場合Γ点を含むメッシュになっています。 12 12 1 1 1 0 なら、メッシュの中央にΓ点がきます。 * pseudopotential .. literalinclude:: files/graphene_band/graphene.scf.in :lines: 28 炭素に対して pseudopotential のファイル C.pz-van_ak.UPF を使うことを意味しています。 これは、quantum ESPRESSO の `website `_ からダウンロードできるファイルの一つです。 pseudopotential の内容については、このファイルの先頭に簡単な記述がありますが、C.pz-van_ak.UPF は Perdew-Zunger のパラメータを用いたLDAを使った ultrasoft型の pseudopotential になっています。 12.0107 は炭素の質量ですが、この計算では利用しません。 .. note:: ここでは、炭素のpseudopotentialとしてこのファイルを用いますが、炭素に対してこの pseudopotential の使用を推奨しているわけではないのでご注意下さい。 pseudopotential の選択（作成）は計算しようとしている物理量や計算量にも関係してくるため、その都度吟味する必要があります。 以上の内容は第一原理計算を行っている論文なら大抵記載されており、これでほぼ計算の内容は決まります。 各パラメータの意味の詳細は `INPUT_PW.html `_ を参照して下さい。 .. note:: ここで設定したパラメータ以外にも非常に多くのパラメータが設定できます。 この計算では default の値で問題ないので省略していますが、本格的に計算する場合は INPUT_PW.html を一読しておくことをお勧めします。 計算を流すためにはさらに pseudopotentialのファイル C.pz-van_ak.upf を pseudo_dirで設定しているディレクトリに置く必要があります。 今の場合、pseudo_dir = './' ですので current directory に `C.pz-van_ak.UPF `_ を置いて下さい。 この状態で以下のコマンドを実行すれば scf計算が行われます。 :: % pw.x < graphene.scf.in > graphene.scf.out 結果は、graphene.scf.out に出力されます。このファイルをみて最後の方に .. literalinclude:: files/graphene_band/graphene.scf.out :language: none :lines: 277-291 とあればうまく計算ができています。 なお、この計算によって得られるグラフェンのエネルギーは :: % grep -a "^\!" graphene.scf.out とすればコマンドラインからも確認できます。 band計算 -------- 次に、scf で計算した charge density を用いて任意のk点に対するKohn-Sham軌道のエネルギーを計算します。 scf計算に利用した graphene.scf.in から以下のハイライトした行だけ変更した `graphene.nscf.in `_ というファイルを作成してください。 .. literalinclude:: files/graphene_band/graphene.nscf.in :language: none :emphasize-lines: 2,22,33-38 nbnd の行は計算するバンド数の指定で、指定しなくても構いません。(デフォルトの値については INPUT_PW.html を参照してください。) 最後の５行で、Γ(0,0,0)-K(2/3,0,0)-M(1/2,1/2sqrt(3),0)-Γのライン上に各30点ずつ計算するよう指定しています。 tpiba_b はこの座標指定に 2pi/a を単位として用いることを表しています。 代わりにcrystal_b を指定すれば以下のように逆格子ベクトル単位で指定することもできます。 :: K_POINTS {crystal_b} 4 0.00000000 0.00000000 0.00000000 30 0.66666667 -0.33333333 0.00000000 30 0.50000000 0.00000000 0.00000000 30 0.00000000 0.00000000 0.00000000 30 また、ibrav を0以外にしている場合は :: K_POINTS {crystal_b} 4 gG 30 K 30 M 30 gG 30 のようにk点の位置をシンボルで表すこともできます。(今の場合、ibrav = 4 を指定しているので、 gG はΓ点、KはK点と解釈されます。) この表記方法の詳細はespressoを展開したディレクトリの中にある、Doc/brillouin_zones.pdfを参照してください。 このファイルを scf計算した directory と同じところにおいて :: % pw.x < graphene.nscf.in > graphene.nscf.out とやればバンド計算終了です。 あとは、計算したKohn-Shamエネルギーをファイルにまとめます。それ専用のプログラム bands.x があるので、以下のようなインプットファイル `graphene.band.in `_ を用意します。 .. literalinclude:: files/graphene_band/graphene.band.in このファイルを用いて以下のコマンドを実行します。 :: % bands.x < graphene.band.in > graphene.band.out これで、graphene.band および graphene.band.gnu というファイルができるはずです。 .. note:: Mac でgfortranを使ってコンパイルした場合、エラーで落ちてしまうかもしれません。 そのときはこちらの `インプット `_ を試してみて下さい。 graphene.band というファイルには、以下のように各k点ごとにKohn-Shamエネルギーのリストが出力されています。 .. literalinclude:: files/graphene_band/graphene.band :language: none :lines: 1-21 :append: ... また、graphene.band.gnu というファイルには gnuplot でプロットしやすいよう加工されたデータが出力されます。 これを gnuplot で `band.plt `_ というファイルを使ってプロットします。（この辺は人によってそれぞれ流儀があるところだと思います。あくまで参考程度に。） :: % gnuplot band.plt これで、次のようなファイルが生成されます。 .. image:: files/graphene_band/band.* .. note:: このプロットでは、E=0がFermi energyになるようにエネルギーをシフトしています。 (もともとの出力におけるエネルギー原点には意味がないので、普通、Fermi energy あるいは valence top などをエネルギー原点に なるようシフトしてプロットします。) Fermi energy は graphene.scf.out に出力されていますが、あくまで離散的なk点を使って推測したものになります。 今の場合、厳密に K点のエネルギーを用いればいいと分かっているので、そちらを優先しています。 DOS計算 ------- 次に状態密度を求めてみましょう。 バンド計算したディレクトリとは別のディレクトリで、scf 計算を行って下さい。 :: % pw.x < graphene.scf.in > graphene.scf.out その後、scf計算に利用した graphene.scf.in から以下のハイライトした行だけ変更した `graphene.nscf.in `_ というファイルを作成してください。 .. literalinclude:: files/graphene_dos/graphene.nscf.in :language: none :emphasize-lines: 2,19,30-31 最後の２行で、36x36x1 のメッシュ状のｋ点を計算するよう指定しています。 実際には対称性を用いて計算するｋ点の数を減らしています。 また、occupationsとしてtetrahedron法を用いるように指定しています。 tetrahedron法の詳細については、 `tetrahedron法について `_ をご覧下さい。 第一原理計算に限らず、限られたk点メッシュの情報からk空間での積分の近似計算を行う上で非常に有用な手法です。 このファイルを用いて nscf 計算を行います。 :: % pw.x < graphene.nscf.in > graphene.nscf.out あとは、計算したエネルギーをまとめてDOSの計算を行います。 専用のプログラム dos.x のインプットとして、以下のようなファイル `graphene.dos.in `_ を用意します。 .. literalinclude:: files/graphene_dos/graphene.dos.in このファイルを用いて dos.x を実行すると、 :: % dos.x < graphene.dos.in > graphene.dos.out graphene.dos というファイルができるはずです。 これはそのままプロットすれば以下のような図が得られます。 .. image:: files/graphene_dos/dos.* E=0はFermi energyで、Dirac点に対応しており、DOSが0になります(紫線)。 一方、DOSの積分(青線)は電子数に対応しており、E=0の値が8になっているのはEF以下に8個電子があることを意味しています。 これは、単位胞に2つ炭素原子があり、それぞれ4つずつ電子を考えていることに対応しています。(1sの電子は擬ポテンシャルの中に入っているため、炭素あたり2s,2pに4つ電子がいると考えている)