はじめに

この記事で書くこと

Pymatgenには、Quantum EspressoやABINITをはじめとする、各種計算用の入力ファイルを作成する機能が付いているが、使い方についての情報がほとんどなく、公式ドキュメントも、「もともと使える人にはわかる」レベルで、導入がなかなか難しかった。そこで、この文章を書いておくことにした。入力ファイルを作成する部分だけのつもりでいたが、せっかくなので、cifファイルのコレクションをもとに、自動で計算を実行してもらうことにした。実際、世に出回っているcifファイルの質はまちまちで、「とりあえず雑に最適化を一回噛ませたほうが使いやすい」というシチュエーションはよくある。とりあえず拾い集めてきたcifファイルをまとめてフォルダに保存しておいて、サクッと一回計算を回しておく、くらいの使い方ができればいいだろう。ここではABINITを使うが、Quantum Espressoでも、ほとんど同じ使い方ができる。

使用するバージョン

pymatgen 2022.0.8
2020.x.xくらいのバージョンだと関数の呼び出し方が違っていてエラーが起きるかもしれないので注意。
ABINIT 9.4.2
操作はMacで行っている。Windowsだとsubprocess内のコマンドを変えなければいけない、、、と思うが、WindowsでABINITをコンパイルできる人にとっては、全く問題にならないだろう。

事前準備

擬ポテンシャルファイルの準備

公式ページから、擬ポテンシャルファイルをダウンロードしておく。"Download the entire dataset table"と書かれたところから、全部ダウンロードできるので、これを持ってきて、"PPs"というフォルダに入れておく。
https://www.abinit.org/psp-tables

フォルダ構成の確認

計算を行なっていくに当たって、必要なものは、
・実行するコード
・計算のもとにするcifファイルの入ったフォルダ
・擬ポテンシャルファイルの入ったフォルダ
で、これらをフォルダの同じ階層に入れておく。要するに、↓のような状態。（9_calc_input.ipynbに実行するコードが書いてある。このファイルの名前はなんでもいい。）

実行するコード

肝になる関数

ここでは、Pymatgenの、BasicAbinitInputという関数を使う。基本的な使い方は、以下の通り。

BasicAbinitInput(structure型のデータ, 擬ポテンシャルファイルへのパスのリスト, 計算用キーワードを格納した辞書)

例えば、適用なcifファイルを読み込んできて、最低限のパラメータ（計算でエラーが出ない、の意味）だけを入力したいなら、以下のようなコードを書けばよい。とりあえず書いておいて、後から部分ごとに解説することにする。なお、今回はformat()を多用するので、簡単な使い方を補遺に記した。

# このコードを実行する前に、適当なフォルダを作成し、そのフォルダまでの相対パスをdir_nameという変数に格納しておく。
# また、そのフォルダの中に、PPsというフォルダを作っておく。

from pymatgen.io.cif import CifParser
from pymatgen.io.abinit.inputs import BasicAbinitInput
import glob
import subprocess

# cifファイルを読み込む
cif_path = glob.glob('./cifs/*.cif')
cif_file = cif_path[0]

# 構造を読み込む
parser = CifParser(cif_file)
mat0 = parser.get_structures()[0]

# 擬ポテンシャルファイルをコピーしてくる
for el in mat0.composition.elements:
    el = str(el)
    cmd_PPs = 'cp ./PPs/{0}.psp8 ./{1}/PPs'.format(el, dir_name)
    subprocess.run(cmd_PPs, shell=True)

# 入力ファイルを作成する
PPs = ['./{0}/PPs/{1}.psp8'.format(dir_name, el) for el in mat0.composition.elements]
param_dict = {'ecut': 5,
              'tolvrs' : 1.000000E-5,
              'pseudos': '"{}"'.format(', '.join(PPs))}
input_file = BasicAbinitInput(mat0, PPs, abi_kwargs=param_dict).to_string()

構造を読み込むところまでは、すでに過去の記事で述べた通り。先のスクリーンショットにあるように、"cifs"というフォルダにcifが入っているので、その中から読み込んでいる。
擬ポテンシャルファイルをコピーしてくるところについては、以下のポイントを押さえればよい。
・擬ポテンシャルファイルのファイル名は、"（元素記号）.psp8"で統一されている。
・必要な擬ポテンシャルファイルは、PPsというフォルダ（相対パスで"./PPs"）に入っている。
・必要な擬ポテンシャルファイルは、組成式に含まれる元素（mat0.composition.elementsで取得可能）。
これらを前提として、subprocessでファイルをコピーしてくればいい。ただし、mat0.composition.elementsが返すリストの要素はElements型なので、これをstr(Elements)として、文字列に直す必要がある。

続いて、入力ファイルを直接作成する関数、BasicAbinitInput()の使い方は、以下の通りである。
・引数は最低3つ必要で、前から順にStructure型、擬ポテンシャルへのパスのリスト、キーワードを格納した辞書
・Structure型は、cifから読み込んだものをそのまま。
・擬ポテンシャルファイルへのパスは、「実行するパイソンファイルを基準とした相対パス」である。
・計算用キーワード（ecut、tolvrsなど）は、引数abi_kwardsに辞書型でkeyにキーワード、valueに値を入れる。
・BasicAbinitInputの返り値はBasicAbinitInput型なので、これをto_string()として、文字列に直す必要がある。
ここで気をつけたいのは、BasicAbinitInput()は擬ポテンシャルファイルへのパスが必須の引数となっているのに、自動ではpseudosに値を割り当ててくれない点（なんで？）。なので、abi_kwardsの辞書には、擬ポテンシャルファイルのパスも含めなければならない。また、ここでは最終的にABINITの入力ファイルに記載されるパスを渡す必要があるので、「ABINITのインプットファイルを保存するディレクトリを基準とした相対パス」を入力しなければならない。

計算を自動で実行するためのコード

ディレクトリの作成

ABINITに限らず、計算用のソフトウェアで計算を行うと、大量の出力ファイルが生成する。これが一緒くたになっていると混乱するので、計算を行う構造ごとに、別々のフォルダを作っておくのがいいだろう。ただし、これを手作業でやるのはめんどくさいので、自動で行えるようにする。ディレクトリ名はなんでもいいが、一意性とわかりやすさを両立するためには、「通し番号_組成式」くらいにしておくのがいいだろう。ここでは、Materials Projectのcifを使っているので、cifファイル名からMaterial IDを持ってきて、これを組成式と合わせてディレクトリ名にしている。

# 計算用のディレクトリを作る
import os
dir_name = '{0}_{1}'.format(cif_file[7:-4], mat0.formula.replace(' ', ''))
if not os.path.exists(dir_name):
    os.makedirs(dir_name)
    os.makedirs(dir_name+'/PPs')

ここは、すでに以前の記事で取り扱ったもので、ディレクトリが存在しなければ、作成するようにしている。
ディレクトリ名の行は、globでとってきたcifファイルへのパスが"./PPs/(materials ID).cif"となっているので、7文字目からスタート、お尻から数えて4文字目より前をとってくれば、Material IDが手に入る。また、structure型のデータ(mat0)から、.formulaで組成式を取り出し、.replace('置換前の文字列', '置換後の文字列')とすることで、空白を削除している。

入力ファイルの保存

ここでは、先ほど作成したディレクトリに入力ファイルを保存している。これは基本的なPythonの使い方なので、特に説明は不要と思う。

input_file_path = './{}/basic_calc.abi'.format(dir_name)
with open(input_file_path, mode='w') as f:
    f.write(input_file)

計算の実行とアウトプットファイルの移動

subprocessコマンドは、実行したコマンドが終了するのを待ってくれるため、普通にABINITの実行コマンドを実行すればよい。微妙に気をつけたいのが、ABINITのチュートリアルでは、実行コマンドを"abinit hoge.abi >& log &"としている場合があるが、こうすると、おそらくバックグランドで計算を実行したまま、コードが進んでしまうので、おそらく大変な目に遭う。最後の&は付けないこと。また、各種アウトプットファイルは「計算を実行した時のカレントディレクトリ」に保存されるため、自動では先に作成したディレクトリに保存してもらえない。手動でmvコマンドを使って移動する。

# 計算を実行する
cmd_calc = 'abinit ./{}/basic_calc.abi >& log'.format(dir_name)
subprocess.run(cmd_calc, shell=True)

# ファイルをディレクトリに保存する
cmd_mv = 'mv basic_calc* ' + dir_name
subprocess.run(cmd_mv, shell=True)

計算前の確認コード

上記が終われば、基本的には計算を自動実行することができる。しかし、存在しない擬ポテンシャルファイルを指定してしまうと、そこでコードの実行がエラーで止まってしまうので、先に以下のコードを実行して、計算のもとになるcifのリストを作成し直しておくと無難。特に、ランタノイド系列の擬ポテンシャルファイルは、ABINITの公式ページからは取得できないので、注意が必要。

cif_path = glob.glob('./cifs/*.cif')
PP_files = glob.glob('./PPs/*.psp8')
PP_els = [element[6:-5] for element in PP_files]
calc_path = []

for i, cif_file in enumerate(cif_path):
    parser = CifParser(cif_file)
    mat0 = parser.get_structures()[0]
    mat0_el = [str(el) for el in mat0.species]
    if all([el in PP_els for el in mat0_el]):
        calc_path.append(cif_file)
    else:    
        print(cif_file[7:-4])
        print(mat0.formula)

やっていることとしては、
・PP_elsに、PPsに入っている擬ポテンシャルファイルの拡張子以前の部分（つまり元素記号の部分）を格納する。
（相対パス"./PPs/*.psp8"なので、7文字目から最後の5文字目より前をとってくればよい。）
・mat0に含まれる元素をmat0.speciesで取得してリスト(mat0_el)に格納
・mat0_elに含まれる文字列が、全てPP_elsに含まれるかを確認。
・全てが含まれていれば、calc_pathに格納。
・もしも、含まれない文字列があれば、printで出力しておく。
となっている。

コードの全体像

最後に、まとめて計算を実行するためのコードを記す。

for cif_file in calc_path:
    #構造を読み込む
    parser = CifParser(cif_file)
    mat0 = parser.get_structures()[0]

    #計算用のディレクトリを作る
    dir_name = '{0}_{1}'.format(cif_file[7:-4], mat0.formula.replace(' ', ''))
    if not os.path.exists(dir_name):
        os.makedirs(dir_name)
        os.makedirs(dir_name+'/PPs')

    #擬ポテンシャルファイルをコピーしてくる
    for el in mat0.composition.elements:
        el = str(el)
        cmd_PPs = 'cp ./PPs/{0}.psp8 ./{1}/PPs'.format(el, dir_name)
        subprocess.run(cmd_PPs, shell=True)

    #入力ファイルを作成する
    PPs = ['./{0}/PPs/{1}.psp8'.format(dir_name, el) for el in mat0.composition.elements]
    param_dict = {'ecut': 5,
                  'tolvrs' : 1.000000E-5,
                  'pseudos': '"{}"'.format(', '.join(PPs))}
    input_file = BasicAbinitInput(mat0, PPs, abi_kwargs=param_dict).to_string()

    #入力ファイルを保存する
    input_file_path = './{}/basic_calc.abi'.format(dir_name)
    with open(input_file_path, mode='w') as f:
        f.write(input_file)

    #計算を実行する
    cmd_calc = 'abinit ./{}/basic_calc.abi >& log'.format(dir_name)
    subprocess.run(cmd_calc, shell=True)
    
    #ファイルをディレクトリに保存する
    cmd_mv = 'mv basic_calc* ' + dir_name
    subprocess.run(cmd_mv, shell=True)

これで、擬ポテンシャルが保存されている構造については、計算を行なってくれる。今回は以上。

補遺

formatの使い方

文字列の操作を行う時に、なんでもかんでも'+'で結合させてもいいけれど、それでは最終的に出来上がる文字列の雰囲気がわかりにくいし、引用符が入り混じってわかりにくくなる（特に、最終的に得たい文字列にも引用符が含まれている場合）ので、エラーの温床になりやすい。ということで、formatを使うことにした。実際のところ、formatは初心者にはわかりにくい面がある。僕の場合は、''や""で囲まれた文字列は、不可触というイメージがあり、その中をいじれるという感覚がよくわからなかったためだ。なんにせよ、formatは、以下の構文で使うことができる。

'文字列xxx{}文字列yyy'.format(置換後の文字列zzz)
# 文字列xxx置換後の文字列zzz文字列yyy

複数の部分に代入したい場合は、

'文字列xxx{0}文字列yyy{1}'.format(置換後の文字列zzz, 置換後の文字列aaa)
# 文字列xxx置換後の文字列zzz文字列yyy置換後の文字列aaa

となる。