This is an old revision of the document!
Siを例に、フォノンの簡単な計算を実行してみる。 事前にpw.xによる電子状態計算を行う。
&control calculation='scf', prefix='si' pseudo_dir = './', outdir='./' / &system ibrav=2, celldm(1)=10.26, nat=2, ntyp=1, ecutwfc=18.0 / &electrons / ATOMIC_SPECIES Si 28.086 Si.pz-vbc.UPF ATOMIC_POSITIONS (alat) Si 0.00 0.00 0.00 Si 0.25 0.25 0.25 K_POINTS automatic 8 8 8 1 1 1
$ pw.x < si.scf.in > si.scf.out
電子状態の計算が終われば、フォノンの振動数を計算できるようになる。 まず、k=(0,0,0)のΓ点におけるフォノンを計算してみる。 以下のファイルを用意する。
phonons of Si at Gamma &inputph tr2_ph=1.0d-14, prefix='Si', epsil=.true., outdir='./', fildyn='Si.dynG', / 0.0 0.0 0.0
第1行はコメント行とみなされ無視されます。 また、最後の行は計算する波数(デカルト座標で、単位は格子定数aを使って2π/a)を表します。
各引数の意味は次のとおりです。
変数 | 初期値 | 説明 |
---|---|---|
tr2_ph | 1.0d-12 | フォノンの計算の収束条件。経験的に1.0d-14を使ったほうがよい。 |
prefix | pwscf | si.scf.inで使用したプレフィックスと同じにする。 |
epsil | .false. | 誘電率を計算する。Γ点の計算に必要。ただし.true.にする場合は、半導体かつΓ点の計算でないとエラーが出る。 |
outdir | ./ | 出力ファイルの場所。 |
fildyn | matdyn | 結果であるdynamical matrixを出力するファイル名。 |
金属の場合は、epsil=.false.として、計算するk点を 0.01 0.0 0.0 のようにΓ点から少しずらしてやるとうまくいきます。 実行は次のようにします。
$ ph.x < Si.phG.in > Si.phG.out
結果はfildynで指定したもの(今の場合はSi.dynG)に出力されます。
(中略) ************************************************************************** freq ( 1) = -0.170404 [THz] = -5.684069 [cm-1] ( 0.274931 0.000000 -0.504256 0.000000 -0.412479 0.000000 ) ( 0.274931 0.000000 -0.504256 0.000000 -0.412479 0.000000 ) freq ( 2) = -0.170404 [THz] = -5.684069 [cm-1] ( -0.161796 0.000000 0.380828 0.000000 -0.573404 0.000000 ) ( -0.161796 0.000000 0.380828 0.000000 -0.573404 0.000000 ) freq ( 3) = -0.170404 [THz] = -5.684069 [cm-1] ( 0.631059 0.000000 0.317327 0.000000 0.032690 0.000000 ) ( 0.631059 0.000000 0.317327 0.000000 0.032690 0.000000 ) freq ( 4) = 15.296893 [THz] = 510.249417 [cm-1] ( -0.463984 0.000000 0.270574 0.000000 0.459900 0.000000 ) ( 0.463984 0.000000 -0.270574 0.000000 -0.459900 0.000000 ) freq ( 5) = 15.296893 [THz] = 510.249417 [cm-1] ( 0.527053 0.000000 0.327503 0.000000 0.339053 0.000000 ) ( -0.527053 0.000000 -0.327503 0.000000 -0.339053 0.000000 ) freq ( 6) = 15.296893 [THz] = 510.249417 [cm-1] ( -0.083269 0.000000 0.565271 0.000000 -0.416575 0.000000 ) ( 0.083269 0.000000 -0.565271 0.000000 0.416575 0.000000 ) **************************************************************************
フォノンの振動に対応する変位ベクトルはXcrysDenを使って見ることができます。 次のファイルを用意します。
&inputfil dyn ='Si.dynG', asr ='simple', lperm = .true., q(1)=1.0, q(2)=0.0, q(3)=0.0 /
Γ点のフォノンは q→0 の極限で定義されます。lpermは q=(q(1), q(2), q(3)) に沿って q→0 の極限をとることを指定します。
実行は次のようにします。
$ dynmat.x < Si.dynmat.in > Si.dynmat.out
すると dynmat.axsf というファイルができるので、以下の手順で変位ベクトルを表示します。
今度はX点(2π/a,0,0)のフォノンを調べる。 手順はさきほどと同じであるが、誘電率が計算できないことに注意する。
phonons of Si at X &inputph tr2_ph=1.0d-14, prefix='si', amass(1)=28.08, outdir='./', fildyn='si.dynX', / 1.0 0.0 0.0
$ ph.x < si.phX.in > si.phX.out
結果はsi.dynXに出力される。
(中略) ************************************************************************** freq ( 1) = 4.207901 [THz] = 140.360482 [cm-1] ( -0.000000 0.000000 -0.023972 0.000000 0.706700 0.000000 ) ( -0.000000 0.000000 0.706700 0.000000 -0.023972 0.000000 ) freq ( 2) = 4.207901 [THz] = 140.360482 [cm-1] ( -0.000000 0.000000 -0.706700 0.000000 -0.023972 0.000000 ) ( -0.000000 0.000000 -0.023972 0.000000 -0.706700 0.000000 ) freq ( 3) = 12.236850 [THz] = 408.177391 [cm-1] ( -0.999882 0.000000 0.000000 0.000000 0.000000 0.000000 ) ( -0.015351 0.000000 -0.000000 0.000000 0.000000 0.000000 ) freq ( 4) = 12.236850 [THz] = 408.177391 [cm-1] ( -0.015351 0.000000 0.000000 0.000000 0.000000 0.000000 ) ( 0.999882 0.000000 -0.000000 0.000000 0.000000 0.000000 ) freq ( 5) = 13.746712 [THz] = 458.540957 [cm-1] ( 0.000000 0.000000 -0.000000 0.000000 0.707107 0.000000 ) ( -0.000000 0.000000 -0.707107 0.000000 0.000000 0.000000 ) freq ( 6) = 13.746712 [THz] = 458.540957 [cm-1] ( 0.000000 0.000000 0.707107 0.000000 0.000000 0.000000 ) ( 0.000000 0.000000 0.000000 0.000000 -0.707107 0.000000 ) **************************************************************************
2重に縮退したものが3つあることがわかる。