##はじめに
###この記事で書くこと
チュートリアル②では、化合物を検索して、cifファイルを保存することができた。

そこで、この保存したcifファイルから、いろいろとデータを取り出してみたい。いろいろというのは、
・組成式
・格子の情報（空間群や格子定数）
・配位環境（結合長や配位サイト）
あたりを指す。もともとpymatgenには、構造を解析するためのツールがいろいろと含まれており、例えば構造の類似性を評価したりできるらしいが、ここではそういった解析手法は扱わない。材料化学の人間がcifを手に入れて最初にやろうとすることは、上に列挙したような事項だからだ。
(2021/7/25 追記)
この記事を公開したあと、pymatgenには配位構造を扱うための機能が備わっていることに気づいた。今のところ使い方がよくわからないので、また使い方がわかったら使うかもしれない。

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

##実際の操作
###cifの読み込み
まずは、cifを読み込むためのモジュールを読み込まないと、話が始まらない。

#CifParser：cifを読み込むためのモジュール
from pymatgen.io.cif import CifParser

早速cifを読み込みたい...が、その前に、どんなcifを所有しているか、一覧する。というわけで、次のコードを使う。

#glob：ファイルを読み込むときに便利なモジュール
import glob
cif_files = glob.glob('cif/*.cif')
display(cif_files)
#出力結果：
#['cif/mp-975179.cif',
# 'cif/mp-1183719.cif',
# 'cif/mp-23456.cif',
# 'cif/mp-865424.cif',
# 'cif/mp-554270.cif',
# 'cif/mp-1187557.cif',
# 'cif/mp-1209329.cif',
# 'cif/mp-19053.cif',
# 'cif/mvc-13900.cif',
# （以下略）

このように、glob(相対パス名)と書くと、指定した相対パス内にあるファイルの一覧（より正確にいうと、文字列型をとるファイルへのパスの一覧）をリスト型で返してくれる。上では、'*.cif'と指定することで、拡張子として「.cif」を持っているファイルを指定した。「 * 」（アスタリスク）はワイルドカード。（ちなみに、僕のフォルダには、前回の記事で読み込んだよりも、大量のcifが入っている。）
というわけで、所有しているcifファイルの一覧が手に入った。これを使って、cifを読み込もう。まずは、なんでもいいので1つ読み込んでみる。

#とりあえず最初のパスを指定した。
#CifParser(パス名)でcifファイルを読み込むことができる。
parser = CifParser(cif_files[0])

#読み込んだものに、.get_structure()[0]をつけることで、構造を読み込む
mat0 = parser.get_structures()[0]
print(type(mat0))
display(mat0)

単に読み込むだけなのに、ちょっとごちゃついている感があるが、とりあえず読み込むことができた。これを実行すると、以下の結果が得られる。（所有しているcifによって、内容自体は異なる。）

出力結果
<class 'pymatgen.core.structure.Structure'>
Structure Summary
Lattice
    abc : 4.465736 4.465736 4.465736
 angles : 90.0 90.0 90.0
 volume : 89.05927106621658
      A : 4.465736 0.0 2.734474649118552e-16
      B : -2.734474649118552e-16 4.465736 2.734474649118552e-16
      C : 0.0 0.0 4.465736
PeriodicSite: Rb (0.0000, 0.0000, 0.0000) [0.0000, 0.0000, 0.0000]
PeriodicSite: Tl (2.2329, 2.2329, 2.2329) [0.5000, 0.5000, 0.5000]
PeriodicSite: O (2.2329, 2.2329, 0.0000) [0.5000, 0.5000, 0.0000]
PeriodicSite: O (2.2329, 0.0000, 2.2329) [0.5000, 0.0000, 0.5000]
PeriodicSite: O (-0.0000, 2.2329, 2.2329) [0.0000, 0.5000, 0.5000]

また冒頭でtypeを判定しているが、ここで読み込まれた構造は、'pymatgen.core.structure.Structure'という種類のデータであることがわかる。見ての通り、格子定数、格子体積、原子座標が含まれている。A、B、Cというのがわかりにくいが、各格子軸方向の基底ベクトルだと思う（多分）。Structure型について、公式は以下のページ。

###構造情報の取得
####組成と格子の情報
ここでは、冒頭で述べた通り、組成式、格子の情報（空間群や格子定数）、配位環境（結合長や配位サイト）の取得の仕方について、記していく。組成式、格子の情報は、簡単に手に入る。

print(mat0.formula)
#出力結果：Rb1 Tl1 O3
print(mat0.get_space_group_info())
#出力結果：('Pm-3m', 221)
print(mat0.lattice)
#出力結果：
#4.465736 0.000000 0.000000
#-0.000000 4.465736 0.000000
#0.000000 0.000000 4.465736

そのまんまと言えばそのまんまだが、formulaで組成、get_space_group_info()で空間群の情報、latticeで基底ベクトルが手に入る。いちいち覚えるのめんどくさいなあ、という場合は、as_dict()を使うと、以下のように、一通りの情報が辞書型で手に入る。

display(mat0.as_dict())
#出力結果：
#{'@module': 'pymatgen.core.structure',
# '@class': 'Structure',
# 'charge': None,
# 'lattice': {'matrix': [[4.465736, 0.0, 2.734474649118552e-16],
#   [-2.734474649118552e-16, 4.465736, 2.734474649118552e-16],
#   [0.0, 0.0, 4.465736]],
#  'a': 4.465736,
#  'b': 4.465736,
#  'c': 4.465736,
#  'alpha': 90.0,
#  'beta': 90.0,
#  'gamma': 90.0,
#  'volume': 89.05927106621658},
# (以下略)

辞書型なので、keyを指定すれば、もちろん以下のように情報を取り出すことができる。スマートなやり方ではないけれど、うんうん悩んでるよりは、こうやった方が早いことも、ままある。

display(mat0.as_dict()['lattice'])
#出力結果：
#{'matrix': [[4.465736, 0.0, 2.734474649118552e-16],
#  [-2.734474649118552e-16, 4.465736, 2.734474649118552e-16],
#  [0.0, 0.0, 4.465736]],
# 'a': 4.465736,
# 'b': 4.465736,
# 'c': 4.465736,
# 'alpha': 90.0,
# 'beta': 90.0,
# 'gamma': 90.0,
# 'volume': 89.05927106621658}

####配位環境
さて、次は配位環境を手に入れたい。配位構造を手に入れるには、
・金属サイト（元素種と座標）の一覧を得る
・金属サイト周りの非金属サイトの一覧を得る
・得た一覧を元に、何らかの幾何構造を調べる
という手順を踏めば良い。とりあえず、金属か非金属かに限らず、サイトの一覧を手に入れたい。これは、sitesを使えば、簡単に手に入る。

print(type(mat0.sites))
#出力結果：<class 'list'>

display(mat0.sites)
#出力結果：
#[PeriodicSite: Rb0+ (0.0000, 0.0000, 0.0000) [0.0000, 0.0000, 0.0000],
# PeriodicSite: Tl0+ (2.2329, 2.2329, 2.2329) [0.5000, 0.5000, 0.5000],
# PeriodicSite: O0+ (2.2329, 2.2329, 0.0000) [0.5000, 0.5000, 0.0000],
# PeriodicSite: O0+ (2.2329, 0.0000, 2.2329) [0.5000, 0.0000, 0.5000],
# PeriodicSite: O0+ (-0.0000, 2.2329, 2.2329) [0.0000, 0.5000, 0.5000]]

display(mat0.sites[0])
#出力結果：PeriodicSite: Rb0+ (0.0000, 0.0000, 0.0000) [0.0000, 0.0000, 0.0000]

つまり、化合物に含まれるサイトが、リスト型で手に入ったことになる。リスト型なのだから、[]で番号を指定してやれば、一個取り出すことができる。リストの各要素はPeriodicSite型で、()内に、cartesian coordinate、[]内にfractional coordinateが書かれている。ちなみに、ここではRb0+とか、O0+とか、チャージが大変座りの悪いことになっているが、これをなんとかする方法は、本筋から外れるので、最後に補遺1として記載した。
さて、ここから金属サイトを取り出すには、チュートリアル①に書いた通り、is_metalを使えばいい。ただし、PeriodicSite型を別の型に変換しなければいけないので、ほんのわずかだけ書くことが増える。いきなり全サイトを処理する前に、1サイト取り出して、金属かどうかを判断してみる。

#specie：PeriodicSite型からSpecie型に変換する
print(type(mat0.sites[0].specie))
print(mat0.sites[0].specie.is_metal)
#出力結果：
#True
#<class 'pymatgen.core.periodic_table.Species'>

見ての通り、pymatgen.core.periodic_table.Species型に変換されている。Element型じゃないのか、と言われるかもしれないが、これらの2つの型の違いは、そのまんま元素と化学種の違いである、と思う。2つは異なるものだが、いずれにせよ金属か非金属かは判断できる。そして、Rb0+は金属なので、is_metalで判定を行った結果、Trueが返り値として得られている。
それでは、全サイトをまとめて処理したい。この文章は、あまりコーディングが得意でない人を意識して書いているので、何種類か書いてみた。

import numpy as np
mask = [site.specie.is_metal for site in mat0.sites]
metal_sites = np.array(mat0.sites)[mask]
display(metal_sites)

#出力結果：
#array([PeriodicSite: Rb0+ (0.0000, 0.0000, 0.0000) [0.0000, 0.0000, 0.0000],
#       PeriodicSite: Tl0+ (2.2329, 2.2329, 2.2329) [0.5000, 0.5000, 0.5000]],
#      dtype=object)

metal_sites = []
for site in mat0.sites:
    if site.specie.is_metal:
       metal_sites.append(site) 
display(metal_sites)

#出力結果：
#[PeriodicSite: Rb0+ (0.0000, 0.0000, 0.0000) [0.0000, 0.0000, 0.0000],
# PeriodicSite: Tl0+ (2.2329, 2.2329, 2.2329) [0.5000, 0.5000, 0.5000]]

metal_sites = [site for site in mat0.sites if site.specie.is_metal]
display(metal_sites)

#出力結果：
#[PeriodicSite: Rb0+ (0.0000, 0.0000, 0.0000) [0.0000, 0.0000, 0.0000],
# PeriodicSite: Tl0+ (2.2329, 2.2329, 2.2329) [0.5000, 0.5000, 0.5000]]

どれを使うかは、好みや個人の習熟度、一緒に働く人の理解度などを基準に決めればいいと思う。最初の例だけ出力結果がnp.arrayになっていることには注意してほしい。一番最初のやり方は、後でも使うので、そっちでの使い方も合わせて、補遺2に詳細を記した。
さて、これで金属サイトの一覧を手に入れることができた。続いて、金属周りの非金属サイト（配位子、ただし酸素も含む）を手に入れたい。ここで頭をコーディングから化学に切り替える必要がある。今は酸化物を扱っているので、例えば、「ある金属から最近接の金属よりも内側にある酸素」を取得すればいいだろう。そしたら、また頭をコーディングに戻して、「(i)まずは最近接の金属サイトまでの距離を取得し、(ii)それよりも結合距離が短い酸素サイトの一覧を得る」という手順をコーディングする。例によって、まずは金属サイト一個について、配位構造を取得してみる。このために、get_neighborsという関数を使う。

#後がごちゃごちゃするので、metal1という変数に、一個目のサイト（Rb）を入れておく
metal1 = metal_sites[0]
#get_neighbors(サイト、距離)：あるサイトから、一定の距離以内のサイトをリストとして取得する
neighbors = mat0.get_neighbors(metal1, 5)
display(neighbors)

#出力結果：
#[PeriodicSite: Tl0+ (-2.2329, -2.2329, -2.2329) [-0.5000, -0.5000, -0.5000],
# PeriodicSite: O0+ (-2.2329, -2.2329, -0.0000) [-0.5000, -0.5000, 0.0000],
# PeriodicSite: Tl0+ (-2.2329, -2.2329, 2.2329) [-0.5000, -0.5000, 0.5000],
# PeriodicSite: O0+ (-2.2329, 0.0000, -2.2329) [-0.5000, 0.0000, -0.5000],
# PeriodicSite: Tl0+ (-2.2329, 2.2329, -2.2329) [-0.5000, 0.5000, -0.5000],
# PeriodicSite: O0+ (-2.2329, 0.0000, 2.2329) [-0.5000, 0.0000, 0.5000],
#(以下省略)

つまり、structure.get_neighbors(サイト、距離(Å))という書き方をすれば、指定したサイトから指定した距離よりも内側にあるサイトをリストとして返してくれる。ちなみに、ここでは処理が重くなることよりも、取りこぼしを恐れているので、かなり長い距離(5 Å)を指定している。そしたら、このリストから、金属だけを取り出す。このやり方は、すでに見た通りで、以下のように行う。個人的な好みで、内包表記を使っている。

neighbors_metal = [site for site in neighbors if site.specie.is_metal]
display(neighbors_metal)
#出力結果：
#[PeriodicSite: Tl0+ (-2.2329, -2.2329, -2.2329) [-0.5000, -0.5000, -0.5000],
# PeriodicSite: Tl0+ (-2.2329, -2.2329, 2.2329) [-0.5000, -0.5000, 0.5000],
# PeriodicSite: Tl0+ (-2.2329, 2.2329, -2.2329) [-0.5000, 0.5000, -0.5000],
# PeriodicSite: Tl0+ (-2.2329, 2.2329, 2.2329) [-0.5000, 0.5000, 0.5000],
#（以下省略）

そしたら、このリストに格納されているサイトと、元のサイト(metal1)との距離を、やはりリストで得て、そのリスト中の最小値を得る、ということを行う。

#distance：2つのサイト間の距離を求める。
#round(少数点数、整数)：少数点数を、整数桁目で四捨五入した値を得る。
mm_dist = [round(metal1.distance(neighbor), 4) for neighbor in neighbors_metal]

#numpyのnonzeroを使って、0以外の値を取り出す。
mm_dist_nonzero = np.array(mm_dist)[np.nonzero(mm_dist)]

#minを使って、最小値を得る。
mm_min = min(mm_dist_nonzero)
print(mm_min)
#出力結果：3.8674

これで、無事Rb周りの金属-金属間距離の最小値を得ることができた。また、ここで使っているnonzeroについては、補遺2に説明を加えた。
続いて、これ以下の結合距離でRb-O結合を作るような酸素の一覧を得ることを目指す。これはすでに一回使ったget_neighbors()を使えばよい。つまり、以下の通り。

neighbors_O = mat0.get_neighbors(metal1, mm_min-0.001)
display(neighbors_O)
#出力結果：
#[PeriodicSite: O0+ (-2.2329, -2.2329, -0.0000) [-0.5000, -0.5000, 0.0000],
# PeriodicSite: O0+ (-2.2329, 0.0000, -2.2329) [-0.5000, 0.0000, -0.5000],
# PeriodicSite: O0+ (-2.2329, 0.0000, 2.2329) [-0.5000, 0.0000, 0.5000],
# PeriodicSite: O0+ (-2.2329, 2.2329, 0.0000) [-0.5000, 0.5000, 0.0000],
# PeriodicSite: O0+ (0.0000, -2.2329, -2.2329) [0.0000, -0.5000, -0.5000],
#（以下省略）

print(len(neighbors_O))
#出力結果：12

無事、酸素だけを取り出すことができた。ちなみに、最近接の金属サイトまでの距離をそのまま使うと、数字の丸め具合によっては金属サイトが入ってしまうので、0.001をmm_minから引いている。また、最後にlenを使って調べた通り、得られた酸素のリストは要素が12個からなる。もともとABO3型の組成を持つ化合物を扱っていたので、Rbサイトの配位数が12と出たら、なるほどペロブスカイト、と判断できて安心。最後に、元の金属サイトと、得られた各酸素サイト、及びその間の距離をひとまとめにしておきたい。以下のように辞書型でまとめておけば、使い勝手が良さそうだ。

{サイト名：[[0サイト1、結合距離]、[0サイト2、結合距離]、[0サイト3、結合距離]...]}

要するに、サイト名をkeyとし、これに対して、入れ子のリストになったvalueが割り当てられている、という構造にする。ただし注意したいのは、今の化合物では金属サイトが1化学種に対して1サイトだったが、一般的にはそうではない。ということで、金属サイトに名前をつける、という点も踏まえてコーディングする必要がある。

#まずは、入れ子になったリストを準備する
coord_dist = [[O_site, metal1.distance(O_site)] for O_site in neighbors_O]
#specieをstr()の中に入れて文字列に変換し、'_1'を足すことで、単純な化学種名からサイト名にしている。
coord_env = {str(metal1.specie)+'_1':coord_dist}
display(coord_env)
#出力結果：
#{'Rb0+_1': [[PeriodicSite: O0+ (-2.2329, -2.2329, -0.0000) [-0.5000, -0.5000, 0.0000],
#   3.157752208588888],
#  [PeriodicSite: O0+ (-2.2329, 0.0000, -2.2329) [-0.5000, 0.0000, -0.5000],
#   3.1577522085888874],
#  [PeriodicSite: O0+ (-2.2329, 0.0000, 2.2329) [-0.5000, 0.0000, 0.5000],
#   3.1577522085888874],
#  [PeriodicSite: O0+ (-2.2329, 2.2329, 0.0000) [-0.5000, 0.5000, 0.0000],
#   3.157752208588888],
#（以下省略）

ここで、'Rb0+_1'の'0+'がうっとおしい、と感じる場合には、次のようにすればいい。

import re
coord_dist = [[O_site, metal1.distance(O_site)] for O_site in neighbors_O]
coord_env = {re.sub(r"[^a-zA-Z]", "", str(metal1.specie))+'_1':coord_dist}
display(coord_env)

出力結果は省略したが、こうすると、keyが'Rb_1'になる。これについては、完全に好みの範囲なので、詳しい説明は、ここでは行わない。気になる場合は、「文字列からアルファベットを取り出す python」などで検索すると、いろいろ出てくる。

これで、Rbサイト周りの非金属サイト（今は酸素のみ）を持ってくることができた。あとは、構造に含まれる金属サイト全てに対し同様の操作を行う、というコードを書こう。金属サイトのリストは手に入っている(metal_sites)ので、このリストの要素に関して、forのループを回せばいいだけだ。せっかくなので、以下には、cifを読み込むところから、最後の辞書を得るところまで、全て一括で書いた。

#特定のフォルダ内にあるcifファイル（へのパス）のリストを得る
cif_files = glob.glob('perovskite_cif/*.cif')

#一個目のcifファイルを読み込んで、mat0という名前を付ける
parser = CifParser(cif_files[0])
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)
    #配位酸素サイトのリストを作成し、辞書に登録する
    coord_dist = [[O_site, metal_site.distance(O_site)] for O_site in neighbors_O]
    coord_env[re.sub(r"[^a-zA-Z]", "", str(metal_site.specie))+'_'+str(i)] = coord_dist

display(coord_env)

長いので出力結果は省略したが、「あるcifファイルに含まれる金属サイト周りの配位環境の一覧」を取り出すことができた。最初のcif_files[0]の部分を[1]とか[2]に変えて、結果が変わるのを確認してほしい。実際に使うときは、これの一部あるいは全体を関数として定義すると使いやすそう。本当は六配位擬八面体の歪み評価（Jahn-Teller歪みの有無など）なども書きたかったが、それを始めると長くなりすぎるので、本記事はここまでとする。以下は補遺。

###補遺
####補遺1：各サイトへのチャージの割り振り
デフォルトでは、全てのサイトのチャージが0になってしまっている。これが気持ち悪いので、各サイトにチャージを割り振りたい...と言いたいところだが、よく考えたら、この目的にRbTl03を使うのは、非常に都合が悪い(素直にチャージを割り振れない。うっかりしながらここまで書いてしまった。)ので、mat0を別の化合物(CrRhO3)に置き換えることにする。この化合物には、以下のようなサイトが含まれている。

display(mat0.sites)
#出力結果：
#[PeriodicSite: Cr (0.0000, 0.0000, 0.0000) [0.0000, 0.0000, 0.0000],
# PeriodicSite: Rh (1.9585, 1.9585, 1.9585) [0.5000, 0.5000, 0.5000],
# PeriodicSite: O (1.9585, 1.9585, 0.0000) [0.5000, 0.5000, 0.0000],
# PeriodicSite: O (1.9585, 0.0000, 1.9585) [0.5000, 0.0000, 0.5000],
# PeriodicSite: O (-0.0000, 1.9585, 1.9585) [0.0000, 0.5000, 0.5000]]

これにチャージを割り振るのは、半自動でできる。Composition型にoxi_state_guesses()を作用させると、タプルに辞書を格納した形で推測したチャージを返してくれる。

ox_guess = mat0.composition.oxi_state_guesses()
display(ox_guess)
#出力結果：
#({'Cr': 3.0, 'Rh': 3.0, 'O': -2.0}, {'Cr': 2.0, 'Rh': 4.0, 'O': -2.0})

ここで得られた辞書を使って、Structure型の構造に対してadd_oxidation_state_by_element(辞書)を作用させると、各サイトにチャージが割り振られる。

mat0.add_oxidation_state_by_element(ox_guess[0])
display(mat0.sites)
#出力結果：
#[PeriodicSite: Cr3+ (0.0000, 0.0000, 0.0000) [0.0000, 0.0000, 0.0000],
# PeriodicSite: Rh3+ (1.9585, 1.9585, 1.9585) [0.5000, 0.5000, 0.5000],
# PeriodicSite: O2- (1.9585, 1.9585, 0.0000) [0.5000, 0.5000, 0.0000],
# PeriodicSite: O2- (1.9585, 0.0000, 1.9585) [0.5000, 0.0000, 0.5000],
# PeriodicSite: O2- (-0.0000, 1.9585, 1.9585) [0.0000, 0.5000, 0.5000]]

もちろん、ここで入力に使う辞書を自分で指定すれば、別のチャージを割り当てることができる。手動でチャージを決めた時は、ちゃんとバランスが取れているか、確かめた方が良い。これをするには、チャージを割り当てたあと、chargeを作用させる。以下の例では、わざとおかしな割り当てを行い（O-）、全体のチャージを調べている。確かに、チャージバランスが崩れていることがわかる。

charge_false = {'Cr': 3.0, 'Rh': 3.0, 'O': -1.0}
mat0.add_oxidation_state_by_element(charge_false)
print(mat0.charge)
#出力結果：3.0

####補遺2：numpyでのフィルタリングの仕方
この記事では、numpyを使って、2種類のフィルタリングを行っている。1つ目は、bool(True / False)を使ってフィルタリングする方法で、
・n個の要素を持つnp.arrayから要素を抜き出すために、
・n個の要素を持つbool型の値を持つリストを用いる方法で、
・i番目の要素がTrueであれば、元のnp.arrayのi番目の要素は出力結果に値が残留し、
・i番目の要素がFalseであれば、元のnp.arrayのi番目の要素は出力結果から除外される。
日本語で書くとよくわからないが、以下のようになる。

test0 = np.array([1,2,3,4,5])
filter = [True, True, False, False, True]
test0[filter]
#出力結果：array([1, 2, 5])

よくわからない、という人はtest0の値と、filterのTrue / Falseの位置、最後の出力結果を、よく見比べてほしい。

もう1つのフィルターは、np.nonzeroを使うもので、読んで字の如し、0かそうでないかを判定する。これを、使ってフィルターになるnp.arrayを作る。続いて、np.arrayは、取り出したい番号のリストやnp.arrayを使って、複数の要素を取り出すことができることを利用し、0でない値を取り出す。これも日本語で書くとよくわからないので、以下の例を示す。

test1 = np.array([1,0,1,0,0,2,3,0.1])
print(np.nonzero(test1))
#出力結果：(array([0, 2, 5, 6, 7]),)

filter = (np.nonzero(test1))
print(test1[filter])
#出力結果：[1.  1.  2.  3.  0.1]

今回の記事は、以上とする。