Nano particle energy¶
これまで、結晶や表面構造をモデリングして様々な特性を見てきました。 現実で触媒を用いる際には、触媒の材料として貴金属が使われることが多くその材料費の削減や触媒部分の表面積向上のために、担体とよばれる基板上に触媒を微粒子の形で散布することにより作成されることがあります。
そのため、まっすぐな表面上ではなく、微粒子(Nano particle)状態での性質をシミュレーションして解析・機構解明したいことがあります。 本節ではこのようなNano particle (クラスター構造)のモデリング・解析について説明します。
[1]:
import pfp_api_client
from pfp_api_client.pfp.calculators.ase_calculator import ASECalculator
from pfp_api_client.pfp.estimator import Estimator, EstimatorCalcMode
from pfcc_extras.visualize.view import view_ngl
print(f"pfp_api_client: {pfp_api_client.__version__}")
estimator = Estimator(calc_mode=EstimatorCalcMode.CRYSTAL_U0, model_version="v3.0.0")
calculator = ASECalculator(estimator)
pfp_api_client: 1.18.0
Cluster の生成方法¶
代表的な構造は、ASE の組み込み関数をもちいて、作成することも可能です。 以下、いくつかの構造を作成してみます。
[2]:
from ase.cluster import Decahedron, Icosahedron, Octahedron, wulff_construction
Decahedron¶
10面体構造の作成ができます。
[3]:
view_ngl(Decahedron("Pt", p=3, q=1, r=0), representations=["ball+stick"])
Icosahedron¶
20面体構造の作成ができます。
[4]:
view_ngl(Icosahedron("Pt", 2), representations=["ball+stick"])
Octahedron¶
8面体構造の作成ができます。
[5]:
view_ngl(Octahedron("Pt", 3), representations=["ball+stick"])
Wulff construction¶
それぞれのミラー指数で定義される表面 surfaces
と、その表面エネルギー esurf
が与えられた際に、 その表面エネルギーが最小となるようなNano particle構造を計算、生成する関数です。 size
で目標となる原子数を指定できますが、実際に生成される原子の数はずれます。
[6]:
from ase.cluster import wulff_construction
surfaces = [(1, 0, 0), (1, 1, 0), (1, 1, 1)]
esurf = [1.0, 1.1, 0.9] # Surface energies.
lc = 3.61000
size = 50 # Number of atoms
atoms = wulff_construction(
"Cu", surfaces, esurf,
size, "fcc",
rounding="above", latticeconstant=lc)
print(f"# atoms target={size}, actual {len(atoms)}")
view_ngl(atoms, representations=["ball+stick"])
# atoms target=50, actual 55
Clusterクラス¶
ここで作られるNano particleはASE上ではCluster
クラスとして与えられます。 これは、ASE Atoms
をWrapしたクラスで、いくつかの追加methodを提供しています。
[7]:
from ase import Atoms
from ase.cluster import Cluster
cluster = Octahedron("Pt", 3)
print("Is Atoms? :", isinstance(cluster, Atoms))
print("Is Cluster?:", isinstance(cluster, Cluster))
print(f"diameter : {cluster.get_diameter()}")
print(f"layers : {cluster.get_layers()}")
print(f"surfaces : {cluster.get_surfaces()}")
Is Atoms? : True
Is Cluster?: True
diameter : 8.17556412847783
layers : [1 2 1 1 1 1 1 1 1 2 2 2 2 2]
surfaces : [[ 1 1 1]
[ 1 0 0]
[-1 -1 -1]
[ 1 -1 -1]
[-1 1 -1]
[ 1 1 -1]
[-1 -1 1]
[ 1 -1 1]
[-1 1 1]
[ 0 0 -1]
[ 0 -1 0]
[-1 0 0]
[ 0 1 0]
[ 0 0 1]]
以下のmethodを使うことで、Cluster
をAtoms
に変換し戻す事もできます。
[8]:
def cluster2atoms(cluster: Cluster) -> Atoms:
return Atoms(cluster.symbols, cluster.positions, cell=cluster.cell)
[9]:
cluster2atoms(cluster)
[9]:
Atoms(symbols='Pt19', pbc=False)
Cohesive energy¶
前節で結晶構造で説明した凝集エネルギーは、Nano particleに対しても計算できます。
以下の例では、Decahedron, Icosahedronに対して、凝集エネルギーを計算しています。
結晶の場合と異なり、nanoparticleは原子数(サイズ)によって、凝集エネルギーも異なることに注意が必要です。 サイズを大きくしていくと最終的には結晶構造の凝集エネルギーに漸近していきます。
今回は以下の参考論文をもとにして、Ru nanoparticleの凝集エネルギーを求めてみます。
参考論文:
Structural Stability of Ruthenium Nanoparticles: A Density Functional Theory Study
Calculations of Real-System Nanoparticles Using Universal Neural Network Potential PFP
まずは孤立した単分子のエネルギーを求めます。
[10]:
from ase import Atoms
symbol = "Ru"
atoms_iso = Atoms(symbol)
atoms_iso.calc = calculator
E_iso = atoms_iso.get_potential_energy()
print(f"E_iso={E_iso:.2f} eV")
E_iso=0.00 eV
次に、 Decahedron
, Icosahedron
に対して用いる格子定数 latticeconstant
を得るのために、結晶構造での最適化を行います。
[11]:
from ase.build import bulk
from ase.constraints import ExpCellFilter, StrainFilter
from ase.optimize import LBFGS
import numpy as np
atoms_bulk = bulk("Ru", "fcc", a=3.8, cubic=True)
atoms_bulk.calc = calculator
opt = LBFGS(ExpCellFilter(atoms_bulk, hydrostatic_strain=True))
opt.run()
a = np.mean(np.diag(atoms_bulk.cell))
print(f"cell={atoms_bulk.cell}, a={a}")
#E_bulk = atoms_bulk.get_potential_energy()
#print(f"E_iso={E_iso:.2f} eV")
Step Time Energy fmax
LBFGS: 0 14:23:12 -26.683119 0.782793
LBFGS: 1 14:23:12 -26.652927 2.537247
LBFGS: 2 14:23:12 -26.686066 0.034455
cell=Cell([3.8100192704244367, 3.8100192704244367, 3.8100192704244367]), a=3.8100192704244367
/tmp/ipykernel_2557/67693790.py:8: FutureWarning: Import ExpCellFilter from ase.filters
opt = LBFGS(ExpCellFilter(atoms_bulk, hydrostatic_strain=True))
得られた格子定数をもとにnanoparticleのエネルギー \(E_{\rm{cluster}}\)を求め、凝集エネルギーを得ることができます。
[12]:
def get_opt_energy(atoms, fmax=0.001, opt_mode: str = "normal", trajectory=None):
if opt_mode == "scale":
opt = LBFGS(StrainFilter(atoms, mask=[1, 1, 1, 0, 0, 0]), logfile=None, trajectory=trajectory)
elif opt_mode == "all":
opt = LBFGS(UnitCellFilter(atoms), logfile=None, trajectory=trajectory)
else:
opt = LBFGS(atoms, logfile=None, trajectory=trajectory)
opt.run(fmax=fmax)
return atoms.get_total_energy()
[13]:
result_list = []
for p in [3, 5, 7, 9, 11]:
print("p", p)
atoms_deca = Decahedron(symbol, p=p, q=1, r=0, latticeconstant=a)
atoms_deca.calc = calculator
E_deca = get_opt_energy(atoms_deca, trajectory=f"output/deca_p{p}.traj")
E_coh = (E_deca / len(atoms_deca) - E_iso)
result_list.append({"structure": "deca", "n_atoms": len(atoms_deca), "E": E_deca, "E_coh": E_coh})
p 3
p 5
p 7
p 9
p 11
[14]:
for p in [3, 5, 7, 9]:
print("p", p)
atoms_icosa = Icosahedron(symbol, p, latticeconstant=a)
atoms_icosa.calc = calculator
E_icosa = get_opt_energy(atoms_icosa, trajectory=f"output/icosa_{p}.traj")
E_coh = (E_icosa / len(atoms_icosa) - E_iso)
result_list.append({"structure": "icosa", "n_atoms": len(atoms_icosa), "E": E_icosa, "E_coh": E_coh})
p 3
p 5
p 7
p 9
構造最適化の途中経過を見てみましょう。
[15]:
from ase.io import Trajectory
view_ngl(Trajectory("output/icosa_5.traj"))
#view_ngl(Trajectory("output/icosa_9.traj"))
#view_ngl(Trajectory("output/deca_p3.traj"))
[16]:
import pandas as pd
df = pd.DataFrame(result_list)
df
[16]:
structure | n_atoms | E | E_coh | |
---|---|---|---|---|
0 | deca | 23 | -108.696635 | -4.727695 |
1 | deca | 105 | -580.091960 | -5.526439 |
2 | deca | 287 | -1681.898443 | -5.862027 |
3 | deca | 609 | -3679.872864 | -6.044238 |
4 | deca | 1111 | -6841.069468 | -6.159332 |
5 | icosa | 55 | -289.233109 | -5.260538 |
6 | icosa | 309 | -1830.022425 | -5.924157 |
7 | icosa | 923 | -5687.111619 | -6.163305 |
8 | icosa | 2057 | -12924.201634 | -6.284789 |
このように、各構造にたいして原子数が変わると凝集エネルギー(eV/atom)がどのくらい変わるかを計算することができました。
論文では、凝集エネルギーは原子数 \(N\)に対して、 \(N^{-1/3}\)と線形比例するような関係になることが図示されています。 以下、確認してみましょう。
[17]:
# df["E_coh"] = df["E"] / df["n_atoms"] - E_iso
df["N^-1/3"] = np.power(df["n_atoms"], -1/3)
[18]:
import plotly.express as px
px.scatter(df, x="N^-1/3", y="E_coh", color="structure", hover_data=["n_atoms"])