This shows you the differences between two versions of the page.
Both sides previous revision Previous revision Next revision | Previous revision | ||
quantumespresso:構造緩和 [2018/02/22 20:33] koudai [注意] |
quantumespresso:構造緩和 [2024/03/30 23:58] (current) koudai |
||
---|---|---|---|
Line 2: | Line 2: | ||
構造緩和計算により、結晶の安定な構造を求めます。 | 構造緩和計算により、結晶の安定な構造を求めます。 | ||
- | |||
構造緩和には次の2通りの方法があります。 | 構造緩和には次の2通りの方法があります。 | ||
- | - 格子定数はそのままで、単位胞内の原子の位置についての構造緩和 | + | - 格子定数はそのままで、単位胞内の原子の内部座標の最適化 |
- | - 単位胞内の原子位置はそのままで、格子定数についての構造緩和 | + | - 内部座標とともに格子定数も変化させる構造緩和 |
通常は2つを繰り返すことで、安定な結晶構造を探します。 | 通常は2つを繰り返すことで、安定な結晶構造を探します。 | ||
以下ではケイ素を例に構造緩和計算を行います。 | 以下ではケイ素を例に構造緩和計算を行います。 | ||
+ | |||
+ | |||
+ | ===== 注意 ===== | ||
+ | |||
+ | * 現在(2024年3月)のところ、noncolin=.true.のときは原子に加わる力の計算ができません | ||
+ | |||
====== 原子位置の構造緩和 ====== | ====== 原子位置の構造緩和 ====== | ||
Line 16: | Line 21: | ||
格子定数はそのままに、単位胞内の原子の位置を緩和させます。 Siの単位胞内における原子の位置は、格子定数を単位として安定な構造で(0, | 格子定数はそのままに、単位胞内の原子の位置を緩和させます。 Siの単位胞内における原子の位置は、格子定数を単位として安定な構造で(0, | ||
- | <file - si.relax.in> | + | <file - Si.relax.in> |
& | & | ||
| | ||
- | | + | |
- | | + | |
- | | + | |
| | ||
| | ||
Line 31: | Line 36: | ||
/ | / | ||
& | & | ||
+ | | ||
/ | / | ||
&ions | &ions | ||
Line 44: | Line 50: | ||
構造緩和計算では、calculation=' | 構造緩和計算では、calculation=' | ||
- | etot_conv_thr 構造緩和の各ステップで、エネルギーの変化がこの値よりも小さくなったときに収束したと判断する。デフォルトは1.d-4だが、経験的にはもう少し厳しいほうが良い。 | ||
- | forc_conv_thr 構造緩和の各ステップで、原子に加わる力の変化がこの値よりも小さくなったときに収束したと判断する。デフォルトは1.d-3だが、経験的にはもう少し厳しいほうが良い。 | ||
- | エネルギーのカットオフecutwfcは普段よりも多めにとっておきます。 | + | ^変数^初期値^説明^ |
+ | |etot_conv_thr|1.d-4|構造緩和の各ステップで、エネルギーの変化がこの値よりも小さくなったときに収束したと判断する| | ||
+ | |forc_conv_thr|1.d-3|構造緩和の各ステップで、原子に加わる力の変化がこの値よりも小さくなったときに収束したと判断する| | ||
+ | DFTで求められる程度の原子位置の精度なので、多くの場合は収束の厳しさはここで指定した値で十分です。 | ||
また、relax計算の際は& | また、relax計算の際は& | ||
- | 原子位置について、原点にあるSiは動いてほしくないので、座標の最後に0 0 0を付け加えます。 | + | 原子位置について、原点にあるSiは動いてほしくないので、座標の最後に0 0 0を付け加えます。 |
+ | 第3成分だけ動いて欲しい場合は0 0 1とします。 | ||
+ | 何も指定していなければ1 1 1となります(すべての方向に動ける)。 | ||
実行は次のようにします。 | 実行は次のようにします。 | ||
- | $ 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 | Begin final coordinates | ||
Line 75: | Line 85: | ||
====== 格子定数の構造緩和 ====== | ====== 格子定数の構造緩和 ====== | ||
- | 今度は格子定数の構造緩和を行います。 | + | 今度は格子定数の構造緩和を行います。 |
+ | 先に原子位置の緩和を行ってください。 | ||
- | <file - si.vc-relax.in> | + | Siの格子定数はボーア半径単位で10.2程度ですが、ここから少しずらして構造緩和をさせてみます。 |
+ | Si.relax.inをコピーしてSi.vc-relax.inを作成します。 | ||
+ | 以下の変更を行います | ||
+ | * calculation=' | ||
+ | * < | ||
+ | * (任意)press_conv_thrを指定 | ||
+ | |||
+ | <file - Si.vc-relax.in> | ||
& | & | ||
| | ||
- | | + | |
- | | + | |
- | | + | |
| | ||
| | ||
Line 91: | Line 109: | ||
/ | / | ||
& | & | ||
+ | | ||
/ | / | ||
&ions | &ions | ||
/ | / | ||
&cell | &cell | ||
+ | press_conv_thr = 0.1 | ||
/ | / | ||
ATOMIC_SPECIES | ATOMIC_SPECIES | ||
Line 105: | Line 125: | ||
</ | </ | ||
- | 今度はcalculation=' | + | ^変数^初期値^説明^ |
+ | |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> |
(略) | (略) | ||
| | ||
Line 135: | Line 156: | ||
</ | </ | ||
- | もともと基本並進ベクトルが格子定数を単位としてa(1)=(-0.5, | + | もともと基本並進ベクトルが格子定数を単位としてa(1)=(-0.5, |
+ | そのため、構造緩和によって得られた格子定数はボーア半径を単位として 10.0 * 0.510132378 / 0.5 = 10.20... となります。 | ||
===== 注意 ===== | ===== 注意 ===== | ||
- | * 実はvc-relax計算は原子位置の緩和も同時に行っていますが、ほとんどの場合はうまくいかないのでrelax計算と組み合わせる必要があります。 | + | * vc-relax計算は原子位置の緩和も同時に行っています。同時緩和はうまくいかないことがあるので、先にrelax計算を行った上でvc-relaxを実行しましょう。 |
- | * vc-relax計算で構造緩和を行うとき、計算に必要なパラメータは最初の格子定数と原子位置の計算から出したものを使い続けます。計算が正しいか確認するため、入力ファイルを得られたパラメータでもう一度vc-relax計算をすると良いでしょう。 | + | * vc-relax計算で構造緩和を行うとき、計算に必要なパラメータは最初の格子定数と原子位置の計算から出したものを使い続けます。結果が変わらなくなるまでvc-relax計算を繰り返してください。 |
- | * 構造緩和の結果、結晶が別の対称性になってしまう場合、プログラムが以下のエラーを出して止まってしまいます。その場合は& | + | * ecutwfcやk点数は多めに必要ですので、これらのパラメタについて収束を確認するようにします。 |
+ | * conv_thrが大きいと、出力ファイルに" | ||
+ | * 構造緩和の結果、結晶が別の対称性になってしまう場合、プログラムが以下のエラーを出して止まってしまいます。その場合は& | ||
+ | < | ||
+ | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | ||
+ | Error in routine checkallsym (2): | ||
+ | not orthogonal operation | ||
+ | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | ||
+ | </ | ||
- | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | ||
- | Error in routine checkallsym (2): | ||
- | not orthogonal operation | ||
- | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | ||
Line 152: | Line 178: | ||
====== 金属の場合 ====== | ====== 金属の場合 ====== | ||
- | 金属の場合はFermi面の情報が必要ですので、 | + | 金属の場合はフェルミ準位の情報が必要ですので、& |
+ | |||
+ | occupations=' | ||
+ | |||
+ | を追加してください。 | ||
+ | degaussの値はk点数に応じて適切な値を入れてください。 | ||
+ | |||
+ | また、k点数が十分たくさんとれる場合は | ||
occupations=' | occupations=' | ||
- | を使用してください。 | + | を使用することもできます。 |
- | このときのk点は、原点を含む | + | |
- | K_POINTS automatic | + | |
- | 12 12 12 1 1 1 | + | |
- | などとします。 | + | |
* http:// | * http:// | ||
+ | |||
+ | |||
+ | ====== 擬ポテンシャルについて ====== | ||
+ | |||
+ | * 一般に構造最適化で得られた格子定数は、PZ型などのLDA計算では過小評価され、PBE型などのGGA計算では過大評価されます。 | ||
+ | * 最近ではPBEを改善したPBESOLという交換相関ポテンシャルが実験の格子定数を比較的よく再現します。 | ||
+ | * 格子定数を求めるのはいろいろな困難があるので、実際の運用では格子定数は測定値を使用し、内部座標についてのみ最適化を行うことも多いです。 |