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, model_version="latest")
estimator = Estimator(calc_mode=EstimatorCalcMode.CRYSTAL_U0, model_version="v3.0.0")
calculator = ASECalculator(estimator)
pfp_api_client: 1.5.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を使うことで、ClusterAtomsに変換し戻す事もできます。

[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
*Force-consistent energies used in optimization.
LBFGS:    0 06:51:35      -26.683124*       0.7828
LBFGS:    1 06:51:35      -26.652447*       2.5549
LBFGS:    2 06:51:35      -26.686065*       0.0313
cell=Cell([3.809979586761215, 3.809979586761215, 3.809979586761215]), a=3.8099795867612145

得られた格子定数をもとに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.696628 -4.727695
1 deca 105 -580.091957 -5.526440
2 deca 287 -1681.898353 -5.862027
3 deca 609 -3679.872592 -6.044238
4 deca 1111 -6841.069454 -6.159333
5 icosa 55 -289.233260 -5.260541
6 icosa 309 -1830.022256 -5.924156
7 icosa 923 -5687.111161 -6.163305
8 icosa 2057 -12924.199909 -6.284788

このように、各構造にたいして原子数が変わると凝集エネルギー(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"])