Exploring Cu Adatom Diffusion with LAMMPSs

Introduction

As a part of my graduate studies, I have been learning how to perform basic molecular dynamics simulations using the LAMMPS software package. This is in service of generating dynamic trajectory data for use in machine learning applications.

In this post, I will be discussing a simulation I performed of a single Cu adatom diffusing on a Cu(100) surface. It is a minor adaptation from a similar simulation written by Eric N. Hahn and is intended to be a learning exercise for myself.

In follow-up posts, I will also show how one can substitute the Cu adatom for a different element, and how one can simulate a different surface. You’ll observe differences in the behavior, especially for heavier atoms like Pt.

Motivation

What is an adatom? Adatom is short for adsorbed atom. An adsorbed atom is an atom that is bound to a surface, but not in the same way that a bulk atom is bound to the surface. An adsorbed atom is bound to the surface by a single bond, whereas a bulk atom is bound to the surface by multiple bonds. This is illustrated in the figure below.

Why simulate adatom diffusion? Well, adatoms are a common feature of surfaces. They are atoms that are not bound to the surface, but rather are free to move around on the surface. They are often the result of a surface reconstruction, or they can be introduced by the experimenter. They are also a common feature of catalytic surfaces, where they can be the active site for a reaction.

In this simulation, we will be looking at a single Cu adatom diffusing on a Cu(100) surface. This is a common surface in the literature, and it is a good starting point for learning how to simulate adatom diffusion.

Simulation Summary

We’ll break down our simulation line-by-line shortly, but first let’s summarize what we’re doing:

• We are simulating a single Cu adatom diffusing on a Cu(100) surface.
• We simulate an NVT ensemble.
• We are using the embedded atom method (EAM) potential for Cu.
• We are using a time step of 5 fs.
• We are running the simulation for 500 ps.
• Atoms are initialized at 600 K, while the lowest atoms are constantly velocity rescaled to 850 K.

Here is the full script, which we’ll spend the rest of the post breaking down:

### Original Created by Eric N. Hahn  ###
### [email protected] ###

### Modifications by Hunter Heidenreich, CSE lab (Harvard, 2023)
### [email protected]
### 2023-09-01

### Simulating adatoms ###
### Version 0.2 ###


units metal
dimension 3
boundary p p s
atom_style atomic

lattice fcc 3.614
variable cubel equal 4
variable fixer1 equal "v_cubel+2"
variable fixer2 equal "v_cubel+1.49"
region  box block -${cubel} ${cubel} -${cubel} ${cubel} -${fixer1} 1 units lattice
region cbox block -${cubel} ${cubel} -${cubel} ${cubel} -${fixer1} 0 units lattice
create_box 1 box
create_atoms 1 region cbox
create_atoms 1 single -0.5 0 0.5 units lattice
region hold block INF INF INF INF -${fixer1} -${fixer2} units lattice
region temp block INF INF INF INF -${fixer2} -${cubel} units lattice
group hold region hold
group temp region temp

pair_style eam/alloy
pair_coeff * * Cu01.eam.alloy Cu

timestep        0.005
compute         new all temp
velocity        temp create 600 12345
fix heater temp temp/rescale 1 850 850 5 1
fix nve all nve
fix freeze hold setforce 0 0 0

variable e     equal pe
variable k     equal ke
variable t     equal etotal
variable T     equal temp
fix energy all ave/time 1 50 50 v_k v_e v_t v_T file energy_avg.txt

minimize 1.0e-4 1.0e-6 1000 10000

dump eve all custom 5 dump.lammpstrj id type xu yu zu   # fx fy fz  # uncomment for forces
dump_modify eve sort id

thermo 50
run 100000  # 100_000 * 5 fs = 500 ps

Line-by-Line Breakdown

Simulation Setup

Units

units metal

This line sets the units for the simulation to metal units. This is a common choice for simulations of metals, and the full list of units can be found in the LAMMPS documentation. For our purposes, we can emphasize this choice by noting that the units for length are Å, the units for energy are eV, and the units for time are ps.

Dimension and Boundary Conditions

dimension 3

This line sets the dimensionality of the simulation to 3D.

boundary p p s

This line sets the boundary conditions for the simulation. The first two letters indicate periodic boundary conditions in the x and y directions, respectively. This means that we are simulating an infinite surface in the x and y directions. Additionally, it means that atoms will interact with their periodic images in the x and y directions.

The third letter indicates a non-periodic boundary condition in the z direction, called a “shrink-wrapped” boundary condition. This means that the surface is finite in the z direction, and that atoms will not interact with their periodic images in the z direction. It is important to use a shrink-wrapped boundary condition and not a fixed boundary condition, because we want the adatom to be able to diffuse off the surface. Where a fixed boundary condition will delete atoms that move out of the simulation box, a shrink-wrapped boundary condition will allow atoms to move out of the simulation box.

Atom Style

atom_style atomic

This line sets the atom style to atomic. This means that atoms are represented as points in space, and that they do not have any internal structure. This is the most common atom style for simulations of metals.

Lattice

lattice fcc 3.614

This line sets the lattice type and lattice constant for the simulation. The lattice type is face-centered cubic (fcc), and the lattice constant is 3.614 Å. This is the lattice constant for Cu.

Variables

variable cubel equal 4
variable fixer1 equal "v_cubel+2"
variable fixer2 equal "v_cubel+1.49"

These lines define variables that we will use later in the script. The first line defines a variable called cubel that is equal to 4. The second line defines a variable called fixer1 that is equal to cubel + 2. The third line defines a variable called fixer2 that is equal to cubel + 1.49. These variables are used to define the size of the simulation box, and they are used to define the regions that we will use to hold the adatom in place.

Regions

region  box block -${cubel} ${cubel} -${cubel} ${cubel} -${fixer1} 1 units lattice
region cbox block -${cubel} ${cubel} -${cubel} ${cubel} -${fixer1} 0 units lattice

These lines define regions that we will use to hold the adatom in place.

The first line defines a region called box that is a block of atoms. The block is defined by the minimum and maximum x, y, and z coordinates. The minimum x and y coordinates are -cubel, and the maximum x and y coordinates are cubel; hence, the block is a square with side length 2 * cubel, centered at the origin. The minimum z coordinate is -fixer1, and the maximum z coordinate is 1. Note, the second line defines a region called cbox that is a block of atoms.

Create Atoms

create_box 1 box
create_atoms 1 region cbox
create_atoms 1 single -0.5 0 0.5 units lattice

These lines create the simulation box and the atoms in the simulation box. The first line creates the simulation box. The second line creates atoms in the simulation box, specifically the bulk atoms. The third line creates an atom at a specific location in the simulation box, specifically the adatom.

Groups of Atoms

region hold block INF INF INF INF -${fixer1} -${fixer2} units lattice
region temp block INF INF INF INF -${fixer2} -${cubel} units lattice
group hold region hold
group temp region temp

Here, we define two critical regions of the simulation box. The first region is called hold, and it defines the bottom of the simulation box where atoms are held in place. The second region is called temp, and it defines the lowest layer of the simulation box where atoms are allowed to move and are velocity rescaled to maintain a constant temperature. The last two lines define groups of atoms based on the regions that we just defined.

Pair Style and Coefficients

pair_style eam/alloy
pair_coeff * * Cu01.eam.alloy Cu

These lines set the pair style and pair coefficients for the simulation. The first line sets the pair style to embedded atom method (EAM) with alloy support.

The second line sets the pair coefficients for all pairs of atoms in the simulation to the Cu EAM potential. This force-field originates from the paper by Mishin et al., and can be downloaded in a LAMMPS-compatible format from the NIST Interatomic Potentials Repository.

Time Step

timestep        0.005

This line sets the time step for the simulation to 0.005 ps. This is a common choice for simulations of metals, and it is small enough to ensure that the simulation is stable.

Velocity

velocity        temp create 600 12345

This line sets the initial velocities for the atoms in the simulation. Here, we’re initializing the velocities to a temperature of 600 K, and we’re using a random seed of 12345.

Fixes

fix heater temp temp/rescale 1 850 850 5 1
fix nve all nve
fix freeze hold setforce 0 0 0

These lines define a few of the fixes that we will use in the simulation.

• The first line defines a fix called heater that rescales the velocities of the atoms in the temp group to maintain a temperature of 850 K.
• The second line defines a fix called nve that integrates the equations of motion for the atoms in the simulation, using the velocity Verlet algorithm.
• The third line defines a fix called freeze that sets the forces on the atoms in the hold group to zero.

variable e     equal pe
variable k     equal ke
variable t     equal etotal
variable T     equal temp
fix energy all ave/time 1 50 50 v_k v_e v_t v_T file energy_avg.txt

These lines define a fix called energy that calculates the average kinetic energy, potential energy, total energy, and temperature of the simulation every 50 time steps. The average values are written to the file energy_avg.txt, where we average every 50 values over 50 time steps.

Run

Relaxation

minimize 1.0e-4 1.0e-6 1000 10000

While we could run the simulation without relaxing the system, it is generally a good idea to relax a system before running a full simulation. You should see this minimize within a step or two, indicating that the system is already fairly relaxed.

Dumping Atoms

dump eve all custom 5 dump.lammpstrj id type xu yu zu   # fx fy fz  # uncomment for forces
dump_modify eve sort id

Before we run the simulation, we define a dump file that will be used to write the positions of the atoms in the simulation. The first line defines a dump file called eve that will be written every 5 time steps. The second line sorts the atoms in the dump file by their ID, which is useful for visualization/post-processing.

If you wish to dump the forces on the atoms in the simulation, you can uncomment the # fx fy fz in the first line.

Run

thermo 50
run 100000  # 100_000 * 5 fs = 500 ps

These lines run the simulation for 100,000 time steps, which corresponds to 500 ps. The first line prints the thermo output every 50 time steps. This prints the temperature, pressure, and other useful information about the simulation.

Visualization

After running the simulation, we can visualize the results using Ovito. Ovito is a free and open-source visualization and analysis tool for atomistic simulations. We won’t go into detail about how to use Ovito here, but we’ll summarize the steps that we used to visualize the results:

• Open the dump file in Ovito.
• Color the atoms by their z coordinate.
• Restrict the height range to 0 to 2.
• Animate the simulation.

Analysis

Energy

We can plot the average kinetic energy, potential energy, total energy, and temperature of the simulation over time.

We cut the first 30 time steps to truncate the relaxation period.

The script to generate this plot is available here:

# Hunter Heidenreich, 2023
#
# Plots the energy of a simulation over time.

import matplotlib.pyplot as plt

from argparse import ArgumentParser


if __name__ == '__main__':
    parser = ArgumentParser()
    parser.add_argument('--input', type=str, required=True)
    parser.add_argument('--output', type=str, required=True)
    parser.add_argument('--skip', type=int, default=1)

    args = parser.parse_args()

    ts = []
    kes = []
    pes = []
    tes = []
    Ts = []

    with open(args.input, 'r') as f:
        for line in f.readlines():
            if line.startswith('#'):
                continue

            line = line.strip()
            if len(line) == 0:
                continue

            # TimeStep v_k v_e v_t v_T
            t, v_k, v_e, v_t, v_T = list(map(float, line.split()))
            ts.append(t)
            kes.append(v_k)
            pes.append(v_e)
            tes.append(v_t)
            Ts.append(v_T)

    # prune data
    ts = ts[args.skip:]
    kes = kes[args.skip:]
    pes = pes[args.skip:]
    tes = tes[args.skip:]
    Ts = Ts[args.skip:]

    # plot
    fig, axs = plt.subplots(2, 2, figsize=(8 * 2, 6 * 2))

    axs[0, 0].plot(ts, kes)
    axs[0, 0].set_xlabel('TimeStep')
    axs[0, 0].set_ylabel('Kinetic Energy')
    axs[0, 0].set_title('Kinetic Energy')

    axs[0, 1].plot(ts, pes)
    axs[0, 1].set_xlabel('TimeStep')
    axs[0, 1].set_ylabel('Potential Energy')
    axs[0, 1].set_title('Potential Energy')

    axs[1, 0].plot(ts, tes)
    axs[1, 0].set_xlabel('TimeStep')
    axs[1, 0].set_ylabel('Total Energy')
    axs[1, 0].set_title('Total Energy')

    axs[1, 1].plot(ts, Ts)
    axs[1, 1].set_xlabel('TimeStep')
    axs[1, 1].set_ylabel('Temperature')
    axs[1, 1].set_title('Temperature')

    plt.tight_layout()
    plt.savefig(args.output)
    plt.close()

Diffusion

We can plot the x and y coordinates of the adatom over time.

We can also plot the z coordinate of the adatom over time.

This plot especially demonstrates the exchange event of our initial adatom with an atom in the surface.

The script to generate these plots is available here:

# Hunter Heidenreich, 2023
#
# Plots the coordinates of the adatom.

import matplotlib.pyplot as plt

from argparse import ArgumentParser


if __name__ == '__main__':
    parser = ArgumentParser()
    parser.add_argument('--input', type=str, required=True)
    parser.add_argument('--output', type=str, required=True)
    parser.add_argument('--id', type=int, default=1665)
    parser.add_argument('--do_z', action='store_true')
    args = parser.parse_args()

    xs = []
    ys = []
    zs = []
    with open(args.input, 'r') as f:
        for line in f:
            if line.startswith(f'{args.id} '):
                x, y, z = line.split()[2:5]
                x, y, z = float(x), float(y), float(z)
                xs.append(x)
                ys.append(y)
                zs.append(z)

    if args.do_z:
        plt.plot(list(range(len(zs))), zs)
        plt.xlabel('step')
        plt.ylabel('z')
        plt.title(f'z vs. step for adatom {args.id}')
    else:
        plt.scatter(xs, ys, s=1)
        plt.xlabel('x')
        plt.ylabel('y')
        plt.title(f'xy coordinates for adatom {args.id}')
    plt.savefig(args.output)

Conclusion

In this post, we’ve simulated the diffusion of an adatom on a Cu(100) surface. Ideally, this script could be useful to someone working in machine learning for atomistic simulations.

If you have any thoughts or questions, don’t hesitate to reach out. Thanks for reading!

References

• LAMMPS
• Ovito
• NIST Interatomic Potentials Repository
• Mishin et al.