Poisoning of Pt(111) by CO: From Atomistic to Mesoscopic Modeling

Note

This example script uses computational engines from the Amsterdam Modeling Suite, and you will need a license to run it. Contact license@scm.com for further questions.

In order to run the example, the AMSBIN environment variable should be properly set. You can test this by typing $AMSBIN/plams -h in a terminal: this will print the PLAMS help message. If this is not the case (e.g. you get ‘No such file or directory’), you need to set up the environmental variable $AMSBIN (see the PLAMS scripting guide for more details).

This example illustrates a procedure for simulating molecular phenomena on catalytic surfaces, starting from an atomic-level description up to a mesoscopic regime in an automated way and with the lowest human supervision possible. We follow the strategy based on the intercommunication and cooperation of three packages: AMS Driver (hereafter AMS for short), EON, and Zacros. The full workflow is carried out in Python.

On a technical level, EON is fully integrated with AMS, being directly accessible through the familiar AMS interface. The PLAMS library is used to access AMS through Python. The Zacros code is coupled through pyZacros, which will generate the required Zacros inputs by directly reading AMS/EON output. pyZacros is also used to manage the Zacros simulations, then parse the output files for subsequent post-processing within the Python environment.

Regarding the chemistry, AMS is used to explore the energy landscape of the system by using specialized atomistic algorithms. This energy landscape is processed in order to determine the binding sites and their connectivity. pyZacros translates this information to clusters, mechanisms, and binding-sites lattices: The building blocks of a Zacros calculation. It then runs the kMC (Kinetic Monte-Carlo) simulations to perform dynamic modeling of adsorption, desorption, surface diffusion, and/or reaction processes. Everything from a simple python script!

The example we show here is a toy model system for the adsorption and diffusion of carbon monoxide on the Pt(111) surface. For this tutorial, we will not be focusing on obtaining an accurate description of the system itself. The CO diffusion will primarily be used as a compact example system for illustrating the automation process. As a simplification, we will not consider any lateral interaction energy corrections among CO molecules at this stage. The simulation described here basically shows the poisoning process of Pt(111) by CO.

This example shows how to conduct a kMC simulation of CO interacting with a Pt(111) surface, starting from its atomic representation. To that aim, we use a 3x3 Pt(111) surface to avoid artificial lateral interactions between the CO and its periodic images. Here, it is essential to point out that both the adsorption sites and the reaction mechanisms will be automatically obtained from the results of the AMS calculation and translated appropriately to Zacros. No pre-defined knowledge of the system is required by the scripts.

The expected mechanisms are sketched in the following figure:

The only necessary information from the system is an initial guess for its geometry. We have used the AMS GUI to generate our CO/Pt(111) system (If you do not have access to AMS, a link to the generated XYZ file is provided below).

We generated a 3x3 Pt(111) surface, put a CO molecule on top of it, and optimized the geometry while keeping the Pt(111) surface frozen. We created two regions for the “adsorbate” and “surface” respectively. The optimized geometry gives a threefold absorption for the CO molecule.

Link to download: CO_ads+Pt111.xyz

Note

The PES exploration tools used in this tutorial will handle the optimization of the initial geometry. Even if the initial coordinates differ slightly, the script will still generate the same landscape.

If you prepare the initial geometry yourself, remember to include the regions. This is required for the script to distinguish the adsorbate. It is also recommended to orient the Pt surface in order to maximize its symmetry. In the provided XYZ file, the Pt surface orientation belongs to the P-3m1 (164) space group. AMS is able to exploit this symmetry to reduce the required number of calculations.

The full example script can be obtained through the following link: CO+Pt111.py

If you have AMS installed, you can use the included Python environment to run the script:

$ amspython CO+Pt111.py

Hereafter, we briefly explain the different sections of the script.

The script starts as follows:

import scm.plams
import scm.pyzacros

mol = scm.plams.Molecule( 'CO_ads+Pt111.xyz' )

scm.plams.init()

Firstly we load the required Python libraries: PLAMS and pyZacros (lines 1-2). Then, we create a PLAMS molecule using the XYZ geometry file we provided above (line 4). Take note that this molecule automatically includes the information about regions that are described in the XYZ file. Finally, we start the PLAMS environment (line 6).

It is convenient to divide our script into four sections for clarity:

• In the first section (Exploring the Energy Landscape), we will obtain the symmetry-irreducible energy landscape for this system, which will indirectly allow us to define the associated reaction mechanisms and the cluster expansion Hamiltonian.

• In the second section (Constructing the kMC Lattice), we will get the kMC lattice, which requires applying the symmetry operators of the Pt surface.

• In the third section (Generating the pyZacros Objects), we will use this information to set up the pyZacros simulation.

• In the fourth section (Running the pyZacros Simulation), we run the kMC simulation itself.

Exploring the Energy Landscape

This section aims to get the energy landscape of the system. By exploiting the symmetry of the system, we are able to significantly reduce the computational effort of the calculation and simplify the analysis of the obtained results. This is achieved using the PESExploration module in AMS:

engine_sett = scm.plams.Settings()
engine_sett.input.ReaxFF.ForceField = 'CHONSFPtClNi.ff'
engine_sett.input.ReaxFF.Charges.Solver = 'Direct'

sett_ads = scm.plams.Settings()
sett_ads.input.ams.Constraints.FixedRegion = 'surface'
sett_ads.input.ams.Task = 'PESExploration'
sett_ads.input.ams.PESExploration.Job = 'ProcessSearch'
sett_ads.input.ams.PESExploration.RandomSeed = 100
sett_ads.input.ams.PESExploration.NumExpeditions = 30
sett_ads.input.ams.PESExploration.NumExplorers = 4
sett_ads.input.ams.PESExploration.SaddleSearch.MaxEnergy = 2.0
sett_ads.input.ams.PESExploration.DynamicSeedStates = 'T'
sett_ads.input.ams.PESExploration.CalculateFragments = 'T'
sett_ads.input.ams.PESExploration.StatesAlignment.ReferenceRegion = 'surface'
sett_ads.input.ams.PESExploration.StructureComparison.DistanceDifference = 0.2
sett_ads.input.ams.PESExploration.StructureComparison.NeighborCutoff = 2.4
sett_ads.input.ams.PESExploration.StructureComparison.EnergyDifference = 0.05
sett_ads.input.ams.PESExploration.StructureComparison.CheckSymmetry = 'T'
sett_ads.input.ams.PESExploration.BindingSites.Calculate = 'T'
sett_ads.input.ams.PESExploration.BindingSites.DistanceDifference = 0.1

job = scm.plams.AMSJob(name='pes_exploration', molecule=mol, settings=sett_ads+engine_sett)
results_ads = job.run()

energy_landscape = results_ads.get_energy_landscape()
print(energy_landscape)

Lines 8-10 enable the ReaxFF engine. We use the CHONSFPtClNi force field, which has been specially designed to study the surface oxidation of Pt(111).

Lines 12-28 specify the PESExploration settings. This task generates the critical points that compose the energy landscape.

The positions of the Pt surface atoms are frozen (line 13). The ProcessSearch method is used to find the escape mechanisms from the different states (line 15), distributed in 30 expeditions with 4 explorers each (lines 17-18), allowing transition state crossing within a 2 eV energy window (line 19). Any newfound local minimum is used as the origin of a new expedition (line 20). For the definitive set of local minima, a geometry optimization of the corresponding independent fragments (CO and Pt surface) is carried out in order to include the gas-phase configurations in the energy landscape (line 21).

For the structure comparison, we establish that the structures are considered the same if their interatomic distances are less than 0.2 A with energy differences less than 0.05 eV (lines 23-25). Symmetry-equivalent structures are also filtered out (line 25).

We request the calculation of the binding sites (line 27). A distance threshold of 0.1 A is used when comparing sites (line 28). The site labels are based on the number of neighboring atoms within a distance of 2.4 A (line 24), as a lower value for the NeighborCutoff may fail to distinguish fcc and hcp sites.

Finally, we create the AMSJob calculation, which requires both the initial molecule and the settings object as input parameters (line 30-31). This calculation should take only a few minutes. Once it has finished, we print out a summary of the energy landscape (lines 33-34). If everything went well, you should see the following output:

PLAMS working folder: /home/user/pyzacros/examples/CO+Pt111/plams_workdir
[06.02|11:22:44] JOB pes_exploration STARTED
[06.02|11:22:44] JOB pes_exploration RUNNING
[06.02|11:23:38] JOB pes_exploration FINISHED
[06.02|11:23:38] JOB pes_exploration SUCCESSFUL
All stationary points:
======================
State 1: COPt36 local minimum @ -7.65164231 Hartree (found 1 times, results on State-1_MIN)
State 2: COPt36 local minimum @ -7.65157184 Hartree (found 1 times, results on State-2_MIN)
State 3: COPt36 local minimum @ -7.62382298 Hartree (found 1 times, results on State-3_MIN)
State 4: COPt36 transition state @ -7.62254754 Hartree (found 6 times, results on State-4_TS_2-3)
  +- Reactants: State 2: COPt36 local minimum @ -7.65157184 Hartree (found 1 times, results on State-2_MIN)
     Products:  State 3: COPt36 local minimum @ -7.62382298 Hartree (found 1 times, results on State-3_MIN)
     Prefactors: 1.549E+13:2.197E+12
State 5: COPt36 transition state @ -7.62243092 Hartree (found 6 times, results on State-5_TS_1-3)
  +- Reactants: State 1: COPt36 local minimum @ -7.65164231 Hartree (found 1 times, results on State-1_MIN)
     Products:  State 3: COPt36 local minimum @ -7.62382298 Hartree (found 1 times, results on State-3_MIN)
     Prefactors: 1.575E+13:2.200E+12
Fragment 1: CO local minimum @ -0.42445368 Hartree (results on Fragment-1)
Fragment 2: Pt36 local minimum @ -7.15428639 Hartree (results on Fragment-2)
FragmentedState 1: CO+Pt36 local minimum @ -7.57874007 Hartree (fragments [1, 2])
  +- State 1: COPt36 local minimum @ -7.65164231 Hartree (found 1 times, results on State-1_MIN)
  |  Prefactors: 8.051E+06:1.667E+16
  +- State 2: COPt36 local minimum @ -7.65157184 Hartree (found 1 times, results on State-2_MIN)
  |  Prefactors: 8.051E+06:1.642E+16
  +- State 3: COPt36 local minimum @ -7.62382298 Hartree (found 1 times, results on State-3_MIN)
     Prefactors: 8.051E+06:2.329E+15

From this output information, we can see that the calculation took less than a minute (lines 1-5) and that the obtained energy landscape contains three local minima (lines 8-10), two transition states (lines 11-18), and one fragmented state (lines 21-27). Additional information is also available, including the absolute energies, the connections between local minima and transition states, and the pre-exponential factors. To get a more amicable and interactive visualization of the energy landscape, you can use the AMSmovie tool by executing the following command:

$ amsmovie plams_workdir/pes_exploration/ams.rkf

Note

AMSmovie currently only includes non-activated exothermic adsorption (X_gas + * ⟷ X*) and bi-molecular surface reactions (X*+Y* ⟷ Z*).

To visualize the binding sites you can use AMSinput:

$ amsinput plams_workdir/pes_exploration/ams.rkf

Note that AMS detected three binding sites, labeled as A, B, and C. In literature, these are commonly labeled as fcc, bridge, and hcp, respectively. These sites were detected automatically without any preconceived information about the system. We will update the labels shortly.

Constructing the kMC Lattice

In the previous section, we obtained both the energy landscape and the associated binding sites in the irreducible symmetry representation. In this section, we are interested in generating the full kMC lattice by using these results.

sett_bs = sett_ads.copy()
sett_ads.input.ams.PESExploration.Job = 'BindingSites'
sett_bs.input.ams.PESExploration.LoadEnergyLandscape.Path= '../pes_exploration'
sett_bs.input.ams.PESExploration.LoadEnergyLandscape.GenerateSymmetryImages = 'T'
sett_bs.input.ams.PESExploration.CalculateFragments = 'F'
sett_bs.input.ams.PESExploration.StructureComparison.CheckSymmetry = 'F'

job = scm.plams.AMSJob(name='binding_sites', molecule=mol, settings=sett_bs+engine_sett)
results_bs = job.run()

We start from the settings object of the previous calculation (line 36) and load its energy landscape information (line 38). Because we do not want to run a new exploration process, the number of expeditions and the number of explorers are both set to 1.

Instead, we request the generation of the symmetry-related images (lines 39-41). The fragment calculation can be disabled to save computational time (line 40). A new AMSJob is then executed, using the same initial molecule and the new settings object (lines 44-45). This calculation creates the images by applying the symmetry operators from the surface to the adsorbate atoms and optimizing the new geometry afterward. Transition states are optimized using the dimer method. The calculation should finish in less then a minute:

[06.02|11:23:38] JOB binding_sites STARTED
[06.02|11:23:38] JOB binding_sites RUNNING
[06.02|11:23:57] JOB binding_sites FINISHED
[06.02|11:23:57] JOB binding_sites SUCCESSFUL

Using AMSinput to visualize the binding sites:

$ amsinput plams_workdir/binding_sites/ams.rkf

We now have the full kMC lattice corresponding to the 3x3 Pt(111) surface, including site connectivity and periodic boundary conditions.

Generating the pyZacros Objects

In the preceding sections, we have obtained the energy landscape and the binding site lattice. These results have to be post-processed to generate the cluster expansion Hamiltonian, the reaction mechanism, and the extended lattice for the Zacros simulation. pyZacros offers a way to do this through the RKFLoader. This class takes PLAMS output and translates it into the required pyZacros objects: Mechanism, ClusterExpansion, and Lattice. The following figure is a schematic representation of reaction processes as defined in AMS and pyZacros, and how the RKFLoader class translates them from one to the other:

In this figure, red crosses represent the binding sites. A and B are the atoms attached to the binding sites (parent atoms), and R is the remainder of the adsorbed molecule.

The following section of the script shows how to use the RKFLoader object and access the corresponding translated objects in pyZacros. It also shows how to customize the binding site labels:

loader_ads = scm.pyzacros.RKFLoader( results_ads )
loader_ads.replace_site_types( ['N333','N221','N331'], ['fcc','br','hcp'] )
loader_bs = scm.pyzacros.RKFLoader( results_bs )
loader_bs.replace_site_types( ['N333','N221','N331'], ['fcc','br','hcp'] )

print(loader_ads.clusterExpansion)
print(loader_ads.mechanism)
print(loader_bs.lattice)
loader_bs.lattice.plot()

The cluster expansion and the mechanism were taken from the symmetry-irreducible energy landscape (lines 46-47) and the lattice from the calculation of the symmetry-generated images (lines 48-49). Lines 51-53 will print out an overview of these parameters using the Zacros input format.

energetics

cluster CO*fcc
  sites 1
  lattice_state
    1 CO* 1
  site_types fcc
  graph_multiplicity 1
  cluster_eng -2.08204e+02
end_cluster

...
end_energetics

mechanism

reversible_step CO*1hcp*2br<->*1hcpCO*2br;(0,1)
  sites 2
  neighboring 1-2
  initial
    1 CO* 1
    2 * 1
  final
    1 * 1
    2 CO* 1
  site_types hcp br
  pre_expon  1.54880e+13
  pe_ratio  7.05110e+00
  activ_eng  7.89792e-01
end_reversible_step

...
end_mechanism

lattice periodic_cell
  cell_vectors
    8.31557575    0.00000000
    4.15778787    7.20149984
  repeat_cell 1 1
  n_site_types 3
  site_type_names br fcc hcp
  n_cell_sites 45
  site_types hcp fcc hcp hcp fcc hcp fcc hcp fcc hcp fcc hcp fcc hcp ...
  site_coordinates
    0.07278720    0.07705805
    0.18391338    0.18801095
    0.07278720    0.41039138
    0.40612054    0.07705805
...
    9-26  self
    7-24  self
    11-30  self
  end_neighboring_structure
end_lattice

Line 54 is used to visualize the lattice:

kMC simulations typically use a larger lattice compared to DFT in order to avoid symmetry-induced biases in the surface structures and to improve sampling statistics. Here, we choose to increase the lattice size by creating a 10x10 supercell:

loader_bs.lattice.set_repeat_cell( (10,10) )
loader_bs.lattice.plot()

Running the pyZacros Simulation

At this point, we finally have all the ingredients we need for our kMC simulation. The following section of the code specifies the simulation conditions and starts Zacros:

settings = scm.pyzacros.Settings()
settings.random_seed = 10
settings.temperature = 273.15
settings.pressure = 1.01325
settings.molar_fraction.CO = 0.1

dt = 1e-8
settings.max_time = 1000*dt
settings.snapshots = ('logtime', dt, 3.5)
settings.species_numbers = ('time', dt)

job = scm.pyzacros.ZacrosJob(
    name='zacros_job',
    lattice=loader_bs.lattice,
    mechanism=loader_ads.mechanism,
    cluster_expansion=loader_ads.clusterExpansion,
    settings=settings,
)
results_pz = job.run()

We have used standard conditions for temperature (273.15 K; line 61) and pressure (1 atm; line 62) with a molar fraction of 0.1 for the CO in the gas phase. The simulation will run for 10 µs of kMC time (line 66), writing snapshots of the lattice state at 0.01, 0.035, 0.123, 0.429, 1.5, and 5.25 µs (line 67, using the logtime option), and saving information about the number of species every 0.01 µs (line 68). Note that by default, pyZacros will start the simulation with an empty lattice.

Lines 70-74 are used to start the Zacros job using the pyZacros objects loaded earlier. This simple simulation should finish in less than a minute:

[06.02|16:08:46] JOB zacros_job STARTED
[06.02|16:08:46] JOB zacros_job RUNNING
[06.02|16:08:47] JOB zacros_job FINISHED
[06.02|16:08:47] JOB zacros_job SUCCESSFUL

Similar to the preceding tutorials, we can now visualize the output of the simulation before finally closing the PLAMS environment.

if job.ok():
    results_pz.plot_lattice_states( results_pz.lattice_states() )
    results_pz.plot_molecule_numbers( ['CO*'] )

scm.plams.finish()

Line 77 will generate snapshots of the lattice states:

Line 78 will plot the transient number of adsorbed CO molecules:

These results show that the Pt surface gets completely poisoned by CO in around 5 µs. (Note that our lattice has 4500 sites.)

In order to obtain a more realistic picture of the CO/Pt(111) system, we could extend this model by including lateral interactions in the cluster expansion. These lateral interactions would reduce the stability of surface CO species, reducing the total coverage. Extension of the reaction mechanism, including e.g. oxidation reactions, would allow study of more complex surface chemistry.

The purpose of this tutorial was to illustrate the setup of automated workflows for connecting atomistic and mesoscopic modeling. The flexible framework provided by PLAMS and pyZacros allows you to easily modify the scripts to include more complex PESExploration tasks when more extensive mechanisms and energetics are considered.

The following tutorial will further illustrate how to connect the resulting pyZacros output to macroscale models.