====== 概要 ====== * 以下のものを計算できます * 最局在ワニエ関数を使ったトランスファー積分 * 誘電関数から求まる動的物理量（電子エネルギー欠損スペクトル、光学伝導度、反射率） * 制限乱雑位相近似 (constraint random-phase approximation, cRPA) を用いた電子間相互作用の大きさの見積もり * cRPAの理論はこちらのテキストの7章がわかりやすいです * https://www.cond-mat.de/events/correl11/ * 日本語および英語のマニュアルでは使用例とその結果が充実しているので、初心者でも特に困ることはないと思います * 最新版 (20200113) の使用方法は論文を参照してください（いろいろ変更点があるので、昔のマニュアルを使うと混乱が生じるかもしれません） * arXiv: https://arxiv.org/abs/2001.02351 * 出版されたもの: https://doi.org/10.1016/j.cpc.2020.107781 * 計算は非常に重いので、普通のパソコンですとユニットセルに原子が5, 6個くらいの物質が限界かと思われます * f電子系やスピン軌道相互作用がある系には非対応 * 擬ポテンシャルはノルム保存型を使用する必要があります。例えば以下のサイトから入手できます * http://www.pseudo-dojo.org/ ... optimized norm-conserving Vanderbilt PPs * https://www.quantum-espresso.org/pseudopotentials/hartwigesen-goedecker-hutter-pp ... Hartwigesen-Goedecker-Hutter PPs ====== インストール ====== ===== RESPACK 20200113 + Ubuntu 18.04 ===== 事前にMPIライブラリと以下のパッケージをインストールしてください。 $ sudo apt install cmake python-numpy - 公式サイトからソースコードをダウンロードします。 * https://sites.google.com/view/kazuma7k6r - 端末を立ち上げ、ファイルを展開し、その中に入ります $ tar xzvf RESPACK-20200113.tar.gz $ cd RESPACK-20200113 - 作業ディレクトリを作成し、その中に入ります $ mkdir build $ cd build - CMakeによりインストールの設定を行います。下の例では /usr/local/respack/ にインストールします * GCCを使う場合 $ cmake -DCONFIG=gcc -DCMAKE_INSTALL_PREFIX=/usr/local/respack ../ * Intel Parallel Studio XE を使う場合 $ cmake -DCONFIG=intel -DMPI_Fortran_COMPILER=mpiifort -DCMAKE_INSTALL_PREFIX=/usr/local/respack ../ - メイクとインストールをします $ make $ sudo make install - パスをとおします $ echo 'export PATH=$PATH:/usr/local/respack/bin/' >> ~/.bashrc $ source ~/.bashrc ======= 使い方 ======= RESPACKは主に次の流れで計算します - QuantumESPRESSO の pw.x によるSCF計算 - RESPACKの calc_wannier による最局在ワニエ関数の計算。事前にSCF計算が必要 - RESPACKの calc_chiqw による誘電関数の計算。cRPAを使う場合は事前にワニエ関数の計算が必要 - RESPACKの calc_w3d による直接相互作用の計算。事前にcalc_chiqwのcRPA計算が必要 - RESPACKの calc_j3d による交換相互作用の計算。事前にcalc_chiqwのcRPA計算が必要 また、データの整理に便利なPythonのスクリプトが RESPACK-20200113/utils/ 内にありますので、マニュアルとともに参照してください。 ===== インプットファイル ===== 基本的には下のファイルをQuantumESPRESSOでの作業ディレクトリにコピーして計算を行います。 !の後ろがコメントで、!nは必ず指定するもの、!oは任意に指定するもの（下の例では断りがない限り初期値が設定されているので消しても動作に影響がない）です。 ¶m_wannier N_wannier = 9, !n: 求めたいワニエ軌道の数 Lower_energy_window =-10.0d0, !n: エネルギーウィンドウ下限 (eV) Upper_energy_window = 35.0d0, !n: エネルギーウィンドウ上限 (eV) N_initial_guess = 9, !n: 初期ガウシアンの個数 icell = 0, !o: 結晶格子情報 EPS_SPILLAGE = 1.0d-4, !o: こぼれ汎関数最小化における収束閾値 (au) EPS_SPREAD = 1.0d-4, !o: 広がり汎関数最小化における収束閾値 (au) DAMP = 0.1d0, !o: こぼれ汎関数最小化における収束因子 MAX_STEP_LENGTH = 4.0d0, !o: 広がり汎関数最小化におけるステップ長 set_inner_window = F, !o: インナーウィンドウの設定有無 (有:T , 無:F) Lower_inner_window = 10.0d0, !o: エネルギーインナーウィンドウ下限 (eV) Lower_energy_windowよりも大きな値にする Upper_inner_window = 30.0d0, !o: エネルギーインナーウィンドウ上限 (eV) Upper_energy_windowよりも小さな値にする flg_initial_guess_direc = 0, !o: 各原子で局所座標系を指定する場合は1 flg_BMAT = 0, !o: 初期ガウシアンの線形結合係数を読み込む場合は1 electron_number_wannier_space = 0.0 !o: ワニエ空間での電子数（最大2*N_wannier個）。ただし0.0とした場合はSCF計算のフェルミエネルギーが使われる flg_fermisurface = 0, !o: フェルミ面を計算 (しない: 0, する: 1) flg_global_dos = 0/ !o: バンド計算の状態密度を計算 (しない: 0, する: 1) s 0.2d0 0.0d0 0.0d0 0.0d0 !n: 軌道型 , 軌道指数 , ガウシアン位置 (結晶座標で指定) s px 0.2d0 0.0d0 0.0d0 0.0d0 !n: px py 0.2d0 0.0d0 0.0d0 0.0d0 !n: py pz 0.2d0 0.0d0 0.0d0 0.0d0 !n: pz dz2 0.2d0 0.0d0 0.0d0 0.0d0 !n: dz2 dzx 0.2d0 0.0d0 0.0d0 0.0d0 !n: dzx dyz 0.2d0 0.0d0 0.0d0 0.0d0 !n: dyz dx2 0.2d0 0.0d0 0.0d0 0.0d0 !n: dx2-y2 dxy 0.2d0 0.0d0 0.0d0 0.0d0 !n: dx ¶m_interpolation N_sym_points = 5, !n: 計算ラインを構成する対称k点数 Ndiv = 40/ !n: 対称k点間の分割数 0.500 0.500 0.500 !n: 対称k点（逆格子ベクトルの分率座標） 0.000 0.000 0.000 !n: 対称k点 0.500 0.000 0.500 !n: 対称k点 0.500 0.250 0.750 !n: 対称k点 0.500 0.500 0.500 !n: 対称k点 ¶m_visualization flg_vis_wannier = 0, !o: 実空間ワニエ関数を計算 (しない: 0, する: 1) N_write_wannier = 9, !o: 求めたい実空間ワニエ関数の数。指定しない場合はN_wannierの値になる ix_vis_min=-1, !o: ワニエ関数の描画範囲 (a1方向始点格子) ix_vis_max= 1, !o: ワニエ関数の描画範囲 (a1方向終点格子) iy_vis_min=-1, !o: ワニエ関数の描画範囲 (a2方向始点格子) iy_vis_max= 1, !o: ワニエ関数の描画範囲 (a2方向終点格子) iz_vis_min=-1, !o: ワニエ関数の描画範囲 (a3方向始点格子) iz_vis_max= 1/ !o: ワニエ関数の描画範囲 (a3方向終点格子) 1 2 3 4 5 6 7 8 9 !o: N_write_wannierを指定したときの、計算したいワニエ軌道の番号 ¶m_chiqw Ecut_for_eps = 3.0d0, !o: 分極関数のカットオフ (Ry) 記述がなければSCF計算のecutwfcの値の1/10になる Num_freq_grid = 70, !o: 計算周波数の総数（メッシュは log ω に対して切られる） N_CALC_BAND = 50, !o: 分極計算で考慮されるバンドの総数。SCF計算のnbndの数よりも小さくする。記述がなければSCF計算のnbndの値になる MPI_num_proc_per_qcomm = 2, !o: コミュニティ当たりのプロセス数。MPIを使用する際に指定 MPI_num_qcomm = 1, !o: コミュニティの数。MPIを使用する際に指定 MPI_io_rank = 0, !o: 標準出力させる MPI の番号 shift_ef = 0.0d0, !o: 人為的なフェルミ準位のシフト値 (eV) Max_excitation_energy = 200.0d0, !o: 計算で考慮する最大励起エネルギー (eV) Green_func_delt = 0.1d0, !o: テトラへドロン計算でのスメアリング値 (eV) ttrhdrn_dmna = 1.0d-3, !o: テトラへドロン計算での縮退判定パラメータ 1 (eV) ttrhdrn_dmnr = 1.0d-3, !o: テトラへドロン計算での縮退判定パラメータ 2 (eV) flg_cRPA = 0, !o: 制限 RPA オプション (0: 通常RPA, 1: 制限RPA) flg_calc_type = 0/ !o: 0: すべてのq, 1: q=0のみ, 2: 任意のq（マニュアル14章） ¶m_calc_int calc_ifreq = 1, !o: 計算出力させたい周波数を1からNum_freq_gridで指定（デフォルトの1は omega=0） ix_intJ_min = 0, !o: 交換積分の計算範囲 (a1 方向始点格子 )(calc_j3d のみ ) ix_intJ_max = 0, !o: 交換積分の計算範囲 (a1 方向終点格子 )(calc_j3d のみ ) iy_intJ_min = 0, !o: 交換積分の計算範囲 (a2 方向始点格子 )(calc_j3d のみ ) iy_intJ_max = 0, !o: 交換積分の計算範囲 (a2 方向終点格子 )(calc_j3d のみ ) iz_intJ_min = 0, !o: 交換積分の計算範囲 (a3 方向始点格子 )(calc_j3d のみ ) iz_intJ_max = 0/ !o: 交換積分の計算範囲 (a3 方向終点格子 )(calc_j3d のみ ) 以下で詳細な計算手順を説明します ===== QuantumESPRESSOによるSCF計算 ===== RESPACKによる計算の前に、QEによる計算を済ませます。 たとえばPWscfの入力ファイルが prefix.scf.in だった場合、次のようにします。 $ mpirun -n 2 pw.x < prefix.scf.in > prefix.scf.out * 既約k点の情報が必要なので、SCF計算 (ネームリスト&control内でcalculation='scf'とする計算) の直後にRESPACKを使った計算をする必要があります * k点数を変更したり、バンド計算など別の計算を行った後にRESPACKを使用する場合は、ディレクトリ prefix.save を削除してから再度SCF計算を実行してください * k点は automatic により生成し、シフトしていないもの (0 0 0) を使用してください outdir='./tmp/'として、RESPACKの作業ディレクトリ dir-wfn/ 内に計算に必要なファイルを生成します $ qe2respack.py tmp/prefix.save/ 以下ではRESPACKを使った計算に入っていきます ===== ワニエ関数の計算 ===== ¶m_wannier, ¶m_interpolation, ¶m_visualization 内のパラメータが読み込まれます $ export OMP_NUM_THREADS=12 $ export MKL_NUM_THREADS=12 $ calc_wannier < respack.in > LOG.wannier * OpenMPによる並列化のみ対応しています * エネルギーウインドウはワニエ基底に射影するバンドの範囲です。射影しない軌道のバンドが少し混じっていても構いません。 * 特にエネルギーインナーウィンドウの範囲を指定すると、その中にあるすべてのバンドがワニエ基底によるバンドと一致するように計算されます。 ===== 分極関数の計算 ===== ¶m_chiqw 内のパラメータが読み込まれます。 制限RPA法 (flg_crpa=1) によって計算する場合は、事前に calc_wannier による計算を済ませておきます。 $ export OMP_NUM_THREADS=12 $ export MKL_NUM_THREADS=12 $ mpirun -n 2 calc_chiqw < respack.in > LOG.chiqw * OpenMPとMPIのハイブリッド計算を行います。OMP_NUM_THREADSおよびMKL_NUM_THREADSにMPIの1スレッドあたりのOpenMP並列数を指定してください * MPI_num_qcommには、すべてのqで計算するときに、いくつの点を同時に並列計算を行うかを指定します * MPI_num_proc_per_qcommには、qの各点の計算に使用するスレッド数を指定します * MPIのプロセス数は MPI_num_proc_per_qcomm * MPI_num_qcomm に一致させます * 励起状態の計算になるので、多くの非占有バンドを取り入れる必要があります * バンドの数、エネルギーカットオフ、k点数、ユニットセルの体積に比例して非常に多くのメモリを消費します * 観測される物理量としての光学応答を調べたければ通常のRPAを、相互作用パラメータを求めたければ制限RPAを使用します。それぞれ収束に必要なパラメータは異なるので注意してください。 * 一般に、通常のRPAの方が計算コストが大きいです ===== 相互作用の計算 ===== ¶m_calc_int 内のパラメータが読み込まれます。 事前に制限RPA法（calc_chiqw で flg_cRPA=1）を使ったすべてのq (flg_calc_type = 0) での計算が必要です。 * 直接相互作用 $ export OMP_NUM_THREADS=12 $ export MKL_NUM_THREADS=12 $ calc_w3d < respack.in > LOG.w3d * 交換相互作用 $ export OMP_NUM_THREADS=12 $ export MKL_NUM_THREADS=12 $ calc_j3d < respack.in > LOG.j3d ====== 計算の収束 ======= 分極関数の計算で収束させなければならないパラメータは次の3つです * k点数（SCF計算） * N_CALC_BAND ... バンドの数 * Ecut_for_eps ... 分極関数のカットオフ (Ry) 有効相互作用を求める場合は、さらに Num_freq_grid（周波数のグリッド数）を増やして有効相互作用の値の収束を確認します。 * 計算が重いので、有効相互作用の有効数値は2桁出すのが精一杯かもしれません ====== プロットファイル ====== ===== 分極関数 ===== ############### # omega range # ############### omega_min=0 omega_max=200 set logscale x ################### # common settings # ################### set terminal postscript eps enhanced color "Helvetica,24" set key right top set xrange [omega_min:omega_max] set style line 1 lw 3 lc "red" set style line 2 lw 3 lc "blue" set style line 3 lw 3 lc "dark-green" set xlabel '{/Symbol w} [eV]' ######## # EELS # ######## set ylabel "-Im {/Symbol e}(0,{/Symbol w})^{-1}" set output "eels.eps" plot 'dat.eels-x' using 1:2 title "x" with lines ls 1,\ 'dat.eels-y' using 1:2 title "y" with lines ls 2,\ 'dat.eels-z' using 1:2 title "z" with lines ls 3 ################################### # Macroscopic dielectric function # ################################### set ylabel "Macroscopic dielectric function" set output "macroscopic_epsilon-x.eps" plot 'dat.macroscopic_epsilon-x' using 1:2 title "Re {/Symbol e}^{-1}" with lines ls 1,\ 'dat.macroscopic_epsilon-x' using 1:3 title "Im {/Symbol e}^{-1}" with lines ls 2 set output "macroscopic_epsilon-y.eps" plot 'dat.macroscopic_epsilon-y' using 1:2 title "Re {/Symbol e}^{-1}" with lines ls 1,\ 'dat.macroscopic_epsilon-y' using 1:3 title "Im {/Symbol e}^{-1}" with lines ls 2 set output "macroscopic_epsilon-z.eps" plot 'dat.macroscopic_epsilon-z' using 1:2 title "Re {/Symbol e}^{-1}" with lines ls 1,\ 'dat.macroscopic_epsilon-z' using 1:3 title "Im {/Symbol e}^{-1}" with lines ls 2 ######################## # Optical conductivity # ######################## set ylabel "Optical conductivity [10^6/{/Symbol W}/m]" set output "optcond-x.eps" plot 'dat.optical_conductivity-x' using 1:2 title "Re {/Symbol s}({/Symbol w})" with lines ls 1,\ 'dat.optical_conductivity-x' using 1:3 title "Im {/Symbol s}({/Symbol w})" with lines ls 2 set output "optcond-y.eps" plot 'dat.optical_conductivity-y' using 1:2 title "Re {/Symbol s}({/Symbol w})" with lines ls 1,\ 'dat.optical_conductivity-y' using 1:3 title "Im {/Symbol s}({/Symbol w})" with lines ls 2 set output "optcond-z.eps" plot 'dat.optical_conductivity-z' using 1:2 title "Re {/Symbol s}({/Symbol w})" with lines ls 1,\ 'dat.optical_conductivity-z' using 1:3 title "Im {/Symbol s}({/Symbol w})" with lines ls 2 ################ # Reflectivity # ################ set ylabel "Reflectivity" set output "reflectivity.eps" plot 'dat.reflectivity-x' using 1:2 title "x" with lines ls 1,\ 'dat.reflectivity-y' using 1:2 title "y" with lines ls 2,\ 'dat.reflectivity-z' using 1:2 title "z" with lines ls 3 exit dir_eps内に設置して、$ gnuplot eps.pltとしてください