cifファイルからQuantum ESPRESSOのインプットファイルを自動生成するよ

Quantum ESPRESSO

結晶構造からエネルギー計算、バンド構造計算、フォノン分散関係計算などが行える量子化学計算アプリケーションが、Quantum ESPRESSOです。Quantum ESPRESSOで計算を行うにあたっては、インプットファイルを作成しなくてはなりません。このインプットファイルに計算条件や結晶構造を書き込むわけですが、有機結晶だと結晶の単位格子中に原子が100個以上含まれていることもザラで、手入力でインプットファイルを作るのは人間業ではありません。そこで、PythonでQuantum ESPRESSOのエネルギー一点計算を行うインプットファイルを自動生成するプログラムを作ってしまおうというのが今回のお話です。Quantum ESPRESSOのインプットファイルの作成に必要なのはcifファイルのみです。

今回はpythonで作成する方法を紹介していますが、Winmostarというアプリを使っても同様のことが可能です。というかその方が一般的な気がします。Quantum ESPRESSOのインプットファイルの作成において、このページを主に参考にしています。

インプットとアウトプット

-インプット
cifファイル

-アウトプット
filename.scf.in

filenameで指定したファイル名のQuantum ESPRESSOのエネルギー一点計算(scf計算)を行うインプットファイルが生成されます。Quantum ESPRESSOが動作する環境において、Terminalなどでfilename.scf.inのあるディレクトリに移動し、以下のコマンドにより計算を実行すると計算結果をfilename.scf.outで得ることができます。

pw.x < filename.scf.in > filename.scf.out

Quantum ESPRESSOの動作環境を整えるのは大変ですが、ソースコードをコンパイル方法¹や、MateriAppsLive!の仮想OSを構築する方法²などがあります。

準備

1. cifファイルを持ってくる
   何はともあれ計算を実行したい構造のcifファイルを手に入れてください。

2. pythoのpymatgenをインストールしておく
   pyhtonへのcifの読み込みはpymatgenを使って行なっています。pymatgenを使えるようにしておいてください³。

3. USPファイルをダウンロードしておく
   USPファイルとは原子の擬ポテンシャルが記述されたファイルのことだそうです。Quantum ESPRESSOはこのポテンシャルを読み込んで計算を行います。USPファイルは原子ごとに指定する必要があるので、たとえばcifの結晶が水素、炭素、酸素でできた有機物の結晶であれば水素、炭素、酸素の三つのUSPファイルが必要です。USPファイルはQuantum ESPRESSOの公式ホームページで公開されています⁴。詳しい選択方法等についてはここを参照してください。

実行

from pymatgen.core.structure import Structure
import numpy as np

## input 

name = "(filename)"
import_file = "(cif file directory)"
usp_files_dir = "(UPS files directory)"
usp_files = "(UPF file names)"

## make in file
# import cif file

file_name = name + ".scf.in"

cif_crystal = Structure.from_file(import_file)

crystal_lattice = cif_crystal.lattice.matrix
crystal_atom_num = np.size(cif_crystal.atomic_numbers)
crystal_atom_posi = cif_crystal.cart_coords
crystal_atom_kind = cif_crystal.atomic_numbers

#control

control_txt = "&control\n"\
" calculation   = 'scf',\n"\
" pseudo_dir = " + str(usp_files_dir) + ",\n"\
"/\n"


# system

system_txt = "&system\n"\
" ibrav         = 0,\n"\
" ntyp          = " + str(len(set(crystal_atom_kind))) + ",\n"\
" nat           = " + str(crystal_atom_num) + ",\n"\
" ecutwfc       = 25.0,\n"\
" ecutrho       = 225.0,\n"\
" nspin         = 1,\n"\
"/\n"

#electrons

electrons_txt = "&electrons\n"\
" conv_thr            = 1d-8,\n"\
" mixing_beta         = 0.7,\n"\
"/\n"

#atomic_species

atomic_species_txt = "ATOMIC_SPECIES\n"\
"" + str(usp_files)+"\n"\
"\n"

#atomic positions

priodic_table = ["H", "He",
                 "Li", "Be", "B", "C", "N", "O", "F", "Ne",
                 "Na", "Mg", "Al", "Si", "P", "S", "Cl", "Ar",
                 "K", "Ca", "So", "Ti", "V", "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn", "Ga", "Ge", "As", "Se", "Br", "Kr",
                 "Rb", "Sr", "Y", "Zr", "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd", "In", "Sn", "Sb", "Te", "I", "Xe"]

atom_posi_txt = ""

for atom_posi in range(crystal_atom_num*3):
    if atom_posi % 3 == 0:
        atom_posi_txt = atom_posi_txt + priodic_table[crystal_atom_kind[atom_posi//3] - 1] + "   " + str("{:.9f}".format(crystal_atom_posi[atom_posi // 3][atom_posi % 3])) + "   "
    elif atom_posi % 3 == 1:
        atom_posi_txt = atom_posi_txt + str("{:.9f}".format(crystal_atom_posi[atom_posi // 3][atom_posi % 3])) + "   "    
    elif atom_posi % 3 == 2:
        atom_posi_txt = atom_posi_txt + str("{:.9f}".format(crystal_atom_posi[atom_posi // 3][atom_posi % 3])) + "   \n"
    else:
        pass

atom_posi_txt = "ATOMIC_POSITIONS angstrom\n" + atom_posi_txt

#k_points

k_point_txt = "K_POINTS {automatic}\n"\
"3 3 3 0 0 0\n"\
"\n"

#crystal_lattice

cryatal_lattice_txt = ""

for cell_parameta in range(9):
    if cell_parameta % 3 == 2:
        cryatal_lattice_txt = cryatal_lattice_txt + str("{:.9f}".format(crystal_lattice[cell_parameta // 3][cell_parameta % 3])) + "   \n "
    else:
        cryatal_lattice_txt = cryatal_lattice_txt + str("{:.9f}".format(crystal_lattice[cell_parameta // 3][cell_parameta % 3])) + "   "
    
cryatal_lattice_txt = "CELL_PARAMETERS angstrom\n " + cryatal_lattice_txt


txt = control_txt + system_txt + electrons_txt + atomic_species_txt + atom_posi_txt + "\n" + k_point_txt + cryatal_lattice_txt + "\n"


with open(file_name, 'wb') as a_file:
    a_file.write(txt.encode('ASCII'))

    
print("finish!")

使い方

本プログラムでは下記の4つの項目を実行前に指定することが必要です。

1. name = "(filename)"
   作成ファイル名の指定です。(filename)を任意のファイル名に書き換えてください。
2. import_file = "(cif file directory)"
   cifの指定です。(cif file directory)をcifファイルが保存してあるディレクトリに書き換えてください。
3. usp_files_dir = "(UPS files directory)"
   UPSファイルの指定です。(UPS files directory)をUPSファイルが保存してあるディレクトリに書き換えてください。
4. usp_files = "(UPF file names)"
   UPSファイルの指定その二です。(UPF file names)をUPSファイル名に書き換えてください。

以上となります。

1. Quantum ESPRESSO公式のダンロード
   https://www.quantum-espresso.org/login/ ↩

2. おそらく東北大学で行われた講習会の資料かと思われます
   http://www.cmpt.phys.tohoku.ac.jp/~koretsune/SATL_qe_tutorial/preparation.html#virtualbox-materiappslive ↩

3. Anacondaでpythonを使っている場合にはこのページに従うとpymatgenをインストールできます。
   https://anaconda.org/conda-forge/pymatgen ↩

4. UPSファイルの公開ページ。
   http://pseudopotentials.quantum-espresso.org/legacy_tables ↩