はじめに

前回の記事（https://qiita.com/ojiya/items/e98a9dd5cb6cd7ad38ca ）に引き続き、PymatgenとOptunaでXRDの解析を目指す。前回の記事では、cifファイルを元に、構造を読み込んで、構造の情報を取り出してみた。実際にXRDの解析を行う場面では、構造、つまり格子や原子座標を最適化しなければならない。最適化の際には、構造を取り出す → 構造を少し変更する → 構造を作り直すというプロセスが必要になるので、ここでは構造を作り直すところについて、解説していく。

実際のコードと説明

Structure型のデータを作るには、Pymatgen.core.structureモジュールの、Structureを使う。これは、以下のコードで読み込むことができる。

from pymatgen.core.structure import Structure

ここで読み込んだStructureは、引数として、格子（lattice）、元素種（species）、座標（coords）をとる。以下では、これらの作り方について、解説していく。

格子（lattice）

格子には、Pymatgen.core.lattice.Lattice型（以下、Lattice型）のデータを入力する。Latticeは、例えば、

from pymatgen.core.lattice import Lattice
lat = Lattice(引数)

のようにすれば、作ることはできる。ただし、この場合の引数はmatrix形式（直交座標空間でabc軸がなすベクトル）なので、第一原理計算などの結果を持ってくるのでもなければ、いささか扱いにくい。そこで、別の方法、つまり格子定数からLattice型のデータを作りたい。これは、基本的には以下のようにできる。

lat = Lattice.from_parameters(4,5,6,90,90,90)
print(type(lat))
#出力：<class 'pymatgen.core.lattice.Lattice'>

ここでは、雑にa = 4、b = 5、c = 6、α = β = γ = 90°の格子を作った。ただし、これは冗長というか、無駄なパラメータが入力されている。つまり、これは直方晶（orthorhombic）なので、それさえ先に知っていれば、後半の角度の情報は要らないはずだ。こういう場合に、最低限のパラメータだけを入力して格子を作る方法はきちんと用意されていて、

lat = Lattice.orthorhombic(4,5,6)

で、ひとつ前に示したコードと同じ格子を作ることができる。もちろん、orthorhombicの他に、立方晶（cubic）など一通り用意されている。使い方はいずれも同じで、Lattice.(晶系の名前)(非自明なパラメータの値)という書き方で格子を作ることができる。ただし、三斜晶（triclinic）だけは用意されていないことに注意。from_parametersと全く同じ作用なのでないのだと思う。
本記事でも、前回の記事で使ったマグネタイトFe3O4を用いて記事を書いていくので、Fe3O4の格子を作ってみる。つまり、以下のコードで表示される格子を作りたい。

from pymatgen.io.cif import CifParser
path = 'Fe3O4_mp-19306_symmetrized.cif'
parser = CifParser(path)
mat = parser.get_structures(primitive=False, symmetrized=False)[0]
lat_Fe3O4 = mat.lattice
display(lat_Fe3O4)
# 出力：
# Lattice
#    abc : 8.53335154 8.53335154 8.53335154
# angles : 90.0 90.0 90.0
# 以下略

見ての通り、a = b = c = 8.533.. Å、α = β = γ = 90°の立方晶格子である。これを作るためのコードは、以下のいずれかである。

lat_from_param = Lattice.from_parameters(8.53335154,8.53335154,8.53335154,90,90,90)
display(lat_from_param)

lat_cubic = Lattice.cubic(8.53335154)
display(lat_cubic)
# 出力は省略

.from_parametersでも.cubicでも同じ格子を作ることができた。

元素種と座標

続いて、元素種と座標を入力する。これら2つは当然セットなので、同時に作った方がいい。辞書型で作っておいて、後でバラすと、変なミスを減らせるだろう。最終的にはそれぞれのリストを作る、ということを念頭におきながら、keyがサイト名、valueが座標のリスト（xyzあるので）になっているような辞書を作っていく。0から作るのは大変なので、前回の記事で扱ったように、既に存在しているcifから構造を読み込んで、これをもとに解説していきたい。まずは、以下のコード（解説は前回記事）で非対称単位を取り出す。

from pymatgen.symmetry.analyzer import SpacegroupAnalyzer

sga = SpacegroupAnalyzer(mat)
mat_sym = sga.get_symmetrized_structure()
eq_sites = mat_sym.equivalent_sites

asym_unit = [site[0] for site in eq_sites]

なぜ非対称単位を取り出すかというと、構造の最適化のためだ。構造を「最適化する」とは、つまり構造を決定するパラメータ（ここでは座標）を少しずつ変えながら、XRDパターンを再現するために最もいいパラメータの値を探すことを言う。各サイトの座標が、非対称単位内のサイトから対称操作によって生成されるということは、換言すれば、空間群（=対称操作のセット）と非対称単位に含まれるサイトの座標が決定すれば、他のサイトの座標は自ずと決まるということになる。従って、最適化するべき座標は、非対称単位に含まれるサイトの座標のみである。

面倒なことに、非対称単位から、単位格子内の元素と座標のセットを出力する機能は、pymatgenには無いらしい。そこで、対称操作によって、これらを作ってやることになる。空間群に属する対称操作のセットを取り出して、座標を変換するのは、既に前回の記事に書いたので、まとめて書くと、

#対称操作を取り出す
sym_ops = sga.get_space_group_operations()
op = sym_ops[1]
rot = op.rotation_matrix
trans = op.translation_vector

# もとになるサイトを取り出す。
site_0 = asym_unit[0]
orig_coord = site_0.frac_coords
print(orig_coord)

# 対称操作を行う
gen_sites = []
for op in sym_ops:
    rot = op.rotation_matrix
    trans = op.translation_vector
    gen_coord = [round(coord-math.floor(coord),4) for coord in np.dot(rot, orig_coord) + trans]
    gen_sites.append(gen_coord)

print(gen_sites[0])
print(gen_sites[21])
# 出力：
# [0.125 0.125 0.625]
# [0.125, 0.875, 0.375]

あとは、サイトおよび対称操作について、for文を回してやればいい、、、というわけには、いかない。異なる対称操作によって、同じ場所にサイトが生成されることがあるからだ。具体例を見た方がわかりやすいと思うので、以下のコードを示す。

gen_sites = []
for op in sym_ops:
    rot = op.rotation_matrix
    trans = op.translation_vector
    gen_coord = [coord-math.floor(coord) for coord in np.dot(rot, orig_coord) + trans]
    gen_sites.append(gen_coord)

print(gen_sites[0])
print(gen_sites[21])
# 出力：
# [0.125, 0.125, 0.625]
# [0.125, 0.125, 0.625]

ちなみに、前回記事のコードとは少し変えて、round（四捨五入する関数）が入っているのは、計算の精度の都合で、「本当は全く同じ座標に来るはずなのに、ちょっとズレた座標に計算されてしまうサイト」があることがあり、これが大変都合が悪いためである。上記のコードでは、生成した座標をgen_sitesに格納しており、確かにgen_sitesの0番目と21番目に同じ座標が入っていることがわかる（他にも同じ座標がいくつかある）。というわけで、生成した座標のリスト（gen_sites）から重複を取り去りたい。ここでは、先人の知恵（ https://note.nkmk.me/python-list-unique-duplicate/ )をお借りして、以下の関数を定義し、用いることにする。

def get_unique_list(seq):
    seen = []
    return [x for x in seq if x not in seen and not seen.append(x)]

これを通すことによって、重複のない座標のリストを手に入れることができた。この部分は、関数にしておくと、後のコードがすっきりするので、関数を定義しておくことにする。つまり、

def gen_sym_sites(orig_coord, set_sym_ops):
    gen_sites = []
    for op in sym_ops:
        rot = op.rotation_matrix
        trans = op.translation_vector
        gen_coord = [round(coord-math.floor(coord),4) for coord in np.dot(rot, orig_coord) + trans]
        gen_sites.append(gen_coord)
    unique_sites = get_unique_list(gen_sites)
    return unique_sites

としておいて、

gen_sites = gen_sym_sites(orig_coord, sym_ops)

とすれば、gen_sitesに、対称操作によって生成した座標のリストが入る。
あとは、非対称単位（asym_unit）に含まれるサイトについてfor文を回してやれば、単位格子内の原子の座標を全て計算することができる。ただし、元素種がなんだったか、という情報が失われないよう、辞書型のデータに計算結果を格納することにして、keyを元素種、valueを座標のセット（gen_sites）とすることにする。コードは以下の通り。

dict_gen_sites = {}
for i, site in enumerate(asym_unit):
    orig_coord = site.frac_coords
    gen_sites = gen_sym_sites(orig_coord, sym_ops)
    gen_key = f'{str(site.species)}_{i}'
    dict_gen_sites[gen_key] = gen_sites
dict_gen_sites

これで、元素の情報を保持したまま、辞書に生成した座標のデータを格納することができた。試しに、1つ座標のセットを取り出してみる。ちなみに、さっき「keyを元素種にする」と書いたのに、通し番号が入るようになっているのは、非対称単位に同じ金属種がいる場合に、辞書の値が上書きされてしまうのを防ぐためである。

key_list = list(dict_gen_sites.keys())
print(key_list)
print(dict_gen_sites[key_list[1]])
#出力は省略

では、Structureの引数を作成すべく、以下のコードを用いる。ここのところは、色々と書き方があるが、入れ子のfor文にならないような書き方にしてみた。ここでポイントになるのは、Structure型のデータを作る際に用いるStructure()の引数のうち、サイトの元素種を指定する部分は、構成要素がComposition型、またはElement型であるようなリストにしなければならない、という点にある。下では、Composition型を使うパターンで書いてみた。実際のところ、置換サイトがある場合にはElement型は使えないので、Composition型の方がエラーが出にくいと思う。Composition型のデータは、pymatgen.core.compositionに含まれるComposition（）を使えば作ることができる。Composition()の引数は、文字列で、組成式をそのまま入れればよい（Fe1とかFe0.5Co0.5とか）。

from pymatgen.core.composition import Composition

list_species = []
list_coord = []
for k in key_list:
    spec = Composition(k.split('_')[0])
    list_gen_sites = dict_gen_sites[k]
    n_gen_site =len(list_gen_sites)
    list_species += [spec] * n_gen_site
    list_coord += list_gen_sites

以上で、無事にStructure()の引数を作ることができた。あとは、以下のようにすれば、Structure型のデータを作ることができる。

from pymatgen.core.structure import Structure
gen_structure = Structure(lat, list_species, list_coord)
display(gen_structure)

今回の記事は、以上。