A theoretical approach to the nucleophilic behavior of benzofused thieno[3,2-b]furans using DFT and HF based reactivity descriptors

Calculations of traditional HF and DFT based reactivity descriptors are reported for the isomeric benzofused thieno[3,2-b ]furans in order to get insight into the factors determining the nature of their interactions with electrophiles. Global reactivity descriptors such as ionization energy, molecular hardness, electrophilicity, frontier molecular orbital energies and shapes, the condensed Fukui functions, total energies were determined and used to identify the differences in the stability and reactivity of benzofused thieno[3,2-b ]furans. Additionally the bond order uniformity analysis, local ionization energy and electrostatic potential energy surfaces revealed structural differences of isomeric thieno[3,2-b ]furans. Calculated values lead to the conclusion that heterocyclic system in thieno[3,2-b ]benzofuran is more aromatic and stable than in isomeric benzothieno[3,2-b ]furan. Theoretical results are in complete agreement with the experimental results and show exceptional reactivity of C(2) atom for both isomers.


Introduction
Benzofused five-membered heterocycles have been the subject of the sustainable interest 1 because they are useful reactants in the organic synthesis.There are many experimental results for benzothieno [3,2-b]furan 1 and thieno [3,2-b]benzofuran 2 (Figure 1) showing their different reactivity and regioselective behavior in the electrophilic substitution reactions. 2,3For example it was reported 2,3 that 2-position of benzofused thieno [3,2-b]furans 1 and 2 are most reactive to the attack of electrophilic reagents in the electrophilic substitution reactions such as chlorination, bromination, formylation, nitration, etc.When the substitution is continued, the 6-position of heterocycles undergoes substitution reaction.The reactivity of heterocyclic compounds 1 and 2 under electrophilic substitution reactions conditions shows that the heterocyclic system in the compound 2 is less reactive than in the compound 1. 2,3  Experimental and theoretical considerations on reaction mechanisms of benzofused heterocycles in the electrophilic substitution reactions point out a dual character in its reactivity. 2- 5One type of the reactions is the electrophilic substitution of benzofused thieno [3,2-b]furans as an aromatic compound, resulting in the substitution of 2-hydrogen via the aromatic electrophilic substitution reaction mechanism (Scheme 1).

Scheme 1
The other course of the electrophilic substitution reaction of benzofused thieno [3,2-b]furans involves an electrophilic addition-elimination mechanism on the double C(2)=C(3) bond (Scheme 2) which was experimentally proved 2 in a bromination reaction of heterocycle 1 by trapping the unstable trans-2,3-dibromo intermediate in the reaction mixture using the 1 H NMR spectroscopy.In case of analogical reactions with heterocycle 2 the appropriate addition intermediate was not detected. 2,3cheme 2 In this contribution we report a study of the benzofused thieno [3,2-b]furans 1 and 2 using computational chemistry methods.The aim of this work is to analyze reactivity features of those molecules using Hartree-Fock method (HF) and Density functional theory (DFT) based reactivity descriptors in order to discover reasons of their different chemical behavior in the electrophilic substitution reactions.Computational chemistry methods offer a unique ability for the synthetic organic chemists to generate optimal geometry structures, and through the structural and electronic properties of reactants and products make decisions as to which of the chemical transformations will occur in reactions.
6][7] It was demonstrated 8,9 that the DFT B3LYP is a reliable method for the calculation of geometries and energies of benzofused heterocycles.1][12] Based on structural uniformity the relative aromaticity of the systems was predicted.The experimental stability of heterocycles was accurately described using the theoretical results.The differences in the stability were explained in terms of aromaticity and delocalization of electron densities on π molecular orbitals.For the stable compounds, a high π molecular orbital delocalization established between two aromatic rings, which may not be presented in the less stable isomers.
It is evident that the aromaticity correlates with the thermodynamic stability of the system. 13he completely filled set of bonding orbitals gives the benzene its thermodynamic and chemical stability.If this concept is applied to a group of aromatic isomers it is clear that isomer having the lowest potential energy is the most thermodynamically stable.Eventually for conjugated cyclic planar ring systems the exceptional thermodynamic and chemical stability was attributed to resonance stabilization.In these cases the electron delocalization enhances the rezonance stabilization energy and the stability and aromaticity of molecules.The more aromatic compound often show greater thermodynamic stability and related properties.
5][16] DFT method provides definitions of important universal concepts of molecular structure stability and reactivity. 17It was developed [18][19][20][21] an approximation for absolute hardness η: (1) In the equation ( 1) I is the vertical ionization energy and A stands for the vertical electron affinity.
According to the Koopman's theorem 22 associated within the framework of HF selfconsistent-field molecular orbital theory the ionization energy and electron affinity can be expressed through HOMO and LUMO orbital energies: The higher HOMO energy corresponds to the more reactive molecule in the reactions with electrophiles, while lower LUMO energy is essential for molecular reactions with nucleophiles. 23hus, the hardness corresponds to the gap between the HOMO and LUMO orbitals.The larger the HOMO-LUMO energy gap the harder molecule. 20 In the past the hardness has been associated with the stability of chemical system. 24This finding reported as the principle of maximum hardness formulated by Parr and Pearson [18][19][20][21] : a rule that "molecules arrange themselves to be as hard as possible".Essentially, as Pearson stated in, 24 hardness measures the resistance to change in the electron distribution in a molecule.The hardness and aromaticity show same relationship.In a number of studies shown 25 that a small HOMO-LUMO gap has been associated with antiaromaticity, and vice versa the larger the HOMO-LUMO energy gap is associated with aromaticity.Moreover Haddon and Fukuhaga 26 showed that a direct relationship exist between the resonance stabilization energies and the HOMO-LUMO gaps in annulenes and demonstrated connection between the thermodynamic stability and kinetic stability (reactivity) of aromatic compounds. 26They presented the following formula for such relation: 24 where RE is the resonance energy and ρ rs the bond order of the r-s bond.
Unlike thermodynamic stability, which is a unique property of ground state, the kinetic stability (reactivity) measures how fast particular reaction goes.The reactivity depends on energies of reactants, reaction transition states and also intermediates with possibility of various subsequent reactions leading to stable products.This illustrates the difficulties of formulating general quantitative reactivity descriptors based on ground state calculations.On the other hand it is well known that the aromatic compounds undergo electrophilic substitution reactions (aromatic substitution) more easily than they do addition reactions.In other words they exhibit tendency to retain their π-electron delocalized structure herewith resonance stabilization energy unchanged.Accordingly the relationship between the change of resonance energy and reaction activation energy exists and it depends on the reaction type. 27Since there is connection between resonance energy and HOMO/LUMO energy separation 26,28 the reactivity can be closely related to the hardness and HOMO/LUMO energies.
So the idea of absolute hardness (half of HOMO/LUMO energies) is commonly used as a criterion of chemical reactivity and stability. 28As a result Aihira et al 29 proposed index using HOMO-LUMO energy separation multiplied by a number of conjugated atoms and successfully applied this index to measure reactivity of policyclic aromatic hydrocarbons. 29This index was found to correlate with chemical reactivity of particular aromatic system.Langenaeker 30 proposed the local hardness reactivity descriptor based on global hardness and demonstrated its superiority in predicting intramolecular reactivity for aromatic electrophilic substitution.Roy et al 31 studied the reactivity of some aromatic aldehides toward acid-catalyzed aromatic exchange reactions with the DFT based reactivity descriptors hardness and local hardness.They interpret the reactivity trends with the trends of aromaticity of aromatic aldehides.They pointed out that in this instance, the aromatic ring influences the reactivity through aromatic π-electron delocalization of positive charge; increasing aromaticity causes the increase of hardness and the decrease of reactivity.
So the presented contributions revealed the fact that high aromaticity and hardness are measures of high stability and low reactivity in the particular aromatic systems.
The electron affinity can also be used in combination with ionization energy to give electronic chemical potential µ defined by Parr and Pearson 21 as the characteristic of electronegativity of molecules : The global electrophilicity index ω was introduced by Parr 32 and calculated using the electronic chemical potential µ and chemical hardness η: According to the definition this index measures the propensity of a species to accept electrons.Under Domingo et al 33 the high nucleophility and electrophility of heterocycles corresponds to opposite extremes of the scale of global reactivity indexes.A good, more reactive, nucleophile is characterized by a lower value of µ, ω; and conversely a good electrophile is characterized by a high value of µ, ω.
5][36] The HSAB principle has been used in a local sense in terms of DFT concepts such as Fukui function f(r). 34Fukui function f(r) is a local reactivity descriptor that indicates the best way to change the number of electrons in a molecule.Hence it indicates the propensity of the electronic density to deform at a given position to accept or donate electrons.  The kui function is defined by Parr and Yang as 34,36 : Where µ is electronic chemical potential defined above, ν is the external potential, ρ corresponds to the electronic density, and N is the total number of electrons of the system.The second formula for f(r), written as [δρ(r)/ δN] ν shows that it is a quantity involving the electron density of the atom or molecule in its frontier valence regions.As ρ(r) is discontinuous function of N, two different types of f(r) can be defined 37 : The problem of reactivity and aromaticity of benzofused heterocycles raises several questions.Surech and Gadre 38 characterized relationship between aromaticity of polycyclic benzenoid hydrocarbons and electrostatic potential topology.The use of molecular surfaces, based on the molecular electron density such as the molecular electrostatic potential (MEP) 39,40 has a long tradition in the qualitative interpretation of chemical reactivity.The molecular electrostatic potential gives a powerful description of molecular properties, such as strong noncovalent interactions, that are predominantly electrostatic in nature.However, much classical chemical reactivity depends on electron donor-acceptor interactions that are not encoded in the MEP.
Another indicator of electrophilic attraction is provided by the local ionization energy potential map, an overlaying of the energy of electron removal (ionization) onto the electron density.Sjoberg et al and Politzer et al 41,42 introduced the local ionization energy potential I(r), defined as : Here ρ i (r) is the electron density of the i-th molecular orbital (MO), and ε i is its energy.Murray and Politzer et al [41][42][43][44] have discussed the properties of the local ionization energy in detail.It is clear that it describes the donor properties of the molecule directly.Results reported by Clark et al 40 suggest that the local ionization energy can represent the visualization of reactivity properties of the aromatic substrate and the regioselectivity of the electrophilic substitution.The absolute reactivity can be judged from the values of the local ionization energy at the π-surface of the aromatic compound.Our goal is to analyze aromaticity of the molecules 1 and 2 and to explain their stability and relative reactivity using MEP, local ionization energy surfaces and bond order characteristics as criteria of their aromaticity.For this purpose, DFT and HF ab-initio calculations were performed on these molecules.In terms of molecular surfaces based on electron density it is possible to explain the aromatic behavior of these compounds.Optimized structures, atomic charges, HOMO-LUMO gaps, Fukui functions, global hardness, electronegativity index are also reported to explain the experimental behavior of these systems.
Since these molecules play a fundamental role in many organic reactions, it is important to make theoretical studies of reactivity descriptors that could help to understand their chemical behavior.Experimentally, the chemical reactivity of those molecules is well known.The purpose of our work was to find reactivity descriptors that explains and confirms the experimental information.In the future for those classes of molecules with unknown reactivity, these parameters could help to understand and predict their behavior.

Geometry parameters and reactivity descriptors
The optimized geometries stability and reactivity descriptors: total energy E, ionization energy I, absolute hardness η, electrophilicity index ω, frontier molecular orbital energy gap Δ H-L, bond length and bond order of isomeric heterocycles 1 and 2 calculated at the HF/6-311G* and DFT B3LYP 6-311G* level of theory are shown in the figures 2 and 3.The computed E for HF/6-311G* and DFT B3LYP 6-311G* methods confirms that thieno[3,2-b]benzofuran 2 is more stable system than benzothieno[3,2-b]furan 1.The energy difference between isomers is 3.6 kcal/mol calculated at the DFT B3LYP 6-311G* level, calculation with HF method shows the same tendency -2.8 kcal/mol difference between isomers.
As shown in figures 2, 3 the hardness and HOMO-LUMO gap as a characteristic of reactivity shows that heterocycle 1 is expected to be more reactive than 2 isomer.The experimental results 1,2 pointed out that heterocycle 1 exhibit high reactivity and antiaromatic behavior in the electrophilic reactions.While heterocycle 2 shows reactivity tendencies typical for aromatics and lower chemical reactivity comparing to 1. 1,2 Those particular results for 1 and 2 confirms the above reported studies that higher aromaticity and hardness correspond to higher stability and lower reactivity for particular aromatic systems.So for more energetically stable and less reactive heterocycle 2 the HOMO-LUMO energy gap and hardness η is larger comparing to izomer 1.
The calculated values of global electrophilicity index ω show the nucleophility power of heterocycles 1 and 2. The obtained ω values for 1 and 2 are similar.However, since heterocycle 1 exhibit a lower ω value up to 0.04 eV comparing to 2, one can expect better propensity of 1 to be involved in the reactions with electrophiles than for heterocycle 2.  The figures 2, 3 also show bond lengths and bond order (italic) values of optimized isomeric heterocycles 1 and 2 using HF/6-311G* and DFT B3LYP 6-311G* level of theory.One can see that optimized geometries of heterocycles 1 and 2 represent planar structures with n-π conjugated bond systems arising due to sulfur and oxygen lone pair electron conjugation with the π system.According to the bond order uniformity approach the ring systems that have the most uniform bond order distributions are the most stable and aromatic ones. 45This can be estimated by the bond order deviations from an average bond order; i.e., for delocalized system of benzene that contains 6 π electrons over 6 carbons average bond order is 1.5.According to our computational study, the structure of heterocycle 2 produces more uniform (more aromatic) ring system.While the less uniform ring system is the heterocycle 1. Aromatic system disarrangement in heterocycle 1 is coursed by weakening of C-O bond in the furan ring of molecule.The C(2)-O(1) bond order 0.85 at HF/6-311G* level and 0.92 at DFT B3LYP 6-311G* level of heterocycle 1 is up to 0.27 and 0.32 weaker comparing to appropriate C(2)-S(1) bond order in the molecule 2. Moreover C(2)-S(1) bond with order values 1.22 and 1.24 ( at HF/6-311G* and DFT B3LYP 6-311G* levels accordingly) is close to aromatic bond.Therefore bond order uniformity study of heterocycles 1 and 2 intimate that heterocycle 1 structurally could be analogues with molecule of aromatic benzothiophene substituted with vinylic moiety -C(2)-C(3) bond, while heterocycle 2 can be considered as a stable aromatic system of thiophene with a joined a phenoxy substituent.
It is worth to mention that two methods HF/6-311G* and DFT B3LYP 6-311G* used in this study gives us the opportunity to compare the performance of both approaches in the interpretation of reactivity descriptors.It has been found 46 that DFT B3LYP method provide a good balance between delocalized and localized bond structures and favour calculations of electron density and reactivity parameters for aromatic structures, while HF ab-initio method tend to favor structures with localized bonds. 47In our calculations, both methods, ab-initio and DFT, provided results very close each other.HF and B3LYP calculated reactivity descriptors: E, I, η, ω, Δ H-L, bond length and bond orders, despite some numerical differences, are qualitatively similar, show very similar reactivity descriptor values, and yield reasonable agreement with the relevant experiment reactivity results.It confirms the suitability of both methods for the interpretation of reactivity tendencies for heterocycles 1 and 2. Hence we may conclude that electron correlation effects are not important for our compounds.This finding is an exception from general rule and should not be extrapolated to other systems.
Further we made an attempt to compare results of bond order uniformity analysis with results of molecular surfaces, based on the molecular electron density analysis.Since the DFT method provides more convenient and accurate way to calculate electron density surfaces and to estimate the ionization energy of a large molecular system than earlier proposed HF method, 48    The presented MEP surface, an overlaying of the electrostatic potential (the attraction or repulsion of a positive charge for a molecule) is valuable for describing overall molecular charge distribution as well as anticipating sites of electrophilic addition.The red color represent negatively charged areas of surface (i.e.those areas where accepting an electrophile is most favorable).
Another indicator of electrophilic attraction is provided by the local ionization potential energy surface, an overlaying of the energy of electron removal (ionization) onto the electron density.The regions with red color represent regions in the molecular surface where electron removal goes (with minimal energy) most easily.The differences in reactivity of heterocycles 1 and 2 can be judged from the values of electrostatic potential and local ionization energy surfaces presented in Figures 4, 5.For heterocycle 1 the lowest local ionization energy values and negatively charged electrostatic potential values are found on the benzothiophene and over the C(2)=C(3) bond of furan ring.This evidently points out why the electrophilic substitution reaction for 1 undergoes via additionelimination mechanism.
On the contrary the electrostatic potential and local ionization energy surfaces for heterocycle 2 shows delocalized π-electron surface that reports stable aromatic system between benzene ring and thiophene heterocycle.This result suggests the presumable possibility of aromatic electrophilic substitution mechanism scenario in the halogenation reactions of heterocycle 2.

Figure 6
While the bond order, I(r) and MEP surfaces values for the of thieno[3,2-b]benzofuran 2 confirms delocalized π-electron surface that reports stable aromatic system between benzene ring and thiophene heterocycle.This result suggests the presumable possibility of aromatic electrophilic substitution mechanism scenario for thieno [3,2-b]benzofuran.

Molecular orbitals
The frontier molecular orbital pictures of the both molecules 1, 2 under study are shown in Figure 7.We present only the HOMO, HOMO-1 and LUMO.The energy difference between the HOMO and HOMO-1 for 1 is smaller than for 2.
For the heterocycles 1 and 2 the π molecular orbital localization exists between benzene, thiophene and benzene, furan rings accordingly as shown in HOMO shapes.For both molecules the greatest extension value of HOMO is observed on C(2) atom.Moreover for 1 the HOMO-1 is delocalized on C(2)=C(3) bond while for 2 the HOMO-1 shape located on C(3) atom.With this molecular orbital analysis the relative reactivity can be explained.The π molecular orbital delocalization agrees well with the reactivity behavior of heterocyclic rings.The greatest extension value of HOMO shape on C(2) atom suggest exceptional reactivity of this atom in the electrophilic reactions.Furthermore the HOMO-1 delocalization on C(2)=C(3) bond in case 1 compatible with addition-elimination mechanism version.

Reactivity parameters for benzofused heterocycles
As outlined in the Introduction it is possible to define atomic reactivity indexes, such as the condensed Fukui functions for given atom in a molecule.For electrophilic addition or substitution reactions that occur with benzofused heterocycles 1 and 2 condensed Fukui functions derived from DFT approach is obtained according to eq.6.
In the Figure 8 the absolute values of condensed Fukui function values for electrophilic attack calculated at HF/6-311G* and DFT B3LYP 6-311G* levels of theory are shown for relevant atoms in heterocyclic compounds 1 and 2. For 1 and 2 the largest of value belong to C(2) atom in the furan and thiophene rings.Fukui values obtained with HF and B3LYP methods shows the same tendencies.This means that C(2) atom should be mostly reactive site towards an electrophilic attack of molecules.The Mulliken charge population calculated at HF/6-311G* and DFT B3LYP 6-311G* levels presented in the Figure 9

Conclusions
A theoretical study of the stability and reactivity was carried out at the density functional theory (DFT) and Hartree-Fock (HF) calculation level for the structures of isomeric benzofused thieno [3,2-b]furans.Reactivity indexes derived from DFT and HF calculations have been successfully applied in understanding of chemical reactivity.Global descriptors such as ionization energy (I), molecular hardness (η), electrophilicity (ω), frontier molecular orbital shapes and energy gaps (ΔH-L), local ionization energy and electrostatic potential energy surfaces were determined and used to identify the differences in the reactivity of heterocycles.
HF/6-311G* and DFT B3LYP 6-311G* calculations were used in this study to compare the performance of both approaches in the interpretation of reactivity descriptors.
HF and DFT, provided results very close each other.It suggests that electron correlation effects are not important for calculated parameters of benzofused thieno [3,2-b]furans.The use of HF method is an exception from general rule applicable only for these compounds.
HF and DFT calculated reactivity descriptors: E, I, η, ω, Δ H-L, bond length and bond orders show very similar reactivity descriptor values, and yield reasonable agreement with the relevant experiment reactivity results.
In general calculated values of E, µ, I, η, ω, ΔH-L lead to the conclusion that thieno[3,2b]benzofuran is more aromatic, more stable and less reactive than isomeric benzothieno [3,2b]furan.This result concurs with the experimental information concerning the different reactivity and stability of both heterocyclic systems.
The study of bond order uniformity, local ionization energy and electrostatic potential energy surfaces analysis revealed structural differences of isomeric thieno [3,2-b]furans that explains their reactivity features: benzothieno [3,2-b]furan could be an analogue of aromatic benzothiophene substituted with C(2)=C(3) vinylic moiety.This evidently points out wherefore the electrophilic substitution reaction for benzothieno [3,2-b]furan goes via addition-elimination mechanism at the C(2).Namely these characteristics for thieno [3,2-b]benzofuran shows delocalized π-electron surface that reports stable aromatic system between benzene ring and thiophene heterocycle.This result suggests the presumable possibility of aromatic electrophilic substitution mechanism scenario for thieno [3,2-b]benzofuran.
Molecular orbital shapes presents the π molecular orbital delocalization between benzene and thiophene, also between benzene and furan rings for the isomeric benzofused thieno [3,2b]furans.Moreover the highest occupied molecular orbital is delocalized on C(2) atom for both isomers suggesting this site to be mostly reactive towards electrophiles.
The calculated values for condensed Fukui function for electrophilic attack shows the largest value belong to C(2) atom for both heterocycles.
Theoretical results from the molecular orbital analysis and Fukui function reactivity indexes are in complete agreement with the observed reactivity of these compounds showing exceptional reactivity of C(2) atom towards electrophilic attack for both isomers.
In general theoretical results are in complete agreement with observed experimental reactivity.

Computational Section
Geometry optimizations for benzothieno [3,2-b]furan 1 and thieno[3,2-b]benzofuran 2 were performed at the HF and DFT level using the GAMESS package. 49It has been found 46,47 that DFT B3LYP method provide a good balance between delocalized and localized bond structures and favour calculations of electron density and reactivity parameters for aromatic structures, while HF ab-initio method tend to favor structures with localized bonds.In our calculations, both methods, ab-initio and DFT, provided results very close each other.So, it seems that electron correlation effects are not important for our compounds, what is an exception from general rule not applicable to other systems.The geometries 1 and 2 were fully optimized at the DFT B3LYP level of theory with a 6-311G* basis set and using ab-initio method with 6-311G* basis set.The structures are minima on potential energy surface and their harmonic vibrational frequencies are positive.Visualization of molecules at their optimized geometries was performed with MOLEKEL 50 program package and images of MEP and ionization potential surfaces was obtained using Spartan. 51
the DFT B3LYP 6-311G* basis set have been used for molecules 1 and 2 to calculate local ionization energy I(r) and molecular electrostatic potential MEP energy surfaces.The visualized results of MEP energy and I(r) surfaces are shown in Figures 4, 5 .
shows the different situation.The negative charge increase in molecule represents the attraction of relevant sites of molecule in reactions with electrophiles.The highest negative charge located on C(3) atom of heterocycle 1.So Mulliken charge and Fukui indexe values for 1 suggest expectation that electrophilic reaction could occur with C(2) and C(3) atoms via three-membered cyclic intermediate.For heterocycle 2 the increase of negative charge and highest Fukui function values were observed on C(2) atom.This is compatible with experimental results showing that C(2) site of molecule proceeds directly by an electrophilic mechanism.Fukui function and Mulliken charge population at HF/6-311G* and DFT B3LYP 6-311G* levels show very similar reactivity descriptor values considering reactivity tendencies.The two methods used in the present work (DFT and HF) for calculation of Fukui function and Mulliken charge population lead to the same qualitatively and quantitatively similar description of the chemistry and reactivity of the heterocycles 1 and 2.