Pymatgenチュートリアル番外① Pymatgenで作った構造をVESTAで表示する

###はじめに
####この記事で書くこと
この記事は、チュートリアル③の次に相当する。

チュートリアル③では、pymatgenを使って、cifファイルを読み込み、金属サイトとその周りの酸素を切り出すということをした。順当な記事の流れから言えば、次は切り出した構造から、さまざまな幾何学的な情報を取り出すという記事を書くのがいいだろう。とは言え、切り出した構造がどのようなものか、目で見て考えたいことはいろいろあるので、まずはこの方法を書いておきたい。構造を可視化するために、VESTA（K. Momma and F. Izumi, J. Appl. Crystallogr., 41, 653 (2008); K. Momma and F. Izumi, J. Appl. Crystallogr., 44, 1272 (2011).）を使う。
今さらVESTAについては説明不要だろうが、実はCUI的な操作が、部分的にではあるが実装されているということは、あまり知られていないように思う。泉先生ご本人によるEvernoteのページに、いくつかのコードが実例として記されている。VESTAの全ての機能が使えるというわけではないが、例えばPythonからの操作で、「古いcifを現代の規格化されたcifに変換するため､単に読み込んでまた保存する」だけでも、自動化できれば便利だ。（これは補遺1に記した。）

####使用環境
この記事で扱うコードは、Macで動作させている。おそらく、subprocessなどでPythonの外に投げるコードの記述がwindowsだと違うので、そこは別の記事を参照してほしい。例えば、「PythonからVESTAを操作してGIF動画を作る」など。他にも、フォルダの置き場所の制限とか、微妙に挙動が違うと思うので、そこは適宜修正してほしい。

pymatgen 2022.0.8
2020.x.xくらいのバージョンだと関数の呼び出し方が違っていてエラーが起きるかもしれないので注意。

###実際の操作
####配位環境を切り出すための関数の定義

まずは、前回の記事を元に、cifを読み込んで金属とその周りの酸素を取り出すための関数を定義する。前回記事の最後の部分を少し書き換えて、cifファイルへのパスを引数として、「金属と配位酸素サイト」からなるリストを作り、それらをvalueとする辞書型（金属サイトの名前がkey）を返す仕組みになっている。また、金属-酸素の結合長は出力しないことにした。

from pymatgen.io.cif import CifParser
import numpy as np
import re

def coenv(file_path):
    #一個目のcifファイルを読み込んで、mat0という名前を付ける
    parser = CifParser(file_path)
    mat0 = parser.get_structures()[0]

    #金属サイトをリストとして取り出す
    metal_sites = [site for site in mat0.sites if site.specie.is_metal]

    #各金属サイト周りの配位環境を取り出す
    coord_env = {}
    for i, metal_site in enumerate(metal_sites):
        #ある金属サイトから5Å以内のサイトを全て持ってくる
        neighbors = mat0.get_neighbors(metal_site, 5)
        #金属か否かでフィルタリングする
        neighbors_metal = [site for site in neighbors if site.specie.is_metal]
        #最近接の金属サイトまでの距離を得る
        mm_dist = [round(metal_site.distance(neighbor), 4) for neighbor in neighbors_metal]
        mm_dist_nonzero = np.array(mm_dist)[np.nonzero(mm_dist)]
        mm_min = min(mm_dist_nonzero)
        #上で得た最近接金属サイトよりも近くにあるサイトを持ってくる
        neighbors_O = mat0.get_neighbors(metal_site, mm_min-0.001)
        #配位酸素サイトのリストを作成し、辞書に登録する
        #リストの1つ目の要素は元の金属サイトにする
        coord_dist = [O_site for O_site in neighbors_O]
        coord_dist.insert(0, metal_site)
        coord_env[re.sub(r"[^a-zA-Z]", "", str(metal_site.specie))+'_'+str(i)] = coord_dist

    return coord_env

####構造の読み込み
さて、この関数を使って、なんでもいいのでcifファイルから配位環境を読み込むことにする。ファイルの読み込みには、例によってglobを使っている。

import glob
cif_list = glob.glob('./perovskite_cif/*.cif')
file_path =cif_list[11]

MOs = coenv(file_path)
print(MOs.keys())
display(MOs)

#出力結果：
# dict_keys(['Nd_0', 'Nd_1', 'Nd_2', 'Nd_3', 'Fe_4', 'Fe_5', 'Fe_6', 'Fe_7'])
#{'Nd_0': [PeriodicSite: Nd (0.0772, 5.3745, 5.8978) [0.0140, 0.9434, 0.7500],
#  PeriodicSite: O (-1.1151, 4.5412, 3.5375) [-0.2027, 0.7971, 0.4499],
#  PeriodicSite: O (-1.1151, 4.5412, 8.2581) [-0.2027, 0.7971, 1.0501],
#  PeriodicSite: O (-0.5209, 3.0118, 5.8978) [-0.0947, 0.5287, 0.7500],
#  PeriodicSite: O (-2.2292, 5.8604, 5.8978) [-0.4053, 1.0287, 0.7500],
#  PeriodicSite: O (1.6350, 4.0044, 4.3262) [0.2973, 0.7029, 0.5501],
#  PeriodicSite: O (1.6350, 4.0044, 7.4694) [0.2973, 0.7029, 0.9499],
#  PeriodicSite: O (1.1151, 6.8530, 4.3262) [0.2027, 1.2029, 0.5501],
#  PeriodicSite: O (1.1151, 6.8530, 7.4694) [0.2027, 1.2029, 0.9499]],
# 'Nd_1': [PeriodicSite: Nd (2.8273, 3.1712, 1.9659) [0.5140, 0.5566, 0.2500],
# （以下省略）

####金属サイトの取り出しとMoleculeへの変換
見ての通り、NdFeO3の構造を読み込み、各金属サイトと、その周りの酸素サイトを、リストで読み込むことに成功した。また、辞書のkeyには'金属_通し番号'が設定されている。
ここから1つの金属とその配位酸素を取り出し、さらに、取り出した構造を元に、pymatgenの定義する'Molecule'という型のデータを作ることで、VESTAで構造を表示する下準備をする。(補遺2にMoleculeの作成の仕方を簡単に記した。)

#Nd0と、周りの配位サイトを取り出す。
Nd0_O = MOs['Nd_0']

#.specie：site型のデータから、元素種を取り出す。
#ここでは内包表記を使って、Nd0_Oに含まれる全ての元素種を取り出している。
atoms = [site.specie for site in Nd0_O]

#.coords：site型のデータから、直交座標系での座標を取り出す。
#同じく内包表記を使って、全てのサイトの座標を取り出している。
cart_coords = [list(site.coords) for site in Nd0_O]

#分子を扱うためのモジュールを読み込む
from pymatgen.core import Molecule

#Moleculeで分子（Molecule型のデータ）を作成する。
#Molecule([元素種のリスト], [直交座標系での座標のリスト])という書き方をする
Nd0_O_mol = Molecule(atoms, cart_coords)

ちなみに、coordsの代わりにsite.frac_coordsのようにすると、fractional coordinatesが得られる。上記のように作ったMoleculeは、以下のように、pymatgen.core.structure.Moleculeで、化学種と座標のデータが含まれている。

print(type(Nd0_O_mol))
display(Nd0_O_mol)

#出力結果
#<class 'pymatgen.core.structure.Molecule'>
#Molecule Summary
#Site: Nd (0.0772, 5.3745, 5.8978)
#Site: O (-1.1151, 4.5412, 3.5375)
#Site: O (-1.1151, 4.5412, 8.2581)
#Site: O (-0.5209, 3.0118, 5.8978)

####xyzファイルへの変換とVESTAでの出力
さて、この分子をファイルに書き出して、VESTAで表示したい。「化学種と座標のデータ」を表示するのに過不足ないデータ形式と言えばxyzファイルで、pymatgenには、もちろんxyzを出力するための機能が備わっている。これは、以下のように使える。XYZ(Molecule型のデータ)という書式は、単にxyzファイルとして出力する前の準備のようなもので、.write_file(ファイル名)と組み合わせてxyzファイルを出力する。

from pymatgen.io.xyz import XYZ

file_name = 'Nd0_O'
XYZ(Nd0_O_mol).write_file(file_name + '.xyz')

!ls *.xyz
#出力結果：Nd0_O.xyz

最後の行は、ターミナルで'ls *.xyz'を実行しろ、という意味で、ワイルドカードを使って、フォルダ内にある拡張子.xyzのファイルの一覧を出力している。見ての通り、Nd0_O.xyzが作れたことがわかる。最後に、これをVESTAで表示しよう。ここで使っているsubprocessは、ターミナルで指定したコマンドを実行するためのもので、補遺3に簡単に使い方を記した。

import subprocess
subprocess.run(['open', '-a', 'VESTA', file_name + '.xyz'])

#xyzファイルは残しておいても邪魔なので、消してしまう。
import time
time.sleep(0.5)
subprocess.run(['rm', file_name + '.xyz'])

これを実行すると、VESTAが起動して、作成したMoleculeに対応した構造が表示される。ついでに、xyzファルは残しておいても邪魔（ファイルが無限に増える）なので、subprocessからrmコマンドを実行して、削除してしまっている。残しておきたい場合は、最後の行をコメントアウトするか、削除すればいい。ここでは、少し処理を待ってやらないといけないようなので、time.sleep(待ち時間)を使って、コードの進行を止めている。ここでは0.5秒止めているが、こんなに止めなくてもいい気はしている。(VESTAでxyzファイルが開かれる間だけ待つことになるので、どの程度必要かは、環境依存。)

####まとめ
最後に、金属サイトを選んでからMoleculeに変換し、VESTAで表示してファイルを削除するところまで、まとめると、以下のようになる。

#金属サイトを選ぶ
Nd0_O = MOs['Nd_0']

#Molecule型に変換する
atoms = [site.specie for site in Nd0_O]
cart_coords = [list(site.coords) for site in Nd0_O]
Nd0_O_mol = Molecule(atoms, cart_coords)

#Molecule型からxyz型に変換し、xyzファイルに変換する
file_name = 'Nd0_O'
XYZ(Nd0_O_mol).write_file(file_name + '.xyz')
subprocess.run(['open', '-a', 'VESTA', file_name + '.xyz'])
time.sleep(0.5)
subprocess.run(['rm', file_name + '.xyz'])

今回の記事のメインは以上。以下は補遺。

###補遺
####補遺1：VESTAの操作
とりあえずコードを記す。ここでは、cifを読み込んで、別のcifとして保存するという操作をする。ちなみに、僕の環境ではカレントディレクトリに読み込むcifがないとエラーが起きてしまったので、cifを保存しているフォルダから持ってくるというコマンドが入っている。環境によっては必要ないと思う。

import time

＃globで作成したcifのリストから、パスを一つ持ってくる
file_path =cif_list[0]

#/という文字でパスを区切ることで、ファイル名だけを取り出す
file_name = cif_list[0].rsplit('/')
#取得したファイル名を.で区切ることで、拡張子なしのファイル名にする
file_name_noex = file_name[-1].rsplit('.')[0]

#新しいcifファイルを保存するフォルダを指定する
#指定したフォルダが存在しない場合は、新しく作成する
path_output_dir = "./converted_cif"
if not os.path.exists(path_output_dir):
    os.makedirs(path_output_dir)

#カレントディレクトリにファイルを持ってくる
subprocess.run(['cp', file_path, '.'])
subprocess.run(['open', '-a', 'VESTA'])

#ここからしばらくは本文で解説する。
com = '''
-open {s}.cif
-save {fold}{s}_cv.cif
-save {fold}{s}.xyz
-close
'''.format(s = file_name_noex, fold = path_output_dir + '/').strip()

with open(file_name_noex + '.vstx', mode = 'w') as f:
  print(com, file = f)

subprocess.run(['open', '-a', 'VESTA', file_name_noex + '.vstx'])

#使ったファイルは消してしまう
time.sleep(0.0001)
subprocess.run(['rm', file_name[-1]])
subprocess.run(['rm', file_name_noex + '.vstx'])

さて、重要な部分は、com=```からしばらくである。.vstxというファイルを実行することでVESTAを操作するので、必然的に、
.vstxファイルの実行 ← .vstxファイルの保存 ← .vstxファイルの中身の作成
という手順を踏むことになる。

中身の作成は、

com = '''
-open {s}.cif
-save {fold}{s}_cv.cif
-save {fold}{s}.xyz
-close
'''.format(s = file_name_noex, fold = path_output_dir + '/').strip()

の部分で行っている。
・各行が1つのコマンドに対応している
・{変数}とし、変数への値の格納は最後の.format以下で行う
の2点を抑えておけばいいと思う。コマンドとして使えるものは、冒頭に示した泉先生のサイトを参照のこと。実際には、'com'という変数に文字列を格納している。
次に、この文字列をファイルに出力する。

with open(file_name_noex + '.vstx', mode = 'w') as f:
  print(com, file = f)

ここはPython入門とかで出てくるお決まりの書き方なので、特に解説はいらないと思う。
最後に、以下の部分で作った.vstxファイルを実行する。subprocessの使い方は、補遺3の通り。

subprocess.run(['open', '-a', 'VESTA', file_name_noex + '.vstx'])

####補遺2：Molecule型の作り方
公式は以下。
https://pymatgen.org/pymatgen.core.structure.html?highlight=molecule#pymatgen.core.structure.Molecule
分子を作るわけなので、最低限、元素と座標が指定されていればいい。つまり、

Molecule([元素のリスト], [座標のリスト])

という書き方をする。さらに情報を加えたい場合は、チャージ（分子全体）、磁気モーメントなどを割り振ることもできる。

####補遺3：!とsubprocessの使い方
すでにネット上に多くの情報があるので、ここでは本記事を読むのに必要な事項についてのみ解説する。!もsubprocessもターミナル（windowsの場合はcmd、Anacondaターミナルの場合もある）でコマンドを実行するために使う。違いは、
・!は、!に続けてコマンドをべた打ちする。変数は使えない。
・subprocessは、コマンドをリストの形で渡す。変数が使える。
である。!の方がサクッと使える、subprocessの方が細かく色々できる、という認識でもいいと思う。subprocessを使う場合は、「ターミナル上で空白で区切る部分ごとに要素に分けたリストを引数で渡す」という使い方をする。日本語で書くとなんのこっちゃという感じなので、以下に簡単な例を示す。

file_name = 'hoge.xyz'
subprocess(['rm', file_name])

これは、ターミナル上で以下のコマンドを実行したのと同じ結果をもたらす。

>> rm hoge.xyz

今回の記事は以上。