Transition state search and characterization of a Ziegler Natta Catalyst

Introduction

Ziegler-Natta catalyzed polymerization

In this tutorial, you will learn how to locate the transition state (TS) for the insertion step of ethylene polymerization catalyzed by a homogeneous Ziegler-Natta catalyst.

/scm-uploads/doc.2023/Tutorials/_images/ZN-Thumbnail.png

Discovered by Karl Ziegler and first applied to polymerization by Giulio Natta, this class of catalysts soon became an industry standard for the production of various polyolefins. Consequently, Ziegler’s and Natta’s work was awarded a Nobel Prize in 1963. Nowadays, the industrial Ziegler-Natta polymerization of polyolefins is among the most important industrial processes.

In the following tutorial, we use the AMS driver and DFTB to locate and characterize the transition state of a Ziegler-Natta reaction and use these results to calculate an approximate rate constant. You are also encouraged to learn how to reuse these results to expedite follow-up ADF calculations for more accurate results.

Strategy

Consider the following Ziegler-Natta reaction:

/scm-uploads/doc.2023/Tutorials/_images/ZNonSingleSite.png

We will use a two-stage approach to locate the transition state of this reaction:

  • Find an initial guess of the TS by scanning the potential energy surface (PES scan) along the bonds newly formed during the reaction (we will use a composite coordinate scan).

  • Optimize and characterize the initial TS guess (transition state search).

Geometry Setup with AMSinput

This section shows how to construct the geometry of the ethene-metallocenium (Cp2Zr+CH3- C2H4) complex in AMSinput.

1. Download the xyz file of the ethene-metallocenium complex
2. In AMSinput select File → Import coordinates and select the file you just downloaded
3. Proceed with the Settings for the Engine section

Settings for the Engine

In the course of this tutorial, the DFTB engine with the GFN-1xTB model and implicit Toluene solvation is used. This particular setup provides a good compromise between computational speed and accuracy well suited for this tutorial. Other engines, such as our DFT code ADF, will yield more accurate results but will require longer computation time.

Here are possible engine settings for this tutorial

For a GFN1-xTB DFTB calculation with implicit Toluene solvation use the following settings:

1. Select the DFTB panel: DFTBPanel
2. Set the Total charge to 1
3. Go to the tab Model → Solvation and select Solvent → Toluene

Settings for the AMS Driver

First, take a look at the reactants and the product in the following figure to identify the two required scan coordinates:

/scm-uploads/doc.2023/Tutorials/_images/ZN-reactants-products.png

During the reaction, two bonds will be newly formed:

  • Zirconium-Carbon (between Zr-Ion and the ethene molecule)

  • Carbon-Carbon (between ethene molecule and methyl group)

These two distances will serve as scan coordinates and these will be scanned starting from their present values down to the distance of a typical Zr-C (measured from Zr-Methyl: ~2.3 Å) and Carbon-Carbon single bond (~1.5 Å).

Setup PES scan

To perform a PES scan, change the task in the main panel of the engine:

1. Task → PES Scan
2. Click on MoreBtn next to the Task to specify the scan settings (or Model → Geometry Constraints and PES Scan)
/scm-uploads/doc.2023/Tutorials/_images/ZN-pes-scan-panel.png

To add a scan coordinate, first select an atom pair

1. Select the Zr (1) and C atom (2). Hold down SHIFT Key for multiple selections.
2. Press the + Button next Zr(..)C(..) distance
3. Enter 2.3 into the second field.
4. Repeat the process for the two Carbon atoms (3) and (4).
5. Enter 1.5 into the second field for the Carbon atoms
/scm-uploads/doc.2023/Tutorials/_images/ZN-scan-coords-before.png

With these settings AMS will scan both coordinates independently, resulting in a 2D projection of the potential energy surface. However, with only one click you can combine the two distances into a single, composite coordinate:

1. Change the SC- index for the C-C distance from 2 to 1
2. Increase the number of scan points to 20

Note

This combination of coordinates implies a concerted formation of bonds. In a situation where this is not known beforehand an independent scan of the two coordinates might yield better results.

/scm-uploads/doc.2023/Tutorials/_images/ZN-scan-coords-after.png

In the last step, loosen the convergence criteria of the geometry optimizer, since a fully converged path is not required at this point, but only a guess close to the transition state geometry.

1. Click on MoreBtn next to Convergence Details
2. Soften the convergence criteria by raising the thresholds to:
- Gradient Convergence 0.005
- Energy convergence 5.0e-5
- Step convergence 5.0e-3
/scm-uploads/doc.2023/Tutorials/_images/ZN-pes-scan-convergence.png

Now you only need to save and run the calculation. The progress can be monitored with AMSmovie.

1. File → Save
2. File → Run
3. SCM → Movie

At runtime, you will see all calculated geometries, including intermediates from the geometry optimizations. Once the PES scan is complete, these intermediates will be hidden. Depending on the AMS version you are using, you might need to close and reopen AMSmovie for this effect to take place. Your plots should roughly look as follows now

/scm-uploads/doc.2023/Tutorials/_images/ZN-pes-scan-result.png

Use the scroll bar to navigate to the highest energy structure or click on the highest point in the energy curve. This geometry will serve as the initial guess for the transition state optimization!

From AMSmovie, save the structure with the highest energy

1. File → Save Geometry and give it the name “TS_initial_guess.xyz”

Results

Energies and barrier height

Once you have optimized a transition state, you can look at the energies to determine the reaction barrier height. You can find the Energy listed in the output file:

1. SCM → Output
2. Type ‘calculation results’ into the search field at the bottom, next to the magnifying glass and hit ENTER
/scm-uploads/doc.2023/Tutorials/_images/ZN-TS-search-bond-energy.png

The Energy for the reactants was already calculated when you ran a geometry optimization of the ethene-metallocenium (Cp2Zr+CH3- C2H4) in the AMSinput Geometry Setup. In the corresponding output file it was found to be Ereactants= -37.6079 Hartree.

The barrier height for this reaction can now be calculated as the difference between these two energies \(E_\text{barrier}= E_\text{TS} - E_\text{reactants} = 0.01937 \text{[Hartree]} = 12.26 \text{[kcal/mol]}\)

Kinetics and Statistical Thermal Analysis

When asked to calculate frequencies, AMS will automatically calculate some useful thermodynamic properties such as entropy, internal energy, constant volume heat capacity, enthalpy, and Gibbs free energy from a statistical mechanics partition function. You can find more information on this in the corresponding section of the AMS manual . For now, we shall use these properties to calculate the rate constant.

Using Gibbs Free Energy of activation we can estimate the rate constant as described by the Eyring-Polanyi equation:

\[k = \kappa (k_{B}T/h) \cdot \exp \left( \frac{-\Delta G^{‡}}{RT} \right)\]

with the transmission coefficient \(\kappa\) assumed to be 1 (By setting \(\kappa = 1\), the resulting rate formula is commonly known as the Transition State Theory rate. This factor corrects for those reactive trajectories that recross the transition state and return without decomposing. It always reduces the reaction rate, so in general \(\kappa \le 1\)), the rate constant for this single step reaction can be calculated from the Gibbs Free Energy of activation ( \(\Delta G^{‡}\)).

To calculate \(\Delta G^{‡}\) we need to know the Gibbs Free Energy of both the transition state and the reactants.

/scm-uploads/doc.2023/Tutorials/_images/ZN-kinetics-Gibbs-Free-E.png

The Free Energy of the transition state ( \(G_{\text{TS}}\) ) can be found in the thermodynamic properties section of the TS optimization output file:

1. SCM → Output
2. Type Gibbs into the search field at the bottom, next to the magnifying glass and hit ENTER
/scm-uploads/doc.2023/Tutorials/_images/ZN-Gibbs-Free-Energy.png
\[G_{\rm{TS}} = -23460.94 \: \textrm{kcal/mol}\]

For the Gibbs Free Energy of the reactants a frequency calculation has to be carried out. Just download the optimized ethene-metallocenium structure and run a single point calculation to obtain the frequencies

1. Use same engine settings as with previous calculations
2. Task -> Single Point
3. Check the frequencies box under the task menu
4. Once done, search for gibbs in the output file

With a resulting Gibbs Free Energy of

\[G_{\rm{reacs}} = -23475.23 \: \textrm{kcal/mol}\]

the Gibbs Free Energy of activation, \(\Delta G^{‡}\), becomes

\[\Delta G^{‡} = G_{\rm{TS}} - G_{\rm{reacs}} = 14.29 \: \textrm{kcal/mol}\]

and the Eyring-Polanyi rate constant at 295.15 K is calculated as

\[k = 1 \cdot (k_{B}T/h) \cdot \exp \left( \frac{-14.29 \: \textrm{kcal/mol}}{RT} \right) = 214 \, \textrm{s}^{-1}\]

Workflows & Scripting

This optional section deals with the workflow automation of this tutorial. It assumes basic knowledge about the command line. If you have not used the command line with AMS before, take a look at our information on how to use the command line .

To execute the workflow:

1. Download optimized ethene-metallocenium complex and place it into an empty directory.
2. Download Ziegler_Natta.py and place into the same directory.
3. Open the command line and type: "$AMSBIN/plams" Ziegler_Natta.py

Some results are printed to the command line. All results are found in the binary ams.rkf and dftb.rkf results files.

KF files can be opened and inspected with the GUI program KF Browser. They should only be viewed in Expert Mode:

1. File → Expert Mode

The PLAMS scripts provided above demonstrate how to read results from these files. A short overview of relevant entries is provided in the tables below.

The PES scan calculation results can be obtained by using the AMSResults.get_pesscan_results() function. It provides gives a dictionary of the most common results from a PESScan. Alternatively, you can find in the KF file *.results/ams.rkf

Table 1 Results in the KF file ams.rkf

Section

Variable

Meaning

PESScan

nScanCoord

Number of independently scanned coordinates

PESScan

nPoints(1-nScanCoord)

Number of points scanned for every scan coordinate

PESScan

RangeStart(1 until nScanCoord)

Starting values for scan coordinate 1 (in Bohr)

PESScan

RangeEnd(1 until nScanCoord)

End values for scan coordinate 1 (in Bohr)

PESScan

PESCoords

Coordinates of the PES (in Bohr)

PESScan

PES

Energies of the PES scan (in Hartree)

PESScan

History Indices

Index of PES point geometries in section History

The frequencies calculated after the transition state optimization are found in the engine specific results file. For DFTB this file is called *.results/dftb.rkf:

Table 2 Results in the KF file dftb.rkf

Section

Variable

Meaning

AMS results

Energy

Final energy of the TS (in Hartree)

Vibrations

Frequencies

Sorted list of final frequencies (in cm-1)

You can also get the frequencies by using the AMSResults.get_frequencies() method.