Quantum ESPRESSOでNaClのDOSをプロットする
Last Update:2025/08/29
以下の手順でNaClの状態密度（DOS）を計算し、プロットします。
-
SCF計算（run_qe.pyを使用）
-
NSCF計算（非自己無撞着計算）
-
dos.x を用いたDOS計算
-
DOSのプロット
1. SCF計算
まずNaとClの擬ポテンシャルを入手します。
wget https://pseudopotentials.quantum-espresso.org/upf_files/Na.pbesol-spn-kjpaw_psl.1.0.0.UPF wget https://pseudopotentials.quantum-espresso.org/upf_files/Cl.pbesol-n-kjpaw_psl.1.0.0.UPF
下の内容のrun_qe.pyを
python3 run_qe.py
Potential energy: -2539.4114330941884 eV
のように実行します。run_qe.pyではaseとpymatgenを使いますので、pipによるインストールが必要です。
また、NaClの結晶構造の情報を得るためにはMaterials ProjectのAPIキーが必要になります。run_qe.pyの”USER API KEY”を個所を書き換えてください。
[run_qe.py]
import os # OS関連の機能を使うためのモジュール
import subprocess # サブプロセスを実行するためのモジュール
from pymatgen.core.structure import Structure # pymatgenの構造を扱うクラス
from pymatgen.ext.matproj import MPRester # Materials ProjectのAPIを使うためのクラス
from ase.io import read # ASEで構造ファイルを読み込むための関数
from ase.calculators.espresso import Espresso # ASEでQuantum ESPRESSOを使うためのクラス
from ase import Atoms # ASEで原子構造を扱うクラス
MY_API_KEY="USER API KEY" # Material ProjectのAPIキー
mpr=MPRester(MY_API_KEY) # MPResterオブジェクトを作成
material_id = "mp-22862" # Material Projectの物質ID (NaCl)
structure = mpr.get_structure_by_material_id(material_id) # 物質IDから構造を取得
# pymatgen StructureオブジェクトをCIFファイルに保存
structure.to(filename="mp-22862.cif") # pymatgenのStructureをCIF形式で保存
# ASEのread関数でCIFファイルを読み込む
atoms = read("mp-22862.cif") # ASEでCIFファイルを読み込み、Atomsオブジェクトを作成
pseudopotentials = {'Na': 'Na.pbesol-spn-kjpaw_psl.1.0.0.UPF', # Naの擬ポテンシャル
'Cl': 'Cl.pbesol-n-kjpaw_psl.1.0.0.UPF'} # Clの擬ポテンシャル
pseudodir = './' # 擬ポテンシャルファイルがあるディレクトリ
input_data = { # Quantum ESPRESSOの入力データを定義
'control': { # 計算制御に関する設定
'calculation': 'scf', # 自己無撞着計算 (SCF)
'prefix': 'nacl', # 計算のprefix (出力ファイル名などに使用)
'pseudo_dir': pseudodir, # 擬ポテンシャルファイルへのパス
'outdir': './', # 出力ディレクトリ
},
'system': { # 系に関する設定
'occupations': 'smearing', # 部分占有を考慮 (金属など)
'smearing': 'gaussian', # スミアリング関数にガウシアンを使用
'degauss': 0.01 # スミアリング幅 (Ry)
},
'electrons': { # 電子に関する設定
'conv_thr': 1.0e-8, # SCF計算の収束閾値
}
}
pw_command = "pw.x" # pw.x (Quantum ESPRESSOの実行ファイル) のパス
calc = Espresso( # Quantum ESPRESSO Calculatorの作成
pseudopotentials=pseudopotentials, # 擬ポテンシャル
kpts=(4, 4, 4), # k点メッシュ
input_data=input_data, # 入力データ
command='mpirun -np 4 ' + pw_command + ' -in espresso.pwi > espresso.pwo' # 実行コマンド
)
#from ase.calculators.espresso import EspressoProfile
#calc = Espresso(
#pseudopotentials=pseudopotentials,
#kpts=(4, 4, 4),
#input_data=input_data,
#profile=EspressoProfile(command='mpirun -np 4 ' + pw_command + ' -in espresso.pwi > espresso.pwo', pseudo_dir=pseudodir)
#)
atoms.set_calculator(calc) # AtomsオブジェクトにCalculatorを設定
energy = atoms.get_potential_energy() # ポテンシャルエネルギーを計算
print(f"Potential energy: {energy} eV") # ポテンシャルエネルギーをeV単位で出力
上の実行によって[espresso.pwi]というpw.x の入力ファイルが作成される。もしcalc = Espressoの箇所でエラー終了してしまったら、calc=Espressoの箇所をコメントアウトし、calc以下にあるコメントアウトを外して試す。
[espresso.pwi]
&CONTROL # 計算の制御に関する設定 calculation = 'scf' # 自己無撞着計算 (SCF) を行う outdir = './' # 出力ファイルを保存するディレクトリ prefix = 'nacl' # 計算の識別子 (ファイル名に使用) pseudo_dir = './' # 擬ポテンシャルファイルが置かれているディレクトリ / &SYSTEM # 計算系に関する設定 occupations = 'smearing' # フェルミ分布で電子占有率を計算する degauss = 0.01 # フェルミ分布の広がり (Ry単位) smearing = 'gaussian' # スミアリング関数としてガウシアン関数を使用 ntyp = 2 # 系の原子の種類数 nat = 2 # 単位セル内の原子の数 ibrav = 0 # ブラベー格子の種類 (0はCELL_PARAMETERSを使用) / &ELECTRONS # 電子に関する設定 conv_thr = 1e-08 # SCF計算の収束判定の閾値 (エネルギーの差) / &IONS # イオンに関する設定（構造緩和計算時に使用） / &CELL # セルに関する設定 (セル形状最適化計算時に使用) / ATOMIC_SPECIES # 原子種の情報 Na 22.98976928 Na.pbesol-spn-kjpaw_psl.1.0.0.UPF # ナトリウムの原子種、質量、擬ポテンシャルファイル名 Cl 35.45 Cl.pbesol-n-kjpaw_psl.1.0.0.UPF # 塩素の原子種、質量、擬ポテンシャルファイル名 K_POINTS automatic # k点の設定 4 4 4 0 0 0 # 4x4x4 のモノクロプケ点メッシュ、原点オフセットなし CELL_PARAMETERS angstrom # 単位胞の格子ベクトル (angstrom単位) 3.95140292000000 0.00000000000000 0.00000000000000 # aベクトル 1.97570124148061 3.42201456944483 0.00000000000000 # bベクトル 1.97570153872327 1.14067218809081 3.22630565117545 # cベクトル ATOMIC_POSITIONS angstrom # 原子座標 (angstrom単位) Na 0.0000000000 0.0000000000 0.0000000000 # ナトリウムの座標 Cl 3.9514028501 2.2813433788 1.6131528256 # 塩素の座標
（注）Quantum ESPRESSOの入力ファイルに全角の文字が入っているとエラーが生じる場合がある。
2. NSCF計算
以下の入力ファイル[nscf.in]を
mpirun -np 4 pw.x -in nscf.in > nscf.out
のように実行してNSCF計算を行います。
[nscf.in]
&control calculation = 'nscf' prefix = 'nacl' pseudo_dir = './' outdir = './' / &system ibrav = 0 nat = 2 ntyp = 2 occupations = 'tetrahedra' / &electrons conv_thr = 1.0d-8 / ATOMIC_SPECIES Na 22.98977 Na.pbesol-spn-kjpaw_psl.1.0.0.UPF Cl 35.45000 Cl.pbesol-n-kjpaw_psl.1.0.0.UPF ATOMIC_POSITIONS crystal Na 0.0 0.0 0.0 Cl 0.5 0.5 0.5 K_POINTS automatic 8 8 8 0 0 0 CELL_PARAMETERS angstrom 3.95140292000000 0.00000000000000 0.00000000000000 1.97570124148061 3.42201456944483 0.00000000000000 1.97570153872327 1.14067218809081 3.22630565117545
3. dos.x を用いたDOS計算
以下の入力ファイル[dos.in]を
mpirun -np 1 dos.x -in dos.in > dos.out
のように実行してdos.x を使用してDOS計算を行います。
無事に実行が終わるとnacl.dosとファイルが出来ます。
[dos.in]
&dos prefix = 'nacl' outdir = './' ngauss = 0 degauss = 0.01 emin = -10.0 ! エネルギー範囲の最小値 (eV) emax = 20.0 ! エネルギー範囲の最大値 (eV) deltae = 0.01 ! エネルギー間隔 (eV) /
4. DOSのプロット
下記の[plot_dos.py] という名前のPythonスクリプトを用います。
[plot_dos.py]
import numpy as np
import matplotlib.pyplot as plt
# データの読み込み
data = np.loadtxt('nacl.dos')
energy = data[:, 0]
dos = data[:, 1]
# プロット
plt.plot(energy, dos)
plt.xlabel('Energy (eV)')
plt.ylabel('DOS (states/eV)')
plt.title('Density of States of NaCl')
plt.grid(True)
plt.xlim([-10, 20]) # エネルギー範囲を調整
plt.savefig('nacl_dos.png') # 図をファイルに保存
plt.show()
[plot_dos.py] を実行します。
python3 plot_dos.py
下図のようなnacl_dos.pngというファイルが出来ます。
