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

[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

Creating cluster structure

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

[2]:
from ase.cluster import Decahedron, Icosahedron, Octahedron, wulff_construction

Decahedron

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

[3]:
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.

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

Octahedron

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

[5]:
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.

[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 class

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

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

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

[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

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:

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

[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

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

[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

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

[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

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

[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

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.

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