Compare shapes

It is sometimes of interest to compare nanoparticles of different structural motifs. This example illustrates how truncated octahedral, decahedral, and icosahedral particles constructed based on the same surface energies can be compared with WulffPack.

Preparations

Before creating the particles, some parameters should be chosen (the values here are arbitrary). For simplicity, we will only consider {111} and {100} facets, but WulffPack can handle an arbitrary number of facets (although the code may take a fairly long time to run). When doing quantitative analyses of the particles, it is recommended to use electronvolts as unit for energy and Ångström as unit for distance.

from wulffpack import (SingleCrystal,
                       Decahedron,
                       Icosahedron)
from ase.build import bulk
import numpy as np

# Specify parameters
prim = bulk('Au',
            crystalstructure='fcc',
            a=4.08)
surface_energies = {(1, 0, 0): 0.066,
                    (1, 1, 1): 0.06}
twin_energy = 1e-3
natoms = 1000

We then create three particles, one of every kind.

oh = SingleCrystal(surface_energies,
                   primitive_structure=prim,
                   natoms=natoms)
dh = Decahedron(surface_energies,
                twin_energy=twin_energy,
                primitive_structure=prim,
                natoms=natoms)
ih = Icosahedron(surface_energies,
                 twin_energy=twin_energy,
                 primitive_structure=prim,
                 natoms=natoms)

Volume

We may now compare the three particles. For peace of mind, we first check that they have the same volume. In WulffPack, the volume is determined from the specified number of atoms and the primitive structure together. It does not take into account that the specified number of atoms may be inadequate to create the specified shape.

print('Volume (sanity check):')
for name, particle in zip(['Oh', 'Dh', 'Ih'], [oh, dh, ih]):
    volume = particle.volume
    print('{}: {:.2f}'.format(name, volume))
print()

This snippet produces the following output:

Volume (sanity check):
Oh: 16979.33
Dh: 16979.33
Ih: 16979.33

Area

We may now turn to something more interesting: the surface area. The syntax is the same as for the volume:

print('Total surface area (excluding twin boundaries):')
for name, particle in zip(['Oh', 'Dh', 'Ih'], [oh, dh, ih]):
    area = particle.area
    print('{}: {:.2f}'.format(name, area))
print()

Total surface area (excluding twin boundaries):
Oh: 3488.01
Dh: 3468.81
Ih: 3400.92

As one might have expected, the surface area of the nearly spherical icosahedron is the smallest.

Fraction of a specified facet

It can also be of interest to compare how large fraction of the area that is, say, {100}:

print('Fraction of {100} facets:')
for name, particle in zip(['Oh', 'Dh', 'Ih'], [oh, dh, ih]):
    fraction = particle.facet_fractions.get((1, 0, 0), 0.)
    print('{}: {:.4f}'.format(name, fraction))
print()

Fraction of {100} facets:
Oh: 0.2775
Dh: 0.2373
Ih: 0.0000

The truncated octahedron has the largest fraction of {100} facets, which is the primary reason decahedra and icosahedra are often energetically favorable for small nanoparticles made of FCC metals, for which the {111} energy is typically lower than {100}.

Surface energy

We may compare the total surface energy of the three particles, this time including the energy of the twin boundaries in the decahedral and icosahedral particles:

print('Total surface energy (including twin energy):')
for name, particle in zip(['Oh', 'Dh', 'Ih'], [oh, dh, ih]):
    energy = particle.surface_energy
    print('{}: {:.2f}'.format(name, energy))
print()

Total surface energy (including twin energy):
Oh: 215.09
Dh: 215.07
Ih: 208.82

The output indicates that the icosahedron is the most stable shape, followed by the decahedron. Note, however, that this model only includes surface energies (including twin boundary energies) and no contributions from, e.g., strain (see below).

Average surface energy

It is sometimes of interest to extract the average surface energy of the particle, i.e. the average of the surface energy per facet \(E_f\) weighted with the area fractions of each facet \(f\),

\[E_\text{average} = \frac{1}{A_\text{tot}} \sum_f A_f E_f,\]

where \(A_f\) is the area of facet \(f\) and \(A_\text{tot}\) is the total area of the particle.

print('Average surface energy (excluding twin energy):')
for name, particle in zip(['Oh', 'Dh', 'Ih'], [oh, dh, ih]):
    energy = particle.average_surface_energy
    print('{}: {:.4f}'.format(name, energy))
print()

Average surface energy (excluding twin energy):
Oh: 0.0617
Dh: 0.0614
Ih: 0.0600

This quantity is of particular interest in the single crystalline (Oh) case, since it provides an estimate for the surface energy that would be measured experimentally.

Strain energy

Decahedral and icosahedral particles can be built from FCC crystalline grain only if they are strained. This strain energy can be estimated based on continuum mechanics if we know the shear modulus and Poisson’s ratio of the material [HowMar84].

shear_modulus = 0.17  # eV / Angstrom^3
poissons_ratio = 0.42  # unitless
print('Estimated strain energies:')
for name, particle in zip(['Dh', 'Ih'], [dh, ih]):
    energy = particle.get_strain_energy(shear_modulus=shear_modulus,
                                        poissons_ratio=poissons_ratio)
    print('{}: {:.2f}'.format(name, energy))
print()

Estimated strain energies:
Dh: 0.52
Ih: 5.94

We note that for this size, the strain energy is sufficient to make the decahedron less stable than the truncated octahedron (the sum of total surface energy and strain energy is larger for the decahedron). Note that these strain energies should only be taken as fairly crude approximations.

Number of corners and total edge length

Finally, it can be interesting to compare how many corners (vertices) the particles have, as well as the total length of their edges. Note that these only includes corners and edges on the surface of the particle, not the 1D or 2D defects inside decahedra and icosahedra.

print('Number of corners and total edge length:')
for name, particle in zip(['Oh', 'Dh', 'Ih'], [oh, dh, ih]):
    ncorners = particle.number_of_corners
    edge_length = particle.edge_length
    print('{}: {} corners, edge length: {:.2f}'.format(name,
                                                       ncorners,
                                                       edge_length))
Number of corners and total edge length:
Oh: 24 corners, edge length: 409.48
Dh: 42 corners, edge length: 610.03
Ih: 72 corners, edge length: 614.70

Note that for both the decahedron and the icosahedron, there are very small “reentrant” surfaces close to where the grains meet at the fivefold axis. These reentrant surfaces increase the number of corners significantly, but are unlikely to manifest in an atomistic representation of these shapes. One can in principle remove these facets by specifying a very small value for the twin energy.

Source code

The complete source code is available in examples/compare_shapes.py

from wulffpack import (SingleCrystal,
                       Decahedron,
                       Icosahedron)
from ase.build import bulk
import numpy as np

# Specify parameters
prim = bulk('Au',
            crystalstructure='fcc',
            a=4.08)
surface_energies = {(1, 0, 0): 0.066,
                    (1, 1, 1): 0.06}
twin_energy = 1e-3
natoms = 1000

# Create one particle per type
oh = SingleCrystal(surface_energies,
                   primitive_structure=prim,
                   natoms=natoms)
dh = Decahedron(surface_energies,
                twin_energy=twin_energy,
                primitive_structure=prim,
                natoms=natoms)
ih = Icosahedron(surface_energies,
                 twin_energy=twin_energy,
                 primitive_structure=prim,
                 natoms=natoms)

# Compare volumes
print('Volume (sanity check):')
for name, particle in zip(['Oh', 'Dh', 'Ih'], [oh, dh, ih]):
    volume = particle.volume
    print('{}: {:.2f}'.format(name, volume))
print()

# Compare areas
print('Total surface area (excluding twin boundaries):')
for name, particle in zip(['Oh', 'Dh', 'Ih'], [oh, dh, ih]):
    area = particle.area
    print('{}: {:.2f}'.format(name, area))
print()

# Compare fractions of {100} facets
print('Fraction of {100} facets:')
for name, particle in zip(['Oh', 'Dh', 'Ih'], [oh, dh, ih]):
    fraction = particle.facet_fractions.get((1, 0, 0), 0.)
    print('{}: {:.4f}'.format(name, fraction))
print()

# Compare total surface energies
print('Total surface energy (including twin energy):')
for name, particle in zip(['Oh', 'Dh', 'Ih'], [oh, dh, ih]):
    energy = particle.surface_energy
    print('{}: {:.2f}'.format(name, energy))
print()

# Compare average surface energies
print('Average surface energy (excluding twin energy):')
for name, particle in zip(['Oh', 'Dh', 'Ih'], [oh, dh, ih]):
    energy = particle.average_surface_energy
    print('{}: {:.4f}'.format(name, energy))
print()

# Compare strain energies
shear_modulus = 0.17  # eV / Angstrom^3
poissons_ratio = 0.42  # unitless
print('Estimated strain energies:')
for name, particle in zip(['Dh', 'Ih'], [dh, ih]):
    energy = particle.get_strain_energy(shear_modulus=shear_modulus,
                                        poissons_ratio=poissons_ratio)
    print('{}: {:.2f}'.format(name, energy))
print()

# Compare edge length and number of corners
print('Number of corners and total edge length:')
for name, particle in zip(['Oh', 'Dh', 'Ih'], [oh, dh, ih]):
    ncorners = particle.number_of_corners
    edge_length = particle.edge_length
    print('{}: {} corners, edge length: {:.2f}'.format(name,
                                                       ncorners,
                                                       edge_length))