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 の組み込み関数をもちいて、作成することも可能です。以下、いくつかの構造を作成してみます。

https://wiki.fysik.dtu.dk/ase/ase/cluster/cluster.html

[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で目標となる原子数を指定できますが、実際に生成される原子の数はずれます。

Wulff construction - Wikipedia

[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を提供しています。

https://wiki.fysik.dtu.dk/ase/ase/cluster/cluster.html#ase.cluster.Cluster

[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に対しても計算できます。

\[E_{\rm{coh}} = E_{\rm{cluster}} - E_{\rm{isolated}}\]

以下の例では、Decahedron, Icosahedronに対して、凝集エネルギーを計算しています。

結晶の場合と異なり、nanoparticleは原子数(サイズ)によって、凝集エネルギーも異なることに注意が必要です。サイズを大きくしていくと最終的には結晶構造の凝集エネルギーに漸近していきます。

今回は以下の参考論文をもとにして、Ru nanoparticleの凝集エネルギーを求めてみます。

参考論文:

まずは孤立した単分子のエネルギーを求めます。

[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"])

確かに、それぞれの構造に対して、x軸を \(N^{-1/3}\)にとり、凝集エネルギーをプロットしてみると直線上にのるような形でプロットされていることがわかります。また、今回の計算では、DecahedralよりもIcosahedralな形のほうがエネルギーが低く、安定であるという結果が得られました。

Ru は、結晶構造ではHCP構造が安定となりますが、上記のDecahadral, IcosahedralはFCC構造をもったNano particleとなっています。上記参考論文では、103原子以上ではHCP構造が最安定だが、それ以下の原子数ではIcosahedral構造の方が安定になりそうであることが示されています。

このように、サイズ依存してどういった構造が安定になるかは異なり、それぞれのエネルギーを計算することで安定構造の分析を行うことができます。今回はHCP構造に関しての計算は行っていませんが、演習問題として取り組んで見ると良いでしょう。