銅結晶の原子配置データ解析チュートリアル¶
このチュートリアルでは、銅（Cu）の結晶構造の原子の配置データを解析します。
データのダウンロードと読み込み¶
ここでは、面心立方（FCC）構造を持つ銅（Cu）の結晶データを解析します。
データは Materials Project から取得します。
銅のFCC結晶には mp-30 というIDが割り当てられており、以下のページで詳細情報を確認できます： https://materialsproject.org/materials/mp-30
このデータは CC-BY ライセンス のもとで提供されており、比較的自由に利用可能です。
ライブラリを読み込みます。
import ase
import ase.io
import homcloud.interface as hc
import numpy as np
import py3Dmol
import io
import homcloud.py3dmolhelper as py3dmolhelper
銅結晶の構造データは、同じディレクトリに保存された CIF 形式のファイルを使用します。
このファイルは、ase.io.read 関数を使って読み込みます。
cu_atoms = ase.io.read("Cu.cif")
構造データを読み込んだ後、原子の位置情報を確認することで、この結晶構造が4つの原子から構成されていることがわかります。
cu_atoms.get_positions()
array([[0. , 0. , 0. ],
[0. , 1.78871534, 1.78871534],
[1.78871534, 0. , 1.78871534],
[1.78871534, 1.78871534, 0. ]])
HomCloudによるパーシステントホモロジー解析を行うには、単位格子内の4原子だけでは小さすぎて解析に適しません。 そのため、結晶構造をXYZ方向にそれぞれ4回ずつ繰り返したスーパーセルを作成して使用します。
cu_atoms_4 = cu_atoms.repeat(4)
作成したスーパーセル構造を、py3Dmol を使って可視化してみましょう。 py3Dmol は xyz フォーマットのデータを読み込む必要があるため、まず ASE の構造データを xyz 形式に変換します。
ここでは、Pythonの StringIO クラスを使って、ファイルを保存せずに直接文字列としてデータを渡します。
# XYZフォーマットへの変換
s = io.StringIO()
ase.io.write(s, cu_atoms_4, format="xyz")
xyz_data = s.getvalue()
# ビューアの準備
view = py3Dmol.view()
# ビューアに原子配置の情報を追加
view.addModel(xyz_data, "xyz")
# 原子データの描画の設定をする
view.setStyle({'stick': {"colorscheme": "default"}})
view.addStyle({"sphere": {"scale": 0.3}})
# 以下で表示
view.zoomTo()
view.show()
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
get_positionsメソッドで原子の座標を取り出します。
atoms_positions = cu_atoms_4.get_positions()
周期構造による問題とノイズの導入¶
HomCloudによるトポロジー解析では、結晶のような周期的構造が原因で、パーシステンス図が正しく計算されない場合があります。特に以下のような幾何的配置が問題となります：
- 3点が同一直線上にある
- 4点が同一平面上にある
- 5点が同一球面上にある
この問題を回避するために、構造に微小なノイズを加えます。
銅原子間の距離は約 2.53 Å であるため、ノイズの大きさを $10^{-4}$ Å 程度に設定すれば、構造の物理的意味を損なうことなく、解析結果への影響もほとんどありません。
atoms_positions += np.random.uniform(-1e-4, 1e-4,size=atoms_positions.shape)
PDを計算して1次元PD，2次元PDを表示します． X軸，Y軸の単位はオングストロームの二乗なのでそれも表示します。
hc.PDList.from_alpha_filtration(atoms_positions, save_boundary_map=True, save_to="cu.pdgm")
PDList(path=cu.pdgm)
pdlist = hc.PDList("cu.pdgm")
pd1 = pdlist.dth_diagram(1)
pd1.histogram((1, 4), 32).plot(colorbar={"type": "log"}, title="1D", unit_name="$\\AA^2$")
pd2 = pdlist.dth_diagram(2)
pd2.histogram((1, 4), 32).plot(colorbar={"type": "log"}, title="2D", unit_name="$\\AA^2$")
1次元PDの解析¶
1次元のパーシステンス図では、対角線から大きく離れた birth-death ペアがすべて (1.6, 2.1) 付近に集中しています。
そこで、この特徴的なペアを抽出して詳しく解析します。
pairs = pd1.pairs_in_rectangle(1.5, 1.8, 2.0, 2.2)
print(len(pairs))
print(pairs[0])
921 Pair(1.5995856591285957, 2.1327148846008974)
ペアの総数は921個であることがわかりました。 この中からランダムに10個を選び、それぞれの stable volume を計算してみましょう。
import random
stable_volumes = [
pair.stable_volume(1e-3)
for pair in random.sample(pairs, 10)
]
計算結果を可視化します。
# ビューアの準備
view = py3Dmol.view()
# ビューアに原子配置の情報を追加
view.addModel(xyz_data, "xyz")
# 原子データの描画の設定をする
view.setStyle({'stick': {"colorscheme": "default"}})
view.addStyle({"sphere": {"scale": 0.3}})
# stable volumeの情報をビューアに追加
for stable_volume in stable_volumes:
py3dmolhelper.add_edges(view, stable_volume.boundary(), "green", 0.4)
# 以下表示
view.zoomTo()
view.show()
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
すると、これらのペアに対応する構造が、面心立方格子における正三角形の配置に対応していることが見えてきます。 このことを確認するために、stable volume で得られたリング構造の各頂点間の距離を計算してみましょう。 距離の計算には、scipy の distance_matrix を使用して距離行列を求めます。
triangle = stable_volumes[0].boundary_points()
print(triangle)
[[5.36620147136236, 2.779296467901551e-05, 8.94358595883734], [3.5775194207006744, 1.3438209613698864e-05, 10.732254114357774], [5.3661543887519185, 1.7887741467357148, 10.732294533350922]]
import scipy as sp
sp.spatial.distance_matrix(triangle, triangle)
array([[0. , 2.52956859, 2.52964264],
[2.52956859, 0. , 2.52960074],
[2.52964264, 2.52960074, 0. ]])
すると、頂点間の距離は約2.529…であり、すべて等しいことがわかります。この値は、既知の原子間距離（2.53 Å）とも一致しています。
三角形を構成する原子だけを抽出してみましょう。boundary_pointsを使えば、頂点の座標を取得できますが、それ以外の情報を取得する方法についても説明します。
boundary_pointsに by="vindexes" を指定すると、各頂点が何番目の原子に対応しているかを調べることができます。なお、番号は0から始まる点に注意してください。
stable_volumes[0].boundary_points(by="vindexes")
[74, 76, 79]
この番号の配列を ASE の Atoms オブジェクトに渡すことで、該当する原子のみを抽出できます。get_positions や get_chemical_symbols を使えば、原子の位置や種類を取得することが可能です。今回はすべての原子が Cu なのであまり面白みはありませんが、より複雑な物質の場合には非常に有用でしょう。
triangle_atoms = cu_atoms_4[stable_volumes[0].boundary_points(by="vindexes")]
triangle_atoms.get_positions()
array([[ 5.36614601, 0. , 8.94357667],
[ 3.57743067, 0. , 10.73229201],
[ 5.36614601, 1.78871534, 10.73229201]])
triangle_atoms.get_chemical_symbols()
['Cu', 'Cu', 'Cu']
get_all_distances を使うと、原子間の距離を計算できます。これは sp.spatial.distance_matrix と同じ結果を返しますが、ASE に慣れている場合はこちらの方が便利かもしれません。
triangle_atoms.get_all_distances()
array([[0. , 2.52962549, 2.52962549],
[2.52962549, 0. , 2.52962549],
[2.52962549, 2.52962549, 0. ]])
三角形のデータを XYZ フォーマットに変換し、ファイルとして保存してみましょう。以下のように実行します。
ase.io.write("cu_triangle_atoms.xyz", triangle_atoms)
cu_triangle_atoms.xyz というファイルに保存されたので、VESTA など、お好みのソフトウェアを使って解析や可視化を行うことができます。これから、XYZファイルの中身を読み込んで表示してみましょう。
with open("cu_triangle_atoms.xyz") as f:
cu_triangle_atoms_xyz_data = f.read()
print(cu_triangle_atoms_xyz_data)
3
Lattice="14.30972268 0.0 0.0 0.0 14.30972268 0.0 0.0 0.0 14.30972268" Properties=species:S:1:pos:R:3:spacegroup_kinds:I:1 spacegroup="P 1" unit_cell=conventional occupancy="_JSON {\"0\": {\"Cu\": 1}, \"1\": {\"Cu\": 1}, \"2\": {\"Cu\": 1}, \"3\": {\"Cu\": 1}}" pbc="T T T"
Cu 5.36614601 0.00000000 8.94357667 2
Cu 3.57743067 0.00000000 10.73229201 0
Cu 5.36614601 1.78871534 10.73229201 3
このデータをpy3Dmolに渡して三角形だけを描画します。
view = py3Dmol.view()
view.addModel(cu_triangle_atoms_xyz_data, "xyz")
view.setStyle({'stick': {"colorscheme": "default"}})
view.addStyle({"sphere": {"scale": 0.3}})
view.zoomTo()
view.show()
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
2次元PDの解析¶
次に、2次元のパーシステンス図（PD）を解析します。
この図では、(2.1, 2.4)付近と (2.1, 3.2)付近の2か所にペアが集中していることが確認できます。
そこで、まずは (2.1, 2.4) 付近に位置するペアを抽出してみましょう。
pairs = pd2.pairs_in_rectangle(2.1, 2.3, 2.3, 2.5)
print(len(pairs)) # Print the number of pairs
print(pairs[0]) # Print the birth and death times
343 Pair(2.132815806696218, 2.3993756154712207)
1次元PDの場合と同様に、ランダムに10個のペアを抽出し、それぞれの stable volume を計算・可視化してみましょう。
# Stable volumes are computed for ten random pairs
stable_volumes = [pair.stable_volume(1e-3) for pair in random.sample(pairs, 10)]
# 3D Visualization
view = py3Dmol.view()
view.addModel(xyz_data, "xyz")
view.setStyle({'stick': {"colorscheme": "default", "opacity": 0.9}}) # 棒の部分を半透明にする
view.addStyle({"sphere": {"scale": 0.3}})
# stable volumeの情報を青色でビューアに追加
for stable_volume in stable_volumes:
py3dmolhelper.add_surface(view, stable_volume.boundary(), "blue")
# 以下表示
view.zoomTo()
view.show()
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
すると、構造が正四面体になっていることがわかります。三角形の場合と同様に、四面体のデータを XYZ 形式に変換し、その部分だけを描画します。 ここでは、以下の2つの工夫も取り入れます：
- 関数化することで再利用しやすくする
- io.StringIO を使って、よりスマートに処理する
# 関数を定義する
def show_atoms(atoms):
f = io.StringIO()
ase.io.write(f, atoms, format="xyz")
view = py3Dmol.view()
view.addModel(f.getvalue(), "xyz")
view.setStyle({'stick': {"colorscheme": "default"}})
view.addStyle({"sphere": {"scale": 0.3}})
view.zoomTo()
view.show()
show_atoms(cu_atoms_4[stable_volumes[0].boundary_points(by="vindexes")])
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
次は、(2.1, 3.2)付近にあるペアを調べてみましょう。
pairs = pd2.pairs_in_rectangle(2.0, 2.2, 3.1, 3.3)
print(len(pairs)) # Print the number of pairs
print(pairs[0]) # Print the birth and death times
108 Pair(2.1330129084253726, 3.199156507490444)
stable_volumes = [pair.stable_volume(1e-3) for pair in random.sample(pairs, 10)]
# 3D Visualization
view = py3Dmol.view()
view.addModel(xyz_data, "xyz")
view.setStyle({'stick': {"colorscheme": "default", "opacity": 0.9}}) # 棒の部分を半透明にする
view.addStyle({"sphere": {"scale": 0.3}})
# stable volumeの情報を青色でビューアに追加
for stable_volume in stable_volumes:
py3dmolhelper.add_surface(view, stable_volume.boundary(), "blue")
# 以下表示
view.zoomTo()
view.show()
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
すると，これらのペアに対応する構造は正八面体のようであることが確認できます。
1つだけ取り出すと、次のようになります。先ほど定義した関数を使用します。
show_atoms(cu_atoms_4[stable_volumes[0].boundary_points(by="vindexes")])
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
つまり、2次元のパーシステント図（PD）に現れる2種類のbirth-deathペアは、それぞれ正四面体と正八面体の空隙に対応していることがわかります。
実際、面心立方格子（FCC格子）には、2種類の空隙、すなわち四面体サイトと八面体サイトが存在し、その比率は2:1であることが知られています。したがって、PDに現れる2種類のペアもこの構造に対応していると考えられます。
しかし、実際に得られたペアの個数比は343対108で、約3:1となっており、理論値とずれがあります。これはなぜでしょうか？
このずれは、データの端の部分で格子の繰り返し構造が考慮されていないことに起因しています。
この問題に対処する方法として、以下の2つが考えられます：
周期境界条件を用いてPDを計算する： これは理想的な解決策ですが、単位格子が直方体である場合にのみ適用可能です。また、周期境界条件を用いると、PHTreesの機能が使用できなくなることがあります。
格子の繰り返し回数を増やす：こちらはより一般的な方法ですが、繰り返し回数を増やすことで計算コストが上昇します。
格子の情報を調べてみましょう．
cu_atoms_4.cell.angles()
array([90., 90., 90.])
格子ベクトルのなす各がすべて90度なので単位胞は直方体です．そのため一つ目の方策を使いましょう。
a, b, c = cu_atoms_4.cell.lengths()
hc.PDList.from_alpha_filtration(atoms_positions, periodicity=((0, a), (0, b), (0, c)), save_to="cu_8_pbc.pdgm")
pd2 = hc.PDList("cu_8_pbc.pdgm").dth_diagram(2)
len(pd2.pairs_in_rectangle(2.1, 2.3, 2.3, 2.5)) / len(pd2.pairs_in_rectangle(2.0, 2.2, 3.1, 3.3))
2.007843137254902
少しずれていますが，ほぼ 2:1 の比になっています．
以上でこのチュートリアルは終わりです。 次はもっと複雑な物質であるアモルファスSiO2を解析していきましょう。
演習問題¶
銅結晶のPDで、対角線付近にあるbirth-deathペアを解析してみてください。