====== 概要 ======

構造緩和計算により、結晶の安定な構造を求めます。
構造緩和には次の2通りの方法があります。

  - 格子定数はそのままで、単位胞内の原子の内部座標の最適化
  - 内部座標とともに格子定数も変化させる構造緩和

通常は2つを繰り返すことで、安定な結晶構造を探します。

以下ではケイ素を例に構造緩和計算を行います。


===== 注意 =====

  * 現在（2024年3月）のところ、noncolin=.true.のときは原子に加わる力の計算ができません


====== 原子位置の構造緩和 ======

格子定数はそのままに、単位胞内の原子の位置を緩和させます。 Siの単位胞内における原子の位置は、格子定数を単位として安定な構造で(0,0,0)と(0.25,0.25,0.25)にありますが、原子を少しずらした点において構造緩和をさせてみます。

<file - Si.relax.in>

&control
   calculation='relax'
   prefix='Si'
   pseudo_dir='./pseudo/'
   outdir = './tmp/'
   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
</file>

構造緩和計算では、calculation='relax'を指定します。 関連するオプションの意味は次の通りです。

^変数^初期値^説明^
|etot_conv_thr|1.d-4|構造緩和の各ステップで、エネルギーの変化がこの値よりも小さくなったときに収束したと判断する|
|forc_conv_thr|1.d-3|構造緩和の各ステップで、原子に加わる力の変化がこの値よりも小さくなったときに収束したと判断する|

DFTで求められる程度の原子位置の精度なので、多くの場合は収束の厳しさはここで指定した値で十分です。
また、relax計算の際は&ionsという項目が必要です。 ここに必要に応じて構造緩和計算のオプションを書き加えますが、今はデフォルトのままとします。

原子位置について、原点にあるSiは動いてほしくないので、座標の最後に0 0 0を付け加えます。
第3成分だけ動いて欲しい場合は0 0 1とします。
何も指定していなければ1 1 1となります（すべての方向に動ける）。


実行は次のようにします。

  $ pw.x < Si.relax.in > Si.relax.out

結果は次のようになります。

<file - 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
（略）
</file>

原子位置が(0.25,0.25,0.25)付近に動いたことが確認できます。


====== 格子定数の構造緩和 ======

今度は格子定数の構造緩和を行います。
先に原子位置の緩和を行ってください。

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'
   pseudo_dir='./pseudo/'
   outdir = './tmp/'
   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
  press_conv_thr = 0.1
/
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
</file>

^変数^初期値^説明^
|press_conv_thr|0.5 (kbar)|構造緩和の各ステップで、単位胞に加わる力がこの値よりも小さくなったときに収束したと判断する|

実行は次のようにします。

  $ pw.x < Si.vc-relax.in > Si.vc-relax.out

結果は次のとおりです。

<file - 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
（略）
</file>

もともと基本並進ベクトルが格子定数を単位として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計算で構造緩和を行うとき、計算に必要なパラメータは最初の格子定数と原子位置の計算から出したものを使い続けます。結果が変わらなくなるまでvc-relax計算を繰り返してください。
  * ecutwfcやk点数は多めに必要ですので、これらのパラメタについて収束を確認するようにします。
  * 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>




====== 金属の場合 ======

金属の場合はフェルミ準位の情報が必要ですので、&systemに

  occupations='smearing', smearing='mp', degauss=0.02

を追加してください。
degaussの値はk点数に応じて適切な値を入れてください。

また、k点数が十分たくさんとれる場合は

  occupations='tetrahedra_opt'

を使用することもできます。


* http://qe-forge.org/pipermail/pw_forum/2017-October/114051.html



====== 擬ポテンシャルについて ======

  * 一般に構造最適化で得られた格子定数は、PZ型などのLDA計算では過小評価され、PBE型などのGGA計算では過大評価されます。
  * 最近ではPBEを改善したPBESOLという交換相関ポテンシャルが実験の格子定数を比較的よく再現します。
    * 格子定数を求めるのはいろいろな困難があるので、実際の運用では格子定数は測定値を使用し、内部座標についてのみ最適化を行うことも多いです。