はじめに

MateriApps LIVE!にプレインストールされたQuantum Espressoで誘電率を計算してみます。環境はVirtualbox上で動かしているMateriApps LIVE! 3.3です。VirtualboxにはCPUを4つ割り当ててみました。

誘電率やフォノン計算については、２０１４年のQuantum Espressoに関するSummer Schoolのページ：
http://www.iiserpune.ac.in/~smr2626/talks-presentations.html
のHands on: Phononsが参考になりそうです。

計算の手順は２ステップあるようです。まず、Quantum Espressoパッケージに含まれているpw.xというプログラムを使い、自己無撞着計算によって基底状態の電子密度やエネルギー、構造を求めます。その結果を使って、ph.xというプログラムで、密度汎関数摂動論（DFPT）に基づいてフォノンモードやボルン有効電荷、誘電率を求めます。

pw.xによる格子定数および内部座標の最適化と電子状態計算

ターミナルを開いて(スタートメニューから”System Tools” > “LXTerminal”)作業ディレクトリを作ります。

$ mkdir qe-GaAs 
$ cd qe-GaAs

次に、計算する元素に対応する擬ポテンシャルの定義ファイルを用意します。Quantum Espressoのホームページでほとんどの元素に対応したファイルが公開されています（http://www.quantum-espresso.org/pseudopotentials）。周期表でGaをクリックすると、様々な方法で作成された擬ポテンシャルファイルの一覧が現れます。今回は一番上の、Ga.pbe-dn-kjpaw_psl.1.0.0.UPFをとりあえず使ってみます。同様に、Asも、一覧の一番上の、As.pbe-n-kjpaw_psl.1.0.0.UPFを使います。ターミナルからは、wgetコマンドでURLを指定してダウンロードすることができます。

$ wget https://pseudopotentials.quantum-espresso.org/upf_files/Ga.pbe-dn-kjpaw_psl.1.0.0.UPF
$ wget https://pseudopotentials.quantum-espresso.org/upf_files/As.pbe-n-kjpaw_psl.1.0.0.UPF

次に、pw.xの入力ファイルを用意します。入力内容についてはreference manualに網羅されており（こちら）、別のレビュー記事でも簡単な解説がなされています。calculation=’vc-relax’を指定して、格子定数と内部座標の同時最適化を行います。

&control
calculation = 'vc-relax' 
restart_mode='from_scratch'
prefix = 'GaAs'
pseudo_dir = './'
wf_collect = .true.
outdir = 'tmp'
/
&system
ibrav = 2
celldm(1) = 10.869
nat = 2
ntyp = 2
ecutwfc = 100
ecutrho = 400
occupations='fixed'
/
&electrons
mixing_mode = 'plain'
mixing_beta = 0.5
conv_thr = 1.0d-10
/
&ions
/
&cell
press_conv_thr = 0.1
/
ATOMIC_SPECIES
Ga 69.723 Ga.pbe-dn-kjpaw_psl.1.0.0.UPF
As 74.921595 As.pbe-n-kjpaw_psl.1.0.0.UPF
ATOMIC_POSITIONS
Ga 0.00 0.00 0.00
As 0.25 0.25 0.25
K_POINTS {automatic}
6 6 6 1 1 1

Linux操作に慣れていましたらemacsなどで入力ファイルを作成します。GUIを好む場合は、

$ mousepad GaAs.scf.in

というコマンドで、GUIエディターのMousepadが立ち上がります。上記の入力ファイル内容をコピー＆ペーストし、File=>Saveで保存できます。これを実行してみます。今回は4CPUをVirtualboxに割り当てているので、mpirunを使って並列計算させてみます。

$ mpirun -np 4 -x OMP_NUM_THREADS=1 pw.x < GaAs.scf.in | tee GaAs.scf.out
....(中略)....
=------------------------------------------------------------------------------=
JOB DONE.
=------------------------------------------------------------------------------=

５分ほどかかりました。1 CPUしか割り当てなかった場合は、

$ pw.x < GaAs.scf.in | tee GaAs.scf.out

で実行できますが、数十分かかってしまいます。平面波基底のカットオフ（ecutwfc）を小さくするなどで計算時間を節約できますが、計算精度とのトレードオフになります。

ph.xを使った密度汎関数摂動論（ＤＦＰＴ）法による誘電率の計算

さて、計算が完了すると、tmpというディレクトリ（GaAs.scf.inの中のoutdirで指定している）ができており、その中に電子密度などの計算結果が出力されています。これをph.xで読み込んでDFPT計算を行います。ph.x用に入力ファイルを別途用意します：

&inputph
outdir='./tmp/',
prefix='GaAs',
fildyn='GaAs.dynmat_G',
tr2_ph=1.0d-14,
amass(1)=69.723,
amass(2)=74.921595,
epsil=.true.
/
0.0 0.0 0.0

ファイル名はGaAs.ph.inにしましょう。mousepadを使う場合は

$ mousepad GaAs.ph.in

でエディターを立ち上げて、上記の内容をコピー＆ペーストし、保存します。それではph.xを実行します。

$ mpirun -np 4 -x OMP_NUM_THREADS=1 ph.x < GaAs.ph.in | tee GaAs.ph.out
....(中略)....
=------------------------------------------------------------------------------=
JOB DONE.
=------------------------------------------------------------------------------=

7分ほどかかりました。並列計算を行わない場合はph.xより前のコマンドは不要ですが、計算時間が数倍かかります。

出力ファイルを調べてみると、誘電率テンソルを見つけることができます。立方晶系なので、対角項のみが有限の値を持ち、15.2という結果が得られました。実験値は12~13程度とされており、良い予測ができています。

Dielectric constant in cartesian axis

( 15.222785223 -0.000000000 -0.000000000 )
( -0.000000000 15.222785223 0.000000000 )
( -0.000000000 0.000000000 15.222785223 )

その他にもボルン有効電荷や

Effective charges (d Force / dE) in cartesian axis

atom 1 Ga 
Ex ( 2.22829 -0.00000 -0.00000 )
Ey ( 0.00000 2.22829 -0.00000 )
Ez ( -0.00000 -0.00000 2.22829 )
atom 2 As 
Ex ( -2.25678 0.00000 0.00000 )
Ey ( -0.00000 -2.25678 0.00000 )
Ez ( 0.00000 0.00000 -2.25678 )

Γ点フォノン

q = ( 0.000000000 0.000000000 0.000000000 )

**************************************************************************
freq ( 1) = -0.305220 [THz] = -10.181035 [cm-1]
freq ( 2) = -0.305220 [THz] = -10.181035 [cm-1]
freq ( 3) = -0.305220 [THz] = -10.181035 [cm-1]
freq ( 4) = 7.580768 [THz] = 252.867194 [cm-1]
freq ( 5) = 7.580768 [THz] = 252.867194 [cm-1]
freq ( 6) = 7.580768 [THz] = 252.867194 [cm-1]
**************************************************************************

が出力されています。freq(1),(2),(3)について、負の値になってしまっています。音響フォノンの振動数はΓ点では0になるのですが、数値精度の問題で有限の値になってしまいます。なお、Γ点以外のフォノンを計算する場合は、こちらのレビューが参考になります。

まとめ

本稿では、MateriApps LIVE!上でQuantum Espressoを使った誘電率の第一原理計算をレビューしました。筆者のパソコンでは、作業時間を含めても1時間程度で実行可能でした。単結晶材料であれば、大規模計算機を使わなくても十分に高精度な計算が可能であることが分かります。