This is an old revision of the document!
構造緩和計算により、結晶の安定な構造を求めます。
構造緩和には次の2通りの方法があります。
通常は2つを繰り返すことで、安定な結晶構造を探します。
以下ではケイ素を例に構造緩和計算を行います。
格子定数はそのままに、単位胞内の原子の位置を緩和させます。 Siの単位胞内における原子の位置は、格子定数を単位として安定な構造で(0,0,0)と(0.25,0.25,0.25)にありますが、原子を少しずらした点において構造緩和をさせてみます。
&control calculation='relax' prefix='si' pseudo_dir='./' outdir = './' etot_conv_thr = 1.d-5 forc_conv_thr = 1.d-4 / &system ibrav= 2, celldm(1) =10.2, nat= 2, ntyp= 1, ecutwfc = 30.0, / &electrons conv_thr = 1.0d-8 / &ions / ATOMIC_SPECIES Si 28.086 Si.pz-vbc.UPF ATOMIC_POSITIONS Si 0.00 0.00 0.00 0 0 0 Si 0.22 0.23 0.24 K_POINTS automatic 6 6 6 1 1 1
構造緩和計算では、calculation='relax'を指定します。 関連するオプションの意味は次の通りです。
変数 | 初期値 | 説明 |
---|---|---|
etot_conv_thr | 1.d-4 | 構造緩和の各ステップで、エネルギーの変化がこの値よりも小さくなったときに収束したと判断する。初期値のままだと粗いので小さくする。 |
forc_conv_thr | 1.d-3 | 構造緩和の各ステップで、原子に加わる力の変化がこの値よりも小さくなったときに収束したと判断する。こちらも初期値だと粗い。 |
また、relax計算の際は&ionsという項目が必要です。 ここに必要に応じて構造緩和計算のオプションを書き加えますが、今はデフォルトのままとします。
原子位置について、原点にあるSiは動いてほしくないので、座標の最後に0 0 0を付け加えます。 第3成分だけ動いて欲しい場合は0 0 1とします。 何も指定していなければ1 1 1となります(すべての方向に動ける)。
実行は次のようにします。
$ pw.x < si.relax.in > si.relax.out
結果は次のようになります。
(略) Begin final coordinates ATOMIC_POSITIONS (alat) Si 0.000000000 0.000000000 0.000000000 0 0 0 Si 0.250003179 0.249999041 0.249996513 End final coordinates (略)
原子位置が(0.25,0.25,0.25)付近に動いたことが確認できます。
今度は格子定数の構造緩和を行います。 先に原子位置の緩和を行ってください。
Siの格子定数はボーア半径単位で10.2程度ですが、ここから少しずらして構造緩和をさせてみます。
&control calculation='vc-relax' prefix='si' pseudo_dir='./' outdir = './' etot_conv_thr = 1.d-5 forc_conv_thr = 1.d-4 / &system ibrav= 2, celldm(1) =10.0, nat= 2, ntyp= 1, ecutwfc = 30.0, / &electrons conv_thr = 1.0d-8 / &ions / &cell / ATOMIC_SPECIES Si 28.086 Si.pz-vbc.UPF ATOMIC_POSITIONS Si 0.00 0.00 0.00 0 0 0 Si 0.25 0.25 0.25 K_POINTS automatic 6 6 6 1 1 1
今度はcalculation='vc-relax'とします。 また、新しい項目&cellが必要です。
実行は次のようにします。
$ pw.x < si.vc-relax.in > si.vc-relax.out
結果は次のとおりです。
(略) crystal axes: (cart. coord. in units of alat) a(1) = ( -0.500000 0.000000 0.500000 ) a(2) = ( 0.000000 0.500000 0.500000 ) a(3) = ( -0.500000 0.500000 0.000000 ) (略) Begin final coordinates new unit-cell volume = 265.50864 a.u.^3 ( 39.34432 Ang^3 ) CELL_PARAMETERS (alat= 10.00000000) -0.510132378 0.000000000 0.510132378 -0.000000000 0.510132378 0.510132378 -0.510132378 0.510132378 -0.000000000 ATOMIC_POSITIONS (alat) Si 0.000000000 0.000000000 0.000000000 0 0 0 Si 0.255066189 0.255066189 0.255066189 End final coordinates (略)
もともと基本並進ベクトルが格子定数を単位としてa(1)=(-0.5, 0, 0.5)だったのが(-0.510132378, 0.000000000, 0.510132378)に拡大しています。 そのため、構造緩和によって得られた格子定数はボーア半径を単位として 10.0 * 0.510132378 / 0.5 = 10.20… となります。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Error in routine checkallsym (2): not orthogonal operation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
金属の場合はFermi面の情報が必要ですので、
occupations='tetrahedra_opt'
を使用してください。 このときのk点は、原点を含む
K_POINTS automatic 12 12 12 0 0 0
などとします。
* http://qe-forge.org/pipermail/pw_forum/2017-October/114051.html