Evaluation of (Z)-5-(Azulen-1-ylmethylene)-2- thioxothiazolidin-4-ones Properties Using Quantum Mechanical Calculations

Derivatives of (Z)-5-(azulen-1-ylmethylene)-2-thioxothiazolidin-4-one are reported as heavy metal (HM) ligands in heterogeneous systems based on chemically modified electrodes. Their ability to coordinate HMs ions has recently been shown to be very selective. In this context, an additional computer-assisted study of their structure was performed using density functional theory (DFT) to achieve a complex structural analysis. Specific molecular descriptors and properties related to their reactivity and electrochemical behaviour were calculated. The correlation between certain quantum parameters associated with the general chemical reactivity and the complexing properties of the modified electrodes based on these ligands was carried out to facilitate the design of molecular sensors. Good linear correlations between DFT-calculated HOMO/LUMO energies and experimental redox potentials were found. A good agreement between the chemical shifts predicted by the DFT method and those determined experimentally from NMR data for these ligands demonstrated the accuracy of the calculations to assess the structural data. Such a computational approach can be used to evaluate other properties, such as electrochemical properties for similar azulene derivatives.


Introduction
Previous studies correlate the electrochemical properties for various organic compounds with structural parameters using density functional theory (DFT) calculations to achieve the rational design of new materials with improved electrochemical properties [1][2][3][4][5][6][7][8]. They are based on the link between the energy levels corresponding to the highest occupied molecular orbital (HOMO) or the lowest unoccupied molecular orbital (LUMO) and the electrochemical oxidation and reduction potentials, respectively [9][10][11]. Strong linear correlations of DFT-calculated HOMO/LUMO energies using B3LYP/6-31G(d) functional [12,13] and experimental redox potentials were found for polycyclic aromatic hydrocarbons by D. Méndez-Hernández and co-workers [14], highlighting the idea that quick, accurate and low-cost predictions using the B3LYP/6-31G(d) functional represent a reliable approach to apply on other molecules to evaluate their electrochemical properties. A comparative computational study using different density functionals on The investigated ligand structures contain a part of rhodanine (III), known for its HMs complexing properties [21]. Azulene derivatives of rhodanine can be used in the precise determination of metals, in the same way that p-dimethylaminophenylenrhodanine was used for the precise determination of Cu, Ni, Fe and Zn ions [22], or in the same way that triarylamine rhodanine derivatives were used as colorimetric sensors for the detection of Ag (I) and Hg (II) ions [23].
The alkyl groups induce a + I effect and increase the electron densities of the molecules. Consequently, the alkyl-substituted compounds are expected to be easier oxidized and harder to be reduced than pattern compounds (T3). However, the steric effect of these groups makes the reductions more difficult to anticipate.
The azulene system is more stable when it is symmetrically substituted, but becomes very reactive by unsymmetrical substitution when the symmetry of the aromatic system is disrupted by the difference in the alkyl groups electronic influence. Therefore, it is expected that T1-unsymmetrically substituted will react faster than T2-symmetrically substituted with 4,6,8-Me 3 despite its higher volume.
Previous studies have highlighted the capacity of azulenes to be immobilized on electrodes through electropolymerization processes if 1,3 positions are free [24], or by π-stacking [25]. For sensor applications, this last way is more advantageous [25], because the number of grafted complexing units on the electrode surface is increased, leading to a better response. T1, T2, and T3 ligands have push-pull structures that can be used to build modified electrodes with complexing properties. These compounds were electrochemically characterized, and their metal-binding properties were reported [16]. Even if their structures are quite similar, different complexing abilities have been found for the modified electrodes based on them. These experimental results led us to look for structural reasons to explain why the behaviour is so different. The results of DFT calculations are reported here and several calculated properties were compared with the experimental ones. Connections were made with the experimental results obtained in the use of the modified electrodes based on these compounds in HM analysis. Table 1 gives several properties of T1-T3, and the experimental features which were found to be different in their use in the modified electrodes preparation/characterization, and utilisation to analyse HMs ions from aqueous solutions. The monomers (mainly characterized in Table 1, lines 1-4) were successfully deposited on glassy carbon electrodes through direct electropolymerization at anodic potentials in millimolar solutions of each ligand. The CMEs were tested then in ferrocene solutions in acetonitrile containing 0.1 M tetrabutylammonium perchlorate (TBAP) as supporting electrolyte (Table 1, lines 5 and 6). The films formations occurred in different ranges of monomer concentrations (line 5). Changes of the ferrocene redox signal of different amplitudes were found (line 6). These materials were tested vs. HMs (Table 1, lines 7-11). They were able to complex HMs ions from aqueous solutions [21,27], but the parameters for their analysis were very different (line 7). T1-modified electrodes have quite high detection limits (10 −6 M for Pb, Cd, Cu and 10 −4 M for Hg). T2-modified electrodes showed an intermediate behaviour, and their detection limits are of about 10 −7 M. T3modified electrodes showed the highest complexing ability leading to attractive values for the lowest detection limits (for instance DL Pb < 10 −8 M).
Lines 5 and 6 in Table 1 collect the characteristics of the film formation and characterization. It is expected that the electrochemical polymerization of azulene ligands on the electrode surface is intensified by the presence of alkyl groups, grafted on the azulene moiety (these effects are seen in the values of first oxidation potential Ea1). However, T1 hardly formed films, while T3 easily formed films. T2 showed intermediate behaviour. The film formation ability seems to be in accord with the analytical performance [17,28].
It is difficult to explain/predict the behaviour only from the analysis of the inductive effects of the substituents on the studied structures [29], because the investigated facts are quite complex involving several processes, as seen in Table 1 (ligand polymerization, CME complexation, etc.). However, the ligand structure is decisive in all these steps, and the rationalization of these structural effects based on quantum mechanical calculations is found as a favourable approach. Similar calculations for other structures are of great use in reducing the number of experiments, which has a major economic impact in finding a way to anticipate the optimal structures.

Computational Procedure Details
The calculations were carried out using Spartan 14 software Wavefunction, Inc. Irvine CA, USA [30]. For the equilibrium geometry at ground state in a vacuum, a series of calculations of molecular properties and quantum chemical parameters was done using Density Functional Theory [31], software algorithm B3LYP method (the Becke's threeparameter hybrid exchange functional with the Lee-Yang-Parr correlation functional [32] and polarization basis set 6-31G (d, p) [33] and wB97XD/6-311++G (d, p) density functional model, stated previously as an appreciable improvement over other empirical dispersioncorrected density functionals [34]. The ab initio calculation of NMR chemical shifts was achieved with a gauge including atomic orbitals (GIAO) [35]. The density functionals models employed by Spartan software used empirically corrected 13 C chemical shifts to reduce the errors in comparison with uncorrected ones [36].

Results
The first step of calculation concerns the generation of 3D structures of (Z)-5-(azulen-1-ylmethylene)-2-thioxothiazolidin-4-one derivatives and their geometry optimization. This step aimed to establish the most stable conformer of each derivative, which has the energy minimum among its conformers. The atomic numbering schemes for the optimized geometries of the analysed structures, arbitrarily assigned by Spartan Software, are illustrated in Figure S1 from Supplementary Material. The calculated bond lengths at ground state (Table S2) and the predicted angles and dihedrals angles for the studied structures (Table S3) confirm the electronic effects of alkyl substituents on the azulene structure for the optimized conformers of T1-T3 molecules. Further, the calculations were made for each structure for these lowest energy conformers.

Predicted and Experimental NMR Chemical Shifts
The ab initio calculation of 1 H-and 13 C-NMR chemical shifts was achieved. Correlations between predicted and experimental (found in [29]) chemical shifts for T1-T3 were performed to check the prediction's accuracy for the investigated structures. They are listed in Tables 2-7. The structures of each compound were introduced in the table's heads using different notations that came from the computational algorithm generated by the software (denoted A1, A2, A3) and from the experimental data [29], according to IUPAC notations (denoted B1, B2, B3) to be in correspondence.         Table 8 lists the results of DFT computations, describing molecular and quantitative structure-activity relationships (QSAR) properties of the investigated ligands using B3LYP/6-31G (d, p) and wB97XD/6-311++G (d, p) density functional models.
The frontier molecular orbitals density distribution obtained by calculation for the studied structures is shown in Figure 2, along with the energy diagram and the gaps between the HOMO and LUMO (∆E). They were calculated using the wB97XD/6-311++G (d,p) density functional model. For B3LYP/6-31G (d,p) model, they are given in Figure S2.  Starting from E HOMO and E LUMO energies given in Table 8, other related quantum descriptors resulted. Table 9 gives the formulas and calculated values for: energy gap (∆E), ionization potential (I), electron affinity (A), electronegativity (χ), global hardness (η), global softness (σ) [37,38], and global electrophilicity index (ω) [39] using B3LYP/6-31G (d,p) and wB97XD/6-311++G (d,p) density functional models. I and A were calculated by applying relationships proposed by Koopmans [40,41], respectively. The softness (σ) and hardness (η) descriptors derived from Pearson's Hard and Soft Acids and Bases Principle (HSAB) [42] and Maximum Hardness Principle (MHP) [43] describe the electronic reactivity, and the response to electronic perturbations, respectively [44]. The global electrophilicity index (ω), as defined by Parr R.G and co-workers [39], is a measure of the reactivity of chemical species in different environments (solvent or biological systems).
The electrostatic potential maps are shown in Figure 3. The colour indicates the value of the electrostatic potential. Red areas suggest negative potentials, colour toward blue designate regions of positive potential. The potential increases in the order: red < orange < yellow < green < blue. Red regions suggest the potential sites for HMs ions complexation, where positive charges are most susceptible to be attracted. Figure 4 gives the maps of local ionization potential obtained using wB97XD/6-311++G (d,p) density functional model. They reveal the regions from which electrons are most easily removed indicating the most susceptible sites to electrophilic attack.
In LUMO maps representation in Figure 5, the absolute value of LUMO is mapped onto an electron density surface (blue colour for large values of LUMO and red colour for small values), allowing to anticipate regions subjected to nucleophilic reactivity.

Correlation between DFT-Calculated Frontier Molecular Orbital's Energies and Experimental Data
The calculated HOMO and LUMO energies (Table 8 lines 6 and 7, respectively) were plotted against experimental oxidation and reduction potentials ( Table 1, lines 2 and 3, respectively), and are shown in Figures S3 and S4. Linear relationships obtained for reduction and oxidation potentials, and their corresponding equations, are given in Table 10, using either B3LYP or wB97XD hybrid functions.

NMR Predictions
The ab initio calculation of NMR chemical shifts listed in Tables 2-7 shows reasonably good results for predicted data in comparison with the experimental ones (taken from [29]). Good agreement between experimental and calculated data concerning 1 H-NMR and 13 C-NMR for T1-T3 structures, except for H4 atom (the hydrogen from -NH of the cycle III) is found. The observed gaps are explained by the proton enolization due to the intermolecular proton transfers from -NH (-SH or -OH). Additionally, small deviations from experimental shifts can occur due to the presence of solvent.

Highlighting Results by Using Different Functional Models
The calculated values given in Table 8

Highlighting Variation of Properties among Investigated Structures
By examining the values of the properties given in Table 8 and taking into account the changes in the compound's structures, the main differences in magnitude are for the calculated dipole moment (line 5), ovality index (line 11), polarizability (line 12), LogP (line 12), and minimum electrostatic potential values (line 16). The dipole moment and the polarizability decrease in magnitude as expected, ranging in the order: T1 > T2 > T3. It is higher for the structure having more alkyl substituents. The ovality index varies in the same order, as the linearity of structures decreases with substitution, being related to molecular surface (line 8) and van der Waals volumes (line 9) which also decrease in the same order. The partition coefficient LogP is the parameter that varies the most for these structures in the same order. The obtained values suggest a hydrophobic character of the structures and confirm their low aqueous solubility ( Table 1, line 4), especially for the ligands substituted with iso-propyl and methyl groups (T3 and T2). However, logP values are not so high due to the contribution of the rhodanine moiety. This moiety is favourable for the complexing properties of all ligands, and also for the CMEs based on such ligands. LogP values could be related to the complexing property which was found for CMEs (Table 1, line 7). The best values are obtained for CMEs based on T3. MinElPot calculated form wB97XD method for the investigated structures decreases in absolute values in the same order. More comments about this parameter are given forward, connected with the discussion of Figure 3.
The other related quantum descriptors for the studied structures resulted in Table 8 such as energy (line 2) in absolute value, energy aq. (line 3) in absolute value, solvation energy (line 4), area (line 8) and volume (line 9), decrease in the order: T1 > T2 > T3, as expected, according to their molecular weight and azulene substituents (methyl and isopropyl). PSA (line 10) value presents insignificant variance for the three investigated structures, suggesting there is no distinction caused by polar substituents or significant disturbing electron distribution grafted on the skeleton. The major contribution in the polar character of the compounds is given by polar functional groups or atoms which are the heteroatoms present in their same complexing unit (the rhodanine). Consequently, PSA is quite similar for the investigated compounds. It is not significantly influenced by the methyl and isopropyl substituents of the azulene moiety, as shown by B3LYP/6-31G (d, p) calculation. However, the values calculated by the wB97XD/6-311++(d,p) model are higher for T1 and T2 than T3, putting in evidence the effect of the alkyl groups on the structure, and confirming the last model is better.
From Figure 2, a quite similar distribution of the frontier molecular orbitals can be observed for the three compounds. The frontier orbitals gap (∆E) which characterizes the chemical reactivity of each molecule, are different for these compounds. The higher value of ∆E gap is found for T2, and it reflects its higher kinetic stability [45,46] expected due to symmetry reasons.
Taking into account our interest in complexing the HMs ions by these ligands, the donor-acceptor interactions were examined. They can occur between the π-electrons of the 2-thioxothiazolidin-4-one and the vacant d-orbital of the metal. As illustrated in Figure 2, HOMO orbitals are preferentially distributed on oxygen and sulphur atoms of the rhodanine cycles. This is in good accordance with their ability to donate their vacant electrons and it is confirmed by the localization of the negative (red and orange regions) electrostatic potential (Figure 3). The complexation unit being the same, 2-thioxothiazolidin-4-one (rhodanine), the ligand's complexing capacity does not differ significantly, which is suggested by the small variance of the electrostatic potential (MinElPot). That observation leads to the assumption that a good complexation can be achieved for all the analysed compounds.
From the graphical representation in Figure 3 it is obvious that for all the investigated structures, the negative areas are mainly localized over the oxygen atoms (red regions). They correspond to the maximum negative values of potential (MinElPot- Table 8

Correlation of Molecular and QSAR Properties
Attempts to find correlations between calculated molecular and QSAR properties (from Table 8) and the electron affinity (A) or ionization potential (I), respectively, are shown in Table 11. Linear relationships were considered for all parameters. The correlation coefficients of the linear dependencies (R 2 ) were of much help to establish correct connections. For instance, in terms of total energy for the correlation with I or A, R 2 is higher for I (0.995) than for A (0.613) when using the wB97XD model. These values are higher than the corresponding ones obtained by the B3LYP method (0.846 and 0.684, respectively). For energy aq. R 2 is 0.613 for the correlation with A and 0.990 with I when using the wB97XD model. For the B3LYP method, the values are lower, being 0.684 and 0.846 respectively. This protocol was followed for all the other properties from Table 8. The best correlations for the parameters obtained by the two methods are collected in Table 12 which gives the most confident dependencies and the method which led to such results. All the best correlations (R 2 over 0.9) were obtained with I, with one exception: the reduction potential which is correlated, as expected, with A. From the 13 best correlations, 11 were obtained using the wB97XD model, which recommends the method as being more performant. Table 11. "y = a + b·x" linear correlations between predicted molecular and QSAR properties and A or I, and their correlation coefficients (R 2 ) for T1-T3 calculated using B3LYP/6-31G (d,p) and wB97XD/6-311++G (d,p) density functional models, respectively, according to  Table 12. "y = a + b·x" linear correlations between predicted molecular and QSAR properties and A or I for T1-T3 calculated using the density functional models which gave us the best correlations coefficients (R 2 ); A and I were expressed in eV.

Correlation of Quantum Chemical Reactivity Parameters
The main quantum parameters of the studied ligands calculated according to B3LYP/6-31G (d, p) and wB97XD/6-311++(d,p) models and given in Table 9 show that I (line 1), χ (line 4) and ω (line 7) range in the order T1 < T2 < T3, regardless the chosen algorithm. This trend is in agreement with experiments that show the same order of complexation for the studied ligands. T3 presents the highest values of these parameters, meaning that it has the highest total energy gain upon saturation with electrons, comparing to the other two ligands and explaining the experimental facts.
A (line 2), ∆E (line 3), η (line 5) and σ (line 6) do not show the same regular variation. The highest value for A is for T3; for T1 and T2, the values for A are close. The highest values for ∆E and η are for T2, being close to those for T3. The highest value for σ is for T1 being close to those of T2 and T3. ∆E of T3 is relatively smaller than T2. This fact indicates T3 is more reactive than T2. Indeed, T3 is found to be more able to interact in the complexation reaction of HMs (Table 1, lines 7-11). According to its lower ∆E, the most reactive compound among all seems to be T1. However, the experiments (Table 1, lines 7-11) show that T1 has the lowest complexation ability. This discrepancy can be attributed to the lowest values for A and I for this ligand. Table 13 shows the linear relationships considered for all quantum chemical reactivity parameters from Table 9, vs. A or I, respectively, as shown in the case of molecular and QSAR properties. The correlation coefficients of the linear dependencies (R 2 ) were used to establish correct connections. For instance, in terms of χ correlation with I or A, R 2 is higher for I (0.995) than for A (0.614) when using the wB97XD model. For ω, a very good correlation was found with A, and the best correlation coefficient is obtained through the B3LYP model. η and σ are not linearly correlated with A and I. The best correlations for the parameters obtained by the two methods are collected in Table 14 which gives the most confident dependencies and the method which led to such results. Table 13. "y = a + b·x" linear correlations between predicted quantum chemical reactivity parameters and A or I (all expressed in eV) for T1-T3 calculated using B3LYP/6-31G (d,p) and wB97XD/6-311++G (d,p) density functional models.  The graphical quantities from Figures 3-5 provide a visual representation of the chemically active sites, and allow the comparison of the local reactivity sites of the analysed structures. The molecular electrostatic potential, previously discussed (Section 4.2.2), is useful to identify the reactive sites for complexing interactions, and to understand the chemical recognition process based on this type of ligands.
The ionization potential map ( Figure 4) is particularly useful to assess chemical reactivity and selectivity, in terms of electrophilic reactions. The blue colour reveals the regions where ionization is relatively difficult. In these regions, localized over -NH group of Cycles III, the values of ionization potential vary as follows: 13.78 ÷ 14.83 (for T1), 13.84 ÷ 14.92 eV (for T2), and 13.84 ÷ 14.90 eV (for T3). The orange areas correspond to the lowest ionization potentials (most accessible to electrophiles), localized over sulphur atoms from Cycles III. These sites present the following values: 7.06 ÷ 7.36 eV (for T1), 7.08 ÷7.44 eV (for T2) and 7.21 ÷ 7.47 eV (for T3). There are no red or orange areas on the local ionization potential maps, indicating there are no clear sites for electrophilic attack. This means the ligands were properly selected for complexation (which involves the nucleophilic attack of HMs ions).
The LUMO map ( Figure 5) indicates nucleophilic reactivity. It can be observed that the colours are toward red, suggesting small values (near zero) of the LUMOs. Consequently, these ligands are not very susceptible to nucleophilic attack.

Correlation between DFT-Calculated Frontier Molecular Orbital's Energies and Experimental
Oxidation and Reduction Potentials E HOMO and E LUMO predicted chemical parameters (Table 8) were correlated with experimental electrochemical properties resulted from the ligand characterization (Table 1) in order to establish the best ligand to be used for the complexation of HMs ions. Linear relationships obtained using B3LYP and wB97XD hybrid functions for oxidation and reduction potentials (Table 10) are illustrated in Figures S3 and S4. R 2 values indicate satisfactory correlations between the calculated and the experimental values. The calculated HOMO orbital energies vary in the same order as the experimental values of the first anodic peak potentials, namely T3 > T2 > T1. The same behaviour was observed for the reduction capacity expressed as Ec1 experimental values, which varies in the opposite direction. So, the redox potential is influenced by the number and position of the alkyl groups, as assumed previously in [21,26,27]. Thus, the evaluation of the oxidation capacity of the investigated azulene systems is in good agreement with previously reported electrochemical data.
In Figure S5 and Table S3, the Mulliken population analysis based on the local electron density reveals differences in charge values on the heteroatoms of the rhodanine cycles. Thus, O sp 2 (red colour) and S sp 2 (orange colour) atoms have negative Mulliken charges, while S sp 3 (yellow colour) shows positive values. Remarkable are changes in the charge of C atoms (blue colour) of the azulene moiety unsubstituted (negative Mulliken charges at T3) and substituted with alkyl groups (at T1 and T2); the substitution leads to positive charges (see differences in charges at C2, C5 and C7 of T1, C4, C5 and C9 of T2, as atom labelling scheme illustrated in Figure S1 (their corresponding Mulliken charges are listed in Table S3).

Conclusions
Quantum chemical calculations for three ligands derivatives of (Z)-5-(azulen-1-ylmeth ylene)-2-thioxothiazolidin-4-one are reported. They are structurally distinguished by the substitution of azulene cycle, respectively, by methyl and isopropyl groups: 3,8-Me 2 -5-iPr (T1), 4,6,8-Me 3 (T2) and H (T3). The chemical calculations resulted in a series of atomic (bond lengths, angles, Mulliken charges) and molecular descriptors particularly valuable in quantitative structure-activity relationships analysis and NMR spectra. These predicted chemical parameters were correlated with the experimental electrochemical characterization, in order to establish the best parameters for a ligand to be used for the complexation of HMs ions. Thus, the results of the evaluation of the oxidation capacity of the investigated azulene systems are in good agreement with previously reported electrochemical data. The calculated HOMO orbital energies vary in the same order as the experimental values of the first anodic peak potentials, namely T3 > T2 > T1. The same behaviour was observed for the reduction capacity which varies in the opposite direction, concluding that the redox potential is influenced by the number and position of the alkyl groups, in accordance with previous assumptions when designing the ligands. Linear correlations of DFT-calculated frontier molecular orbitals' energies and the experimental oxidation and reduction potentials were found. The computer-aided study turned out to be a complex structural approach, being an alternative to find parameters that matter when designing new ligands.