Differences

This shows you the differences between two versions of the page.

--- quantumespresso:構造緩和 [2018/02/25 22:25]
koudai [原子位置の構造緩和]
+++ quantumespresso:構造緩和 [2024/03/30 23:58]
koudai
@@ Line 2: / Line 2: @@
 構造緩和計算により、結晶の安定な構造を求めます。
 構造緩和には次の2通りの方法があります。
-  - 格子定数はそのままで、単位胞内の原子の位置についての構造緩和
+  - 格子定数はそのままで、単位胞内の原子の内部座標の最適化
-  - 単位胞内の原子位置はそのままで、格子定数についての構造緩和
+  - 内部座標とともに格子定数も変化させる構造緩和
 通常は2つを繰り返すことで、安定な結晶構造を探します。
 以下ではケイ素を例に構造緩和計算を行います。
+===== 注意 =====
+  * 現在（2024年3月）のところ、noncolin=.true.のときは原子に加わる力の計算ができません
 ====== 原子位置の構造緩和 ======
@@ Line 16: / Line 21: @@
 格子定数はそのままに、単位胞内の原子の位置を緩和させます。 Siの単位胞内における原子の位置は、格子定数を単位として安定な構造で(0,0,0)と(0.25,0.25,0.25)にありますが、原子を少しずらした点において構造緩和をさせてみます。
-<file - si.relax.in>
+<file - Si.relax.in>
 &control
    calculation='relax'
-   prefix='si'
+   prefix='Si'
-   pseudo_dir='./'
+   pseudo_dir='./pseudo/'
-   outdir = './'
+   outdir = './tmp/'
    etot_conv_thr = 1.d-5
    forc_conv_thr = 1.d-4
@@ Line 31: / Line 36: @@
 /
 &electrons
+   conv_thr = 1.0d-8
 /
 &ions
@@ Line 44: / Line 50: @@
 構造緩和計算では、calculation='relax'を指定します。 関連するオプションの意味は次の通りです。
-etot_conv_thr	構造緩和の各ステップで、エネルギーの変化がこの値よりも小さくなったときに収束したと判断する。デフォルトは1.d-4だが、経験的にはもう少し厳しいほうが良い。
-forc_conv_thr	構造緩和の各ステップで、原子に加わる力の変化がこの値よりも小さくなったときに収束したと判断する。デフォルトは1.d-3だが、経験的にはもう少し厳しいほうが良い。
+^変数^初期値^説明^
+|etot_conv_thr|1.d-4|構造緩和の各ステップで、エネルギーの変化がこの値よりも小さくなったときに収束したと判断する|
+|forc_conv_thr|1.d-3|構造緩和の各ステップで、原子に加わる力の変化がこの値よりも小さくなったときに収束したと判断する|
+DFTで求められる程度の原子位置の精度なので、多くの場合は収束の厳しさはここで指定した値で十分です。
 また、relax計算の際は&ionsという項目が必要です。 ここに必要に応じて構造緩和計算のオプションを書き加えますが、今はデフォルトのままとします。
@@ Line 56: / Line 65: @@
 実行は次のようにします。
-  $ pw.x < si.relax.in > si.relax.out
+  $ pw.x < Si.relax.in > Si.relax.out
 結果は次のようになります。
-<file - si.relax.out>
+<file - Si.relax.out>
 （略）
 Begin final coordinates
@@ Line 76: / Line 85: @@
 ====== 格子定数の構造緩和 ======
-今度は格子定数の構造緩和を行います。 Siの格子定数はボーア半径単位で10.2程度ですが、ここから少しずらして構造緩和をさせてみます。
+今度は格子定数の構造緩和を行います。
+先に原子位置の緩和を行ってください。
-<file - si.vc-relax.in>
+Siの格子定数はボーア半径単位で10.2程度ですが、ここから少しずらして構造緩和をさせてみます。
+Si.relax.inをコピーしてSi.vc-relax.inを作成します。
+以下の変更を行います
+  * calculation='vc-relax'に変更
+  * <nowiki>&cellフィールド</nowiki>を追加
+  * （任意）press_conv_thrを指定
+<file - Si.vc-relax.in>
 &control
    calculation='vc-relax'
-   prefix='si'
+   prefix='Si'
-   pseudo_dir='./'
+   pseudo_dir='./pseudo/'
-   outdir = './'
+   outdir = './tmp/'
    etot_conv_thr = 1.d-5
    forc_conv_thr = 1.d-4
@@ Line 92: / Line 109: @@
 /
 &electrons
+   conv_thr = 1.0d-8
 /
 &ions
 /
 &cell
+  press_conv_thr = 0.1
 /
 ATOMIC_SPECIES
@@ Line 106: / Line 125: @@
 </file>
-今度はcalculation='vc-relax'とします。 また、新しい項目&cellが必要です。
+^変数^初期値^説明^
+|press_conv_thr|0.5 (kbar)|構造緩和の各ステップで、単位胞に加わる力がこの値よりも小さくなったときに収束したと判断する|
 実行は次のようにします。
-  $ pw.x < si.vc-relax.in > si.vc-relax.out
+  $ pw.x < Si.vc-relax.in > Si.vc-relax.out
 結果は次のとおりです。
-<file - si.vc-relax.out>
+<file - Si.vc-relax.out>
 （略）
      crystal axes: (cart. coord. in units of alat)
@@ Line 136: / Line 156: @@
 </file>
-もともと基本並進ベクトルが格子定数を単位としてa(1)=(-0.5, 0, 0.5)だったのが(-0.510132378, 0.000000000, 0.510132378)に拡大しています。 そのため、構造緩和によって得られた格子定数はボーア半径を単位として10.0*0.510132378/0.5=10.20・・・となります。
+もともと基本並進ベクトルが格子定数を単位としてa(1)=(-0.5, 0, 0.5)だったのが(-0.510132378, 0.000000000, 0.510132378)に拡大しています。
+そのため、構造緩和によって得られた格子定数はボーア半径を単位として 10.0 * 0.510132378 / 0.5 = 10.20... となります。
 ===== 注意 =====
-  * 実はvc-relax計算は原子位置の緩和も同時に行っています。ほとんどの場合はうまくいかないので、先にrelax計算を行ってからvc-relaxを実行しましょう。
+  * vc-relax計算は原子位置の緩和も同時に行っています。同時緩和はうまくいかないことがあるので、先にrelax計算を行った上でvc-relaxを実行しましょう。
-  * vc-relax計算で構造緩和を行うとき、計算に必要なパラメータは最初の格子定数と原子位置の計算から出したものを使い続けます。計算が正しいか確認するため、入力ファイルを得られたパラメータでもう一度vc-relax計算をすると良いでしょう。
+  * vc-relax計算で構造緩和を行うとき、計算に必要なパラメータは最初の格子定数と原子位置の計算から出したものを使い続けます。結果が変わらなくなるまでvc-relax計算を繰り返してください。
-  * 結果がエネルギーのカットオフやk点数に対してちゃんと収束しているか必ず確認しましょう（得られた構造に対してエネルギーのカットオフやk点数を増やしてもう一度構造緩和してみる）。
+  * ecutwfcやk点数は多めに必要ですので、これらのパラメタについて収束を確認するようにします。
-  * 構造緩和の結果、結晶が別の対称性になってしまう場合、プログラムが以下のエラーを出して止まってしまいます。その場合は&systemに対称性を考慮しないというオプションnosym=.true.を加えてください。
+  * conv_thrが大きいと、出力ファイルに"SCF correction compared to forces is large: reduce conv_thr to get better values"と注意が出てくるので小さくしましょう。
+  * 構造緩和の結果、結晶が別の対称性になってしまう場合、プログラムが以下のエラーを出して止まってしまいます。その場合は&systemに対称性を考慮しないというオプションnosym=.true.を加えてください（ただし計算時間が非常に増大します）。
+<code>
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    Error in routine checkallsym (2):
+    not orthogonal operation
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+</code>
-  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-      Error in routine checkallsym (2):
-      not orthogonal operation
-  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ Line 154: / Line 178: @@
 ====== 金属の場合 ======
-金属の場合はFermi面の情報が必要ですので、
+金属の場合はフェルミ準位の情報が必要ですので、&systemに
+  occupations='smearing', smearing='mp', degauss=0.02
+を追加してください。
+degaussの値はk点数に応じて適切な値を入れてください。
+また、k点数が十分たくさんとれる場合は
   occupations='tetrahedra_opt'
-を使用してください。
+を使用することもできます。
-このときのk点は、原点を含む
-  K_POINTS automatic
-12 12 0 0 0
-などとします。
 * http://qe-forge.org/pipermail/pw_forum/2017-October/114051.html
+====== 擬ポテンシャルについて ======
+  * 一般に構造最適化で得られた格子定数は、PZ型などのLDA計算では過小評価され、PBE型などのGGA計算では過大評価されます。
+  * 最近ではPBEを改善したPBESOLという交換相関ポテンシャルが実験の格子定数を比較的よく再現します。
+    * 格子定数を求めるのはいろいろな困難があるので、実際の運用では格子定数は測定値を使用し、内部座標についてのみ最適化を行うことも多いです。

雑多な記録

User Tools

Site Tools

Differences

Page Tools