Quantum Espressoは第一原理計算を用いて物質のバンド構造などの電子状態を計算してくれるパッケージです。

インストール方法

科学計算用のスパコンなどでは予めコンパイルしたものが入っていると思いますが、自分でやるかり方を紹介します。Linux環境(あるいは仮想環境)を想定しています。

１つ目は直接的なやり方です。Quantum Espressoのホームページからダウンロードしてインストールします：

それを解凍して make allでコンパイルします。ただし、BlasやLApackなどの線形代数計算のパッケージが必要で、また環境に応じて色々とmake.incのパラメーターを設定しなければなりません。

2つ目はMateriAppsLive!を使うやり方です。これは東大の物性研究所が提供している環境構築ツールで、VirtualBoxとをインストールすればそこの仮想環境でQuantum Espressoを動かすことができます。こちらの方がコンパイルのパラメーター設定の手間が省けて楽だと思います。
VirtualBoxはここからダウンロードできます：

また、MateriAppsLive!はここからダウンロードできます：

両方インストールしたら同じディレクトリに置き、MateriAppsLive-(バージョン).ovaをダブルクリックします。
すると、VirtualBoxが起動して、インポートの画面になるので「インポート」をクリックします。
数分で読み込みが終わってマネージャーが起動します。
マネージャー画面でMateriAppsLive...を選択してバーの起動ボタンをクリックすると、ログイン画面が出てくるのでユーザー名:userとパスワード:liveを入力します。

これらの手順を行うと、コンパイル済みのQuantum Espressoが入ったLinux仮想環境が使えるようになります。
Quantum Espressoは/usr/share/にあります。
また、操作は基本的にスタートメニュー -> System Tools -> LXTerminalでターミナルを開いてコマンド入力で行います。

ホストOSと仮想環境のファイル共有

仮想環境は基本的にホストOSとは切り離されていますが、ホストOSのフォルダを共有することができます。

仮想環境を閉じます。
VirtualBoxマネージャーでMateriAppsLiveを選択してバーの設定ボタンをクリック
共有フォルダーのタブを開いて右側の「+」マークをクリック
「フォルダーのパス」の右にある「ｖ」マークをクリックして「その他」を選択
共有するフォルダを選択して「自動マウント」にチェックを入れてＯＫボタンをクリック
仮想環境で/media/sf_***/に指定したフォルダがつくられ、ホストＯＳと共有されます。

計算の下準備

結晶の電子構造を計算するために、入力する結晶構造の情報と、結晶を作る原子の情報の2つが必要となります。

結晶の情報

を参考にしてください。

原子の情報

次に結晶を構成する原子のデータです。これは擬ポテンシャルと言って、単独の原子のイオンのポテンシャル(つまりクーロンポテンシャル)を結晶構造を作った場合に合わせて
変形させたものになります。Quantum Esspressoのホームページからダウンロードできます：

リンク先の周期表から必要な元素をクリックすると、擬ポテンシャルの一覧を見ることができます。必要なUPFファイルをダウンロードしてください。
ただし、様々な種類があるので適切なものを選ぶ必要があります。
擬ポテンシャル選択の着眼点として例えば以下の3つがあります。

Pseudopotential type：PAWの方がUSPPよりも良いです
Functional type:LDA、PBEがあると思いますが、LDAは最低次の近似で、PBEは高次の補正も入っています。
Full-relativistic/Scalar-relativistic：Full-relativisticはスピン軌道相互作用が入っていて、Scalar-relativisticは入っていません。ただしFull-relativisticは複数の元素がある場合、1種類にしか適用できず、他の元素の擬ポテンシャルはScalar-relativisticにしなければいけないので注意が必要です。

他にもUPFファイルの中身を見てみるとどの軌道の電子が自由電子として設定されているかや、推奨されるカットオフエネルギーなどが書かれているのでそうした情報もチェックしておくとよいかもしれません。

計算例

以下、$\mathrm{Bi_2Se_3}$を例に計算してみます。詳しい入力ファイルの作り方は

を参照してください。

ファイル構造

/usr/share/にあったQuantum Espressoのフォルダをホームディレクトリにコピーして、q-e-qe-6.5というフォルダ名にしました(僕の使っていたパッケージのバージョンは6.5でした)。
そして~/q-e-qe-6.5/workと~/q-e-qe-6.5/pseudoの２つのフォルダを作ってpseudoフォルダには擬ポテンシャルのファイルを入れ、基本的にはworkディレクトリで作業をしました。

構造データの取得

の例で得た入力データを使います。

Bi2Se3.in

#************************************************************************************
#*                  Generated by cif2cell 1.2.10 2022-05-19 13:30                   *
#*  T. Bjorkman, Comp. Phys. Commun. 182, 1183-1186 (2011). Please cite generously. *
#*                                                                                  *
#*                                         ()                                       *
#*             Failed to get author information, No journal information             *
#************************************************************************************

&SYSTEM
  ibrav = 0
  A =    4.13400
  nat = 5
  ntyp = 2
/
CELL_PARAMETERS {alat}
  0.577350269189626   0.000000000000000   2.306079664570230 
 -0.288675134594813   0.500000000000000   2.306079664570230 
 -0.288675134594813  -0.500000000000000   2.306079664570230 
ATOMIC_SPECIES
  Bi  208.98000  Bi_PSEUDO
  Se   78.96000  Se_PSEUDO
ATOMIC_POSITIONS {crystal}
Bi   0.400460000000000   0.400460000000000   0.400460000000000 
Bi   0.599540000000000   0.599540000000000   0.599540000000000 
Se   0.209700000000000   0.209700000000000   0.209700000000000 
Se   0.790300000000000   0.790300000000000   0.790300000000000 
Se   0.000000000000000   0.000000000000000   0.000000000000000

今回はBiにBi.rel-pbe-dn-kjpaw_psl.1.0.0.UPF、SeにSe.pbe-dn-kjpaw_psl.1.0.0.UPFを使います。Biにはスピン軌道相互作用入りのものを使いました。このBi2Se3.inを元に各ステップの入力ファイルを作成します。

また、Quantum espressoは基本的に単位をÅからボーア(a.u.)に変える必要があります:

4.13400 Å /0.529177 \simeq 7.81213 \quad a.u.

しかし、モードがibrav=0のときは$\mathrm{A}=\cdots$とオングストローム単位で書いても大丈夫です。

構造最適化

こちらの記事をご覧ください。必ずしもこの過程を踏む必要はありませんが、圧力をかけた結晶構造を扱う場合などには以降の計算の前に必要な処理になります。

電荷密度分布計算

第一原理計算はシステムの電荷密度分布を元に様々なデータを計算するので、まずは大元となる電荷密度分布を計算します。

入力ファイルの例

入力ファイルは以下の通りです。

Bi2Se3.scf.in

&CONTROL
                 calculation = 'scf' ,
                  wf_collect = .true. ,
                      outdir = './output' ,
                  pseudo_dir = '~/q-e-qe-6.5/pseudo/' ,
                      prefix = 'Bi2Se2' ,
                     disk_io = 'low' ,
                   verbosity = 'high' ,
                     tstress = .true. ,
                     tprnfor = .true. ,
  max_seconds = 86000,
  restart_mode = from_scratch,
 /
 &SYSTEM
  ibrav = 0,
  celldm(1) = 7.81213,
  nat = 5,
  ntyp = 2,
                       ecutwfc = 150 ,
                     ecutrho = 1500 ,
                 occupations = 'smearing' ,
                     degauss = 0.00001 ,
		     noncolin = .true. ,
		     lspinorb = .true. ,
  /
 &ELECTRONS
            electron_maxstep = 10000000,
                    conv_thr = 0.0000000001 ,
                 mixing_beta = 0.5 ,
/
CELL_PARAMETERS {alat}
  0.577350269189626   0.000000000000000   2.306079664570230 
 -0.288675134594813   0.500000000000000   2.306079664570230 
 -0.288675134594813  -0.500000000000000   2.306079664570230 
ATOMIC_SPECIES
  Bi  208.98000  Bi.rel-pbe-dn-kjpaw_psl.1.0.0.UPF
  Se   78.96000  Se.pbe-dn-kjpaw_psl.1.0.0.UPF

ATOMIC_POSITIONS {crystal}
Bi   0.400460000000000   0.400460000000000   0.400460000000000 
Bi   0.599540000000000   0.599540000000000   0.599540000000000 
Se   0.209700000000000   0.209700000000000   0.209700000000000 
Se   0.790300000000000   0.790300000000000   0.790300000000000 
Se   0.000000000000000   0.000000000000000   0.000000000000000 

K_POINTS automatic 
  4 4 4   0 0 0

追加した部分でしばしばいじる必要がある部分を解説します。

&controlのパラメーター

outdir：計算結果を保存するディレクトリです。ない場合は自動生成します。
prefix：出力ファイル名のキーワードです。何でもよいのですが僕はいつも物質名にしています。
pseudo_dir：擬ポテンシャルのデータファイルが置いてあるディレクトリです。
max_seconds：計算の時間制限です(単位は秒)。この時間を経過すると計算が途中で打ち切られます。
restart_mode：from_scratchとrestartの2種類があります。restartはmax_secondで打ち切られた場合、その続きから計算するモードです。from_scratchは前回までの計算に関係なく0から計算するモードです。

&SYSTEMのパラメーター

noncolinとlspinorb：スピン軌道相互作用を考えるときにtrueにしておかなければいけないパラメーターです。スピン軌道相互作用がないときは２つともfalseにします。

&ELECTRONのパラメーター

conv_thr：密度関数の最適化の収束についての閾値です。
mixing_beta：割線法の更新パラメーターで、更新前のデータの重みです。スピン軌道相互作用があるときは非凸性が大きく山に引っかかりやすいので、そういうときは更新前のデータをあまり含めないように0.2~0.3くらいにするとよいです。割線法については以下のページを参考にしてください：

K_POINTS：空間のメッシュです。automaticモードでは対称性を元により細かいメッシュが切られるので、密度計算のときは$4\times 4 \times 4$くらいあれば十分です。

次のコマンドで実行できます。標準出力結果をBi2Se3.scf.outに書き出しています。

../bin/pw.x <Bi2Se3.scf.in >Bi2Se3.scf.out

ただ、簡単な物質の計算なら手元のPCで数分くらい回せば終わると思いますが、一般には並列化しないと計算時間が膨大になると予想されるので

mpirun -np 60 ../bin/pw.x <Bi2Se3.scf.in >Bi2Se3.scf.out

のようにMPIなどを使った方がよいと思います。ものにもよりますが並列数は少なくとも数十ノードくらいはあった方がいいと思います。

ワイコフ位置を使った書き方

ワイコフ位置は対称操作で移りあう他の内部座標を省略したときの最小限必要な原子の内部座標です。各構造のワイコフ位置は

などで調べることができます。

Qunatum espressoでは結晶構造の入力の際にワイコフ位置の原子から残りの原子を自動補間してくれるモードもあります。ただ僕はエラーでうまく動かせたことがありません。

出力

出力の最後に次のように終了時刻が表示されたら無事計算が完了しています。

Bi2Se3.scf.out

:
:
:

   This run was terminated on:  15: 2:23  19May2022            

=------------------------------------------------------------------------------=
   JOB DONE.
=------------------------------------------------------------------------------=

他にも出力にはFermi準位の情報なども出力されます。

Bi2Se3.scf.out

:
:
the Fermi energy is     8.7764 ev
:
:

バンド計算

次の記事を参考にしてください。

状態密度計算

以下の記事を参考にしてください。

フェルミ面の描画

こちらの記事をご覧ください。

ジョブシーケンス

僕はこれらの計算の実行をいちいち入力するのが面倒なので、次のようなBashファイルを実行して、まとめてこれらの計算を行うようにしています。ちなみにこれは60並列でMPI計算をしています。MPIの実行コマンドの表記は使っているマシンやパッケージによって異なると思うので、適当なものに書き換えて使ってください。

sub.sh

### input parameters
dir_exe='../..'
prefix='Bi2Se3'

##########################################################################################################

Exe_mpi(){
    #var1: executable file
    #var2: input file
    #var3: output file
    mpirun -np 60 $1 <$2 >$3
}

Sub_qe_seq(){
    #################################
    #var1: prefix
    #var2: executable file directry
    #################################
    echo "scf calculation start"
    Exe_mpi "$2/bin/pw.x"      "${1}.scf.in"       "${1}.scf.out"
    echo "pdos calculation start"
    Exe_mpi "$2/bin/pw.x"      "${1}.nscf.pr.in"   "${1}.nscf.pr.out"
    Exe_mpi "$2/bin/projwfc.x" "${1}.pr.in"        "${1}.pr.out"
    echo "band calculation start"
    Exe_mpi "$2/bin/pw.x"      "${1}.nscf.band.in" "${1}.nscf.band.out"
    Exe_mpi "$2/bin/bands.x"   "${1}.band.in"      "${1}.band.out"
}

Sub_qe_seq    $prefix $dir_exe

Wannier90を使ったタイトバインディングモデルの作成

この物質はBiとSeのp軌道だけを考えれば性質が説明できるのでこれらの軌道だけに限定して、より簡単なタイトバインディング模型を作れば、物質の重要な構造をはっきりさせたり、さらにいろいろな計算が手軽に行ったりできるようになります。そこで、Quantum EspressoにWannier90というパッケージを組み合わせればこのような模型を作ることができます。
詳しくはこちらを参考にしてください。

注意・補足

Quantum Espressoの入力ファイルはタブ空白が入っているとエラーを起こすので気を付けてください。
{prefix}.EXITという空のファイルをoutdirに作ると計算を途中で安全に打ち切ることができます。今の計算例だとBi2Se3.EXITとなります。こうやって打ち切った計算はrestart_mode=restartとして再度実行すれば再開できます。また、作成した.EXITファイルは計算中断の後に自動的に削除されます。ただし、作成してから計算が止まるまでに多少時間がかかります。長いときは１～２時間くらい待たないといけません。
入力ファイルのパラメーター一覧は次のページを参考にしてください

参考ページ、文献