構造緩和計算により、結晶の安定な構造を求めます。

構造緩和には次の2通りの方法があります。

1. 格子定数はそのままで、単位胞内の原子の位置についての構造緩和
2. 単位胞内の原子位置はそのままで、格子定数についての構造緩和

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

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

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

si.relax.in

&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
/
&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だが、経験的にはもう少し厳しいほうが良い。

エネルギーのカットオフecutwfcは普段よりも多めにとっておきます。

また、relax計算の際は&ionsという項目が必要です。 ここに必要に応じて構造緩和計算のオプションを書き加えますが、今はデフォルトのままとします。

原子位置について、原点にあるSiは動いてほしくないので、座標の最後に0 0 0を付け加えます。 他にも動いてほしくない原子があれば、同様に書き加えます。

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

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

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

  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程度ですが、ここから少しずらして構造緩和をさせてみます。

  si.vc-relax.in

&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 / &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

結果は次のとおりです。

  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)

1. 0.510132378 0.000000000 0.510132378
2. 0.000000000 0.510132378 0.510132378
3. 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・・・となります。 ↑ 注意 †

  実はvc-relax計算は原子位置の緩和も同時に行っていますが、ほとんどの場合はうまくいかないのでrelax計算と組み合わせる必要があります。
  vc-relax計算で構造緩和を行うとき、計算に必要なパラメータは最初の格子定数と原子位置の計算から出したものを使い続けます。計算が正しいか確認するため、入力ファイルを得られたパラメータでもう一度vc-relax計算をすると良いでしょう。
  構造緩和の結果、結晶が別の対称性になってしまう場合、プログラムが以下のエラーを出して止まってしまいます。その場合は&systemに対称性を考慮しないというオプションnosym=.true.を加えてください。

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      Error in routine checkallsym (2):
      not orthogonal operation
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%