収束のチェック ============== ここでは、計算がちゃんと収束しているかどうか確認する方法について紹介します。 scfループの収束条件 ------------------- まず、scfのループに対する収束条件を見てみます。 graphene.scf.out を見てみて下さい。 :: iteration # 1 ecut= 30.00 Ry beta= 0.70 Davidson diagonalization with overlap ethr = 1.00E-02, avg # of iterations = 5.7 negative rho (up, down): 1.315E-04 0.000E+00 total cpu time spent up to now is 0.4 secs total energy = -22.81159749 Ry estimated scf accuracy < 0.23064343 Ry のような感じでループ一回ごとに total energy などが計算されていることが分かります。 この total energy が conv_thr よりも高い精度で得られた時点でscf ループが終わります。 一度、conv_thr をいろいろ変えてみていつループが終わるか試してみるとよいでしょう。 conv_thr は求めたいものの精度や計算時間との兼ね合いになりますが、default の conv_thr = 1.d-6 は収束条件として若干あまいかもしれません。 ただし、あまり収束条件をきつくすると収束しなくなるので、注意しましょう。 (普通に使われる倍精度の実数の精度がそもそも16桁程度しかありません。) cutoffエネルギーに対する収束 ---------------------------- 次にカットオフ(ecutwfc)を変えて全エネルギーの変化を見てみます。 カットオフは平面波の基底の数を制限するために用いているので、大きければ大きいほど精度があがることになります。 (ただし、その分計算が重くなっていきます。) どの程度のカットオフを用いればよいかは使うpseudopotenitalに大きく依存します。 大きなカットオフが必要なpseudopotentialのことを堅い(hard)pseudopotentialといい、小さいカットオフで済むpseudopotentialのことを やわらかい(soft)pseudopotentialといいます。 scf計算に必要なファイル(graphene.scf.in)を用意した上で以下のようなスクリプト(`conv_cutoff.sh `_)を実行してみましょう。 .. literalinclude:: files/graphene_conv/conv_cutoff.sh これで、graphene_10.scf.out のような計算結果のファイルができるはずです。 この10の部分が波動関数の対するカットオフで上記のスクリプトでは、10から50まで5通り計算しています。 この結果は以下のコマンドで確認できます。 :: % grep "^\!" *.out graphene_10.scf.out:! total energy = -21.92216798 Ry graphene_20.scf.out:! total energy = -22.79213625 Ry graphene_30.scf.out:! total energy = -22.84944327 Ry graphene_40.scf.out:! total energy = -22.85150423 Ry graphene_50.scf.out:! total energy = -22.85536979 Ry 見てわかる通り、カットオフが 30 Ry ぐらいになると全エネルギーが 0.01 Ry ぐらいの精度がでているようにみえます。 ただし、50Ryでも完全に収束しているわけではないですし、これだけでは実際にどの程度のカットオフを使えばよいのか判断が難しいと思います。 カットオフに関しては、実際に必要な物理量（バンドギャップや有効質量、状態密度など）をカットオフを変えて計算してみるのがよいでしょう。 例えば、異なる相(ダイヤモンドとグラフェンなど)の間のエネルギー差は、同じカットオフを使えば同程度の誤差が入る(同程度エネルギーが高めにでる)ため、その精度は0.01 Ry よりも高くなります。 k点メッシュに対する収束 ----------------------- k点メッシュとは、計算上現れる波数積分を離散的なメッシュ状の点の和で置き換えたときのメッシュのことを表しています。 従って、k点メッシュは細かければ細かい方がよいことになります。 エネルギーのk点メッシュ依存性をみるために、以下のようなスクリプト(`conv_kpoint.sh `_)を実行してみましょう。 .. literalinclude:: files/graphene_conv/conv_kpoint.sh この結果を以下のようにみてみると、 :: % grep "^\!" *.out graphene_k12.scf.out:! total energy = -22.84944327 Ry graphene_k16.scf.out:! total energy = -22.84999504 Ry graphene_k4.scf.out:! total energy = -22.85248079 Ry graphene_k8.scf.out:! total energy = -22.84985642 Ry となります。 上からk点メッシュが 12x12x1, 16x16x1, 4x4x1, 8x8x1 でのエネルギーになります。 8x8x1 ぐらいで0.001 Ry 程度の精度がでていることが分かると思います。 もし、この物質が絶縁体であれば、このようにしてk点メッシュの収束を確認すれば終わりということになりますが、金属では以下のsmearingの幅についても考慮する必要があります。 Smearing について ----------------- Smearing は金属の計算において、Fermi面近傍の状態に対して、どの状態にどれだけ電子を詰めるかを決めるために用います。 絶対零度においては、EF以下の軌道は完全につまっていてEF以上の軌道は完全に空になるわけですが、これを離散化したk点を使ってそのまま実装すると、複数のバンドがEFにかかっている場合などでは EF のちょっとの変化に対して全エネルギーが大きく変わってしまう、ということが起こりえます。 そこで、本来E=EFで不連続になる占有数を鈍らせることにより計算を収束しやすくすることができます。 (Fermi分布関数をT=0のものを使う代わりにT有限のものを使う、ということをイメージしてもらえればよいと思います。) ここで、どのように鈍らせるかを決めるオプションが「smearing」であり、その幅を決めるのが「degauss」になります。 従って、degauss はk点メッシュと密接に関係しています。 一般にsmearing の幅が大きければ、荒いk点メッシュでも収束しますが、求めるべき積分値とは少しずれることが予想されます。 理想的にはk点メッシュを細かく、smearing の幅も小さくすると求めるべき積分値に近づいていきます。 個人的には、smearing='m-p', degauss=0.01-0.02 (Ry) ぐらいを用いることが多いですが、これも必要な物理量とその精度に依存することですので、適宜確認してみて下さい。