はじめに

材料研究とX線回折（X-ray diffraction、XRD）パターンの測定・解析は、切っても切り離せない関係にある。というか、「材料粉体を合成したら、まずはXRD」が基本中の基本と認識している。要するに合成した材料は、全てXRDをとる、というのが、材料化学者の嗜みである。必然的に、手元には大量のXRDデータが得られる。これを効率的に解析することは、業務効率化に欠かせないテクニックだろう。近年、こうしたアイデアは随所で見られ、XRDを自動解析するソフトウェアはいくつか報告されている。中でも、optunaとGSAS-IIを利用して、自動でRietveld解析を行うという論文（およびコード）は、界隈でかなりバズった。これに対して、忖度なく、すごい！と思うと同時に、optunaを使うなら、最初からoptunaで全部最適化すればいいのでは？という気持ちもあった。そういうわけで、簡易的にこれを実装するコード(フォークト関数でパターンをシミュレーションして、格子定数および原子座標、デバイ・ワラー因子を最適化するだけ)を組んでみたので、解説記事を書こうと思う。ただし、かなり長くなったので、いくつかの記事に分けることにした。それから、せっかくなので、用いるPymatgenの各機能の説明も加えることにする。本記事では、Pymatgenを使って構造を扱う部分について、詳しく述べる。

やりたいことの全体像

1.結晶構造の作成
2.1.からXRDパターンの計算
3.実測のパターンとの差分を計算する
4.3.で得た差分が最小化するように結晶構造のパラメータを調節する

この記事で書くこと

・結晶構造と対称性の話
・Pymatgenでの構造の読み込みかた
やりたいことの全体像からすると、1.の前半部分くらいの内容を書く。

使用するライブラリのバージョン

Pymatgen 2022.0.7
2021.x.x以前では、構造の読み込みの際に設定できるパラメータが違うようなので注意。

結晶構造と対称性の話

非材料化学系の人のために、少しだけ結晶構造の説明を加えておく。詳しく知りたい方は、「物質の対称性と群論」を読むこと。(リンク → https://www.amazon.co.jp/物質の対称性と群論-今野-豊彦/dp/4320034090 ) すでにわかってるよ、という方は、読み飛ばしてください。なお、ここでは並進対称性を持った構造についての話だけをすることにする（準結晶には当てはまらないことを書く、という意味。ほとんど確実に並進対称性は満たしているので、特に気にする必要はない。）。
ある構造をなす結晶が、空間中に存在しているとき、当然、その構造を目で見て観測するようなことはできず、多くの場合はX線のような光線との相互作用を観測し、それを生データとして解析をおこなっていくことになる。このとき、観測されたデータの直接の解釈では、実際の結晶の座標のようなもの（長さの次元を持つ、「実空間」での座標）は得られない。得られるものは、長さの逆数の次元を持つ、「逆空間」での情報である。長さの逆数の次元をもつ、ということは、その意味合いは「単位長さあたりの数」に相当しているはずだ。すると、結晶構造には、本質的に「周期性」や「対称性」のようなものが含まれていることがわかる(*)。では、結晶構造を記述するにも、「周期性や対称性」を用いれば良いだろう。では、周期性や対称性を用いて結晶構造を記述する、とは、具体的にどういう行為のことを指しているだろうか。
いきなり答えを書くと、「周期性や対称性を用いて結晶構造を記述する」とは、「結晶の一部に対して対称操作を行うことで、別の部分を生成できるようにする」ことを指す。ここは、具体例を見た方が話が早いと思う。下に雑な図を書いたが、丸1つと四角4つが結合したような分子があって、それが平行に並んで結晶を構成していたとする。図には合計で丸を2つ、四角を8つ書いたが、この構造を記述するのに、合計で10個の座標を使う必要はない。まず、左側の分子を記述するには、丸の座標が1つ、四角の座標が1つあれば良い。残った3つの四角は、すでにある1つの四角を、原点（丸があるところ）周りに90°回転させる操作と、原点に対して対称の点を生成する操作があればいい。従って、構造を記述するためには、座標2つ分と、対称操作2つがあればいい。このような、「これだけ座標を記述してあれば、他の部分は対称操作で生成できますよ」という最小単位のセットを、「非対称単位」と呼ぶ。また、「非対称単位に対して、これだけ対称操作を加えますよ」という対称操作のセットには決まりがあって、これを「空間群」と呼ぶ（230個ある。全部そらで言える人はあまりいないと思う）。また、図中の右の分子は、左の分子を平行移動することにより生成することができる。このような操作は「並進操作」と呼び、「並進操作によって、結晶の他の部分を生成することができるような単位」のことを「単位格子」と呼んでいる。下の図では単位格子を2次元の正方形で描いたが、実際にはもちろん3次元で、一般には立方体ではない（立方体のものもある）。単位格子を構成する辺は「軸」と呼び、a, b, c軸がある。また、各軸のなす角は、α, β, γと呼ばれる。

(*)準結晶は、このような周期性や対称性を持たない。結晶学では、長らく結晶の定義を、ここで述べているように「並進対称性を持つもの」と定めていたが、色々あって、今では準結晶も結晶に含まれるようになった。そういう意味で、ここに書いてあることは少し古い考え方だが、現実問題として、準結晶を扱うよりも、並進対称性を持った構造を扱うことの方が圧倒的に多いので、このように説明している。

Pymatgenでの構造の読み込みかた

cifを読み込んでStructure型のデータを得る

Pymatgenでは、構造を、独自の型である、pymatgen.core.structure.Structure型（以降では単にStructureと書く）で扱う。XRDを計算するにも構造がないと計算できないので、まずは、Structureを簡単に扱えるようになろう。雛形には、以下のマグネタイト（Fe3O4）を用いる。リンク先からsymmerrizedされたcifをダウンロードして、コードと同じフォルダに入れておく。

まずは、下記のコードで構造を読み込んでくる。

#cifファイルの名前はデフォルトにしたとする。変更した場合は、その名前にすること。
path = 'Fe3O4_mp-19306_symmetrized.cif'
parser = CifParser(path)
mat = parser.get_structures(primitive=False, symmetrized=False)[0]
print(type(mat))
#出力：<class 'pymatgen.core.structure.Structure'>

ここのポイントは2つで、
・CifParserのget_structuresはリストを返すこと
・get_structuresの引数
に気を付ける。

1つ目は詳しくは過去の記事(https://qiita.com/ojiya/items/9ba8bfadab96cdddb0b9) に書いたことで、CifParserにcifファイルへのパスを入力（CifParser(path)）し、それを元に構造を取得（parser.get_structures()[0]）する。一般に、cifには複数の構造が入っていることがあるため、リストで取得する構造になっている。そこで、[0]をつけて、最初の構造を取り出す必要がある。
2つ目のポイントは、引数であるprimitiveとsymmetrizeについて。get_structuresは引数を2つとることができて、primitiveは読み込むときにprimitive cellで読み込むか、conventional cellで読み込むか、を指定する。上で述べた、単位格子はconventional cellである（prmitive cellが何か、は気にしなくて良い。）。conventionalで読み込んでおいた方が、後で格子を作るときのパラメータと食い合わせが良いので、primitive=False（デフォルトではTrue）とする。また、symmetrizedは、2022.0.7ではうまく動いていないようなので、Falseにしておく（デフォルトでもFalse）。本当は、非対称単位を取り出すことができる機能のはず。
最後にprint(type(mat))をつけたので、matの型（'pymatgen.core.structure.Structure'、以下、単に「Structure型」と呼ぶ）が表示されたと思う。以降では、このStructure型からの情報の引き出し方を記載していく。

Structure型からのデータの取り出し

単純に、入っている情報を大雑把に引き出したいなら、

#ここでは、上記で読み込んだ'mat'をそのまま使うことにする
dict_mat = mat.as_dict()
print(dict_mat.keys())
print(dict_mat)

のように、.as_dict()を使えばよい。Structure型に入っている情報のうち、主要なものを、辞書型のデータで扱うことができる。この辞書のkeyには'at_module'、'at_class'、'charge'、'lattice'、'sites'（at_はアットマーク）があり、たとえば、以下のようにlattice（結晶格子）の情報を引き出してみると、

mat_lattice = dict_mat['lattice']
print(type(mat_lattice))
print(mat_lattice.keys())
#出力：<class 'dict'>
#dict_keys(['matrix', 'a', 'b', 'c', 'alpha', 'beta', 'gamma', 'volume'])

やはり辞書型のデータで結晶格子の情報が入っている。keyを見れば何が書かれているかは想像がつくと思う。計算科学関連で手を動かしたことがないと、'matrix'はピンとこないかもしれないので補足しておくと、これは格子のa, b, c軸が直交座標系（cartesian座標系）で、どのようなベクトルになっているかを示した行列になっている。「直交座標空間における格子座標空間での基底ベクトル」と言った方が、わかりやすい人にはわかりやすいかもしれない。

さて、大雑把に構造の情報を取り出す方法を書いたが、毎回辞書型に変換するのは面倒なので、頻繁に使うものは、使い方を紹介しておくことにする。

結晶格子

mat_lattice = mat.lattice
print(type(mat_lattice))
print(mat_lattice)
print(mat_lattice.parameters)

実行するとわかるが、.latticeとすることで、pymatgen.core.lattice.Lattice型（以降では単にLattice型）のデータを得ることができる。また、単に「（Lattice型のデータ）を表示せよ」というコードを実行すると、matrixが表示される。この辺は、実験系の人には扱いにくい点ではある。
見慣れた格子状数の状態でデータを得たい場合は、.parametersをつければ、タプル型のデータとして、格子状数を得ることができる。あるいは、

print(mat_lattice.a)
print(mat_lattice.alpha)

のようにして、.aや.alphaというような記法で、各格子状数を取り出すこともできる。もちろん、.b、.c、.beta、.gammaを使えば、他の格子状数も取り出すことができる。

含まれる元素種

ここは微妙にややこしい（前情報なしだと確実にエラーの元になるコードを書いてしまう）ところ。以下には、二種類の方法で元素種の取り出し方を書いた。

mat_spec = mat.species
print(mat_spec)
mat_spec_occ = mat.species_and_occu
print(mat_spec_occ)
# 出力：
# [Element Fe, Element Fe, Element Fe....
# [Comp: Fe1, Comp: Fe1, Comp: Fe1....

.speciesと、.species_and_occという、似たような働きを持った2つの手法があり、いずれも構造中の各サイトにどのような化学種がいるのかを、リスト型で返す。というのはいいとして、リストの中に入っているのは、.speciesではElement型（元素を扱うための型）で、一方、.species_and_occではComposition型（組成式を扱うための型）になっている。これが何でエラーの元になるかと言うと、Element型では、たとえば'Fe0.8In0.2'のようなサイトを扱うことができない。要するに、どこかのサイトが置換されている場合、.speciesを使うと、エラーが出てくる。また、欠損サイトがあるとき、つまり、占有率が1より小さい値を取る場合を表現することができない。どちらがいいい、ということはないので、場合によって使い分ける必要がある。

サイトごとの化学種および座標

続いて、構造中に含まれる化学種と、その座標を取得するには、

mat_sites = mat.sites
print(type(mat_sites))
# 出力：<class 'list'>

のように、（構造データ）.sitesとする。この結果、リスト型が返されて、リストの各要素が構造に含まれるサイトの化学種および座標になる。リスト型なので、個々の要素を見るためには[i]とすればよい。

print(type(mat_sites[0]))
print(mat_sites[0])
# 出力：<class 'pymatgen.core.sites.PeriodicSite'>
# [1.06666894 1.06666894 5.33334471] Fe

とりあえず、最初の要素（[0]）を取り出してみた。出力を見ればわかるように、'pymatgen.core.sites.PeriodicSite'型（以下、PeriodicSite型）のデータになっており、どうやら、そのサイトの化学種と座標のデータが含まれているらしい。
さらに細かくデータを取り出したい時は、

print(mat_sites[0].species)
print(mat_sites[0].coords)
print(mat_sites[0].frac_coords)
# Fe1
# [1.06666894 1.06666894 5.33334471]
# [0.125 0.125 0.625]

のようにする。ここで引っかかりやすいのは、Structure型に対するのとは打って変わって、.speciesを用いると、Composition型を返すところで、非常に混乱しやすい。PeriodicSite型からElement型を取り出したい時は、.specieを使う（sが無い）。また、単に.coordsとすると、cartesian座標（直交座標系での座標）が返ってくるので、fractional座標（格子座標系での座標）が必要な時は、.frac_coordsを用いることにも注意が必要。

対称性の情報

構造から空間群の情報を取り出したい時は、以下のように.get_space_group_info()を使う。

mat_sg = mat.get_space_group_info()
print(mat_sg)
# 出力：('Fd-3m', 227)

タプル型で、1つ目の要素にヘルマン・モーガン記号、2つ目に空間群番号が入っている。ヘルマン・モーガン記号の方は文字列型、空間群番号の方は整数型のデータになっているので、これ以上の情報を、このデータから取り出すことはできない。そこで、これ以上の情報が欲しい場合は、Pymatgenの別の機能を使う。

from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
sga = SpacegroupAnalyzer(mat)
print(sga)
# 出力：<pymatgen.symmetry.analyzer.SpacegroupAnalyzer object at 0x7fe481a1d250>

SpaceGroupAnalyzerを使うには、まずは引数として構造のデータを入れ、それをもとに、色々と対称性に関する情報を引き出す。ここでは、非対称単位と対称操作を取り出すことにする。まずは、非対称単位を取り出す。

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

print(type(eq_sites))
print(type(eq_sites[0]))
print(eq_sites[0])
# 出力：
# <class 'list'>
# <class 'list'>
# [PeriodicSite: Fe (1.0667, 1.0667, 5.3333) [0.1250, 0.1250, 0.6250],
# PeriodicSite: Fe (1.0667, 7.4667, 3.2000) [0.1250, 0.8750, 0.3750], 
# 以下略

.equivalent_sitesの返り値は入れ子のリストで、「『互いに結晶学的に等価なサイトを要素とするリスト』を要素とするリスト」になっている。上記のコードでは、「互いに結晶学的に等価なサイトが入ったリスト」のうち、最初の1つ（[0]）を取り出してきている。このリストの要素の型はPeriodicSiteになっている。続いて、以下の要領で、非対称単位を取り出す。これは、互いに結晶学的に等価なサイトのうち1つずつを取り出せばいいので、以下の通り、内包表記で簡単に取り出せる。今回使っているFe3O4では、非対称単位中に、Feが2サイト、Oが1サイトある。

asym_unit = [site[0] for site in eq_sites]
print(asym_unit)
# 出力：
# [PeriodicSite: Fe (1.0667, 1.0667, 5.3333) [0.1250, 0.1250, 0.6250],
# PeriodicSite: Fe (0.0000, 0.0000, 0.0000) [0.0000, 0.0000, 0.0000],
# PeriodicSite: O (1.0296, 1.0296, 3.2371) [0.1207, 0.1207, 0.3793]]

これで、asym_unitには、matの非対称単位が入っていることになる。
続いて、対称操作を取り出す。

sym_ops = sga.get_symmetry_operations()
print(type(sym_ops))
print(type(sym_ops[0]))
#出力：
# <class 'list'>
# <class 'pymatgen.core.operations.SymmOp'>

.get_symmetry_operationsは、リストを返し、その要素は'pymatgen.core.operations.SymmOp'型（以下、SymmOp型）となっている。さらに、SymmOpの中身はというと、

print(sym_ops[1].rotation_matrix)
print(sym_ops[1].translation_vector)
# 出力：
# <class 'numpy.ndarray'>
# [[ 0. -1.  0.]
#  [ 1.  0.  0.]
#  [ 0.  0.  1.]]
# <class 'numpy.ndarray'>
# [0.25 0.75 0.75]

のように、.rotation_matrixと.translation_vectorで取り出すことができる。見ての通り、回転操作と並進操作を行うための行列が手に入る。なお、sym_ops[0]でなく[1]にしたのは、[0]が恒等操作で面白みが無いため。ちなみに.affine_matrixを使えば、掛け算のみで対称操作を行うことができる行列を手に入れることができるが、ちょっと説明がややこしくなってしまうので、ここでは言及するに留める。
さて、ここまでのコードで非対称単位と対称操作を取り出すことができているので、これらを組み合わせれば、単位格子中の原子を全て生成することができる。座標変換は行列演算になるので、以下のように計算してやればよい。

\mathbf{X'=RX+T}

ただし、Xは元の座標、X'は変換後の座標、Rは回転を表す行列、Tは並進を表すベクトルである。上記でみたような、非対称単位中のFeの片方の座標（0.1250, 0.1250, 0.6250）に対してsym_ops[1]の対称操作を行うとすると、以下の行列を計算することになる。

\begin{pmatrix}
0 & -1 & 0 \\
1 & 0  & 0 \\
0 & 1  & 1 \\
\end{pmatrix}
\begin{pmatrix}
0.125\\
0.125\\
0.625\\
\end{pmatrix}
+
\begin{pmatrix}
0.25\\
0.25\\
0.75\\
\end{pmatrix}

というわけで、まずは以下のように行列を定義する。

op = sym_ops[1]
rot = op.rotation_matrix
trans = op.translation_vector
print(rot)
print(trans)
# 出力：
# [[ 0. -1.  0.]
#  [ 1.  0.  0.]
#  [ 0.  0.  1.]]
# [0.25 0.75 0.75]

それから、numpyの行列演算の機能（numpy.dot）を用いて計算を行うと、以下のごとく、対称操作によって生成された座標が得られる。

frac_coord = asym_unit[0].frac_coords
print(frac_coord)

sym_generated_raw = np.dot(rot, frac_coord) + trans
print(sym_generated_raw)
＃出力：
# [0.125 0.125 0.625]
# [0.125 0.875 1.375]

ただし、ここで得られた座標は、fractionalなのに座標が1を越えており、お行儀が良くないので、これを修正したい。fractionalでの座標を0~1に収めるには、「その値自体を越えない最大の整数」で引けば良く、これはpythonモジュール「math」の機能であるfloorが、まさに「その値自体を越えない最大の整数」を得るための機能になっている。というわけで、以下のようにすればよい。

import math
sym_generated = [coord - math.floor(coord) for coord in sym_generated_raw]
print(sym_generated)
#出力：
# [0.125, 0.875, 0.375]

これで、非対称ユニット内のあるサイトに対称性を加えて、セル内の別のサイトを作ることができるようになった。これで、対称操作のリスト（sym_ops）と、非対称ユニット内のサイトのリスト（asym_unit）に関してfor文を回せば、セル内のサイトは全て生成することができる。ただし、生成するサイトと、対称操作は、一対一に対応しているわけではない。複数の操作が、空間的に同じ位置のサイトを生成することがある。たとえば、座標(0.5, 0.5, 0.5)のサイトに対して原点を中心とする180°の回転操作を行うと、(-0.5, -0.5, -0.5)となり、元の位置に戻ってくることになる。そこで、for文で回してやった後に、重複する要素を除外する、という作業が必要になる。ただし、この記事でそれを書き始めると、ちょっと長くなりすぎるので、今後の記事に記載することにする。

以上、本記事では、結晶構造と対称性について、結晶をあまり扱ったことない人向けの説明に加え、Pymatgenを用いた結晶構造の読み込み方を記し、読み込んだデータから、格子の情報、サイトの情報、対称性の情報を引き出す方法について説明した。次の記事では、Pymatgenを用いてStructure型のデータを作る方法について、記す。