The compound shown below in figure 1 is listed by Aldrich as 4,5,6,7-tetrahydroindole, but is it really?
To a hardcore organic chemist it is clear that this is not an indole but a pyrrole because the lack of aromaticity in the fused ring gives this molecule the same reactivity as 2,3-diethyl pyrrole. If you search the ChemSpider database for ‘tetrahydroindole’ the search returns the following compound with the identical chemical formula C8H11N but with a different hydrogenation pattern: 2,3,3a,4-Tetrahydro-1H-indole
The real indole, upon an electrophilic attack, behaves as a free enamine yielding the product shown in figure 3 in which the substitution occurs in position 3. This compound cannot undergo an Aromatic Electrophilic Susbstitution since that would imply the formation of a sigma complex which would disrupt the aromaticity.
On the contrary, the corresponding pyrrole is substituted in position 2
These differences in reactivity towards electrophiles are easily rationalized when we plot their HOMO orbitals (calculated at the M062X/def2TZVP level of theory):
If we calculate the Fukui indexes at the same level of theory we get the highest value for susceptibility towards an electrophilic attack as follows: 0.20 for C(3) in indole and 0.25 for C(2) in pyrrole, consistent with the previous reaction schemes.
So, why is it listed as an indole? why would anyone search for it under that name? Nobody thinks about cyclohexane as 1,3,5-trihydrobenzene. According to my good friend and colleague Dr. Moisés Romero most names for heterocyles are kept even after such dramatic chemical changes due to historical and mnemonic reasons even when the reactivity is entirely different. This is only a nomenclature issue that we have inherited from the times of Hantzsch more than a century ago. We’ve become used to keeping the trivial (or should I say arbitrary) names and further use them as derivations but this could pose an epistemological problem if students cannot recognize which heterocycle presents which reactivity.
So, in a nutshell:
Chemistry makes the chemical and not the structure.
A thing we all know but sometimes is overlooked for the sake of simplicity.
I found it surprising that the trichloromethyl group could be chemically reduced into a methyl group quite rapidly in the presence of thiophenol, but once again a failed reaction in the lab gave us the opportunity to learn some nuances about the chemical reactivity of organic compounds. Even more surprising was the fact that this reduction occured through a mechanism in which chlorine atoms behave as electrophiles and not as nucleophiles.
We proposed the mechanism shown in figure 1 to be consistent with the 1H-NMR kinetic experiment (Figure 2) which shows the presence of the intermediary sulfides and leads to the observed phenyl-disulfide as the only isolable byproduct. The proposed mechanism invokes the presence of σ-holes on chlorine atoms to justify the attack of thiophenolate towards the chlorine atom leaving a carbanion behind during the first step. The NMR spectra were recorded at 195K which implies that the energy barriers had to be very low; the first step has a ~3kcal/mol energy barrier at this temperature.
To calculate these energy barriers we employed the BMK functional as implemented in Gaussian09. This functional came highly recommended to this purpose and I gotta say it delivered! The optimized geometries of all transition states and intermediaries were then taken to an MP2 single point upon which the maximum electrostatic potential on each atom (Vmax) was calculated with MultiWFN. In figure 3 we can observe the position and Vmax value of σ-holes on chlorine atoms as suggested by the mapping of electrostatic potential on the electron density of various compounds.
We later ran the same MP2 calculations on other CCl3 groups and found that the binding to an electron withdrawing group is necessary for a σ-hole to be present. (This fact was already present in the literature, of course, but reproducing it served us to validate our methodology.)
We are pleased to have this work published in PhysChemChemPhys. Thanks to Dr. Moisés Romero for letting us into his laboratory and to Guillermo Caballero for his hard work both in the lab and behind the computer; Guillermo is now bound to Cambridge to get his PhD, we wish him every success possible in his new job and hope to see him again in a few years, I’m sure he will make a good job at his new laboratory.
Literature in synthetic chemistry is full of reactions that do occur but very little or no attention is payed to those that do not proceed. The question here is what can we learn from reactions that are not taking place even when our chemical intuition tells us they’re feasible? Is there valuable knowledge that can be acquired by studying the ‘anti-driving force’ that inhibits a reaction? This is the focus of a new manuscript recently published by our research group in Tetrahedron (DOI: 10.1016/j.tet.2016.05.058) which was the basis of Guillermo Caballero’s BSc thesis.
It is well known in organic chemistry that if a molecular structure has the possibility to be aromatic it can somehow undergo an aromatization process to achieve this more stable state. During some experimental efforts Guillermo Caballero found two compounds that could be easily regarded as non-aromatic tautomers of a substituted pyridine but which were not transformed into the aromatic compound by any means explored; whether by treatment with strong bases, or through thermal or photochemical reaction conditions.
These results led us to investigate the causes that inhibits these aromatization reactions to occur and here is where computational chemistry took over. As a first approach we proposed two plausible reaction mechanisms for the aromatization process and evaluated them with DFT transition state calculations at the M05-2x/6-31+G(d,p)//B3LYP/6-31+G(d,p) levels of theory. The results showed that despite the aromatic tautomers are indeed more stable than their corresponding non-aromatic ones, a high activation free energy is needed to reach the transition states. Thus, the barrier heights are the first reason why aromatization is being inhibited; there just isn’t enough thermal energy in the environment for the transformation to occur.
But this is only the proximal cause, we went then to search for the distal causes (i.e. the reasons behind the high energy of the barriers). The second part of the work was then the calculation of the delocalization energies and frontier molecular orbitals for the non-aromatic tautomers at the HF/cc-pVQZ level of theory to get insights for the large barrier heights. The energies showed a strong electron delocalization of the nitrogen’s lone pair to the oxygen atom in the carbonyl group. Such delocalization promoted the formation of an electron corridor formed with frontier and close-to-frontier molecular orbitals, resembling an extended push-pull effect. The hydrogen atoms that could promote the aromatization process are shown to be chemically inaccessible.
Further calculations for a series of analogous compounds showed that the dimethyl amino moiety plays a crucial role avoiding the aromatization process to occur. When this group was changed for a nitro group, theoretical calculations yielded a decrease in the barrier high, enough for the reaction to proceed. Electronically, the bonding electron corridor is interrupted due to a pull-pull effect that was assessed through the delocalization energies.
The identity of the compounds under study was assessed through 1H, 13C-NMR and 2D NMR experiments HMBC, HMQC so we had to dive head long into experimental techniques to back our calculations.
Last week we posted some insights on finding Transitions States in Gaussian 09 in order to evaluate a given reaction mechanism. A stepwise methodology is tried to achieve and this time we’ll wrap the post with two flow charts trying to synthesize the information given. It must be stressed that knowledge about the chemistry of the reaction is of paramount importance since G09 cannot guess the structure connecting two minima on its own but rather needs our help from our chemical intuition. So, without further ado here is the remainder of Guillermo’s post.
METHOD 3. QST3. For this method, you provide the coordinates of your reagents, products and TS (in that order) and G09 uses the QST3 method to find the first order saddle point. As for QST2 the numbering scheme must match for all the atoms in your three sets of coordinates, again, use the connection editor to verify it. Here is an example of the input file.
link 0 --blank line-- #p b3lyp/6-31G(d,p) opt=(qst3,calcfc) geom=connectivity freq=noraman --blank line-- Charge Multiplicity Coordinates of reagents --blank line— Charge Multiplicity Coordinates of products --blank line-- Charge Multiplicity Coordinates of TS --blank line---
As I previously mentioned, it happens that you find a first order saddle point but does not correspond to the TS you want, you find an imaginary vibration that is not the one for the bond you are forming or breaking. For these cases, I suggest you to take that TS structure and manually modify the region that is causing you trouble, then use method 2.
METHOD 4. When the previous methods fail to yield your desired TS, the brute force way is to acquire the potential energy surface (PES) and visually locate your possible TS. The task is to perform a rigid PES scan, for this, the molecular structure must be defined using z-matrix. Here is an example of the input file.
link 0 --blank line-- #p b3lyp/6-31G(d,p) scan test geom=connectivity --blank line-- Charge Multiplicity Z-matrix of reagents (or products) --blank line--
In the Z-matrix section you must specify which variables (B, A or D) you want to modify. First, locate the variables you want to modify (distance B, angle A, or dihedral angle D). Then modify those lines within the Z-matrix, here is an example.
B1 1.41 3 0.05 A1 104.5 2 1.0
What you are specifying with this is that the variable B1 (a distance) is going to be stepped 3 times by 0.05. Then variable A1 (an angle) is going to be stepped 2 times by 1.0. Thus, a total of 12 energy evaluations will be performed. At the end of the calculation open the .log file in gaussview and in Results choose the Scan… option. This will open a 3D surface where you should locate the saddle point, this is an educated guess, so take the structure you think corresponds to your TS and use it for method 2.
I have not fully explored this method so I encourage you to go to Gaussian.com and thoroughly review it.
Once you have found your TS structure and via the imaginary vibration confirmed that is the one you are looking for the next step is to verify that your TS connects both your reagents and products in the potential energy surface. For this, an Intrinsic Reaction Coordinate (IRC) calculation must be performed. Here is an example of the input file for the IRC.
link 0 --blank line-- #p b3lyp/6-31G(d,p) irc=calcfc geom=connectivity --blank line-- Charge Multiplicity Coordinates of TS --blank line--
With this input, you ask for an IRC calculation, the default numbers of steps are 20 for each side of your TS in the PES; you must specify the coordinates of your TS or take them from the .chk file of your optimization. In addition, an initial force constant calculation must be made. It often occurs that the calculation fails in the correction step, thus, for complicated cases I hardly suggest to use irc=calcall, this will consume very long time (even days) but there is a 95% guaranty. If the number of points is insufficient you can put more within the route section, here is such an example for a complicated case.
link 0 --blank line-- #p b3lyp/6-31G(d,p) irc=(calcall,maxpoints=80) geom=connectivity --blank line-- Charge Multiplicity Coordinates of TS --blank line--
With this route section, you are asking to perform an IRC calculation with 80 points on each side of the PES, calculating the force constants at every point. For an even complicated case try adding the scf=qc keyword in the route section, quadratic convergence often works better for IRC calculations.
Guillermo Caballero, a graduate student from this lab, has written this two-part post on the nuances to be considered when searching for transition states in the theoretical assessment of reaction mechanisms. He’s been quite successful in getting beautiful energy profiles for organic reaction mechanisms, some of which have even explained why some reactions do not occur! A paper in Tetrahedron has just been accepted but we’ll talk about it in another post. I wanted Guillermo to share his insight into this hard practice of computational chemistry so he wrote the following post. Enjoy!
Yes, finding a transition state (TS) can be one of the most challenging tasks in computational chemistry, it requires both a good choice of keywords in your route section and all of your chemical intuition as well. Herein I give you some good tricks when you have to find a transition state using Gaussian 09 Rev. D1
METHOD 1. The first option you should try is to use the opt=qst2 keyword. With this method you provide the structures of your reagents and your products, then the program uses the quadratic synchronous transit algorithm to find a possible transition state structure and then optimize it to a first order saddle point. Here is an example of the input file.
link 0 --blank line-- #p b3lyp/6-31G(d,p) opt=qst2 geom=connectivity freq=noraman --blank line-- Charge Multiplicity Coordinates of reagents --blank line-- Charge Multiplicity Coordinates of products --blank line---
It is mandatory that the numbering must be the same in the reagents and the products otherwise the calculation will crash. To verify that the label for a given atom is the same in reagents and products you can go to Edit, then Connection. This opens a new window were you can manually modify the numbering scheme. I suggest you to work in a split window in gaussview so you can see at the same time your reagents and products.
The keyword freq=noraman is used to calculate the frequencies for your optimized structure, it is important because for a TS you must only observe one imaginary frequency, if not, then that is not a TS and you have to use another method. It also occurs that despite you find a first order saddle point, the imaginary frequency does not correspond to the bond forming or bond breaking in your TS, thus, you should use another method. I will give you advice later in the text for when this happens. When you use the noraman in this keyword you are not calculating the Raman frequencies, which for the purpose of a TS is unnecessary and saves computing time. Frequency analysis MUST be performed AT THE VERY SAME LEVEL OF THEORY at which the optimization is performed.
The main advantage for using the qst2 option is that if your calculation is going to crash, it generally crashes at the beginning, in the moment of guessing your transition state structure. Once the program have a guess, it starts the optimization. I suggest you to ask the algorithm to calculate the force constants once, this generally improves on the convergence, it will take slightly more time depending on the size of your structure but it pays off. The keyword in the route section is opt=(qst2,calcfc). Indeed, I hardly encourage you to use the calcfc keyword in any optimization you want to run.
METHOD 2. If method 1 does not work, my next advice is to use the opt=ts keyword. For this method, the coordinates in your input file are those for the TS structure. Here is an example of the input file.
link 0 --blank line-- #p b3lyp/6-31G(d,p) opt=ts geom=connectivity freq=noraman --blank line-- Charge Multiplicity Coordinates of TS --blank line--
The question that arises here is how should I get the coordinates for my TS? Well, honestly this is not a trivial task, here is where you use all the chemistry you know. For example, you can start with the coordinates of your reagents and manually get them closer. If you are forming a bond whose length is to be 1.5Å, then I suggest you to have that length in 1.6Å in your TS. Sometimes this becomes trial and error but the most accurate your TS structure is, based on your chemical knowledge, the easiest to find your TS will be. As another example, if you want to find a TS for a [1,5]-sigmatropic reaction a good TS structure will be putting the hydrogen atom that migrates in the middle point through the way. I have to insist, this method hardly depends on your imagination to elucidate a TS and on your chemistry background.
Most of the time when you use the opt=ts keyword the calculations crashes because of an error in the number of eigenvalues, you can avoid it adding noeigen to the route section; here is an example of the input file, I encourage you to use this method.
link 0 --blank line-- #p b3lyp/6-31G(d,p) opt=(ts,noeigen,calcfc) geom=connectivity freq=noraman --blank line-- Charge Multiplicity Coordinates of TS --blank line--
If you have problems in the optimization steps I suggest you to ask the algorithm to calculate the force constants in every step of the optimization opt=(ts,noeigen,calcall) this is quite a harsh method because will consume long computing time but works well for small molecules and for complicated TSs to find.
Another ‘tricky’ way to get your coordinates for your TS is to run the qst2 calculation, then if it fails, take the second- or the third-step coordinates and used them as a ‘pre-optimized’ set of coordinates for this method.
By the way, here is another useful trick. If you are evaluating a group of TSs, let’s say, if you are varying a functional group among the group, focus on finding the TS for the simplest case, then use this optimized TS as a template where you add the moieties and use this this method. This works pretty well.
For this post we’ll leave it up to here and post the rest of Guillermo’s tricks and advice on finding TS structures next week when we’ll also discuss the use of IRC calculations and some considerations on energy corrections when plotting the full energy profile. In the mean time please take the time to rate, like and share this and other posts.
Thanks for reading!
As part of an ongoing collaboration with the University of Arizona (UA) and the Center for Advanced Research and Studies (CINVESTAV – Saltillo), we are looking into the use of calix[n]arenes for bio-remediation agents capable to extract Arsenic (V) and (III) species from water. Water contamination by arsenic is a pressing issue in northern Mexico and the southern US, therefore any efforts aiming to their elimination has strong social and health repercussions.
As in previous studies, all calixarenes were optimized along with their corresponding guests within the cavity, namely H3AsO4, H2AsO4– and HAsO42- at the DFT level with the so-called Minnesota functionals by Truhlar and Cao, M06-2X/6-31G(d,p) level of theory. Interaction energies were calculated through the NBODel procedure. Calixarenes with R = SO3H and PO3H are the most promising leads. This study is now publishes in the Journal of Inclusion Phenomena and Macrocyclic Chemistry (DOI 10.1007/s10847-016-0617-0) as an online first article.
This article is also the first to be published by our undergraduate (and almost grad student in a month) Gustavo Mondragón who took this project on a side to his own research on photosynthesis.
Now my colleagues in Arizona and Saltillo, Prof. Reyes Sierra and Dr. Eddie López, respectively, will work on the experimental side of the project. Further calculations are being undertaken to extend this study to As(III) and to the use of other potential extracting materials such as metallic nanoparticles to which calixarenes could be covalently linked.
Last week at the congress of the Mexican Society of Chemistry I presented some of our results in the study of photosynthesis. Below I embeded the talk. Unfortunately for the wider audience of this blog, the talk is in Spanish (if anyone out there is willing to make subtitles for it I’ll hire you on the spot!)
The slides are also in Spanish although they should be easier to follow for non-Spanish speakers and they are uploaded in SlideShare at this link.
A big thank you to Maria Eugenia “Maru” Sandoval for all the hard work and time invested in this project!
Thanks for clicking!