Computational study of the chemical reactivity properties of bis ( trimethyl tetrathiafulvalenyl ) thiophene

The chemical reactivity of four bis (trimethyltetrathiafulvalenyl) thiophene is determined by its potential (electronic) energy (hyper) surface. All the quantum chemical calculations have been carried out using DFT level of theory, B3LYP functional and 6-31G(d,p) as basis set. Molecular electrostatic potential (MEP) and HOMO-LUMO energy levels have been performed. The local reactivity descriptor such as Fukui function is also performed to determine the reactive sites within the title molecules. The chemometric methods PCA and HCA were employed to find the subset of variables that could correctly classify the compounds according to their reactivity. Indexing terms/


INTRODUCTION
Density functional theory (DFT) is a powerful, formally exact reformulation of quantum mechanics [1][2][3]. It is distinct from quantum chemical methods because in its revolutionary perspective the electronic density, rather than the many-electron wave function, plays the central role. However, chemical reactivity of a molecule cannot be described by its electron density alone, because the course of a reaction is rather determined by its response toward different perturbations caused by an approaching reagent. Sensitivities of an electron density toward structural modifications and its responses to changes in external potential and conditions are actually more important in reflecting the reactivity of the corresponding system, than its absolute values. Within the framework of the Density Functional Theory, it is possible to define universal concepts of molecular structure stability and reactivity through global and local reactivity parameters. Note that even when these reactivity parameters were derived from DFT, they may be employed to analyze the reactivity at other levels of theory [4][5][6]. These global and local reactivity indices, such as chemical potential (µ), hardness (η), and Fukui function are defined as the first or second derivative of electronic energy and electron density. In order to have the convenient evaluations and synthetic usefulness for a number of compounds with range of properties for variable applications and vast exploitation, there is need to evaluate his structural and his reactivity. Using quantum chemical calculations these aspects have been evaluated. Molecular electrostatic potential (MEP) and electronic descriptors are used for the chemical reactivity descriptions of molecule. Evaluation of chemical reactivity finds their suitability to react and give new class of compounds suiting pharmacological activities and material applications.
In this work we have analyzed the molecular reactivity of series of bis (trimethyltetrathiafulvalenyl) thiophene through global and local reactivity descriptors derived from the Density Functional Theory (DFT). We consider that this kind of study will contribute to get a better understanding of the chemical behavior of the title compounds.

MATERIALS AND METHODS
All computational calculations have been performed on personal computer using the Gaussian 09W program packages developed by Frisch and coworkers [7]. The Becke's three parameter hybrid functional using the LYP correlation functional (B3LYP), one of the most robust functional of the hybrid family, was herein used for all the calculations, with 6.31G(d,p) basis set [8,9]. Gaussian output files were visualized by means of GAUSSIAN VIEW 05 software [10]. Principal

Molecular geometry
It is generally accepted that the accuracy of theoretical bond lengths and angles largely depends on the functional and basis set used in the calculations. Therefore, DFT calculations using B3LYP functional 6-31G(d,p) as basis sets were selected to simulate the title compounds to obtain accurate molecular structure. The optimized geometry corresponds to the most stable conformation revealed by the lack of imaginary frequency in the vibration mode calculation (figure 1). The calculated bond lengths and angles for the title compounds were collected in tables 1-4, respectively.  I S S N 2 3 2 1 -807X V o l u m e 1 3 N u m b e r 1 J o u r n a l o f A d v a n c e s i n c h e m i s t r y 5940 | P a g e J a n u a r y 2 0 1 7 w w w . c i r w o r l d . c o m

Molecular electrostatic potential (MEP) diagrams
The molecular electrostatic potential, V(r) is related to the electronic density and is a powerful descriptor in determining sites for electrophilic attack, nucleophilic attack and hydrogen-bonding interaction. MEP values can be expressed as the following equation [16]: V o l u m e 1 3 N u m b e r 1 J o u r n a l o f A d v a n c e s i n c h e m i s t r y 5941 | P a g e J a n u a r y 2 0 1 7 where ZA is the charge of nucleus A located at RA , ρ(r') is the electronic density function of the molecule, and r' is the dummy integration variable. Potential decreases in the order blue> green> yellow> orange> red. Electrophilic regions are represented by red, nucleophilic by blue and green indicates neutral electrostatic potential. Negative electrostatic potential corresponds to attraction of the proton by the concentrated electron density in the molecules (from lone pairs, π-bond, etc). Positive electrostatic potential corresponds to repulsion of the proton by the atomic nuclei in regions where low electron density exists and the nuclear charge is incompletely shielded.
As shown in figure 2, large positive region is localized on the C-H bonds, indicating a possible site for nucleophilic attack (blue color, strongest attraction) and large negative region is localized on sulfur atoms indicating electrophilic attack (red color, strongest repulsion). These sites give information about the region from where the compound can have intramolecular (C-H•••S) interactions.

Global reactivity descriptors
The chemical potential (µ), electronegativity (ϰ), hardness (η), and softness (S) are defined as follows [17]: where E is the energy, N the number of electrons, V the external potential of the molecular system under consideration. Evaluation of µ and η faces a practical difficulty due to the discontinuity of the energy E with respect to the variation of N [18,19]. The implication of this discontinuity in conceptual DFT has also been pointed out [20]. One generally makes a finite difference approximation to calculate these quantities. The global hardness is an indicator of the overall stability of the system. A maximum hardness principle (MHP) relating hardness to stability at a constant chemical potential has been proposed by Pearson [21,22]. The proof of the same has been provided by Parr and Chattaraj [23,24]. The chemical potential, µ, is identified as the negative of the electronegativity (ϰ) by Iczkowski and Margrave [25]. Although, these quantities have been defined based on the Kohn -Sham (KS) energy, the same definitions have been routinely used in the canonical molecular orbital theory. The most popularly used formulas for the computation of these quantities use the three point finite difference approximation and express µ and η through the ionization potential (I) and the electron affinity (A).
Another DFT approximation originating from the Koopmans' theorem [3] is I ≈ -EHOMO, A ≈ -ELUMO where εHOMO and εLUMO are the KS one electron eigenvalues associated with the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO), respectively, from approximate DFT calculations on neutral I S S N 2 3 2 1 -807X V o l u m e 1 3 N u m b e r 1 J o u r n a l o f A d v a n c e s i n c h e m i s t r y 5942 | P a g e J a n u a r y 2 0 1 7 w w w . c i r w o r l d . c o m molecules. There have been, in fact, numerous analytical and numerical evidences that softness has a closer link with those properties, especially with polarizability [26][27][28][29] and can thus be used for the hard-soft classification of molecules in a quantitative way.

Figure 3. Highest occupied molecular orbitals and lowest unoccupied molecular orbitals of compounds 1 and 4
As shown in table 5, the molecule 3 is the most molecule has the ability to accept electrons (EHOMO = -4.478eV) while 2 has the highest HOMO energy (EHOMO = -4.470eV) that allows him to be the best electron donor molecule. The high value of the energy gap indicates that the molecule shows high chemical stability, while a small HOMO-LUMO gap means small excitation energies to the manifold of excited states, table 5 shows that compound 2 is the most stable. High ionization energy indicates high stability and chemical inertness and small ionization energy indicates high reactivity of the atoms and molecules. Compound 2 has the lowest ionization potential value (I = 4.470eV) which indicate that it is the best electron donor. The electronic affinity (A) is defined as the energy released when an electron is added to a neutral molecule. A molecule with high (A) values tend to take electrons easily. From table 5 it is clear that Compound 3 is the best reactive. According to the parameters mentioned above, the chemical reactivity varies with the structural of molecules. Chemical hardness (softness) value of compound 3 is lesser (greater) among all the molecules. Thus, compound 3 is found to be more reactive than all the molecules. Also, compound 3 possesses higher electronegativity value than all compounds so; it is the best electron acceptor. The values of ω for compounds 1-4 indicate that they are two series (3,4) and (1,2). The first one has the high value of electrophilicity index which, shows that the compounds of this group are a strong electrophiles.

Local reactivity descriptors
In order to understand the detailed reaction mechanism such as the region-selectivity, apart from the global properties, local reactivity parameters are necessary for differentiating the reactive behavior of atoms forming a molecule. The Fukui function [30] (f) and local softness [31] (s) are two of the most commonly used local reactivity parameters. The Fukui function is primarily associated with the response of the density function of a system to a change in number of electrons I S S N 2 3 2 1 -807X V o l u m e 1 3 N u m b e r 1 J o u r n a l o f A d v a n c e s i n c h e m i s t r y 5943 | P a g e J a n u a r y 2 0 1 7 w w w . c i r w o r l d . c o m (N) under the constraint of a constant external potential [v(r)]. To probe the more global reactivity, indicators in the grand canonical ensemble are often obtained by replacing derivatives with respect to N, by derivatives with respect to the chemical potential µ. As a consequence, in the grand canonical ensemble, the local softness s(r) replaces the Fukui function f(r). Both quantities are thus mutually related and can be written as follows: Accordingly, Fukui function also represents the response of the chemical potential of a system to a change in external potential. As the chemical potential is a measure of the intrinsic acidic or base strength, and the local softness incorporates the global reactivity, both parameters provide us with a pair of indices to probe for example, the specific sites of interaction between two reagents. Generally, it is demonstrated that the larger the value of the Fukui function, the greater the reactivity of the corresponding site. Once again, due to the discontinuity of the electron density with respect to N, finite difference approximation leads to three types of Fukui function for a system, namely (1) ) (r f  for nucleophilic attack measured by the electron density change following addition of an electron, (2) ) (r f  for electrophilic attack measured by the electron density change upon removal of an electron, and (3) ) ( 0 r f for radical attack approximated as the average of both previous terms. They are defined as follows: Using one-electron orbital picture, Fukui functions can be approximately defined as: These relations highlight the fact that the formalism of DFT-based chemical reactivity built by Parr and coworkers, captures the essence of the pre DFT formulation of reactivity under frontier molecular orbital theory (FMO). Berkowitz showed that similar to FMO, DFT could also explain the orientation or stereo-selectivity of a reaction [32]. In addition, DFTbased reactivity parameters are augmented by more global terms expressed in the softness. For studying reactivity at the atomic level, however, a more convenient way of calculating the f(r) functions at atomic resolution is used. The condensed-to-atom Fukui functions for an atom k in a molecule are expressed as [33]: Where K q is the electronic population of atom k in the molecule under consideration.

Principal Component Analysis (PCA)
In this work, we auto scaled all calculated variables in order to compare them in the same scale. Afterwards, PCA (principal component analysis) was used to reduce the number of variables and select the most relevant ones, i.e. those responsible for the tetrathiafulvalenes derivatives reactivity. After performing many tests, a good separation is obtained between more active and less active tetrathiafulvalenes compounds using ten variables: I, A, χ , ɳ, s, µ, ω, EHOMO , ELUMO , ΔEgap (see tables 5 and 6).
We can observe from PCA results that the first three principal components (PC1, PC2 and PC3) describe 99.92% of the overall variance as follows: PC1 = 98.89%, PC2 = 1.01% and PC3 = 0.02%. The score plot of the variances is a reliable representation of the spatial distribution of the points for the data set studied after explaining almost all of the variances by the first two PCs. The most informative score plot is presented in figure 4 (PC1 versus PC2) and we can see that PC1 alone is responsible for the separation between more active (3 and 4) and less active compounds (1 and 2) where PC1>0 for the more active compounds and PC1<0 for the less active ones. The same results follow in the case of global reactivity trend based on ω.

Figure 4. Score plot for bis (trimethyltetrathiafulvalenyl) thiophenes 1-4 in gas phase
The loading vectors for the first two principal components (PC1 and PC2) are displayed in figure 5. We can see that more active compounds (PC1 ˃ 0) can be obtained when we have higher A, I, S, χ, ω, values. In this way, some important features on the more active compounds can be observed.  Figure 6 shows HCA analysis of the current study. The horizontal lines represent the compounds and the vertical lines the similarity values between pairs of compounds, a compound and a group of compounds and among groups of compounds. We can note that HCA results are very similar to those obtained with the PCA analysis, i.e. the compounds studied were grouped into two categories: more actives compounds 2 and 3 and less active compounds 1 and 2.

CONCLUSION
From the whole of the results presented in this contribution it has been clearly demonstrated that the sites of interaction of the title compounds can be predicted by using DFT-based reactivity descriptors such as the hardness, softness, and electrophilicity, as well as Fukui-function calculations. Compound 3 is found to be more reactive than all the molecules and the parameters of local reactivity descriptors show that 4C is more reactive site for nucleophilic and free radical attacks and 11C for electrophilic attack.