Methanol adsorption & Vibrational frequencies

10. Methanol adsorption onto zeolite structure

As explained breifly in section 9, zeolites are used industrialy as acid catalysts in a variety of different reactions, and these often include the adsorption of methanol. More information concerning the calculations, geometries and reactivity of such reactions can be read about in the following article:

Methanol Adsorption in Zeolites - A First-Principles Study

During the adsorption process, methanol can undergo conversion reactions to products such as gasoline or dimethyl ether.

This tutorial will describe how methanol is adsorbed onto a zeolite surface, and the adsorption-energy calculated, computationally. To find the adsorption energy, 3 individual energies must be calculated:

  • The energy of an optimized zeolite system.

  • The energy of an optimized methanol molecule.

  • The energy of an optimized methanol-bound zeolite system.

All of the above can be completed by running the input file methanol_adsorption.py found in the directory zeol_covalent/10_methanol_adsorb

To begin with, a zsm5 cluster fragment is created. The cluster is optimized using a basic optimization script. This is the same script we have seen repeatedly throughout the preceeding tutorials.

The value of the optimized energy is saved in the variable ecalc_zsm5_optimized.

ecalc_zsm5_optimized = opt.result.energy

Next, a methanol fragment is created in the same way, with coordinates read from methanol.pun:

methanol = Fragment(coords="methanol.pun")

The methanol molecule is in the qm_region only. To denote this, all of the atom labels have a ‘1’ suffix too. It is also worth noting that the coordinates of the methanol molecule are set so that the methanol molecule adsorbs onto the zsm5 structure in the correct position and orientation.

At this point, we would usually define the qm_theory and mm_theory to produce the combined qmmm object. However, as the methanol molecule is part of the zeolite qm_region only, we only define the qm_theory, with the fragment methanol set as a parameter.

qm_theory = NWChem(frag=methanol, method='scf', basis='3-21g')

In the opt object, the theory is set to qm_theory as opposed to qmmm as usual. Again, because there is are no mm atoms in the structure.

opt = Opt(theory=qm_theory,
          algorithm='lbfgs', tolerance=0.0045, trust_radius='energy',
          maxcycle=200, maxene=100, maxstep=1.0)

The value of the optimized energy is saved in the variable ecalc_methanol_optimized.

ecalc_methanol_optimized = opt.result.energy

Finally, a methanol-adsorbed zeolite cluster fragment is produced. We redefine both constituent fragments used before, and combine them as such:

zsm5.append(methanol)

The above code illustrates the use of the append() function. It adds the content of the methanol fragment to the zsm5 cluster fragment. As such, this fragment now contains the atoms, charges and connectivity of a zeolite, and a methanol molecule.

Because of these changes to the zsm5 cluster fragment, we have to redefine any qm data that will be input to the qmmm object.

qm_theory = NWChem(method='scf', basis='3-21g', maxiter=100) 

qm_region = zsm5.getRegion(1)

The atoms frozen during optimization must also be updated.

frozen = getNotQM(zsm5.natoms, qm_region)

The value of the optimized energy is saved in the variable ecalc_zsm5_optimized.

ecalc_adsorbed_optimized = opt.result.energy

Now that we have the energy of the three structures, we can compare the variable’s ecalc_zsm5_optimized, ecalc_methanol_optimized and ecalc_adsorbed_optimized to find the adsorption energy. The equation is:

\[E_{adsorption} = E_{methanol-zsm5} - (E_{methanol} + E_{zsm5})\]

The energy is calculated, with an appropriate a.u. to eV conversion factor, and displayed to screen using the following code:

adsorption_energy = ecalc_adsorbed_optimized - (ecalc_zsm5_optimized + ecalc_methanol_optimized)
print("Using these values, the energy of adsorption = ", adsorption_energy*27.21138602, " eV")

11. Vibrational frequencies of free and adsorbed methanol

The vibrational frequency of the free and adsorbed methanol molecule can be calculated by running the input file methanol_vib.py found in the directory zeol_covalent/11_methanol_vib.

When finding the vibrational frequency of the free methanol molecule, the code is very similar to a traditional optimization script explained in the previous tutorials. To begin, the fragment methanol is created using coordinates read from methanol.pun.

As all of the methanol atoms are in the qm_region, the mm_theory can be ignored.

The only difference to a usual optimization object is the thermal=True parameter. Here, DL_FIND calculates a finite-difference Hessian and thermal corrections to the enthalpy and entropy. In this mode there is no optimisation stage. If there are no frozen atoms, as in the case of the free methanol example, the rotational and translation modes are projected out before determining the vibrational frequencies. If frozen atoms are found this projection is not applicable but the softest modes can be selectively ignored in the thermochemical analysis.

Note that thermal calculations work in mass-weighted coordinates throughout. The displacements used and resulting Hessian will therefore not agree exactly with the output of a force module calculation (even if the same masses and del value are specified), where the mass-weighted Hessian is generated from an initial Cartesian Hessian.

When finding the vibrational frequency of the adsorbed methanol molecule, we must include also some atoms from the zsm5 fragment. This follows in a similar fashion to that seen in section 10.

Towards the end of the output from this run, DL_FIND will report the thermochemical analysis, which covers the vibrational frequency modes, the vibrational energy correction to E_electronic, the total zero point energy (ZPE), total vibrational enthalpy contribution (E vib) and the total vibrational entropy contibution (S vib).