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

メタン分子CH4を題材に、GW近似で準粒子エネルギーがどのように変化するかを見る。


====== 計算の手順 ======

GW近似計算の手順は、次の3ステップである。

  - pw.xによる自己無撞着計算
  - pw4gww.xによるGW近似計算に必要なデータの生成
  - gww.xによる実際のGW近似計算

GW近似計算が終われば、あとは通常の手続きでポストプロセスの計算（バンドや状態密度など）を行う。


===== pw.xによる自己無撞着計算 =====

PWscfは平面波展開によって基底状態を計算するので、周期的な配置のものしか扱うことができない。
メタンは分子なので、立方格子の各頂点にメタン分子をおいて格子定数を非常に大きくすることで、メタン一分子を考えることと等価になる。

擬ポテンシャルとしてH.pz-vbc.UPFとC.pz-vbc.UPFを使う。
GWLではノルム保存型の擬ポテンシャルのみが使えることに注意する。

電子配置は水素が1s1、炭素が[He]2s22p2なので、価電子数は全部で8である。
つまり、スピンの縮退を考慮すればフェルミ準位より下にエネルギー準位は4つできるはずだが、フェルミ準位よりも上の状態も見るためにnbnd=5とする。

なお、分子なので分散はないから、k点はガンマ点だけでよい（k点を指定しない）。

入力ファイルは次のようになる。

<file - methane.scf.in>
&control
   calculation = 'scf',
   restart_mode='from_scratch',
   prefix='ch4',
   tprnfor = .true.,
   pseudo_dir = './',
   outdir='./tmp/'
/
&system
   ibrav = 1,
   celldm(1) = 15.0,
   nat = 5,
   ntyp= 2,
   ecutwfc = 40.0,
   nbnd = 5
/
&electrons
   diagonalization = 'cg'
   mixing_beta = 0.5,
   conv_thr = 1.0d-8
/
ATOMIC_SPECIES
 H  1.0   H.pz-vbc.UPF
 C  12.0  C.pz-vbc.UPF
ATOMIC_POSITIONS {bohr}
 H        1.198204546   1.198204546   1.198204546
 H       -1.198204546  -1.198204546   1.198204546
 H        1.198204546  -1.198204546  -1.198204546
 H       -1.198204546   1.198204546  -1.198204546
 C        0.000000000   0.000000000   0.000000000
</file>

PWscfを実行する。

  $ pw.x < methane.scf.in > methane.scf.out


===== pw4gww.xによるGW近似計算に必要なデータの生成 =====

次に、GW近似に必要なデータの生成を行う。 ここでは必要最低限、指定しなければならないもののみを指定した。

入力ファイルを作成する。

<file - methane.pw4gww.in>
&inputpw4gww    
   prefix = 'ch4'
   num_nbndv(1) = 4
   num_nbnds = 5
   l_truncated_coulomb = .true.
   truncation_radius = 7.5d0
   numw_prod = 100
   pseudo_dir = './',
   outdir = './tmp/'
/
</file>

ここで指定しているパラメータは以下のとおりである。

^変数^説明^
|prefix|プレフィックスはpw.xで計算したものと同じにする|
|num_nbndv(1)|価電子数。電子が全て詰まっているバンドの数と部分的に詰まっているバンドの数の合計を指定する。今の場合は4である。スピン分極があって、アップスピンとダウンスピンで価電子数が異なる場合はnum_nbndv(2)も指定する。|
|num_nbnds|バンドの数。pw.xの計算で指定した数と同じにする。|
|l_truncated_coulomb|クーロン相互作用を長距離で切り捨てるかどうかの指定。もし.true.であればtruncation_radiusも指定する。分子の場合はこのように切り捨てて構わないが、結晶を考える場合は.false.を指定し、事前にhead.xを別途実行する必要がある（マニュアル参照）。|
|truncation_radius|クーロン相互作用を切り捨てるときの距離。単位はBohr。クーロン相互作用は長距離力なので、分子を並べて計算するときにお互いが相互作用しあってしまっては困るので、このようなカットオフを入れる必要がある。|
|numw_prod|polarizability basisの次元を指定。大きいほど結果は正確になるが、計算に時間がかかる。詳しくは下リンク先のチュートリアル参照。|

詳しくはマニュアルを参照のこと。

  http://www.gwl-code.org/tutorial.html

入力ファイルができたら、pw4gww.xの実行を行う。

  $ pw4gww.x < methane.pw4gww.in > methane.pw4gww.out


===== gww.xによる実際のGW近似計算 =====

最後に、実際のGW近似計算を行う。

次の入力ファイルを用意する。

<file - methane.gww.in>
&inputgww
ggwin%prefix='ch4'
ggwin%max_i=5
ggwin%i_min=1
ggwin%i_max=5
ggwin%omega=20
ggwin%n=118
ggwin%grid_freq=5
ggwin%second_grid_i=3
ggwin%second_grid_n=10
ggwin%omega_fit=20
ggwin%n_grid_fit=240
ggwin%n_fit=120
ggwin%l_truncated_coulomb=.true.
ggwin%outdir='./tmp/'
/
</file>

ここで指定しているパラメータは以下のとおりである。

^変数^説明^
|ggwin%max_i              |pw.x計算のときに指定したnbndと同じにする。|
|ggwin%i_min, ggwin%i_max |バンドのi_min番目からi_max番目を計算。普通はi_minを1に、i_maxをmax_iと同じにする。|
|ggwin%omega              |虚数周波数の幅。単位はRy。|
|ggwin%n                  |omegaの刻み数。|
|ggwin%grid_freq          |omegaの刻み方。3なら等間隔、5ならomega=0付近に細かくとる。|
|ggwin%second_grid_i      |ggwin%grid_freq=5のとき、omegaの刻みの1番目からggwin%second_grid_i番目までを細かく刻むことにする。|
|ggwin%second_grid_n      |ggwin%grid_freq=5のときの、細かく刻む部分でもともとの刻みをさらに何分割にするか。分割数はggwin%second_grid_nで指定した数の2倍。|
|ggwin%omega_fit          |自己エネルギーを計算するときの虚数周波数の幅。単位はRy。|
|ggwin%n_grid_fit         |自己エネルギーを計算するときの虚数周波数の刻み数。|
|ggwin%n_fit              |多極展開するときの刻み数。|
|ggwin%l_truncated_coulomb|クーロン相互作用を切り捨てるか。pw4gww.xで計算するときと同じにする|

詳しくはマニュアルを参照のこと。

gww.xの実行は次のようにする。

  $ gww.x < methane.gww.in > methane.gww.out


===== エネルギー準位の計算 =====

部分状態密度を計算することで、各軌道のエネルギー準位を確認する。

tmp/ch4-bands.datを、作業ディレクトリにbands.datという名前でコピーする。

  $ cp tmp/ch4-bands.dat bands.dat

次に、以下の入力ファイルmethane.pdos.inを作成する。

<file - methane.pdos.in>
&projwfc
  outdir='./tmp/'
  prefix='ch4'
  lgww=.true.
  Emin=-30.0, Emax=5.0, DeltaE=0.05,
  ngauss=0, degauss=0.01559
/
</file>

オプションlgww=.true.によって、bands.datを読み込んでGW近似における部分状態密度を計算する。
部分状態密度を計算するため、projwfc.xを実行する。

  $ projwfc.x < methane.pdos.in > methane.pdos.out

すると各軌道の部分状態密度が得られる。
GWLは tmp/ch4.save/ の中身を変更しないので、lgww=.true.を削除してprojwfc.xを実行すると、GW近似をする前の部分状態密度を得られる。