• Siを例にフォノンの分散と状態密度を計算します。
• 計算には超格子の作成が必要です。
• 事前に cif2cell をインストールしてください。
  • CIF2Cell

SiのCIFファイルと擬ポテンシャルを入手します。 今回は結晶構造は Crystallography Open Database から、擬ポテンシャルはQEのウェブサイトから入手しました。

• CIFファイル … http://www.crystallography.net/cod/9008565.html
• 擬ポテンシャル … https://www.quantum-espresso.org/pseudopotentials/ps-library/si

擬ポテンシャルは Si.pbe-n-kjpaw_psl.1.0.0.UPF を使用します

CIFファイル名を9008565.cifとします。 cif2structを使って、普通の入力ファイルと、超格子の入力ファイルを作成します。

$ cif2cell -p pwscf -f 9008565.cif -o Si.pw.in
$ cif2cell --supercell=[3,3,3] -p pwscf -f 9008565.cif -o Si333.pw.in

• 超格子のサイズは最初は2×2×2など小さいものから、結果が収束するまで徐々に大きくしていきます

超格子のインプットファイルに足りない部分を書き加えます。 cif2cellのデフォルトで格子定数はA（オングストローム単位）で指定されますが、celldm(1)（ボーア半径単位）に直しておきます。 また、原子に加わる力の大きさを出力するために tprnfor=.true. を忘れないようにします。

Si333.pw.in

&control
   calculation='scf',
   prefix='Si',
   pseudo_dir = './pseudo/',
   outdir='./tmp/'
   tprnfor=.true.
/

&SYSTEM
  ibrav = 0
  celldm(1) = 10.2625
  nat = 54
  ntyp = 1
  ecutwfc = 40
  ecutrho = 320
/

&electrons
  conv_thr = 1.0d-8
/

K_POINTS automatic
 4 4 4 1 1 1

CELL_PARAMETERS {alat}
  1.500000000000000   1.500000000000000   0.000000000000000 
  1.500000000000000   0.000000000000000   1.500000000000000 
  0.000000000000000   1.500000000000000   1.500000000000000 
ATOMIC_SPECIES
  Si   28.08500  Si.pbe-n-kjpaw_psl.1.0.0.UPF
ATOMIC_POSITIONS {crystal}
Si   0.000000000000000   0.000000000000000   0.000000000000000 
Si   0.083333333333333   0.083333333333333   0.083333333333333 
（略）
Si   0.750000000000000   0.416666666666667   0.416666666666667

フォノンの計算に必要な原子間力定数(Interatomic force constant)の計算を行います

超格子の原子位置を少しだけずらしたQEの入力ファイルを作成します。 &cellフィールドと&positionフィールドには、さきほどの超格子のインプットファイルで得られた基本並進ベクトルと原子位置をはりつけてください。

Si.alm0.in

&general
  PREFIX = Si333
  MODE = suggest
  NAT = 54; NKD = 1
  KD = Si
/

&interaction
  NORDER = 1
/

&cell
  10.2625
  1.500000000000000   1.500000000000000   0.000000000000000 
  1.500000000000000   0.000000000000000   1.500000000000000 
  0.000000000000000   1.500000000000000   1.500000000000000 
/

&cutoff 
  Si-Si None
/

&position
1   0.000000000000000   0.000000000000000   0.000000000000000 
1   0.083333333333333   0.083333333333333   0.083333333333333 
（略）
1   0.750000000000000   0.416666666666667   0.416666666666667 
/

できたらALMを実行します。

$ alm Si.alm0.in > Si.alm0.log

すると Si333.pattern_HARMONIC というファイルができるので、ALAMODEに付属のPythonスクリプトを使ってQEの入力ファイルを作成します。

$ python -m displace --QE Si333.pw.in --mag 0.01 --prefix disp -pf Si333.pattern_HARMONIC

• mag … 原子をどの程度ずらすかで、0.01から0.04程度の値（単位はオングストローム）が推奨されていますので、いろいろ試してください
• prefix … 生成されるQEの入力ファイル名に使われますので、自由に設定してください

今回の場合はQEの入力ファイルが1つだけ(disp1.pw.in)できます。 出てきた入力ファイルをすべてQEで実行します

$ pw.x < disp1.pw.in > disp1.pw.out

• 超格子の計算なので時間がかかります。並列計算による実行をお勧めします。

disp1.pw.outに各原子にはたらく力が出力されていますので、これをALAMODEに付属のPythonスクリプトで整理します。

$ python -m extract --QE Si333.pw.in *.pw.out > DFSET_harmonic

最後にALMを使って原子間力定数を計算します。 Si.alm0.inをコピーしてSi.alm1.inを生成し、次の点を変更します

• &generalフィールドでMODE=optimizeに変更
• &optimizeフィールドを作成して、中に DFSET = DFSET_harmonic と記述

Si.alm1.in

&general
  PREFIX = Si333
  MODE = optimize
  NAT = 54; NKD = 1
  KD = Si
/

&interaction
  NORDER = 1
/

&optimize
  DFSET = DFSET_harmonic
/

&cell
  10.2625
  1.500000000000000   1.500000000000000   0.000000000000000 
  1.500000000000000   0.000000000000000   1.500000000000000 
  0.000000000000000   1.500000000000000   1.500000000000000 
/

&cutoff 
  Si-Si None
/

&position
1   0.000000000000000   0.000000000000000   0.000000000000000 
1   0.083333333333333   0.083333333333333   0.083333333333333 
（略） 
1   0.750000000000000   0.416666666666667   0.416666666666667 
/

できたらALMを実行します

$ alm Si.alm1.in > Si.alm1.log

Si333.fcsとSi333.xmlが出力されたら成功です。

Si_phband.in

&general
  PREFIX = Si333
  MODE   = phonons
  FCSXML = Si333.xml

  NKD = 1; KD = Si
  MASS = 28.0855
/

&cell
  10.2625
  0.500000000000000   0.500000000000000   0.000000000000000 
  0.500000000000000   0.000000000000000   0.500000000000000 
  0.000000000000000   0.500000000000000   0.500000000000000 
/

&kpoint
  1
  G 0.0 0.0 0.0 X 0.5 0.5 0.0 51
  X 0.5 0.5 1.0 G 0.0 0.0 0.0 51
  G 0.0 0.0 0.0 L 0.5 0.5 0.5 51
/

• &cellフィールドは、元の結晶の基本並進ベクトルです。Si.pw.inのものをコピーしてください
• &kpointフィールドの1行目の 1 は、バンド分散の計算を指定しています。2行目以降は2π/aを単位とする分率座標で指定した波数の経路です

$ anphon Si_phband.in > Si_phband.out

フォノンの分散は Si333.bands に出力されます

Si_phband.in をコピーして si_phdos.in というファイルを作成し、&kpointフィールドを次のようにかきかえます。

Si_phdos.in

&general
  PREFIX = Si333
  MODE   = phonons
  FCSXML = Si333.xml

  NKD = 1; KD = Si
  MASS = 28.0855
/

&cell
  10.2625
  0.0 0.5 0.5
  0.5 0.0 0.5
  0.5 0.5 0.0
/

&kpoint
  2
  20 20 20
/

• &kpointフィールドの1行目の 2 は、状態密度の計算を指定しています。2行目はメッシュの数です。

$ anphon Si_phdos.in > Si_phdos.out

状態密度は Si333.dos に出力されます