DFT calculations and electrochemical studies on azulene ligands for heavy metal ions detection using chemically modified electrodes

A computational study on three related derivatives of 5-[(azulen-1-yl)methylene]-2thioxoimidazolidin-4-one was conducted using density functional theory by calculating a series of molecular descriptors and properties of their optimized geometries (electrostatic and local ionization potentials, molecular frontier orbitals, etc.). Thermodynamic properties (zero-point energy, enthalpy, constant volume heat capacity, entropy and Gibbs energy) for these derivatives have been calculated and related to ligands electrochemical behavior. Reduction and oxidation potentials have been correlated to their calculated energy levels for LUMO and HOMO orbitals. Chemically modified electrodes based on these derivatives have been tested in view of heavy metal ions recognition, and their detection limits have been correlated to the calculated values of electron affinity.


Introduction
Azulenes have a push-pull structure, with a five-member cyclic moiety (electron-rich) connected to a seven-member cyclic moiety (electron-poor).Azulene derivatives undergo irreversible electrooxidations and irreversible or quasi-reversible reductions when polarized at anodic or cathodic potentials, respectively [1].Due to their low oxidation potential, they can lead to modified electrodes which could be used to build new electrochemical sensors.The sensors based on azulene derivatives are the basis of most recent methods for determination of heavy metals from water samples [2].This method leads to comparable results to the best-known methods known to date, such as those based on iron oxide/graphene composite [3], bismuth nanoparticle-porous carbon paste [4], etc.
The azulenes investigated in this study are derivatives of 5-[(azulen-1-yl)methylene]-2-thioxoimidazolidin-4-one (L1), shown in Fig. 1, which have in their structure complexing imidazolidine group, specific for the recognition of heavy metal ions [5].One of these ligands, L2, has been investigated in molecular recognition for building complexing chemically modified electrodes for heavy metals, and lead to promising results for the heterogeneous recognition of Pb(II) and Cu(II) ions [6].Similar azulenes with the thiourea moieties as complexing monomers have led to modified electrodes with good complexing abilities for heavy metals, with very good responses for Cd(II) and Pb(II) ions [7].
The present work reports on the investigation of L1 -L3 compounds ability to coordinate heavy metals.Chemically modified electrodes have been prepared from these ligands under similar conditions and used for heterogeneous recognition.Their detection limits for some heavy metal ions (Hg, Pb, Cd, Cu) have been evaluated in order to be correlated to their molecular structural descriptors.Density functional theory (DFT) investigations on these azulenes derivatives have been performed to reveal some global reactivity parameters derived from frontier molecular orbitals energy diagram, to be linked with their ability to complex heavy metals, for electrochemistry applications.Predicted structural characteristics are optimized, compared and correlated in terms of NMR chemical shifts, aiming to determine the structure thermodynamically favored by synthesis [5].Then, their computed parameters are obtained and analyzed in connection with some complexing and redox properties determined by electrochemistry.Our interest concerns the use of L1 -L3 as ligands in heterogeneous recognition of heavy metals (Fig. 1).The main parts of their structure are illustrated with different colors in Fig. 1.They are the seven-member azulene ring (I), the five-member azulene ring (II) and the imidazolidine ring (III).

Experimental
Acetonitrile (CH 3 CN) and tetrabutylammonium perchlorate (TBAP) from Fluka were used as received for the solvent and supporting electrolyte.The azulene derivatives were synthesized according to the previous described method [5].
The electrochemical characterizations were carried out using a PGSTAT 12 AUTOLAB potentiostat connected to a three-electrode cell.The working electrode was a glassy carbon disk (3 mm diameter).The active surface was polished before each determination with diamond paste (1 μm) and cleaned with acetonitrile.Ag/AgNO 3 (10 mM) in 0.1 M TBAP, CH 3 CN was used as reference electrode and a platinum wire was used as auxiliary electrode.The electrochemical experiments were performed by cyclic voltammetry (CV), differential pulse voltammetry (DPV) and rotating disk electrode (RDE).CV curves were recorded from 0.1 V/s to 1.0 V/s (for different scan rate).DPV curves were recorded at 0.01 V/s with a pulse height of 0.025 V and a step time of 0.2 s.RDE curves were recorded at 0.01 V/s, for various rotation rates of the disk (between 100 and 2000 rpm).The potential was referred to the potential of the ferrocene/ferrocenium redox couple (Fc/Fc + ) which in our experimental conditions was +0.07 V.The preparation of the modified electrodes has been done by controlled potential electrolysis in millimolar solutions of each ligand, in acetonitrile solutions containing TBAP (0.1 M) as supporting electrolyte.All electrochemical experiments were performed at 25 °C under argon atmosphere.
The experiments for heavy metals ions detection were carried out in a three-electrode cell (transfer cell).The working electrode was a glassy carbon disk (diameter of 3 mm) electrode modified with each ligand, the reference electrode was Ag/AgCl, 3 M KCI, and the auxiliary electrode was a platinum wire.0.1 M buffer acetate (pH 5.5) was used as supporting electrolyte.The recognition experiments were performed at 25 °C and under argon atmosphere.The buffer acetate solution was prepared from 0.2 M acetic acid and 0.2 M sodium acetate solutions.Metal cation solutions were purchased from Sigma Aldrich and Fluka: mercury (II) acetate -Sigma Aldrich, cadmium nitrate tetrahydrate -Sigma Aldrich, copper (II) acetate monohydrate -Fluka, lead (II) nitrate -Sigma Aldrich.Heavy metals solutions with concentrations of 10 −4 and 10 −6 mol/L, were prepared by successive dilutions from 10 -2 M stock solutions of mercury (Hg), cadmium (Cd), copper (Cu) and lead (Pb).
The procedure to modify the glassy carbon electrodes consisted in electropolymerization of each ligand in millimolar solutions in 0.1 M TBAP, CH 3 CN, followed by cleaning in acetonitrile, equilibration of the modified electrode in 0.1 M acetate buffer at pH 5.5, between -0.9 V and + 0.6 V vs. Ag/AgCl, over oxidation in 0.1 M acetate buffer at pH 5.5 between -0.2 V to + 1.2 V.
The accumulation of heavy metal ions has been performed in the assay solutions containing heavy metal ions in different concentrations (starting from 10 -10 to 10 -4 ) for 20 minutes under magnetic stirring.This step has been followed by the stripping of the accumulated metal ions, after a reduction at -1.2 V for 120 s, in CV or DPV anodic scans with a scan rate of 0.01 V/s.The stripping currents for Cd, Pb, Cu, and Hg currents have been evaluated for each concentration.

Computational procedure details
The calculations were carried out using Spartan 14 software Wavefunction, Inc. Irvine CA USA [8] on Intel(R) Core i5 at 3. 2 Ghz CPU PC.First, the 3D CPK models of the studied compounds were generated.Then, a systematic conformational search and analysis was achieved to establish the more stable conformers of each of the three derivatives of 5-[(azulen-1-yl) methylene]-2-thioxoimidazolidin-4-one, presenting the energy minima.The most stable conformer (lowest energy) was established using MMFF (Molecular Mechanics Force Field for organic molecules developed by Merck Pharmaceuticals) model by refining the geometry for each studied compound.For these refined structures, a series of calculations of molecular properties and topological descriptors has been carried out using [9], software algorithm hybrid B3LYP model (the Becke's three parameter hybrid exchange functional with the Lee-Yang-Parr correlation functional) [9][10][11] and polarization basis set 6-31G* [8,12] in vacuum, for equilibrium geometry at ground state.

Quantum mechanical calculations
The atomic numbering scheme given by the Spartan 14 software for structures under study is shown at the beginning of the Table 1.The parameters (bond lengths and dihedral angles) of the compounds L1-L3 have been obtained for the optimized geometries.They are listed in Tables 1-2.The molecular structures comprise two coplanar cycles (I and II) and the cycle III, which can be rotated around C14-C13 bond, showing a small deviation from co-planarity given by the small values of the dihedrals (H7, C14, C13, C12) and (C3, C14, C13, N2).L3 conformer configuration presents also the torsion angles of 119.04° (C18, C15, C7, C4) and 61.03° (C16, C15, C7, C9).In Table 3 are given the results of quantum chemical calculations for a) valuable predicted molecular properties and b) descriptors for L1 -L3 compounds, which provide information to conduct quantitative structure-property relationships (QSPR) for Z isomers.
From Table 3 it can be observed that area and volume for investigated compounds increase in the order: L1< L2 < L3, as expected, correlated with the molecular weight and presence of methyl and isopropyl substituents.The ovality index shows the same trend, suggesting the deviation of molecules from the spherical shape, considering the minimum surface for the spherical shape.This parameter is related to molecular surface area and van der Waals volume and increases with the linearity of the structures in the same order as above-mentioned descriptors (area, volume).The polarizability increases in the order: L1 < L2 < L3, showing increased induction (polarization) interactions, resulting from an ion or a dipole inducing a temporary dipole in an adjacent molecule.The same tendency for the dipole moment was observed, as expected.Hydrogen bond donor (HBD) and acceptor (HBA) counts have the same values for all structures under study, due to the presence of the same thioxoimidazolidine cycle III, containing the same donors (nitrogen) and acceptors (sulphur and oxygen).The octanol-water partition coefficient (log P) and polar surface area (PSA) are also evaluated and listed in Table 3. Log P values are positive, showing their good affinity for hydrophobic phases, and increase in the order: L1 < L2 < L3.L1 had also the smallest PSA (polar surface area) value reinforcing the trend of log P; PSA values for L2 and L3 are quite similar.

HOMO -LUMO analysis
The molecular frontier orbitals (the highest occupied molecular orbital -HOMO and the lowest unoccupied molecular orbital -LUMO) energies of the studied compounds (in Z form) have been calculated (Table 4) and their density distributions are given in the bottom (HOMO) and in the top (LUMO) of Fig. 2, respectively.Analysis of Fig. 2 shows the distribution of HOMO orbitals for all studied compounds is delocalized as follows: over atoms C7 and C8 from cycle I, over bonds C1-C3 and C0-C2 from cycle II, over C13-C14 bond (between cycle II and cycle III) and over heteroatoms N2, O1 and S1 form cycle III.For all compounds, LUMO orbitals (Fig. 2) are entirely delocalized over cycle I and cycle II, while for L2 structure, they are distributed supplementary over functional groups C=O (O1), C=S (S1) and C12-C13 bond from cycle III.The other occupied orbitals energy levels are listed in Table 5, and they can be related with UV Vis allowed transitions.HOMO and LUMO energies from Table 4 have been examined in order to find the tendency of each ligand to donate electrons (to empty molecular orbitals with low energy of convenient molecules), or to accept electrons, respectively.The resulting band gap E (E HOMO -E LUMO ) can be used to provide useful information on the chemical reactivity and kinetic stability of each ligand [13,14].The band gap between the HOMO and LUMO orbitals energy has been found to increase in the order: L3 < L2 < L1 (Table 4).Consequently, according the obtained values for E, L1 presents the lowest reactivity, followed by L2 and L3 (the most reactive).Starting from HOMO and LUMO energy (Fig. 2), other related descriptors have been calculated: ionization potential (I), electron affinity (A), electronegativity (χ), chemical hardness (η), global softness (σ) [15,16], chemical potential (µ), global electrophilicity index (ω), and also ionization potential (defined as I = -E HOMO ) and electron affinity (defined as A = -E LUMO ), by applying Koopmans' theorem [17,18].Derived calculated quantum chemical parameters for the most stable conformer (Z) of each studied compound are listed in Table 4.The parameters from Table 4 are important global chemical reactivity descriptors of the studied molecules in terms of their reactivity and site selectivity.L1 has the highest ionization potential and electron affinity, which reflects his superior capability to interact with heavy metal cations and to accept one electron from a donor.In terms of chemical hardness, large E signifies hard molecules and small E refers to soft molecules.In this context, L1 presents the hardest structure.The values obtained for the electrophilicity index (ω) are not in the same order.It indicates an unexpected global electrophilic nature of these molecules, as measure of energy lowering due to maximal electron flow between donor and acceptor [19].

Molecular electrostatic potential and local ionization potential
Molecular electrostatic potential (MEP) is a valuable descriptor for overall molecular charge distribution and for the way of electrophilic addition assessment [20,21].Graphical quantities resulted from quantum chemical calculations for MEP and local ionization potential (LIP) are given in Fig. 3.These diagrams provide a visual representation of the chemically active site and comparative reactivity of atoms.In MEP graphical representation from Fig. 3, the red surface ( system) indicates a negative potential, suggesting that this region is attracted to a positive charge; the blue surface ( system) indicates a positive potential, thus, this region is attracted to negative charges.Potential increases in the order: red < orange < yellow < green < blue.For all the studied compounds, the negative area is mainly localized over the oxygen (red region) and less on sulphur atoms (orange region).The maximum negative values (over sulphur oxygen) have been read and they vary as follow: -203 kJ/mol to -154 kJ/mol (for L1); -168 kJ/mol to -151 kJ/mol (for L2); -171 kJ/mol to -153 kJ/mol (for L3).The maximum positive regions (blue) are localized on the -NH groups from cycle III, assigned to atoms H4 and H5, respectively, and vary between 154 and 204 kJ/mol for all studied structures.For our further applications (complexation of heavy metals) the most interesting are orange and red regions (negative potential), which are most likely to be involved in the complexation event.
Regarding LIP maps from Fig. 3, the red regions suggest areas where electrophilic attack could happen, thus, atoms from which electron removal (ionization) is easy.By convention, in blue color are illustrated the regions where ionization is relatively difficult to occur.Orange regions (the more susceptible to ionization) present close values between 6.9 and 7.02 eV for all studied compounds, so the results of the complexation event for similarly modified electrodes should be quite similar.

Thermodynamic properties
Thermodynamic properties such as zero-point vibrational energy, enthalpy, constant volume heat capacity, entropy and Gibbs energy have been calculated with DFT method for equilibrium geometry at ground state for Z and E conformers, over a temperature range between 273.15 and 373.15 K. Except for the zero-point energy (ZPE), all depend on temperature.Their calculated values are listed in Tables 6 for 298.15K.The values of enthalpy, constant volume heat capacity (C v ), entropy (S) and Free Gibbs energy (G), increase with temperature, in agreement with the enhancement of the molecular vibrational intensities as the temperature rises [22].Obtained thermodynamic functions are given in Table 7 and they are useful to estimate directions of further chemical reactions, by computing those to other thermodynamic functions and relationships.The dependence between calculated thermodynamic properties (enthalpy, constant volume heat capacity, entropy and Gibbs energy) and temperature were fitted by quadratic equations and the correlation fitting factors (R 2 ) vary between 0.9951 and 1.The obtained values for R 2 show very good correlation of thermodynamic parameters and their fitting equations (Table 7).
The observed predicted and experimental values for chemical shifts from 13 C and 1 H NMR given in the Supporting Information, show good correlation for Z conformer of all compounds.It leads to the assumption that Z conformers are formed by synthesis, as assumed from the experimental procedure [5].

Electrochemical characterization of compounds
The electrochemical characterization of L1 -L3 compounds has been performed and the redox potentials were evaluated from the electrochemical curves of L1 -L3 azulene derivatives recorded on glassy carbon electrode by: cyclic voltammetry (CV), differential pulse voltammetry (DPV), and rotating disk electrode voltammetry (RDE) in 0.1 M TBAP, CH 3 CN.Selected anodic and cathodic CV, DPV, and RDE curves recorded in millimolar solutions (0 -2 mM) of L1 -L3 in 0.1M TBAP in CH 3 CN have been shown in Figs.4-6.When comparing the curves obtained for these compounds, the following conclusions connected to the electrochemical behavior of L1, L2 and L3 can be drawn: • Their oxidation potentials for the first oxidation peak (in V) are, respectively, 0.347, 0.297, and 0.266, being ranged in the order: L3 < L2 < L1, in correlation with HOMO energies (Table 4).This fact is in concordance with the presence of the electron releasing alkyl groups on azulene moiety in case of L3 and L2.The presence of larger and bulky substituents conducts to an easier oxidation of the ligand.• Their reduction potentials for the first reduction peak (in V) are, respectively, -1.598, -1.760, and -1.772, being ranged in the order: L3 < L2 < L1, in correlation with LUMO energies (Table 4).The small differences between the calculated LUMO values and the experimental values of the reduction potentials could be explained by a strong interaction of the ligand with the solvent, which is neglected in quantum calculus.These values also are also connected to the presence of alkyl groups on the azulene structure.• All the electrochemical processes in Figs.4-6 are irreversible.While the vinylazulene moiety is more prone to oxidation, the vinylhydantoin part is more likely to be reduced.• All compounds seem to accumulate on the electrode surface (leading to polymer films) at more positive potentials then the first anodic peak (denoted a1 on each DPV curve), in agreement with the sharp decrease of the oxidation currents in RDE and DPV curves after reaching this potential.

Electrochemical recognition using modified electrodes
The compounds L1 -L3 were used to prepare modified electrodes by electrochemical polymerization in similar conditions, as mentioned in the experimental part.
The recognition experiments performed with L1, L2, and L3 modified electrodes lead to evaluate the detection limits for Pb(II), Cd(II), Cu(II) and Hg(II) ions.In Fig. 7a are given the stripping curves obtained using poly L1 modified electrodes for the recognition of different concentrations of Pb(II), Cd(II), Cu(II) and Hg(II) ions, while in Fig. 7b are given the dependences of the stripping peak currents on concentration for the metal ions, as resulted from Fig. 7a, for the most efficient recognition events that has been found with this modified electrode for Cd(II) and Pb(II).Similar curves have been obtained for the modified electrodes prepared with L2 and L3 ligands, in similar conditions.The detection limits have been estimated for the electrodes modified with L1, L2 and L3 (Table 8).These values show that the best ligand is L1, as by using this ligand, the lowest limits for all tested heavy metals have been obtained.A correlation between ligand molecular descriptors and the obtained detection limits is evident.The best behavior of L1 can be explained by the fact that L1 presents the smallest values for HOMO and LUMO energies, and consequently, the highest ionization potential (5.41 eV) and the highest capability to accept electrons (electron affinity, A = 2.50 eV) among the three studied structures (Table 4).Also, judging by the molecular electrostatic potential, L1 structure shows the most negative potential over the oxygen with the most negative value of -203 kJ/mol, so it is the most capable to attract positive heavy metal ions charges, as previously discussed for Fig. 3.
The close values of detection limits for L2 and L3 are explained by the similar effect of substituents in L2 and L3, seen in their close values of electron affinity (Table 4).Consequently, their complexing capacity is similar, fact which is in a good agreement with the detection limits obtained for heavy metal ions using these ligands.

Conclusions
Chemical reactivity descriptors were predicted from the DFT calculations.From molecular frontier orbitals energy level diagrams and gaps, it can be concluded that among all studied structures, L1 is the most stable compound.The oxidation potentials were found to be ranged in accordance with calculated HOMO energies trend, while the variation of reduction potentials versus LUMO energies were perturbed by the occurrence of steric effects.
Thermodynamic analysis reveals that all thermodynamic calculated parameters, except the zeropoint energy (ZPE), are increasing with temperature.The chosen fitting equations for thermo- Cd Pb dynamic properties are in a good agreement with calculated data.The NMR experimental and theoretical results supporting each other, being consistently supplemented by computed molecular descriptors and thermodynamic properties.The obtained characteristics are in a good agreement with the detection limits obtained for heavy metals using modified electrodes obtained from these ligands.The corresponding detection limits evaluated for Pb, Cd, Cu and Hg are of about 10 -8 for Pb and Cd with L1 and of about 10 -5 for Cu and Hg.These limits are larger for L2 and L3, which are very close due to the similarity of their structures.They are of about 10 -6 M for Pb and Cd, and 10 -4 M for Cu and Hg.This study contributes to the prediction of specific structures able to coordinate heavy metal ions in view of recognition.

Figure 7 .
Figure 7. Stripping curves (a) using L1 modified electrodes for different concentrations of Cd(II), Pb(II), Cu(II) and Hg(II) ions and the corresponding calibration curves (b) for Cd(II) and Pb(II) ions.

Table 4 .
Calculated quantum chemical parameters of the studied compounds

Table 6 .
Predicted thermodynamic properties at 298.15 K for Z and E isomers of L1 -L3

Table 8 .
Estimated detection limits of investigated heavy metal ions for L1 -L3.