slab 関連のモデリング Tips をまとめました。
slab の作成は pfcc-extras の make_surface, 編集・描画には SurfaceEditor が便利です。
目次
make_surface
"make_surface" は bulk 構造をAtoms として与え, 引数でミラー指数、 繰り返し数、 層数、 真空層の厚さを指定することで slab 構造を作成します。
引数はこちらです。与えない場合デフォルト値(以下)が適用されます。
- miller_indices: Tuple[int] = (1, 1, 1),
- layers: int = 4,
- rep: Tuple[int] = (4, 4, 1),
- vacuum: float = 30.0
from ase.build import bulk
from pfcc_extras.structure.surface import makesurface
from pfcc_extras.visualize.view import view_ngl
atoms = bulk("Si", "diamond", a=5.43)
slab = makesurface(atoms, miller_indices=(1,1,0), rep=(6,6,1), layers=6, vacuum=20)
view_ngl(slab, representations=["ball+stick"])
-
step構造の作成(make_surface)
引数として与えるミラー指数を工夫することで表面に step 構造を作成できます。
from pfcc_extras.structure.surface import makesurface
from pfcc_extras.visualize.view import view_ngl
slab = makesurface(atoms, miller_indices=(5,4,4), rep=(4,4,1), layers=48)
view_ngl(slab)
SurfaceEditor
SurfaceEditor の利用にはcalculator のセットが必要です。以下のように使用します。
from pfp_api_client import Estimator, ASECalculator
from pfcc_extras.visualize.surface_editor import SurfaceEditor
calc_mode="CRYSTAL_PLUS_D3" # Matlantis 力場の計算モード "CRYSTAL_PLUS_D3", "CRYSTAL", "MOLECULE", etc...
model_version="v7.0.0" # Matlantis 力場のバージョン "v4.0.0", "v3.0.0", etc...
estimator = Estimator(calc_mode=calc_mode, model_version=model_version)
calculator = ASECalculator(estimator)
slab.calc = calculator
se = SurfaceEditor(slab)
se.display()
atom index の表示:
上記、se.display() を実行後、別のセルで以下のように実行すると元素種と原子のインデックスを原子上に表示できます。
se.vh.view.add_label(
color="black", labelType="text",
labelText=[slab[i].symbol + str(i) for i in range(slab.get_global_number_of_atoms())],
zOffset=1.5, attachment='middle_center', radius=1.0)
直方体 slab に変換
pymatgen を利用して直方体セルに変換します。
例:
from ase.build import bulk
from pfcc_extras.structure.surface import makesurface
from pfcc_extras.visualize.view import view_ngl
atoms = bulk("Si", "diamond", a=5.43)
slab100 = makesurface(atoms, miller_indices=(1,0,0), rep=(6,6,1), layers=6, vacuum=20)
view_ngl(slab100, representations=["ball+stick"])
from pfcc_extras.structure.surface import transform_rectangular_slab
slab_transfromed = transform_rectangular_slab(slab100)
view_ngl(slab_transfromed, representations=["ball+stick"])
slabの x, y, z 座標を確認
以下のようにして XYZ 座標を確認できます。
import matplotlib.pyplot as plt
import pandas as pd
# Check z_position of atoms
atoms = slab.copy()
df = pd.DataFrame({
"x": atoms.positions[:, 0],
"y": atoms.positions[:, 1],
"z": atoms.positions[:, 2],
"symbol": atoms.symbols,
})
#display(df)
fig, axs = plt.subplots(1, 3, figsize=(15, 5))
# ax.scatter(df.index, df["z"])
# ax.grid(True)
# ax.set_xlabel("atom_index")
# ax.set_ylabel("z_position")
# plt.show(fig)
# x座標のプロット
axs[0].scatter(df.index, df["x"])
axs[0].set_title('X Coordinates')
axs[0].set_xlabel('Index')
axs[0].set_ylabel('X Value')
# y座標のプロット
axs[1].scatter(df.index, df["y"])
axs[1].set_title('Y Coordinates')
axs[1].set_xlabel('Index')
axs[1].set_ylabel('Y Value')
# z座標のプロット
axs[2].scatter(df.index, df["z"])
axs[2].set_title('Z Coordinates')
axs[2].set_xlabel('Index')
axs[2].set_ylabel('Z Value')
# グラフのレイアウトを調整
plt.tight_layout()
plt.show()
- plotly を使う場合
Z座標をプロット
import numpy as np
import plotly.graph_objects as go
x, y = np.array(range(len(slab))), sorted(slab.get_positions()[:,2])
fig = go.Figure(data=[
go.Scatter(x=x, y=y, name="sin"),
])
fig.update_xaxes(title="sorted_atom_index") # X軸タイトルを指定
fig.update_yaxes(title="height") # Y軸タイトルを指定
fig.update_layout(width=800, height=600) # 図の高さを幅を指定
fig.show()
下端層を固定(Fixatom)
Z座標を利用して下端層を固定します。
from ase.constraints import FixAtoms
# Fix atoms under z=1A
c = FixAtoms(indices=[atom.index for atom in slab if atom.position[2] <= 1])
slab.set_constraint(c)
SurfaceEditor で 固定した原子に ”Fixed” と表示させて確認します。
from pfp_api_client import Estimator, ASECalculator
calc_mode="CRYSTAL_PLUS_D3" # Matlantis 力場の計算モード "CRYSTAL_PLUS_D3", "CRYSTAL", "MOLECULE", etc...
model_version="v7.0.0" # Matlantis 力場のバージョン "v4.0.0", "v3.0.0", etc...
estimator = Estimator(calc_mode=calc_mode, model_version=model_version)
calculator = ASECalculator(estimator)
from pfcc_extras.visualize.surface_editor import SurfaceEditor
slab.calc = calculator
se = SurfaceEditor(slab)
se.display()
se.vh.view.add_label(
color="black", labelType="text",
labelText=[("Fixed" if i in slab.constraints[0].index else "") for i in range(slab.get_global_number_of_atoms())],
zOffset=1.5, attachment='middle_center', radius=1.0)
表面を水素終端
取得したZ座標を利用して水素終端します。
まずはZ座標から最表面の原子の index を取得します。
surfaceO_index = [atom.index for atom in slab if atom.position[2] >= 14]
surfaceO_index
取得した原子の上に水素を配置します(posi)
O_source = slab[surfaceO_index].copy()
O_source.symbols = "H"+str(len(surfaceO_index)) #配置する原子を指定
O_source.positions += [0,0,1] #水素を配置する位置
slab_H_terminated = slab+O_source
view_ngl(slab_H_terminated, representations=["ball+stick"])
表面欠陥構造を作成
取得した X, Y, Z座標を利用して、表面構造を削除し、表面の欠陥構造を作成します。
from ase import Atoms
from ase.io import read, write
slab_step = slab.copy()
# 削除する x,y z 座標と座標の閾値を設定します
x_threshold = 10.0
y_threshold = 10.0
z_threshold = 12.0
# x, y, z 座標がそれぞれの閾値を超える原子インデックスを取得します
indices_to_delete = [
atom.index for atom in slab_step
if (atom.x >= x_threshold or atom.y >= y_threshold) and atom.z >= z_threshold
]
# 原子を削除します
del slab_step[indices_to_delete]
view_ngl(slab_step*(2,2,1), representations=["ball+stick"])