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