Help us understand the problem. What is going on with this article?

Quantum ESPRESSO入力ファイル作成手順5.格子ベクトルの最適化

クロスアビリティ Winmostarサポートチームです。
WinmostarはQuautum ESPRESSOを簡単に利用するための統合GUI環境となっています。無料トライアルを配布していますので、是非 WinmostarのWebサイト をご覧ください。

0. はじめに

第一原理計算ソフトQuantunm ESPRESSO https://www.quantum-espresso.org の入力ファイルの書き方について、複数回にわたって解説していきます。今回は第5回目です。格子ベクトル最適化について解説したいと思います。

この記事では、初歩的な内容しか紹介しません。
入力ファイルの詳細について知りたい方は、マニュアル http://www.quantum-espresso.org/Doc/INPUT_PW.html をご参照下さい。

1. 入力ファイル

格子ベクトル最適化計算を実施する場合、入力ファイルをこのように書きます。

&CONTROL
    calculation   = "vc-relax"
    forc_conv_thr = 1.0e-03
    nstep         = 100
/

&SYSTEM
    ibrav       = 2
    a           = 5.46873
    ntyp        = 1
    nat         = 2

    ecutwfc     = 40
    ecutrho     = 160
    occupations = "fixed"
/

&ELECTRONS
    conv_thr = 1.0e-06
/

&IONS
    ion_dynamics = "bfgs"
/

&CELL
    press_conv_thr = 0.5
    cell_dynamics  = "bfgs"
/

K_POINTS {automatic}
 4  4  4  0 0 0

ATOMIC_SPECIES
Si     28.08550  Si.pbe-rrkj.UPF

ATOMIC_POSITIONS {alat}
Si      0.000000   0.000000   0.000000
Si      0.250000   0.250000   0.250000

上記は、Si結晶の例です。必要最小限の内容を記したものです。入力ファイルの各項目について、以降の章で解説します。

2. 基本設定

先ずは、calculationという変数に"vc-relax"と設定します。これを設定すると、格子ベクトル最適化計算が実施されます。必須の内容ですので、忘れずに設定して下さい。

&CONTROL
    calculation = "vc-relax"
/

3. 収束条件

格子ベクトル最適化計算では、ストレステンソルに従って格子ベクトルを少しずつ変形させます。そして、ストレステンソルが十分小さくなった時に、計算が収束したと判定します。収束判定の閾値は、press_conv_thrにて設定します。単位は、kbarです。デフォルト値は0.5kbarで、通常の目的であればデフォルト値のままで問題ありません。

&CELL
    press_conv_thr = 0.5
/

また、最適化における繰り返し計算の最大回数は、nstepにて設定します。nstepには数十~数百程度を設定していれば、大抵の計算は収束に至ります。

&CONTROL
    nstep = 100
/

格子ベクトルの変形に伴い、原子座標も変化します。このため、格子ベクトル最適化計算の際には、原子座標も同時に最適化される必要があります。原子に働く力の閾値は、forc_conv_thrにて設定します(単位は、Ry/bohr)。

&CONTROL
    forc_conv_thr = 1.0e-03
/

4. 最適化アルゴリズム

最適化のアルゴリズムには、いくつかのものが用意されています。しかしながら、格子ベクトル最適化計算においては、多くの場合に BFGSのみで事足ります。以下のように設定します。

&CELL
    cell_dynamics = "bfgs"
/

格子ベクトルをBFGSで最適化する場合には、原子座標の最適化アルゴリズムにもBFGSを適用する必要があるので留意してください。以下のように設定します。

&IONS
    ion_dynamics = "bfgs"
/

5. ストレステンソル

原子座標を最適化する場合、原子に働く力が0となる構造にてエネルギーが最小となります。これは、変分原理が成り立つためです。すなわち、波動関数を展開する平面波基底と原子座標の直積にて表現された変分空間にて、エネルギーが最小となる点を探索しているわけです。しかしながら、格子ベクトル最適化計算では状況が異なります。ストレステンソルが0となる格子ベクトルにおいても、エネルギーは最小とはなりません。格子ベクトルが変形するとフーリエ変換のメッシュが変形するため、平面波基底の数が変化します。この平面波基底の増減が、ストレステンソルとエネルギーの矛盾をもたらします。つまり、平面波基底の成す変分空間は格子ベクトルの値に依存するため、格子ベクトルの変分空間との直積空間では変分原理が成り立たないわけです。

しかしながら、変分原理が成り立たないと最適化のアルゴリズムが適応できません。そこで、格子ベクトル最適化の最中では暫定的に変分原理が成り立つように、平面波基底の数を固定します。そうすると、ストレステンソルは0へと収束し、”エネルギー”は最小となります。ただし、この”エネルギー”は固定された平面波基底にて計算された偽のエネルギーです。最適化後の格子ベクトルと入力ファイルにて指定されたカットオフエネルギーにて再定義された平面波基底にて、再度SCF計算を実施すると真のエネルギーが得られます。実際に、Quantum ESPRESSOにてvc-relaxを実行すると、格子ベクトルを最適化した後に、最適化構造にて真のエネルギーが計算されます。多くの場合、真のエネルギーは偽のエネルギーよりも大きくなります。偽のエネルギーが変分的に最小値であるのに対して、真のエネルギーは最小値ではないためです。

真のエネルギー ≠ 偽のエネルギー という問題は、数値計算上の都合によります。つまり、波動関数を平面波で展開しているということが原因なわけです。仮に、Kohn-Sham方程式の厳密解が常に得られるのであれば、このような問題は存在しません。平面波は完全基底であり、カットオフエネルギーを大きくすることでその精度を任意に高めることが出来ます。したがって、波動関数のカットオフエネルギー(ecutwfc)を無限に大きくすれば、真のエネルギー = 偽のエネルギー となります。実際には、ecutwfcを通常の1.5~2.0倍程度とすることで、概ね 真のエネルギー ≒ 偽のエネルギー が達成できます。

6. 計算結果の見方

vc-relaxを実行すると、最適化された格子ベクトルなどの計算結果は標準出力に表示されます。残念ながら出力内容が少々読みづらいので、ここで解説します。計算を実施すると、最適化計算の各ステップにて更新された格子ベクトルの値が下記のように出力されます。

CELL_PARAMETERS (alat= 10.33440199)
  -0.501531417  -0.000000000   0.501531417
   0.000000000   0.501531417   0.501531417
  -0.501531417   0.501531417  -0.000000000

1行目の alat= 10.33440199 というのが、格子定数(bohr単位)です。2~4行目が格子ベクトルです。この格子ベクトルは、格子定数(alat)を単位とした無次元量です。bohr単位の格子ベクトルは、無次元量の格子ベクトルに格子定数(alat)を乗じることで得られます。この例では、bohr単位の格子ベクトルは、

 -5.1830273  0.0000000  5.1830273
  0.0000000  5.1830273  5.1830273
 -5.1830273  5.1830273  0.0000000

となります。また、Quantum ESPRESSOのvc-relaxでは、格子定数を固定して、無次元量の格子ベクトルを最適化します。真のエネルギーを計算するSCFの直前に出力されたCELL_PARAMETERSが、格子ベクトルの最適化値です。

7. 最適化された格子ベクトルのチェック

波動関数のカットオフエネルギーを十分に大きくすることで、ストレステンソルが0となる格子ベクトルにてエネルギーをほぼ最小値とすることが出来ます。しかしながら、真のエネルギー = 偽のエネルギー が完全に成り立つわけではないので、念の為に、最適化された格子ベクトルの妥当性をチェックすることが望ましいです。具体的には、最適化された格子ベクトルを用いて入力ファイルを作成し、再度、格子ベクトル最適化計算を実施します。再計算にて格子ベクトルが変化しなければ、妥当な結果と言えます。格子ベクトルが変化した場合には、そのまま計算を続行して最適化された格子ベクトルを再度チェックします。


次回は、バンド構造計算か状態密度計算について解説したいと思います。
ご意見、ご要望ありましたら、コメントください。

Quantum ESPRESSO入力ファイル作成手順シリーズ
1.結晶構造の作成
2.SCF計算の設定
3.擬ポテンシャルファイルの選択方法
4.原子座標のみの最適化
5.格子ベクトルの最適化
6.状態密度の計算
7.局所状態密度の計算
8.バンド構造の計算
9.バンド数の設定
10.van der Waals相互作用
11.LDA+U法
12.NEB法
13.Phonon計算(特定q点)
14.Phonon計算(バンド構造)
15.擬ポテンシャルの作成
16.擬ポテンシャルのテスト
17.SCF計算の入力ファイルの実例

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away