The Contribution of Density Functional Theory to the Atomistic Knowledge of Electrochromic Processes

In this review, we provide a brief overview of the contribution that computational studies can offer to the elucidation of the electronic mechanisms responsible for the electrochromism phenomenon, through the use of the density functional theory (DFT) and its time-dependent formulation (TDDFT). Although computational studies on electrochromic systems are not as numerous as those for other physico-chemical processes, we will show their reliability and ability to predict structures, excitation energies, and redox potentials. The results confirm that these methods not only help in the interpretation of experimental data but can also be used for the rational design of molecules with interesting electrochromic properties to be initiated for synthesis and experimental characterization.


Introduction
Electrochromism refers to the phenomenon in which a chemical system changes color following the application of an electric field of suitable value, returning to the original color once the applied field is eliminated. The change in the absorption band is due to the triggering of a redox reaction induced by the electric field. Compounds possessing more than two redox states are subject to different color changes as the applied electric field reaches the given redox potential. The whole process is described in Scheme 1.

Introduction
Electrochromism refers to the phenomenon in which a chemical system changes color following the application of an electric field of suitable value, returning to the original color once the applied field is eliminated. The change in the absorption band is due to the triggering of a redox reaction induced by the electric field. Compounds possessing more than two redox states are subject to different color changes as the applied electric field reaches the given redox potential. The whole process is described in Scheme 1.  Although the term electrochromic originally referred to a color change from a transparent (absorption in the UV region) to a colored species, today it is extended to all shifts of the absorption band in the entire electromagnetic spectrum caused by the application of an electric field.
Despite this amount of work, accurate information at the atomistic level of electrochromic processes is somehow lacking, and modern theoretical methods can help to shed light on these aspects. Among these, those based on the Density Functional Theory (DFT) have proved to be very useful since they are able to provide accurate results even for medium to large systems.
In this short review, we summarize the contribution of DFT methods in improving the current knowledge of the electrochromism phenomenon at the atomistic level. Particular emphasis will be devoted to those studies in which they computed and analyzed properties often difficult to obtain experimentally. Moreover, the outcomes of such investigations are not only useful to explain the material properties but also to build a rational design strategy to synthesize new molecules with given properties.
The computed properties we are interested in include structural determination, absorption behaviors and their origins, redox potentials, and electronic features.
The structures of the considered systems are reported in Scheme 2.

Structural Determination
In order to show the accuracy of DFT methods in determining geometric parameters with accuracy, we report the calculated structures for quadruply benzannulated nonplanar borepin (1) [60] and for 4-
In particular, for diphenyl- [1,2] dithiol [4,3-c] [1,2] we obtain C-C, C-S, and S-S bond length values that are very similar for the three functionals used and in the range of the experimental errors. The same agreement is found for the C-C-C, C-C-S, and C-S-S valence angles for which the error, with respect to the experimental counterpart [71], is about one degree.
Concerning the dihedral angles φ1 and φ2, a higher deviation from the experiment is noted for all the employed functionals. This can be explained by considering that, in the solid phase, the value of this parameter is influenced by the other molecules present in the The values of the main geometrical parameters reported in Table 1 clearly show very good agreement with the previously reported X-ray crystallographic data [71].
In particular, for diphenyl- [1,2] dithiol [4,3-c] [1,2] we obtain C-C, C-S, and S-S bond length values that are very similar for the three functionals used and in the range of the experimental errors. The same agreement is found for the C-C-C, C-C-S, and C-S-S valence angles for which the error, with respect to the experimental counterpart [71], is about one degree.
Concerning the dihedral angles ϕ1 and ϕ2, a higher deviation from the experiment is noted for all the employed functionals. This can be explained by considering that, in the solid phase, the value of this parameter is influenced by the other molecules present in the crystalline structure, while our result is related to a single molecule.
For borepin molecules, we note that the maximum deviation from the experimental B-C and C-C bond lengths [41] is found in the case of the M06 functional, which does, however, provide values in the range of experimental error (Table 1).
Theoretical values obtained for the C-B-C valence angles and the C-C-C-C torsional angle obtained, using the three functionals, are in excellent agreement with the experimental X-ray results [41].
As previously mentioned, the experimental structures relative to the species with different oxidation states obtained upon the application of the electric field are often missing and difficult to measure. The use of calculations based on quantum mechanics can be useful to overcome this lack of information. In Figure 1, we report the main geometrical parameters obtained for viologen 3 as a function of the oxidation state.
The values of the main geometrical parameters reported in Table 1 clearly show very good agreement with the previously reported X-ray crystallographic data [71].
In particular, for diphenyl- [1,2] dithiol [4,3-c] [1,2] we obtain C-C, C-S, and S-S bond length values that are very similar for the three functionals used and in the range of the experimental errors. The same agreement is found for the C-C-C, C-C-S, and C-S-S valence angles for which the error, with respect to the experimental counterpart [71], is about one degree.
Concerning the dihedral angles φ1 and φ2, a higher deviation from the experiment is noted for all the employed functionals. This can be explained by considering that, in the solid phase, the value of this parameter is influenced by the other molecules present in the crystalline structure, while our result is related to a single molecule.
For borepin molecules, we note that the maximum deviation from the experimental B-C and C-C bond lengths [41] is found in the case of the M06 functional, which does, however, provide values in the range of experimental error (Table 1).
Theoretical values obtained for the C-B-C valence angles and the C-C-C-C torsional angle obtained, using the three functionals, are in excellent agreement with the experimental X-ray results [41].
As previously mentioned, the experimental structures relative to the species with different oxidation states obtained upon the application of the electric field are often missing and difficult to measure. The use of calculations based on quantum mechanics can be useful to overcome this lack of information. In Figure 1, we report the main geometrical parameters obtained for viologen 3 as a function of the oxidation state. The bond between the carbon atom of the CH 3 group and the nitrogen of the pyridine ring (C1-N) undergoes a decrease, going from +2 to +1 radical to neutral charge states while the N-C2 increases following the consecutive addition of electrons. The C2-C3 and C4-C5 bond lengths involved in the conjugation of the ring decrease.
A more significant structural change is evidenced by analyzing the inter-ring dihedral angle that assumes a value of 39 • in the case of the more-stable 32+ oxidation state but becomes planar in the case of 3+1 (0.13 • ) and 3 (0.02 • ) states.
These two examples account well for the importance of determining the structural variations induced by the change in the oxidation state on electrochromic systems, achievable by using reliable theoretical tools.

Excitation Energies
To assess the performance of DFT-based methods, we report the computed excitation energies for 3 [57] and squaraine 4 [70] in their different oxidation states by using the B3LYP, M06, and wB97XD exchange correlation functionals. The range-separated hybrid wB97XD functional has been chosen since it also takes into account weak interactions. The results are reported in Table 2 together with the corresponding experimental UV-Vis values [70,72]. Table 2. Excitation energies (nm) and oscillator strength of three and four molecules in different oxidation states computed in dichlorometane solvent by using different exchange-correlation potentials and the 6-31+G* basis set. The % of the transition is that obtained at the B3LYP level. Experimental data are taken from refs [70,72]. For system 3, the maximum deviation with respect to the experiment [72] (113 nm) is found in the case of the wB97XD functional that also shows an average error of about 57 nm considering all the oxidation states. The other two functionals used give similar results with average errors of 24 and 20 nm for M06 and B3LYP, respectively. The higher average error was obtained with the functional WB97XD, although it falls within an acceptable error range (0.1-0.3 eV), which could be due to the fact that at short range, it contains 22% of Hartree-Fock exchange and in the intermediate region is smoothly described by an error function parameter of 0.2.

System
In discussing the spectral behavior following the increase of the applied electric field, we note that they are in close agreement with the experiment. In fact, the dication 3 2+ species has a measured maximum absorption band of 337 that is well reproduced by all the employed functionals.
The addition of an electron generating the 3 + cation radical slightly red-shifts this band and two new bands in the Vis/NIR region appear (447 and 1100 nm). Both M06 and B3LYP functionals reproduce these spectral changes well. A further injection of an electron causes the 3 species to change the spectral feature that now presents two bands at 558 and 712 nm. Again, the obtained excitation energies for this oxidation state give results very close to the experimental counterparts. The analysis of the excitation behavior allows one to individuate the origin of the transitions. Table 2 shows the absorption at 368 nm (B3LYP) is mainly HOMO-LUMO (H→L) in nature (70 %), while it becomes βHOMO-1→LUMO in the radical cation. The other two peaks present in the Vis/NIR spectrum are due to the βH→L and αH→L transitions, respectively. Likewise, the computed excitation energies for the 3 neutral molecule allows one to assign the two recorded absorptions to the transitions HOMO→ LUMO + 3 and H→L, respectively (see Table 2).
The squaraines containing thiophene 4 are interesting since their structures contain donor-acceptor-donor moieties, which make them useful for a variety of applications. Molecule 4 is also important because it also has, in addition to the neutral one, positive (+1, +2) and negative (−1) oxidation states that are generated in Vis/NIR spectroelectrochemistry experiments. Table 2 shows how for these systems, the excitation energies obtained are also close to the experimental counterparts especially when the M06 and B3LYP functionals are used. A major deviation is found in the neutral species, but it is due to the possible intermolecular aggregation as demonstrated recently [70]. The average error obtained with these methods (about 25 nm) confirms their reliability in the treatment of electronic properties of complex systems. The analysis of the HOMO and LUMO orbitals involved in the transitions ( Figure S1) allows one to establish that the band in the 450-500 nm range can be attributed to a charge-transfer between olate donors to the carbon atoms of the central cyclobutenediylium fragment and the thienyl moiety, for the 4 2+ species [70], while the main transition computed for molecule 4 is a π−π* band, as also demonstrated by the computed density difference plot with the first excited state considered ( Figure S2).

Redox Potentials
The possibility to obtain reliable redox potentials at the theoretical level allows one to examine the energetics of the redox processes of a chemical species. Furthermore, it is useful to set up the spectroelectrochemical and cyclovoltametric experiments as we have indications on the potentials of the electric fields to apply.
Theoretically, the redox potentials can be obtained following the Nernst equation where the Faraday constant (F) and the number of electrons transferred in the process (n) are present. The free energy change, ∆(∆Gs) redox , is computed on the basis of the Born-Haber thermodynamic cycle [73]. For a given system (A), the ∆(∆Gs) redox can be obtained by computing the solvation energies (∆ Gs) for the neutral and cationic state by using the following equation: where ∆Gg redox , ∆Gs(Red), and ∆Gs(Ox) are the change of standard Gibbs free energy for the gas phase reaction and solvation energies of reduced and oxidized species, respectively. In order to compare the results with the experimental data, the computed redox potentials are provided relative to that of ferrocene (E(Fc+/Fc).
In Table 3, we report the computed reduction potentials (V) for a series of electrochromic materials obtained by using the M06 exchange correlation potential, that in our experience is the most suitable for these determinations [57,60]. Comparison with the relative experimental data [41,43,72,74,75] is, generally, quite good. Table 3. Reduction potentials (V) for studied compounds computed by using the M06/6-31+G* level of theory. The used theoretical and experimental solvents are 5 and 6 in tetrahydrofuran; 3 and 7 in dichloromethane. a from refs [41,72,74,75]. A more detailed comparison between computed and experimental data reveals that the difference varies from 0 to 0.43 and the average error turns out to be 0.15 and 0.7 V for the first and second reduction potentials, respectively.

Electronic Properties
During the redox process, radical species can be obtained. As an example, in viologenlike systems, starting from the 3 2+ , after the first reduction step, the 3 + radical species is formed. It is well known that the stability of the radicals strongly depends on their capability to delocalize the unpaired electron: Higher delocalization means great stability.
Theoretical methods allow one to compute the spin distribution over the entire set of molecules as reported in Figure 2 for the 3 + radical. the first and second reduction potentials, respectively.

Electronic Properties
During the redox process, radical species can be obtained. As an example, in viologen-like systems, starting from the 3 2+ , after the first reduction step, the 3 + radical species is formed. It is well known that the stability of the radicals strongly depends on their capability to delocalize the unpaired electron: Higher delocalization means great stability. Theoretical methods allow one to compute the spin distribution over the entire set of molecules as reported in Figure 2 for the 3 + radical. The examination of this figure indicates that the unpaired electron is not only localized on the electron-rich nitrogen, but it is distributed on all the conjugated pyridine rings and on the C-C connecting bond with resulting system stabilization.
Furthermore, the computation of the molecular electrostatic potential allows one to follow the charge flow during the reduction steps. The visualization of what happens after the electron addition for the 3 2+ + e→3 + . + e→3 reduction behavior is given in Figure 3. The examination of this figure indicates that the unpaired electron is not only localized on the electron-rich nitrogen, but it is distributed on all the conjugated pyridine rings and on the C-C connecting bond with resulting system stabilization.
Furthermore, the computation of the molecular electrostatic potential allows one to follow the charge flow during the reduction steps. The visualization of what happens after the electron addition for the 3 2+ + e→3 +. + e→3 reduction behavior is given in Figure 3.

Electronic Properties
During the redox process, radical species can be obtained. As an example, in viologen-like systems, starting from the 3 2+ , after the first reduction step, the 3 + radical species is formed. It is well known that the stability of the radicals strongly depends on their capability to delocalize the unpaired electron: Higher delocalization means great stability. Theoretical methods allow one to compute the spin distribution over the entire set of molecules as reported in Figure 2 for the 3 + radical. The examination of this figure indicates that the unpaired electron is not only localized on the electron-rich nitrogen, but it is distributed on all the conjugated pyridine rings and on the C-C connecting bond with resulting system stabilization.
Furthermore, the computation of the molecular electrostatic potential allows one to follow the charge flow during the reduction steps. The visualization of what happens after the electron addition for the 3 2+ + e→3 + . + e→3 reduction behavior is given in Figure 3. The analysis indicates how the positive charge, although reduced, continues to be distributed over the entire molecular surface upon the first reduction step. The second reduction process is characterized by different charge distribution with negative charges concentrated in the aromatic rings and between them.
The knowledge of the frontier molecular orbitals and their composition can give further information on the electronic changes occurring during the redox processes. In Figure 4, we report the HOMO-LUMO energies and the orbital composition for the three oxidation states of 3. It is interesting to note that both HOMO and LUMO orbitals are significantly destabilized in going from the +2 to 0 oxidation states. The HOMO-LUMO energy gap changes during the reduction process, assuming the highest value in 3 2+ (5.27 eV) and the lowest in the neutral (3) oxidation state. The changes in the molecular orbital composition evidence the effect of the addition of electrons in going from 3 2+ to 3, especially for the HOMO orbital. oxidation states of 3. It is interesting to note that both HOMO and LUMO orbitals are significantly destabilized in going from the +2 to 0 oxidation states. The HOMO-LUMO energy gap changes during the reduction process, assuming the highest value in 3 2+ (5.27 eV) and the lowest in the neutral (3) oxidation state. The changes in the molecular orbital composition evidence the effect of the addition of electrons in going from 3 2+ to 3, especially for the HOMO orbital.

Materials and Methods
Full geometry optimizations in solvent media have been performed using the 6-31G* basis set in conjunction with the hybrid B3LYP [76,77], meta-hybrid M06 [78], and the range-separated hybrid wB97XD [79,80] exchange-correlation functionals. After the geometry optimization step, frequency calculations have been performed, at the same level of theory, to verify if the obtained structures are real minima. The unrestricted Kohn-Sham formalism has been used in the case of radical species. On these species, the analysis of the 〈S 2 〉 value (in all cases, reported at about 0.77) ensures the lack of spin contamination. The time-dependent density functional linear response theory (TDDFT) has been used to obtain the excitation energies starting from the previous optimization structures and employing the 6-31+G* basis set. Solvent effects have been taken into account by using the non-equilibrium polarizable continuum model [81]. All the computations have been performed with the Gaussian 09 code [82].

Conclusions
In this contribution, we have summarized the great potential of density functional theory in explaining the atomistic processes involved in the electrochromism of some

Materials and Methods
Full geometry optimizations in solvent media have been performed using the 6-31G* basis set in conjunction with the hybrid B3LYP [76,77], meta-hybrid M06 [78], and the rangeseparated hybrid wB97XD [79,80] exchange-correlation functionals. After the geometry optimization step, frequency calculations have been performed, at the same level of theory, to verify if the obtained structures are real minima. The unrestricted Kohn-Sham formalism has been used in the case of radical species. On these species, the analysis of the S 2 value (in all cases, reported at about 0.77) ensures the lack of spin contamination. The timedependent density functional linear response theory (TDDFT) has been used to obtain the excitation energies starting from the previous optimization structures and employing the 6-31+G* basis set. Solvent effects have been taken into account by using the nonequilibrium polarizable continuum model [81]. All the computations have been performed with the Gaussian 09 code [82].

Conclusions
In this contribution, we have summarized the great potential of density functional theory in explaining the atomistic processes involved in the electrochromism of some important molecular materials. The main benefits of the theoretical computational investigation can be briefly summarized: − A reliable prediction of the molecular structures with errors comparable to that coming from solid-state X-ray measurements. In addition, DFT computations allow one to determine the structures of systems with unstable oxidation states that are, consequently, difficult to measure experimentally in laboratory conditions; − The reproduction of the excitation energies in all the molecular oxidation states is close to the experimental counterparts, with an error of about 0.2-0.3 eV. It is possible to assign the transitions and analyze the orbital origin and possible charge transfer processes; − The redox potentials can be easily obtained with results very close to the experimental data; − Spin charge distribution, molecular electrostatic potentials, charge distribution, and orbital picture substantially contribute to the interpretation of experimental data and shed light on their atomistic mechanisms.
We believe that the performances highlighted herein are sufficiently positive to encourage future theoretical studies with the aim of deepening knowledge on electrochromism in organometallic systems with multiple oxidation states, ionic materials, and molecules Molecules 2021, 26, 5793 9 of 12 inserted in a polymeric or liquid crystal matrix, and at the same time, to rationally design new molecules with particular electrochromic and chemical properties to be suggested for syntheses and future applications.
Supplementary Materials: The following are available online, Figure S1: Plots of the molecular orbitals involved in the transitions for 4 and 4 2+ , Figure S2