Adding an Adsorbate

The next step in this tutorial will go through how to add an adsorbate to a surface. Here we will consider the adsorption of CO₂. To work out the energy of adsorption (E_ads), we need to first know the energy of the surface (MgO), the energy of the adsorbate (CO₂) and the energy of the combined system (MgO–CO₂):

\[E_{ads} = E_{MgO_{suf}+CO_2} - E_{MgO_{surf}} - E_{CO_2}\]

We have already calculated the optimised energy of MgO in the previous examples. We now need to do the same for both CO₂ and MgO–CO₂.

1. Optimised Energy of Adsorbate

Note

The scripts for this step of the tutorial can be found in mgo_surface/4_co2_opt/.

Unlike in the previous example, there is no need to consider MM methods for CO₂ as it will be part of the QM region when the adsorbed system is calculated. Instead the script co2_opt_energy.py contains a simple QM-only calculation:

# Import initial CO2 geometry
co2 = Fragment(coords="co2.xyz")

# QM settings
qm_theory = NWChem(frag=co2,
                   method='dft',
                   functional='b3lyp',
                   basis='co2_nwchem.basis',
                   path="mpirun -np 4 nwchem")

# Run QM optimisation
opt = Opt(theory=qm_theory,
          algorithm='lbfgs', trust_radius='energy',
          maxcycle=200, maxene=300) 
opt.run()

For consistency, B3LYP density functional calculations are used as in the MgO calculation, and the basis set used for CO₂ is the same Ahlrichs pVDZ basis used for oxygen in MgO.

2. Optimised Energy of the Adsorbed State

Note

The scripts for this step of the tutorial can be found in mgo_surface/5_mgo+co2_opt/.

The combined MgO–CO₂ system requires a full QM/MM calculation of a similar nature to the bare MgO surface.

In the script mgoco2.py the combined system is first created by appending a CO2 molecule to the optimised MgO cluster. This is the starting structure for the optimisation and is saved to a new punch file mgoco2start.pun.

# Merge MgO surface cluster and CO2 geometries
mgo = Fragment(coords='mgo_shells_opt.pun', connect_mode=None)
co2 = Fragment(coords='co2.xyz', connect_mode=None)
mgo.append(co2)
mgo.save('mgoco2start.pun','pun')

We now have a total of 2373 atoms in the system (with indices 0-2372). The last three of these are the atoms of the adsorbed CO₂ molecule. The qm_region now contains both region 1 of the MgO cluster and CO₂:

# QM region contains MgO cluster region 1 and CO2
co2_atoms = [2370, 2371, 2372]
qm_region = np_append(mgo.getRegion(1), co2_atoms)

Similarly, the active optimisation region contains MgO regions 1-3 plus CO₂:

# Active atoms for optimisation are regions 1-3 plus CO2
active_region = np_append(mgo.getRegion(1,2,3), co2_atoms)

The QM/MM setup is otherwise very similar to the MgO optimisation step, except that the basis mgoco2_nwchem.basis passed to the NWChem interface contains basis functions for both MgO and CO₂.

With NWChem running on 4 cores, the geometry optimisation should complete in approximately 90 minutes.

As in previous stages, the resulting optimised coordinates are saved to a new punch file mgoco2end.pun, and the result is validated against a pre-calculated optimised energy.

You can also visualise the optimised geometry of the combined system using the file mgoco2end.xyz. You should see that the CO₂ has been bent away from the initial geometry on adsorption to MgO.

Note

Once all three optimised energies have been found, the energy of adsorption can be calculated using the formula above. You should obtain a result of -0.0245 a.u. (-15.4 kcal/mol).