クロスアビリティ Winmostarサポートチームです。
WinmostarはQuautum ESPRESSOを簡単かつ高度に利用するための統合GUI環境となっています。学生は無料で利用でき、学生以外も無料トライアルを入手することができますので、是非 WinmostarのWebサイトをご覧ください。

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

この記事では、初歩的な内容しか紹介しません。
入力ファイルの詳細について知りたい方は、マニュアル 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. ストレステンソル
原子座標のみを最適化する場合、全ての原子に働く力が閾値以下になれば、エネルギーは極小値になり、収束した結果をそのまま使うことができます。原子座標のみの最適化では格子ベクトルは固定で、格子ベクトルに依存する平面波基底の数は一定であるためです。しかしながら、格子ベクトルの最適化も含む計算では状況が異なります。Quantum ESPRESSOの格子ベクトルの最適化では、格子ベクトルが変化しても計算開始時に決めた平面波基底の数のまま計算を進め、ストレステンソルが閾値以下になるまで行います。そのため、構造が収束した時のエネルギーは、初期格子ベクトルから決められた平面波基底の数で算出されます。

そこで、最適化された格子ベクトルを基に平面波基底の数を決め直して、あらためてエネルギーを1度だけ計算します。そうすることで、入力ファイルで指定するカットオフエネルギーと格子ベクトルから決める平面波基底の数の対応付けを正確に行うことができます。平面波基底の数が変わるため、最後のエネルギーは1つ前の構造収束時のエネルギーに比べて上がる場合もあれば、下がる場合もあります。初めから十分な平面波基底の数で計算する、つまりカットオフエネルギーを大きくして計算すると、このエネルギーの違いを小さくすることができます。

#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

１行目の alat= 10.33440199 というのが、格子定数(bohr単位)です。２～４行目が格子ベクトルです。この格子ベクトルは、格子定数(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となる格子ベクトルにてエネルギーをほぼ最小値とすることが出来ます。しかしながら、真のエネルギー = 偽のエネルギーが完全に成り立つわけではないので、念の為に、最適化された格子ベクトルの妥当性をチェックすることが望ましいです。具体的には、最適化された格子ベクトルを用いて入力ファイルを作成し、再度、格子ベクトル最適化計算を実施します。再計算にて格子ベクトルが変化しなければ、妥当な結果と言えます。格子ベクトルが変化した場合には、そのまま計算を続行して最適化された格子ベクトルを再度チェックします。

次回は、バンド構造計算か状態密度計算について解説したいと思います。
なお、Winmostarを利用するとQuautum ESPRESSOをGUI上から簡単かつ高度に活用することができます。学生は無料で利用でき、学生以外も無料トライアルを入手することができますので、是非 WinmostarのWebサイトをご覧ください。

Quantum ESPRESSO入力ファイル作成手順シリーズ
１．結晶構造の作成
 ２．SCF計算の設定
 ３．擬ポテンシャルファイルの選択方法
 ４．原子座標のみの最適化
 ５．格子ベクトル及び格子内部の原子座標の最適化
 ６．状態密度の計算
 ７．局所状態密度の計算
 ８．バンド構造の計算
 ９．バンド数の設定
 １０．van der Waals相互作用
 １１．LDA+U法
 １２．NEB法
 １３．Phonon計算（特定q点）
１４．Phonon計算（バンド構造）
１５．擬ポテンシャルの作成
 １６．擬ポテンシャルのテスト
 １７．SCF計算の入力ファイルの実例