====== 概要 ======
* GaAsを例に、ガンマ点におけるフォノンの縦光学(LO)フォノンと横光学フォノン(TO)の分裂を扱います
* 極性半導体ではLO-TO分裂が起こるので、この計算は必須です。
====== 計算手順 ======
===== 準備 =====
GaAsのCIFファイルは以下のサイトのものを使い、格子定数はvc-relax計算で得られたものに変更しました
* https://materials.springer.com/isp/crystallographic/docs/sd_0311662
* vc-relax計算で得られた格子定数は 5.663 Å = 10.702 Bohr
擬ポテンシャルはPS Libraryからウルトラソフト型を利用しました。
* 擬ポテンシャル ... https://www.quantum-espresso.org/pseudopotentials/ps-library/
* As.pbesol-n-rrkjus_psl.1.0.0.UPF
* Ga.pbesol-dnl-rrkjus_psl.1.0.0.UPF
第一原理計算にはQEを使用します。
入力ファイルは以下のように作成します。
$ cif2cell -p pwscf -f GaAs.cif -o GaAs.scf.in
===== ボルン有効電荷の計算 =====
DFPTを用いてボルン有効電荷を計算します。
SCF計算で基底状態を求めたあと、ph.xを使って計算します
&control
prefix = 'GaAs',
calculation = 'scf',
pseudo_dir = './pseudo/',
outdir = './tmp/',
/
&SYSTEM
ibrav = 0
nat = 2
ntyp = 2
ecutwfc = 40.0
ecutrho = 240.0
/
&electrons
conv_thr = 1.0d-10
/
CELL_PARAMETERS {angstrom}
2.8314 2.8314 0.0000
2.8314 0.0000 2.8314
0.0000 2.8314 2.8314
ATOMIC_SPECIES
As 74.92100 As.pbesol-n-rrkjus_psl.1.0.0.UPF
Ga 69.72300 Ga.pbesol-dnl-rrkjus_psl.1.0.0.UPF
ATOMIC_POSITIONS {crystal}
As 0.2500000000 0.2500000000 0.2500000000
Ga 0.0000000000 0.0000000000 0.0000000000
K_POINTS automatic
8 8 8 1 1 1
phonons of GaAs
&inputph
tr2_ph = 1.0d-14
prefix = 'GaAs'
outdir = './tmp/'
fildyn = 'GaAs.dynG'
epsil = .true.
zeu = .true.
trans = .false.
/
0.0 0.0 0.0
epsilは誘電率を、zeuは原子に加わる力からボルン有効電荷を計算するように指定します。
このほか、分極からボルン有効電荷を計算するzueもありますが、どちらを使ってもほぼ同じ値になります。
また、フォノンの振動数の計算は必要ないので、transは.false.にしておきます。
$ pw.x < GaAs.scf.in > GaAs.scf.out
$ ph.x < GaAs.phG.in > GaAs.phG.out
誘電率およびボルン有効電荷は次のように出力されます。
（略）
Dielectric constant in cartesian axis
( 14.192331204 0.000000000 -0.000000000 )
( 0.000000000 14.192331204 -0.000000000 )
( -0.000000000 -0.000000000 14.192331204 )
Effective charges (d Force / dE) in cartesian axis
atom 1 As
Ex ( -2.13215 -0.00000 0.00000 )
Ey ( -0.00000 -2.13215 -0.00000 )
Ez ( -0.00000 -0.00000 -2.13215 )
atom 2 Ga
Ex ( 2.14310 0.00000 -0.00000 )
Ey ( -0.00000 2.14310 0.00000 )
Ez ( 0.00000 0.00000 2.14310 )
（略）
===== Dynamical matrixの計算 =====
QE用の超格子の入力ファイルを作成します
$ cif2cell --supercell=[3,3,3] -p pwscf -f GaAs.cif -o GaAs_333.pw.in
超格子のインプットファイルに足りない部分を書き加えます。
cif2cellのデフォルトで格子定数はA（オングストローム単位）で指定されますが、celldm(1)（ボーア半径単位）に直しておきます。
また、原子に加わる力の大きさを出力するために tprnfor=.true. を忘れないようにします。
&control
prefix = 'GaAs',
calculation = 'scf',
pseudo_dir = './pseudo/',
outdir = './tmp/',
tprnfor=.true.
/
&SYSTEM
ibrav = 0
celldm(1) = 10.702
nat = 54
ntyp = 2
/
&electrons
conv_thr = 1.0d-8
/
CELL_PARAMETERS {alat}
1.500000000000000 1.500000000000000 0.000000000000000
1.500000000000000 0.000000000000000 1.500000000000000
0.000000000000000 1.500000000000000 1.500000000000000
ATOMIC_SPECIES
Ga 69.72300 Ga.pbesol-dnl-rrkjus_psl.1.0.0.UPF
As 74.92100 As.pbesol-n-rrkjus_psl.1.0.0.UPF
K_POINTS automatic
4 4 4 1 1 1
ATOMIC_POSITIONS {crystal}
As 0.083333333333333 0.083333333333333 0.083333333333333
As 0.416666666666667 0.083333333333333 0.083333333333333
（略）
Ga 0.666666666666667 0.333333333333333 0.333333333333333
また、ALAMODE用の原子間力定数の計算の入力ファイルを作成します
&general
PREFIX = GaAs_333
MODE = suggest
NAT = 54; NKD = 2
KD = As Ga
/
&interaction
NORDER = 1
/
&cell
10.702
1.500000000000000 1.500000000000000 0.000000000000000
1.500000000000000 0.000000000000000 1.500000000000000
0.000000000000000 1.500000000000000 1.500000000000000
/
&cutoff
*-* None
/
&position
1 0.083333333333333 0.083333333333333 0.083333333333333
1 0.416666666666667 0.083333333333333 0.083333333333333
（略）
2 0.666666666666667 0.333333333333333 0.333333333333333
/
できたらALMを実行します。
$ alm GaAs_333.harmonic.in > GaAs_333.harmonic.log
GaAs_333.pattern_HARMONIC というファイルができるので、ALAMODEに付属のPythonスクリプトを使ってQEの入力ファイルを作成します。
$ python -m displace --QE GaAs_333.pw.in --mag 0.01 --prefix disp -pf GaAs_333.pattern_HARMONIC
今回の場合はQEの入力ファイルが2つできます。
出てきた入力ファイルをすべてQEで実行します
$ pw.x < disp1.pw.in > disp1.pw.out
$ pw.x < disp2.pw.in > disp2.pw.out
* 計算は重いのでOpenMPやMPIによる並列計算をおすすめします
これをALAMODEに付属のPythonスクリプトで整理します。
$ python -m extract --QE GaAs_333.pw.in *.pw.out > DFSET_harmonic
ALMを使って原子間力定数を計算します。
GaAs_333.harminic.inをコピーして GaAs.harmonic_opt.in を生成し、次の点を変更します
* &generalフィールドで MODE=optimize に変更
* &optimizeフィールドを作成して、中に DFSET = DFSET_harmonic と記述
&general
PREFIX = GaAs_333
MODE = optimize
NAT = 54; NKD = 2
KD = As Ga
/
&interaction
NORDER = 1
/
&cell
10.702
1.500000000000000 1.500000000000000 0.000000000000000
1.500000000000000 0.000000000000000 1.500000000000000
0.000000000000000 1.500000000000000 1.500000000000000
/
&cutoff
*-* None
/
&optimize
DFSET = DFSET_harmonic
/
&position
1 0.083333333333333 0.083333333333333 0.083333333333333
1 0.416666666666667 0.083333333333333 0.083333333333333
（略）
2 0.666666666666667 0.333333333333333 0.333333333333333
/
ALMを実行します
$ alm GaAs_333.harmonic_opt.in > GaAs.harmonic_opt.log
GaAs_333.fcs と GaAs_333.xml が出力されたら成功です。
===== フォノンの分散 =====
DFPT計算で得た GaAs.phG.out から誘電率とボルン有効電荷を記録したファイルを作ります
14.200395680 -0.000000000 -0.000000000
-0.000000000 14.200395680 -0.000000000
-0.000000000 -0.000000000 14.200395680
-2.12591 -0.00000 0.00000
-0.00000 -2.12591 -0.00000
0.00000 0.00000 -2.12591
2.14222 0.00000 0.00000
0.00000 2.14222 0.00000
0.00000 0.00000 2.14222
* 誘電率、各原子のボルン有効電荷の順に書きます
* 各原子の有効ボルン電荷の順番は、以下のanphonを実行することで得られる GaAs.phband.log 内のものに合わせてください
（略）
Atomic positions in the primitive cell (fractional):
1: 2.500000e-01 2.500000e-01 2.500000e-01 As
2: 0.000000e+00 0.000000e+00 0.000000e+00 Ga
（略）
ボルン有効電荷なし
&general
PREFIX = GaAs_333
MODE = phonons
FCSXML = GaAs_333.xml
NKD = 2
KD = As Ga
/
&cell
10.702
0.500000000000000 0.500000000000000 0.000000000000000
0.500000000000000 0.000000000000000 0.500000000000000
0.000000000000000 0.500000000000000 0.500000000000000
/
&kpoint
1
G 1.00 1.00 1.00 X 0.50 0.50 1.00 51
X 0.50 0.50 1.00 G 0.00 0.00 0.00 51
G 0.00 0.00 0.00 L 0.50 0.50 0.50 51
L 0.50 0.50 0.50 X 0.50 0.50 0.00 51
X 0.50 0.50 0.00 W 0.75 0.50 0.25 51
W 0.75 0.50 0.25 L 0.50 0.50 0.50 51
/
ボルン有効電荷あり
&general
PREFIX = GaAs_333
MODE = phonons
FCSXML = GaAs_333.xml
NONANALYTIC = 3
BORNINFO = GaAs.born
NKD = 2
KD = As Ga
/
&cell
10.702
0.500000000000000 0.500000000000000 0.000000000000000
0.500000000000000 0.000000000000000 0.500000000000000
0.000000000000000 0.500000000000000 0.500000000000000
/
&kpoint
1
G 1.00 1.00 1.00 X 0.50 0.50 1.00 51
X 0.50 0.50 1.00 G 0.00 0.00 0.00 51
G 0.00 0.00 0.00 L 0.50 0.50 0.50 51
L 0.50 0.50 0.50 X 0.50 0.50 0.00 51
X 0.50 0.50 0.00 W 0.75 0.50 0.25 51
W 0.75 0.50 0.25 L 0.50 0.50 0.50 51
/
$ anphon GaAs.phband.in > GaAs.phband.log
{{ :alamode:gaas_phonon-disp.png?400 |}}
赤い実線がボルン有効電荷あり、青い破線がなしです。
ガンマ点でTOフォノンとLOフォノンが分裂するのがわかります。