クロスアビリティ 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   = "relax"
    forc_conv_thr = 1.0e-03
    nstep         = 100
/

&SYSTEM
    ibrav       = 4
    a           = 5.71191
    c           = 24.6638
    ntyp        = 1
    nat         = 12

    ecutwfc     = 25
    ecutrho     = 100

    occupations = "smearing"
    smearing    = "gaussian"
    degauss     = 0.01
/

&ELECTRONS
    conv_thr = 1.0e-06
/

&IONS
    ion_dynamics = "bfgs"
/

K_POINTS {automatic}
 3  3  1  0 0 0

ATOMIC_SPECIES
Al     26.98154  Al.pbe-rrkj.UPF

ATOMIC_POSITIONS {angstrom}
Al      1.427977   0.824443  14.663754  1 1 1
Al      4.283932   0.824443  14.663754  1 1 1
Al      0.000000   3.297772  14.663754  1 1 1
Al      2.855955   3.297772  14.663754  1 1 1
Al      0.000000   1.648886  12.331877  1 1 1
Al      2.855955   1.648886  12.331877  1 1 1
Al     -1.427977   4.122215  12.331877  1 1 1
Al      1.427977   4.122215  12.331877  1 1 1
Al      0.000000   0.000000  10.000000  0 0 0
Al     -1.427977   2.473329  10.000000  0 0 0
Al      1.427977   2.473329  10.000000  0 0 0
Al      0.000000   4.946658  10.000000  0 0 0

上記は、Al[111]スラブの例です。必要最小限の内容を記したものです。入力ファイルの各項目について、以降の章で解説します。

#2. 基本設定
先ずは、calculationという変数に"relax"と設定します。これを設定すると、構造最適化計算が実施されます。必須の内容ですので、忘れずに設定して下さい。

&CONTROL
    calculation = "relax"
/

#3. 収束条件
構造最適化計算の収束判定には、原子に働く力を用います。すべての原子に働く力の絶対値が、forc_conv_thrに設定された値より小さくなった時に収束したと判定されます。forc_conv_thrの単位は、Ry/bohrです。通常の目的であれば、1.0e-3 Ry/bohr程度を設定していれば、十分な精度が得られます。また、構造最適化における繰り返し計算の最大回数は、nstepにて設定します。nstepには数十～数百程度を設定していれば、大抵の計算は収束に至ります。

&CONTROL
    forc_conv_thr = 1.0e-03
    nstep         = 100
/

#4. 最適化アルゴリズム
構造最適化計算のアルゴリズムには、BFGSとDamped-MDの２つが使えます。BFGSを使う場合には、以下のように設定します。

&IONS
    ion_dynamics = "bfgs"
/

構造最適化計算による原子座標の変化が小さい場合には、BFGSが適しています。例えば、データベースなどから結晶構造を取得して、これを初期構造とする場合には原子座標がほとんど動かないことが予想されます。

Damped-MDを使う場合には、以下のように設定します。

&IONS
    ion_dynamics = "damp"
/

原子座標の変化が大きい場合には、Damped-MDが有効です。データベースやツールなどを使わずに、手作業で原子座標を作成した場合には、構造最適化計算にて座標が大きく動くことが予測されます。そのような場合には、Damped-MDを使うことをお勧めします。

#5. 原子座標の設定
原子座標は、入力ファイル内のATOMIC_POSITIONSにて設定します。

ATOMIC_POSITIONS {angstrom}
Al      1.427977   0.824443  14.663754  1 1 1
Al      4.283932   0.824443  14.663754  1 1 1
Al      0.000000   3.297772  14.663754  1 1 1
Al      2.855955   3.297772  14.663754  1 1 1
Al      0.000000   1.648886  12.331877  1 1 1
Al      2.855955   1.648886  12.331877  1 1 1
Al     -1.427977   4.122215  12.331877  1 1 1
Al      1.427977   4.122215  12.331877  1 1 1
Al      0.000000   0.000000  10.000000  0 0 0
Al     -1.427977   2.473329  10.000000  0 0 0
Al      1.427977   2.473329  10.000000  0 0 0
Al      0.000000   4.946658  10.000000  0 0 0

ATOMIC_POSITIONS以下の各行が、原子座標を意味します。1列目が元素名、2~4列目が原子のXYZ座標です。5~7列目は、原子のXYZ座標を可動とするか否かを指示します。1の場合には該当の座標を可動として、0の場合には固定します。つまり、上記の入力ファイルの例では、1~8番目の原子が可動であり、構造最適化計算にてその座標が最適化されます。9~12番目の原子は固定であり、構造最適化計算の最中その座標は一定値となります。原子座標の特定の成分のみを可動とすることも可能です。例えば、

Al      1.427977   0.824443  14.663754  0 0 1

とすれば、XY座標が固定され、Z座標が最適化されます。また、XYZの3成分をすべて可動とする場合には、行末の「1 1 1」は省略可能です。例えば、

Al      1.427977   0.824443  14.663754

と書けば、当該原子は可動となります。

#6. 計算時間短縮のテックニック

上述の一連の設定を行えば、正しく構造最適化計算が実施されます。しかしながら、構造最適化計算は、その内部にてSCF計算を何度も行うため、比較的に計算時間が長くなります。そこで、計算コストが低減して、計算時間を短縮するためのテクニックを紹介します。基本的な考え方としては、構造最適化計算の序盤ではSCF計算の精度を低下させ、終盤ではSCF計算の精度を向上させるというものです。つまり、序盤においては、原子座標が十分に適切なものとは言えないため、当該の座標にて電子状態を厳密に扱うのはコストパフォーマンスが悪いというわけです。原子座標の精度が悪いときは、SCF計算の精度が悪くても問題ない、というアイデアです。具体的な計算手順としては、先ず、SCF計算の精度を落として構造最適化計算を実施します。計算完了後、出力ファイルに記載された最終構造を入力ファイルの初期構造に設定して、SCF計算の精度を向上させて再度構造最適化計算を実施します。このように構造最適化計算を2段階とすることで、序盤での計算時間を短縮することが可能です。原子座標の動きが大きい場合には、このテクニックはかなり効果的で、全計算時間を1/10近く短縮できる場合もあります。

SCF計算の精度を落とす場合、具体的には、電子状態の収束条件(conv_thr)、波動関数カットオフ(ecutwfc)、k点サンプリング(K_POINTS)を下表のように調整すると良いです。

設定項目	設定値
conv_thr	1.0e-4 ~ 1.0e-3
ecutwfc	16 Ry ~ 20 Ry
K_POINTS	Γ点のみ ~ 高精度計算の半値

また、序盤の低精度計算ではDamped-MDを使うのも良いです。低精度での構造最適化計算では、構造を完全に収束させる必要はないので、適当なところで計算を中断して後続の高精度計算に移っても問題ありません。ここでは2段階の計算について言及しましたが、初期構造が複雑な場合には 3段階や4段階にするとさらに効率が良くなるケースもあります。

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

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