原子配置データの解析¶
このチュートリアルでは銅の結晶とアモルファスシリカの原子配置データを解析します．最初は銅のデータです．
銅の結晶データの解析¶
まずは面心立方格子を持つ銅の結晶を解析しましょう． Materials Projectからダウンロードしたデータを用います．
銅の面心立方格子結晶はmp-30というidが付いており，
https://materialsproject.org/materials/mp-30 に詳しい情報があります．
このデータの使用条件については https://materialsproject.org/about/terms を見てください． 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
ase.io.readでCIFファイルを読み込みます．ファイルは同じディレクトリに置いてあります．
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で可視化しましょう．py3Domlは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を表示します．
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")
pd2 = pdlist.dth_diagram(2)
pd2.histogram((1, 4), 32).plot(colorbar={"type": "log"}, title="2D")
1次元PDの解析¶
1次元PDには対角線から離れたbirth-death pairは(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.5995711542379136, 2.132745949527214)
ペアの個数は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)
[[3.577337367319096, 1.7887448036458313, 8.943594348243618], [3.5775067714130175, 3.577388081868643, 10.732353815511095], [5.3661249218916005, 1.78879255860797, 10.732234592479681]]
import scipy as sp
sp.spatial.distance_matrix(triangle, triangle)
array([[0. , 2.52960575, 2.52962346],
[2.52960575, 0. , 2.52947205],
[2.52962346, 2.52947205, 0. ]])
すると頂点間の距離は2.529…で等しいことがわかります． この値は既知の原子間距離(2.53Å)とも一致しています．
三角形をなす原子だけを取り出してみましょう。boundary_pointsを使えば頂点の座標は取り出せますが，それ以外の情報を取り出す方法も説明します。
boundary_pointsにby="vindexesを渡すと何番目の原子なのかを調べることができます。番号は0からスタートすることに注意してください。
stable_volumes[0].boundary_points(by="vindexes")
[73, 92, 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([[ 3.57743067, 1.78871534, 8.94357667],
[ 3.57743067, 3.57743067, 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 3.57743067 1.78871534 8.94357667 1
Cu 3.57743067 3.57743067 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.1328447328590237, 2.3993872246457144)
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.1329399796726025, 3.19916509445435)
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.
すると今度は正八面体のようであることが見てとれます．
一個だけ取り出すと次のようになります。さきほど定義した関数を使います。
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機能が使えなくなったりします．二つめは繰り返し回数を増やす必要があって計算コストが上昇します．
とりあえず格子の情報を調べてみましょう．
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 の比になっています．
このずれは，正八面体のほうの個数が1個少ないことによります。以下のようにして分母に1を足してやるとちょうど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)) + 1)
2.0
正八面体のほうの個数が1個少ない理由は，3次元周期境界条件のもとで空洞の個数を数えることに関する数学的な問題に由来します。
アモルファスシリカの原子配置¶
次にアモルファスシリカの原子配置をPHで解析していきましょう． データは https://isaacs.sourceforge.net/ex.html からダウンロードできるものを使います．
データのダウンロード¶
次のセルを実行するとデータがダウンロードできます．
import urllib
data = urllib.request.urlopen("https://isaacs.sourceforge.net/tests/sio2-cart.chem3d").read()
with open("sio2-cart.chem3d", mode="wb") as f:
f.write(data)
ファイルフォーマットの変換¶
ASEやpy3dMolはchem3dというファイルフォーマットを取り扱えないので，xyzに変換します．
import re
with open("sio2-cart.chem3d") as infile:
with open("sio2-cart.xyz", "w") as outfile:
outfile.write(infile.readline())
outfile.write("SiO2\n")
for line in infile:
atom, index, x, y, z = re.split(r"\s+", line.strip())
outfile.write(f"{atom} {x} {y} {z}\n")
ファイルの読み込み¶
変換後のデータをASEで読み込みます
sio2 = ase.io.read("sio2-cart.xyz")
可視化¶
py3DMolでアモルファスシリカの立体構造を可視化します．
with open("sio2-cart.xyz") as f:
sio2_xyz = f.read()
view = py3Dmol.view()
view.addModel(sio2_xyz, "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.
sio2
Atoms(symbols='O2000Si1000', pbc=False)
PDの計算¶
PDを計算します．ここで，sio2.get_chemical_symbols()とすると原子の名前('O'や'Si')が得られるので，これを使って各点に原子の種類の情報を覚えさせておきます．
hc.PDList.from_alpha_filtration(sio2.get_positions(),
vertex_symbols=sio2.get_chemical_symbols(),
save_boundary_map=True,
save_to="sio2.pdgm")
PDList(path=sio2.pdgm)
1次元PDの解析¶
銅の場合と同様に、1次元PDを計算・出力してみましょう。
pd1 = hc.PDList("sio2.pdgm").dth_diagram(1)
pd1.histogram().plot(colorbar={"type": "log"})
拡大してみます．
pd1.histogram((0, 8), 256).plot(colorbar={"type": "log", "max": 10})
birth が 0.7 あたりにペアが垂直に集まっています．ここから1点取って調べてみましょう．
pair = pd1.nearest_pair_to(0.7, 4.0)
pair
Pair(0.6844174999999998, 3.9994768570248276)
stable_volume = pair.stable_volume(0.01)
以下の点からなるリング構造が検出されました。
stable_volume.boundary_points()
[[-9.603, -9.559, 1.597], [-10.398, -14.357, 1.84], [-8.603, -5.624, 0.501], [-5.639, -4.191, 7.279], [-7.506, -3.719, 4.804], [-10.06, -8.238, 2.43], [-7.814, -8.827, -0.738], [-5.83, -5.47, 9.373], [-8.239, -12.941, -1.155], [-7.0, -5.139, 5.33], [-8.946, -7.346, 7.504], [-7.953, -4.238, 1.024], [-7.123, -8.344, 0.673], [-7.385, -7.487, 9.515], [-7.628, -6.841, 0.923], [-8.222, -3.702, 3.441], [-6.028, -6.958, 8.796], [-10.073, -11.659, 3.81], [-8.811, -7.828, 9.084], [-5.895, -3.961, 8.813], [-9.057, -8.223, 6.139], [-9.508, -10.811, 2.543], [-10.167, -13.282, 4.091], [-9.609, -7.204, 3.434], [-8.925, -7.408, 4.76], [-9.633, -13.041, -0.301], [-9.838, -14.651, 3.234], [-10.892, -13.612, 0.54], [-8.686, -10.033, -1.214], [-7.848, -11.389, -1.05], [-5.691, -5.448, 6.213], [-7.933, -3.018, 2.003]]
boundary_points_symbolsを使うと，さきほど点に与えた原子名の情報を取り出すことができます．
stable_volume.boundary_points_symbols()
['Si', 'O', 'Si', 'O', 'Si', 'O', 'O', 'O', 'Si', 'O', 'O', 'O', 'Si', 'O', 'O', 'O', 'Si', 'Si', 'Si', 'Si', 'Si', 'O', 'O', 'Si', 'O', 'O', 'Si', 'Si', 'Si', 'O', 'Si', 'Si']
Si原子とO原子が混ざったリング構造であることがわかります．
このリング上での原子の並びを調べましょう。 boundary_loop_symbolsで得られます。
実はこのリング構造は（Si-O-Si-O-... という具合に）Si原子とO原子が交互に繋がってできています。
stable_volume.boundary_loop_symbols()
['Si', 'O', 'Si', 'O', 'Si', 'O', 'Si', 'O', 'Si', 'O', 'Si', 'O', 'Si', 'O', 'Si', 'O', 'Si', 'O', 'Si', 'O', 'Si', 'O', 'Si', 'O', 'Si', 'O', 'Si', 'O', 'Si', 'O', 'Si', 'O']
リング構造を可視化してみましょう．
view = py3Dmol.view()
view.addModel(sio2_xyz, "xyz")
view.setStyle({'stick': {"colorscheme": "default", "opacity": 0.8}}) # 棒の部分は半透明で
view.addStyle({"sphere": {"scale": 0.3, "opacity": 0.8}}) # ボールの部分も半透明で
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.
これはSiとOの間の化学結合によって実現されているリングであるようです．他にもこの縦の分布の点をいくつか取り出して同様に調べてみると良いでしょう．
show_atoms(sio2[stable_volume.boundary_loop(by="vindexes")])
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
次にペアが集中している 2 ≦ birth, death ≦ 3 の範囲を調べてみます．
pd1.histogram((2, 3), 32).plot(colorbar={"type": "log", "max": 10})
この集中している部分は何でしょうか？ (2.6, 2.64)の所からペアを一つ取り出します．
pair = pd1.nearest_pair_to(2.7, 2.75)
stable_volume = pair.stable_volume(0.005)
stable_volume.boundary_points()
[[-4.716, -12.817, -0.02], [-3.266, -12.775, -0.712], [-4.159, -9.9, -1.433]]
stable_volume.boundary_points_symbols()
['O', 'Si', 'O']
Si原子1つとO原子2つからなる三角形のようです．
view = py3Dmol.view()
view.addModel(sio2_xyz, "xyz")
view.setStyle({'stick': {"colorscheme": "default", "opacity": 0.8}}) # 棒の部分は半透明で
view.addStyle({"sphere": {"scale": 0.3, "opacity": 0.8}}) # ボールの部分も半透明で
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.
show_atoms(sio2[stable_volume.boundary_loop(by="vindexes")])
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
O-Siという結合を持つペアと，それとは少し離れたOが抽出されていることがわかります。
OとSiの間にしか化学結合がないにもかかわらずリング構造が抽出されている，というのは興味深い事実です． PDの計算では化学結合の情報は無視して半径を膨らませていことでリング構造を抽出するわけなので 別におかしな話ではないのですが，結合情報を見ているだけでは見えてこない構造が見えてくるのは PH解析の利点の一つと言えるでしょう．
ここの分布を見てみると，実は色々な構造が混ざっていることがわかります。より詳しく調べたければ逆解析結果に対して統計的処理をしていくとよいでしょう。
演習問題¶
- 銅結晶のPDで、対角線付近にあるbirth-deathペアを解析してみてください。
- アモルファスシリカの例題で，"Si0001", "Si0002", ...というように原子を番号付けしてみてください。
- お好きな原子配置データをPH解析してみてください。
発展的話題¶
原子の初期半径を決める¶
原子の種類をPH解析に反映させるため，各原子の初期半径を指定することができます．次の例ではOの半径を1.275，Siの半径を0.375として解析します． この半径は https://www.pnas.org/content/113/26/7035.full で使われているもので，アモルファスシリカのPH解析に実際に使われた値です． 初期半径の決めかたについてはリンク先の論文を読んください．
また，半径の値は二乗する必要があることに注意してください．
weights = np.array([1.275**2 if atom == "O" else 0.375**2 for atom in sio2.get_chemical_symbols()])
pd1 = hc.PDList.from_alpha_filtration(
sio2.get_positions(),
weight=weights,
vertex_symbols=sio2.get_chemical_symbols(),
save_boundary_map=True
).dth_diagram(1)
pd1.histogram().plot(colorbar={"type": "log"})
pd1.histogram((-0.2, 2), 256).plot(colorbar={"type": "log"})
pair = pd1.nearest_pair_to(0.00, 0.2)
stable_volume = pair.stable_volume(0.005)
stable_volume.boundary_points_symbols()
['O', 'O', 'Si']
view = py3Dmol.view()
view.addModel(sio2_xyz, "xyz")
view.setStyle({'stick': {"colorscheme": "default", "opacity": 0.8}}) # 棒の部分は半透明で
view.addStyle({"sphere": {"scale": 0.3, "opacity": 0.8}}) # ボールの部分も半透明で
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.
show_atoms(sio2[stable_volume.boundary_loop(by="vindexes")])
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
左下のペアが集中している部分はO-Si-Oが作る三角形に対応しているようです。
pair = pd1.nearest_pair_to(0.71, 0.85)
stable_volume = pair.stable_volume(0.005)
stable_volume.boundary_points_symbols()
['O', 'O', 'O']
view = py3Dmol.view()
view.addModel(sio2_xyz, "xyz")
view.setStyle({'stick': {"colorscheme": "default", "opacity": 0.8}}) # 棒の部分は半透明で
view.addStyle({"sphere": {"scale": 0.3, "opacity": 0.8}}) # ボールの部分も半透明で
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.
show_atoms(sio2[stable_volume.boundary_loop(by="vindexes")])
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
斜めに延びる分布は3つのO原子が作る三角形に対応しているようです．
birth = 0.0 の所で縦に伸びているペアはSi-O-Si-…という結合によるループなのですが，逆解析を適用してみると少し違った様子が見えます。
pair = pd1.nearest_pair_to(0.00, 1.50)
stable_volume = pair.stable_volume(0.005)
stable_volume.boundary_points_symbols()
['O', 'Si', 'O', 'Si', 'O', 'O', 'O', 'O', 'O']
O-Si-Oと並んでいる部分がありますが，O-Oと並んでいる場所もあります。 これは
- 隣接するSiとOのボールが初期半径でだいたい接する
- 最近接の2つのOのボールが初期半径でだいたい接する
ように初期半径を決めているので，O-Oが優先されるかSi-Oが優先されるかが微妙な原子の間隔のゆらぎで決まってしまうからです。
このように重みを変えると今まで見えにくかった構造が見えると同時に，別の構造が見えにくくなることもよくあります。
一部の原子だけを取り出してPDを計算する¶
複数種類の原子からなる物質を解析する場合には，上のように原子ごとに重み(初期半径)を変える方法の他に，特定の種類の原子だけを取り出してPDを計算する方法があります。 ここでは，Siを取り除いてOだけでPDを計算してみましょう。
すると，O原子だけの番号を取り出してやる必要があります。
get_chemical_symbolsで原子の名前が見れたわけなので，以下のようにしてPythonのリスト内包表現を使って番号のリストを作ります。
o_atom_indexes = [i for (i, symbol) in enumerate(sio2.get_chemical_symbols()) if symbol == "O"]
o_atom_indexes
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260, 261, 262, 263, 264, 265, 266, 267, 268, 269, 270, 271, 272, 273, 274, 275, 276, 277, 278, 279, 280, 281, 282, 283, 284, 285, 286, 287, 288, 289, 290, 291, 292, 293, 294, 295, 296, 297, 298, 299, 300, 301, 302, 303, 304, 305, 306, 307, 308, 309, 310, 311, 312, 313, 314, 315, 316, 317, 318, 319, 320, 321, 322, 323, 324, 325, 326, 327, 328, 329, 330, 331, 332, 333, 334, 335, 336, 337, 338, 339, 340, 341, 342, 343, 344, 345, 346, 347, 348, 349, 350, 351, 352, 353, 354, 355, 356, 357, 358, 359, 360, 361, 362, 363, 364, 365, 366, 367, 368, 369, 370, 371, 372, 373, 374, 375, 376, 377, 378, 379, 380, 381, 382, 383, 384, 385, 386, 387, 388, 389, 390, 391, 392, 393, 394, 395, 396, 397, 398, 399, 400, 401, 402, 403, 404, 405, 406, 407, 408, 409, 410, 411, 412, 413, 414, 415, 416, 417, 418, 419, 420, 421, 422, 423, 424, 425, 426, 427, 428, 429, 430, 431, 432, 433, 434, 435, 436, 437, 438, 439, 440, 441, 442, 443, 444, 445, 446, 447, 448, 449, 450, 451, 452, 453, 454, 455, 456, 457, 458, 459, 460, 461, 462, 463, 464, 465, 466, 467, 468, 469, 470, 471, 472, 473, 474, 475, 476, 477, 478, 479, 480, 481, 482, 483, 484, 485, 486, 487, 488, 489, 490, 491, 492, 493, 494, 495, 496, 497, 498, 499, 500, 501, 502, 503, 504, 505, 506, 507, 508, 509, 510, 511, 512, 513, 514, 515, 516, 517, 518, 519, 520, 521, 522, 523, 524, 525, 526, 527, 528, 529, 530, 531, 532, 533, 534, 535, 536, 537, 538, 539, 540, 541, 542, 543, 544, 545, 546, 547, 548, 549, 550, 551, 552, 553, 554, 555, 556, 557, 558, 559, 560, 561, 562, 563, 564, 565, 566, 567, 568, 569, 570, 571, 572, 573, 574, 575, 576, 577, 578, 579, 580, 581, 582, 583, 584, 585, 586, 587, 588, 589, 590, 591, 592, 593, 594, 595, 596, 597, 598, 599, 600, 601, 602, 603, 604, 605, 606, 607, 608, 609, 610, 611, 612, 613, 614, 615, 616, 617, 618, 619, 620, 621, 622, 623, 624, 625, 626, 627, 628, 629, 630, 631, 632, 633, 634, 635, 636, 637, 638, 639, 640, 641, 642, 643, 644, 645, 646, 647, 648, 649, 650, 651, 652, 653, 654, 655, 656, 657, 658, 659, 660, 661, 662, 663, 664, 665, 666, 667, 668, 669, 670, 671, 672, 673, 674, 675, 676, 677, 678, 679, 680, 681, 682, 683, 684, 685, 686, 687, 688, 689, 690, 691, 692, 693, 694, 695, 696, 697, 698, 699, 700, 701, 702, 703, 704, 705, 706, 707, 708, 709, 710, 711, 712, 713, 714, 715, 716, 717, 718, 719, 720, 721, 722, 723, 724, 725, 726, 727, 728, 729, 730, 731, 732, 733, 734, 735, 736, 737, 738, 739, 740, 741, 742, 743, 744, 745, 746, 747, 748, 749, 750, 751, 752, 753, 754, 755, 756, 757, 758, 759, 760, 761, 762, 763, 764, 765, 766, 767, 768, 769, 770, 771, 772, 773, 774, 775, 776, 777, 778, 779, 780, 781, 782, 783, 784, 785, 786, 787, 788, 789, 790, 791, 792, 793, 794, 795, 796, 797, 798, 799, 800, 801, 802, 803, 804, 805, 806, 807, 808, 809, 810, 811, 812, 813, 814, 815, 816, 817, 818, 819, 820, 821, 822, 823, 824, 825, 826, 827, 828, 829, 830, 831, 832, 833, 834, 835, 836, 837, 838, 839, 840, 841, 842, 843, 844, 845, 846, 847, 848, 849, 850, 851, 852, 853, 854, 855, 856, 857, 858, 859, 860, 861, 862, 863, 864, 865, 866, 867, 868, 869, 870, 871, 872, 873, 874, 875, 876, 877, 878, 879, 880, 881, 882, 883, 884, 885, 886, 887, 888, 889, 890, 891, 892, 893, 894, 895, 896, 897, 898, 899, 900, 901, 902, 903, 904, 905, 906, 907, 908, 909, 910, 911, 912, 913, 914, 915, 916, 917, 918, 919, 920, 921, 922, 923, 924, 925, 926, 927, 928, 929, 930, 931, 932, 933, 934, 935, 936, 937, 938, 939, 940, 941, 942, 943, 944, 945, 946, 947, 948, 949, 950, 951, 952, 953, 954, 955, 956, 957, 958, 959, 960, 961, 962, 963, 964, 965, 966, 967, 968, 969, 970, 971, 972, 973, 974, 975, 976, 977, 978, 979, 980, 981, 982, 983, 984, 985, 986, 987, 988, 989, 990, 991, 992, 993, 994, 995, 996, 997, 998, 999, ...]
以下のようにして酸素だけからなるデータを作ります。
sio2_only_o = sio2[o_atom_indexes]
show_atoms(sio2_only_o)
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
ここから先はこれまでと同様です。
hc.PDList.from_alpha_filtration(sio2_only_o.get_positions(),
vertex_symbols=sio2_only_o.get_chemical_symbols(),
save_boundary_map=True,
save_to="sio2_only_o.pdgm")
PDList(path=sio2_only_o.pdgm)
pd1 = hc.PDList("sio2_only_o.pdgm").dth_diagram(1)
pd1.histogram().plot(colorbar={"type": "log"})
pd1.histogram((1.0, 7.5), 128).plot(colorbar={"type": "log"})
必要に応じてこれまでと同様逆解析を適用してみたりするとよいでしょう。