銅の結晶の原子配置データの解析¶

このチュートリアルでは銅の結晶の原子配置データを解析します．

準備¶

解析には ASE，py3Dmolが必要です．次のようにして pip でインストールします．

pip3 install ase py3Dmol

データのダウンロードと読み込み¶

面心立方格子を持つ銅の結晶を解析しましょう． Materials Projectからダウンロードしたデータを用います．

銅の面心立方格子結晶はmp-30というidが付いており， https://materialsproject.org/materials/mp-30 に詳しい情報があります．

このデータの使用条件については https://materialsproject.org/about/terms を見てください． CC-BY なので比較的自由に利用可能なデータです．

ライブラリを読み込みます．

In [1]:
import ase
import ase.io
import homcloud.interface as hc
import numpy as np
import py3Dmol
import io
import homcloud.py3dmolhelper as py3dmolhelper

ase.io.readでCIFファイルを読み込みます．ファイルは同じディレクトリに置いてあります．

In [2]:
cu_atoms = ase.io.read("Cu.cif")

次のようにして原子の位置を見ると，4つの原子からなることがわかります．

In [3]:
cu_atoms.get_positions()
Out[3]:
array([[0.        , 0.        , 0.        ],
       [0.        , 1.78871534, 1.78871534],
       [1.78871534, 0.        , 1.78871534],
       [1.78871534, 1.78871534, 0.        ]])

HomCLoudによる解析は4原子だと小さすぎて都合が悪いので，XYZ方向に4回ずつ繰り返したデータを使うことにします．

In [4]:
cu_atoms_4 = cu_atoms.repeat(4)

このデータをpy3Dmolで可視化しましょう．py3Domlはxyzフォーマットに変換して読み込む必要があります．ここではデータの変換にPythonのStringIOクラスを利用します．

In [5]:
# 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メソッドで座標を取り出します．

In [6]:
atoms_positions = cu_atoms_4.get_positions()

これはHomCloudによる解析において，結晶のような整然とした構造は不都合です． 以下のような構造がある場合には，計算結果がおかしくなります．

  1. 3点が同一直線上にある
  2. 4点が同一平面上にある
  3. 5点が同一球面上にある

そこで小さなノイズを加えてこの問題を解決します． 原子間の距離は2.53Åなので$10^{-4}$程度の大きさのノイズにしておけば解析結果には影響を与えません．

In [7]:
atoms_positions += np.random.uniform(-1e-4, 1e-4,size=atoms_positions.shape)

PDを計算して1次元PD，2次元PDを表示します．

In [8]:
hc.PDList.from_alpha_filtration(atoms_positions, save_boundary_map=True, save_to="cu.pdgm")
Out[8]:
PDList(path=cu.pdgm)
In [9]:
pdlist = hc.PDList("cu.pdgm")
pd1 = pdlist.dth_diagram(1)
pd1.histogram((1, 4), 32).plot(colorbar={"type": "log"}, title="1D")
pd2 = pdlist.dth_diagram(2)
pd2.histogram((1, 4), 32).plot(colorbar={"type": "log"}, title="2D")
No description has been provided for this image
No description has been provided for this image

1次元PDの解析¶

1次元PDには対角線から離れたbirth-death pairは(1.6, 2.1)あたりに全て集中しています．そこでこのペアを取り出します．

In [10]:
pairs = pd1.pairs_in_rectangle(1.5, 1.8, 2.0, 2.2)
In [11]:
print(len(pairs))
print(pairs[0])

921
Pair(1.599613313402053, 2.132723280439032)

ペアの個数は921ということがわかります．ここから10個ランダムに取ってきてstable volumeを計算してみましょう．

In [12]:
import random
stable_volumes = [
  pair.stable_volume(1e-3)
  for pair in random.sample(pairs, 10)
]

計算結果を可視化します．

In [13]:
# ビューアの準備
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 を使って距離行列を計算します．

In [14]:
triangle = stable_volumes[0].boundary_points()
print(triangle)

[[10.732324543027627, 8.943571738108423, 1.7886897402816557], [10.732274710368953, 10.732239693008147, 3.5775216970677057], [12.520973637710009, 8.943581432388873, 3.5774936901518357]]
In [15]:
import scipy as sp
sp.spatial.distance_matrix(triangle, triangle)
Out[15]:
array([[0.        , 2.52967445, 2.52964131],
       [2.52967445, 0.        , 2.52957353],
       [2.52964131, 2.52957353, 0.        ]])

すると頂点間の距離は2.529…で等しいことがわかります． この値は既知の原子間距離(2.53Å)とも一致しています．

三角形をなす原子だけを取り出してみましょう。boundary_pointsを使えば頂点の座標は取り出せますが，それ以外の情報を取り出す方法も説明します。 boundary_pointsにby="vindexesを渡すと何番目の原子なのかを調べることができます。番号は0からスタートすることに注意してください。

In [16]:
stable_volumes[0].boundary_points(by="vindexes")
Out[16]:
[225, 244, 231]

この番号の配列をASEのAtomsオブジェクトに渡すと，その原子だけを取り出すことができます。get_positionsやget_chemical_symbolsで原子の位置や種類を取り出すことができます。今回は原子は全てCuなのであまり面白くありませんが，もっと複雑な物質であれば役立つでしょう。

In [17]:
triangle_atoms = cu_atoms_4[stable_volumes[0].boundary_points(by="vindexes")]
In [18]:
triangle_atoms.get_positions()
Out[18]:
array([[10.73229201,  8.94357667,  1.78871534],
       [10.73229201, 10.73229201,  3.57743067],
       [12.52100735,  8.94357667,  3.57743067]])
In [19]:
triangle_atoms.get_chemical_symbols()
Out[19]:
['Cu', 'Cu', 'Cu']

get_all_distancesを使うと原子間距離が計算できます。sp.spatial.distance_matrixと結果は同じです。ASEに慣れている場合はこちらのほうが便利かもしれません。

In [20]:
triangle_atoms.get_all_distances()
Out[20]:
array([[0.        , 2.52962549, 2.52962549],
       [2.52962549, 0.        , 2.52962549],
       [2.52962549, 2.52962549, 0.        ]])

三角形のデータをXYZフォーマットに変換してファイルに保存しましょう。以下のようにします。

In [21]:
ase.io.write("cu_triangle_atoms.xyz", triangle_atoms)

cu_triangle_atoms.xyzというファイルに保存されたので，VESTAなど自分の好きなソフトウェアで解析や可視化することもできます。XYZのファイルの中身を読み込んで表示します。

In [22]:
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      10.73229201       8.94357667       1.78871534        1
Cu      10.73229201      10.73229201       3.57743067        0
Cu      12.52100735       8.94357667       3.57743067        3

以下ではこのデータをpy3Dmolに渡して三角形だけを描画しましょう。

In [23]:
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)あたりのペアを取り出します．

In [24]:
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.1328264661314558, 2.3993758240272927)

1次元PDと同様，ランダムに10個取り出してstable volumeを計算，可視化しましょう．

In [25]:
# 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を使ってちょっとかっこよくやる
In [26]:
# 関数を定義する
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()
In [27]:
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)あたりのペアを調べましょう．

In [28]:
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.132953930258449, 3.1991644960714676)
In [29]:
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.

すると今度は正八面体のようであることが見てとれます．

一個だけ取り出すと次のようになります。さきほど定義した関数を使います。

In [30]:
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 pairは正四面体と正八面体に対応しているということがわかります．

実は面心立方格子というのは2種類の隙間(四面体サイトと八面体サイトと呼ばれる)が2:1の比率であることが知られています． するとこれに対応していると考えられます．

しかし2種類のペアの比率は 343 / 108 でおよそ3:1でちょっと個数にずれがあります．これはどういうことでしょうか．

実はこれはデータの端の部分で繰り返し構造が考慮されていないせいです．

これに対応する方法は2つあります．

  • 周期境界条件でPDを計算する
  • 繰り返し回数を増やす

一つ目は理想的解決策ですが，単位格子が直方体の場合にしか適用できません．また周期境界ではPHTrees機能が使えなくなったりします．二つめは繰り返し回数を増やす必要があって計算コストが上昇します．

とりあえず格子の情報を調べてみましょう．

In [31]:
cu_atoms_4.cell.angles()
Out[31]:
array([90., 90., 90.])

格子ベクトルのなす各がすべて90度なので直方体です．そのため一つ目の方策を使いましょう。

In [32]:
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))
Out[32]:
2.007843137254902

少しずれていますが，ほぼ 2:1 の比になっています．

このずれは，正八面体のほうの個数が1個少ないことによります。以下のようにして分母に1を足してやるとちょうど2になります。

In [33]:
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)) + 1)
Out[33]:
2.0

正八面体のほうの個数が1個少ない理由は，3次元周期境界条件のもとで空洞の個数を数えることに関する数学的な問題に由来します。

以上でこのチュートリアルは終わりです。 次はもっと複雑な物質であるアモルファスSiO2を解析していきましょう。

演習問題¶

銅結晶のPDで、対角線付近にあるbirth-deathペアを解析してみてください。