概要
第一原理計算では、できるだけ外部で調節するパラメータを減らして電子状態を計算する。 しかし、どうしても自分の手で設定しないといけないものがある。 それは
- k点数
- 波動関数についての運動エネルギーのカットオフecutwfc
- 電荷密度についての運動エネルギーのカットオフecutrho
- degauss(金属のみ)
である。
以下ではecutrhoを除くこれらの値の決定方法について簡単に解説する。
理想としては、できるだけ大きいk点数、できるだけ大きいカットオフで計算を行うことである。 ところが計算資源や時間というものは有限である上に、ある程度以上の精度をあげた計算はそもそも意味をなさないということもあるので、どこかで打ち切る基準が必要である。
まず最初に、出発点となる値を決めなければならない。 経験的には、次の値から出発するとよい。
- k点数・・・単位胞内の原子数が1から3なら6*6*6、4から10なら4*4*4、10から20なら2*2*2、それ以上は1*1*1。物質の次元性に合わせて各方向にとるk点数を変えてもよい。
- ecutwfc・・・まずは20, 30, 40, 60 Ryで計算して様子を見る。ウルトラソフト型の擬ポテンシャルを使う場合はその半分。
- ecutrho・・・ecutwfcの4倍。ただし擬ポテンシャルがウルトラソフト型の場合はそれ以上(およそ8~12倍が必要と言われているが、とりあえず最初は12倍しておく)。
- degauss・・・0.02程度。occupations='smearing'とsmearing='mp'を同時に指定。
なお、occupations='tetrahedra'を使うと金属でもdegaussを指定する必要がなくなるが、smearingの場合よりも多くのk点数を要求されること、低次元物質ではうまくいかないことから、SCF計算では推奨されない。 degaussというフェルミ準位の値に影響する任意の値が入ってしまわない利点もあるので、状態密度計算などSCF計算後にk点を多く取った上で利用すると便利である。
QE6.1からテトラへドロン法を改良した occupations='tetrahedra_opt' が追加実装されており、これを使ってSCF計算しても問題ないと見る向きもある。 ただしsmearingを使う場合に比べて多くのk点数が必要であることに変わりない。
最近はecutwfcやecutrhoの大きさの目安が擬ポテンシャルファイルに書いてあるので、それを使ってk点数とdegaussの収束のみを調べるのが便利である。
各パラメータの収束判定方法
ecutwfc
ecutwfcを変えながらSCF計算を行い、total energyの収束を見る(金属の場合はFermi energyも)ことでecutwfcの値を決定する。 ecutwfcに対するtotal energyのグラフを作成してみるとわかりやすい。 ほとんどの場合、ecutwfcが大きくなるにしたがってtotal energyはある値に向かって収束していくはずである。 total energyに関して欲しい有効数値の桁数が得られる最小のecutwfcを探す。 ただし、ecutwfcを増やせば増やすほど精度が増すかというとそうでもなく、配布されている擬ポテンシャルの有効桁数から数値計算の精度の限界は1mRy程度である。
ちなみに、ecutwfcが小さいと必要な状態が足りずに収束にかかるiterationの回数が大きくなることを利用して、ecutwfcを決定する方法もある。
ecutrho
ecutrhoについては、擬ポテンシャルの種類に対して
- ノルム保存型 ・・・ ecutwfcの4倍
- PAW型 ・・・ ecutwfcの4倍以上
- ウルトラソフト型 ・・・ ecutwfcの8から12倍程度
が推奨されている。 ノルム保存型の場合は ecutrho = 4 * ecutwfc が正確に成り立つ。 ウルトラソフト型に関してはecutrhoの収束も見ないといけない。
入力ファイルでecutrhoを指定しなければ自動的にecutwfcの4倍になる。
k点数
ecutwfcが決まったら、k点数を変えて計算する。 最初が2*2*2であれば、次は3*3*3、4*4*4と増やしていけばよい。 物質の1次元性や2次元性が強いのであれば、そちらの方向を優先して増やしていく。 total energyとFermi energy(金属の場合)の変化を見て、収束に必要なk点数を探す。 なおFermi energyの誤差を10meV程度まで追求するのはあまり意味がない。
degauss
金属の場合はdegaussの設定が必要であるが、これが大変むずかしい。 k点は離散的にしかとれないので、scf計算に使った少しのk点から準位をブロードニングによってひろげてフェルミ面がどこにあるか判断する。 基本的にk点数が多いほどフェルミ準位を正確に見積もることができるが、当然、計算コストは増大する。 そのため、k点数を減らしてその分をdgaussを大きくすることで補うことになる。
degaussに対して0.001~0.1の範囲で全エネルギーとフェルミ準位をプロットしてみて、もっともらしいと思われるdegaussを探すのが、地道だが確実な方法である。 smearingが'gaussian'や'mv'の場合は結果がdegaussに強く依存するが、mpは比較的依存性が弱い。
参考:https://www.researchgate.net/post/Which_degauss_value_should_I_use_for_QUANTUM_ESPRESSO_input_file
smearingは通常はmpで良いが、負の電荷密度を許した計算になっているので、計算でエラーが出ることがある。 そのような場合はmvに変更するとよい。 ただしdegauss依存性がmpに比べて大きいので注意する。
時間はかかるものの、mpとmvで計算したフェルミ準位をdegaussに関してプロットして、それらが交差する点をもってdegaussを決めるという方法もある。 また、k点数を上で収束させた値よりも多めにとってみて、occupations='tetrahedra'(あるいは'tetrahedra_opt')の結果と比較してみるのも良い。
まとめ
以上の流れをまとめると次のようになる。
- k点数、ecutwfc、degauss(金属のみ)について出発点となる適当な値を決める。ウルトラソフト型擬ポテンシャルの場合はecutrhoをecutwfcの12倍にする。
- ecutwfcを変えてtotal energyとFermi energy(金属のみ)を収束させる。この段階での収束は粗めでよい。
- k点数を変えて収束させる
- 再びecutwfcを変えて収束させる
- (ウルトラソフト型擬ポテンシャルの場合のみ)ecutrhoを収束の影響がない程度まで減らす。
- (金属のみ)degaussを変化させてFermi energyを決定する
以上により、外部から与えるパラメータの値を決定できる。
注意
全エネルギーは比較的収束性が良いので、欲しい物理量によってはより多くのk点やカットオフエネルギーが必要になる場合がある。 必ずカットオフを変えて計算してみて最終的に欲しい値がちゃんと収束しているか確認すること。