Nanoparticle energy

So far, we have modeled crystals and surface structures and calculated several properties. In the real world, precious metals are often used as catalysts, and in order to reduce the material cost and increase the surface area of the catalyst, the catalyst is often dispersed in the form of nanoparticles on a substrate called support.

Therefore, it is sometimes desirable to simulate the properties of the catalyst not as a surface, but in the form of nanoparticles, in order to analyze and clarify the mechanism. This section describes the modeling and analysis of such nanoparticles (also called “cluster” structures).

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

Creating cluster structure

Typical cluster structures can be created using the built-in methods of ASE. Let’s make some structures.

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

from ase.cluster import Decahedron, Icosahedron, Octahedron, wulff_construction

Decahedron

It is possible to create an decaheral (polydedron with 10 faces) structure.

view_ngl(Decahedron("Pt", p=3, q=1, r=0), representations=["ball+stick"])

Icosahedron

It is possible to create an icosahedral (polydedron with 20 faces) structure.

view_ngl(Icosahedron("Pt", 2), representations=["ball+stick"])

Octahedron

It is possible to create an octahedral (polydedron with 8 faces) structure.

view_ngl(Octahedron("Pt", 3), representations=["ball+stick"])

Wulff construction

Given surfaces defined by their respective Miller indices surfaces and their surface energies esurf, the Wulff method calculates and generates a nanoparticle structure that minimizes its surface energy. You can specify the target number of atoms with size, but the actual number of atoms generated will be different.

• Wulff construction - Wikipedia

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 class

The nanoparticle created here is given as a Cluster class on ASE. This is a wrap class for Atoms, providing some additional methods.

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

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]]

The Cluster can be converted back to Atoms by using the following method if necessary.

def cluster2atoms(cluster: Cluster) -> Atoms:
    return Atoms(cluster.symbols, cluster.positions, cell=cluster.cell)

cluster2atoms(cluster)

Atoms(symbols='Pt19', pbc=False)

Cohesive energy

The cohesive energy described for crystal structures in the previous section can also be calculated for nanoparticles.

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

In the following example, cohesive energies are calculated for Decahedron and Icosahedron.

Note that, unlike crystals, nanoparticles have different cohesive energies depending on the number of atoms (size). As the size is increased, the cohesive energy eventually asymptotically approaches the cohesive energy of the crystal structure.

In this tutorial, we will try calculating the cohesive energy of Ru nanoparticles based on the following reference paper.

Reference papers:

• Structural Stability of Ruthenium Nanoparticles: A Density Functional Theory Study

• Calculations of Real-System Nanoparticles Using Universal Neural Network Potential PFP

First, we find the energy of an isolated single molecule.

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

Next, we optimize the crystal structure to obtain the lattice constants latticeconstant, which are used for the Decahedron and Icosahedron.

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))

Based on the obtained lattice constants, the energy of the nanoparticle \(E_{\rm{cluster}}\) can be determined to obtain the cohesive energy.

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()

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

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

Let’s take a look at the process of structural optimization.

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

import pandas as pd

df = pd.DataFrame(result_list)
df

	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

Thus, we were able to calculate how much the cohesive energy (eV/atom) changes as the number of atoms changes for each structure.

The paper illustrates that the cohesive energy is linearly proportional to \(N^{-1/3}\) with respect to the number of atoms \(N\). Let’s check it out below.

# df["E_coh"] = df["E"] / df["n_atoms"] - E_iso
df["N^-1/3"] = np.power(df["n_atoms"], -1/3)

import plotly.express as px

px.scatter(df, x="N^-1/3", y="E_coh", color="structure", hover_data=["n_atoms"])

If we take \(N^{-1/3}\) as the x-axis and plot the cohesive energy for each structure, we see that it is plotted on a straight line. In addition, this calculation shows that the icosahedron has lower energy and is more stable than the decahedron.

The HCP structure is stable in Ru bulk, but decahedral and icosahedral nanoparticles above are made from the FCC structure. The reference paper shows that the HCP structure is the most stable with 103 atoms or more, but the icosahedral structure is likely to be more stable with fewer atoms.

Thus, what kind of structure will be stable depends on the size, and by calculating the energy of each, we can analyze the stable structure. Although we have not performed calculations for the HCP structure in this tutorial, you may work on it as an exercise.