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.2020/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.2020/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.

Import geometry from xyz file
Build the structure yourself
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

See also

If you are not yet familiar with the editing tools in AMSinput, take a look at our Introduction to Building structures.

In AMSinput
1. Select the ferrocene structure in the structure tool: StructTool → Metal complexes → Sandwiches → Ferrocene Sandwich
2. Replace Fe with a Zr: Right-click on the Fe atom → Change atom type → Zr
3. Click on the Zr atom, select the Carbon tool (CTool) and add a carbon atom
4. Click on Atoms → Add hydrogens

Change the angle between the two benzene rings:

1. Select one of the dummy atoms in the ring center
2. Holding down shift click on the Zr atom and then on the other dummy atoms in the ring center
3. Change the Xx-Zr-Xx angle to 170 (or use the slider)
/scm-uploads/doc.2020/Tutorials/_images/ZN-metallocenium-structure.png

Finally we can add a ethene molecule:

1. In the search box Search type ethene
2. Select Molecules → C2H4: Ethene (ADF)
3. Position the ethane molecule roughly as shown in the picture below
/scm-uploads/doc.2020/Tutorials/_images/ZN_ethene_position.png

We now want to pre-optimize the structure we just built:

1. Switch to the DFTB panel: ADFPanel → DFTBPanel
2. Set the Total charge to 1
3. Go to the tab Model → Solvation and select Solvent → Toluene
4. Go to the Main panel and click on Pre-Optimize

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

DFTB
ADF
MOPAC

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

For a DFT calculation with ADF, try the following settings.

1. Select the ADF panel: ADFPanel
2. Set the Total charge to 1
3. XC functional → GGA-D → PBE-D4(EEQ)
4. Basis set → DZP
5. Go to panel Model → Solvation
6. Select Solvation method → COSMO
7. Select COSMO Solvent → Toluene

The pre-optimized geometries provided throughout this tutorial are optimized with DFTB. Make sure to re-optimize them when you switch to ADF.

Tip

To just confirm a TS with ADF, do not run a full frequency calculation! Instead use the much faster PES point characterization available under Properties -> IR (Frequencies) -> Characterize PES point.

For a semiempirical calculation with MOPAC, we suggest starting with following settings

1. Total charge select 1
2. Method → PM7
3. Go to panel Model → Solvation
4. Check the box Use COSMO
5. 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.2020/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.2020/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.2020/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.2020/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
/scm-uploads/doc.2020/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.2020/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.2020/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.2020/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.2020/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
PLAMS Workflows
Run scripts
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 .

The tutorial splits into two basic steps that can be combined into one but are provided as separate scripts here for the sake of simplicity:

  • PES-scan.py runs a PES scan and extracts the highest energy structure as xyz file
  • TS-optimize.py uses the structure generated by PES-scan.py as initial guess to optimize and characterize the TS.

To execute the workflow:

1. Download optimized ethene-metallocenium complex and place it into an empty directory.
2. Download PES-scan.py and TS-optimize.py and place into the same directory.
3. Open the command line and type: plams PES-scan.py
4. Once the calculation has finished you can run the second script with plams TS-optimize.py

The scripts are meant to be run after another as the first script generates the input structure for the second script. The results are printed to the command line.

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 are found 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)

The calculations can be run directly from a console input by executing the following run scripts:

To run the scripts from the command line use the following commands

1. chmod +x ZN-PES-Scan.run
2. ZN-PES-Scan.run > ZN-PES-Scan.out
3. chmod +x ZN-TS-Opt.run
4. ZN-TS-Opt.run > ZN-TS-Opt.out

The results are then found in the corresponding output files ZN-PES-Scan.out and ZN-PES-TS-Opt.out