Tips and Tricks for Transition State Searches for Click Reactions¶
This tutorial will guide you through more a complicated transition state search where the reactants do not approach from the same plane. The transition state for the Pyrithione tautomerization in the AMS-Tutorials
was rather simple as the hydrogen atom moves along the O…S direction and the transition state (TS) can be expected somewhere in between. A Potential Energy Surface (PES) scan along the H…S distance can help to find a maximum in energy which can serve as input for a TS calculation. You may want to revisit this tutorial if you are not familiar with PES scans.
Some reactions are more complicated and have more complex moves which do not easily translate to a simple PES coordinate scan. So how can we get a good guess geometry?
- Find an Organic Textbook of your choice that discusses reaction mechanisms.
- Collaborate with an Organic Chemist.
- Use a molecule kit and build the structures you need and act out possible mechanisms.
We will use the Diels-Alder reaction to explain the thought processes of a modeller for click reactions. The Diels-Alder reaction mechanism is a click reaction to build six membered rings from scratch using a diene and a dienophile. The simplest six membered ring can be formed with a molecule of butadiene as diene and a molecule of ethene as dienophile which react to form Cyclo-1-hexene.
The reaction above is presented correctly but simplified. If you study the literature you will find out that
- Diene and Dienophile are not approaching each other in the same plane.
- Cyclo-1-hexene is not flat.
As Cyclo-1-hexene is non-aromatic and has only one double bond we expect the ring to be more a distorted chair than a planar ring. Equipped with this information we will use the ADF GUI and start our work from there.
We will work with the ADF and DFTB engines. We will use routines in the AMS driver for geometry optimization, frequency calculations, Potential Energy Surface (PES) scans, and Nudged Elastic Band method (NEB).
You should be proficient to comfortably draw a structure, set up a calculation when given hints, run it and retrieve results. You may want to refresh your skill in the GUI tutorials.
Setting up the Structure¶
There are many ways how diene and dienophile can approach each other – too many to try them all. We are, however, very certain what the result of our reaction looks like, so we will start from there.
- Open AMSinputUse the drawing tools create a molecule of Cyclo-1-hexeneMake sure it’s not planar and pull one of the carbon atoms out of the plane
The algorithm to find a minimum in ADF is programmed in such a way to find the nearest local minimum. By forcing the ring to be planar, indeed, a minimum can be found. If you perform frequency calculations on this planar structure, you would notice that there are 2 small negative frequencies – and you are stuck on some saddle point of the second order. Making the structure aplanar opens a new search space and you will find a true local minimum with no negative frequencies.
If you are uncertain how to manipulate, you can use a start structure from
Setting up the calculations¶
After manipulating our input to an aplanar structure we will perform an optimization to find a minimum for Cyclo-1-hexene.
- From the Main tab choose the following optionsTask → Geometry OptimizationXC-Functional → GGA-D → PBE-D3(BJ)Basis set → TZPFrozen Core → None
These settings will provide you with a good starting point for organic reactions in general.
- Save the file as DA_product.ams and run the calculation
When you look at the output split the molecule in your mind to the reaction partners that must have formed this product.
C2-C1-C6-C5 are carbons that stem from the diene and C3-C4 are carbons that stem from the dienophile. Thus, our next task will be to reverse the reaction and move the reactants apart, similar to unclicking of LEGO bars. This is the intuitive approach to finding a transition state.
Follow your intuition¶
Starting from the optimized structure from the previous section, we will setup a PES scan and we will use a DFTB method to speed up the process.
- Select the C2 and C3 atoms in the GUI using your left mouse button.Press + to add the scan coordinate.Set Number of scan points for coordinate to
10.Set the scan range from
3Å.Save the file as
DA_PES_DFTB.amsand run the calculation.
The calculation may finish with a warning such as:
This is ok for us as we are just after a good intuitive guess, so we can ignore this message.
You will be presented with the corresponding structure and you can recognize already the diene in one plane and the dienophile approaching from above, as described in a textbook for this type of reaction. We will use this structure to perform a frequency calculation to see where on the PES we actually are and to check if this structure is suitable for a transition state (TS) calculation.
- In AMSmovie, using either arrow keys or the slider, select Frame number 8 (or the frame where the maximum is)Click on File → Update Geometry in InputThis will bring AMSinput to the front updated with the geometry of Frame 8In AMSinput use the Bonds → Guess Bonds as the user-drawn bonds do no longer match the geometry
Performing a TS calculation is not simple; we are looking for a special point of the PES, a saddle point of the first order. This is where the curvature of the PES area is positive in all directions but one. In other words, we are looking for one negative frequency that corresponds to the move along our reaction coordinates, i.e. C2-C3 and C5-C6.
Thus, we will have to do a frequency calculation.
- In AMSinput, go to the Main panelSelect Task → Single PointSelect Frequencies
We can now run the Frequency calculation, which will compute the Hessian and the normal modes:
The first mode should have an imaginary frequency (which in the table is shown as a negative frequency). Click on the line in the table corresponding to the imaginary frequency to visualize the mode. You will see C2 and C3 moving towards and away from each other as you would expect them to do in the actual reaction.
Now we are ready to let DFTB find the real TS.
- Open the AMSinput window of the job “DA_Freq_DFTB”Select Task → Transition StateMake sure that Frequencies is still checked
You will need the frequency file to assist the TS calculation. The algorithm needs the Hessian as input.
- Go to the panel Details → Geometry OptimizationSelect Initial Hessian → From FileClick on the folder icon next to Initial Hessian From:Select the file dftb.rkf in the folder DA_Freq_DFTB.results
Check if your input still carries the geometry constraint from the previous PES scan calculation, which has to be removed before starting the new calculation:
- Go to the panel Model → Geometry Constraints and PES ScanRemove the C(2) C(3) constraint by clicking on “-“ in front of it.
We can now run the TS search calculation:
There is only ONE negative frequency defining a transition state!
You will see to molecules, clearly identifiable as butadiene and ethene, only one negative frequency, and when you animate the vibration you see the molecules moving towards each other as if they would form the product and away from each other, as if they were splitting up again.
What can possibly go wrong?¶
You could find more than one negative frequency which means you are lost on the PES and have not found the TS. This can happen if the PES scan goes wrong because you did not have a good starting structure. In the case of click reactions, rearrange the two reactants in such a fashion that the resemble the suggested text book mechanism, i.e.
… then repeat the PES.
Animate the negative frequencies and identify the one that represents the reaction coordinate, i.e. the one that looks like a movie of the reaction to happen. Then animate the other negative ones that have nothing to do with the reaction. Some common scenarios can be:
|geometries are forced to be planar, but should not||force a slight aplanarity, repeat TS calculation|
|local geometries are staggered, but should be eclipsed||change the torsion angle, repeat TS calculation|
What are the next steps?¶
- Reuse the DFTB geometry and Hessian to calculate the TS at the DFT level, e.g. PBE-D3(BJ)/TZP. This works well if the DFT and DFTB potential energy surfaces are reasonably similar. Switch to ADF with the updated geometry, select Transition State as Task. In the Details → Geometry Convergence panel click on the folder next to the Initial Hessian field and select the dftb.rkf from the DFTB TS calculation.
- If the TS does not converge, you could retry the above steps taken with DFTB, but now at the DFT level with ADF.
- After you’ve obtained a TS at the DFT level, perform a frequency calculation again to confirm that there is really only one negative frequency, and thus you are confident about your calculation.
- Report geometries and frequency values in your paper (supporting information).
- Perform an IRC search to find minima of reactants and products, see here.
Starting from reactant and product – Nudged Elastic Band (NEB)¶
Another way to locate a TS would be using the NEB method. This requires you to provide the structure of your reactants and your products. These structures you provide have to have a realistic geometry but do not necessarily have to be minima on the PES. In case of our Diels-Alder reaction we do know about the product but we are uncertain how exactly the diene and dienophile look in their local minima before they “started to click”.
The NEB algorithm does not know about chemical reactions; it will interpolate between the geometries of your reactants and products. A preliminary path on the PES can be established which will then be discretized into a finite number of points. Each point corresponds to a geometry which the reactants are adopting on their way to reach the product. Imagine you would thread these images like pearls on a necklace, but using a rubber band. The NEB algorithm will check the slope of the PES and “nudges” the rubber band toward the minimum energy path (MEP). The MEP is defined as the path on the PES where every point is at an energy minimum in all directions perpendicular to the path. It also passes through at least one ﬁrst-order saddle point which is the transition state we are looking for.
To help you select two suitable input structures, let’s revisit your file DA_PES_DFTB.ams.
- Choose File → Save Geometry… and save the file as
As our reactant, we could use frame 9 in our results.
- Click on the Frame 9 of the red curveChoose File → Save Geometry… and save the file as reactants.xyz
We now need to load both files into the ADF GUI.
- Switch to AMSinput by loading the DA_PES_DFTB job to retrieve the right initial settingsSelect the molecule and delete itClick on File → Import coordinates… and select the file reactants.xyzAdd a new molecule: Edit → New MoleculeClick on File → Import coordinates… and select the file product.xyzRename the tabs at the bottom to “reactants” and “product”
You are now ready to setup your calculation. For speed we will use DFTB, but you could also use ADF if you wish:
- Go to the Geometry Optimization panel in the Details sectionIn the Main tab chooseTask → NEB
Select the details for the NEB task:
In the output file you can find details of the TS found by NEB (and that it is indeed a TS):
NEB found a transition state! ------------------------------------------------------------ TS barrier height from the left 0.01354825 Hartree 8.502 kcal/mol 35.571 kJ/mol TS barrier height from the right 0.11569338 Hartree 72.599 kcal/mol 303.753 kJ/mol ------------------------------------------------------------ Transition state geometry -------- Geometry -------- Atoms Index Symbol x (angstrom) y (angstrom) z (angstrom) 1 C -3.63967395 -3.82921998 0.26868933 2 C -2.66160984 -4.29800077 -0.55265218 3 C -0.92086846 -2.96782502 0.13862010 4 C -1.35288245 -1.68672028 0.25903016 5 C -3.46534601 -1.54747433 -0.55164324 6 C -4.04880363 -2.47192783 0.26002929 7 H -3.97590594 -4.44356012 1.09897194 8 H -2.45626585 -3.84733667 -1.51209988 9 H -2.28741349 -5.30617793 -0.43258950 10 H -0.35610869 -3.27806195 -0.72983162 11 H -0.87110841 -3.62272212 0.99600867 12 H -1.10613001 -0.94389870 -0.48594712 13 H -1.68138227 -1.30524209 1.21589823 14 H -3.01840962 -1.81938947 -1.49659877 15 H -3.71555662 -0.49941080 -0.45233597 16 H -4.68469490 -2.13638404 1.07393647
This will show you the NEB path, the highest point is the calculated transition state.
If you play the movie you should see the reaction happen. You can also open the job in AMSspectra, to see the normal mode along the reaction coordinate.
What can possibly go wrong?¶
You could fail with the NEB and find no TS. In this case revisit the structure you are least sure – which would be the reactants in our case. Imagine again each point on the NEB preliminary path corresponds to a geometry which the reactants are adopting on their way to reach the product. Make a movie out of it in your mind’s eye. Does the reaction look smooth? Are atoms forced to penetrate each other to change the geometry from reactant to product? If so, correct the geometries in such a way that this does not happen and start the NEB afresh.
What are the next steps?¶
- Use the TS obtained with DFTB and perform a frequency calculation using ADF with the settings from before, i.e. PBE-D3(BJ) with a TZP basis set to get a guess for the Hessian. (initial guess for TS calculation)
- Use the Hessian to perform a TS calculation. (actual TS calculation)
- Do perform a frequency calculation again to confirm that there is really only one negative frequency, and thus you are confident about your calculation. (confirmation of TS)
- Report geometry and frequency values in your paper.
- Perform an IRC search to find minima of reactants and products, see here.