TD-DFT Prediction of the Intermolecular Charge-Transfer UV-Vis Spectra of Viologen Salts in Solution

The absorption spectrum of viologen salts in a medium or low polar solvent is an essential feature that influences all its “chromic” applications, whether we are considering thermochromic, electrochromic, photochromic or chemochromic devices. The prediction by quantum chemical methods of such absorption bands, typically observed in the visible range and due to charge transfer (CT) phenomena, is a very challenging problem due to strong solvent effects influencing both the geometry and the electronic transitions. Here we present a computational protocol based on DFT to predict with very high accuracy the absorption maxima of the CT bands of a series of viologen salts in solvents of low and medium polarity. The calculations also allow a clear dissection of the solvent effects, direct and indirect, and orbital contributions to the CT band.

One interesting aspect of the photochemistry of viologens is the color of the salts and corresponding solutions in low-or nonpolar solvents due to charge-transfer from the donor anion to the bipyridinium unit. As an example, N-N'-dialkyl-4,4 -bipyridinum salts, in the solid state, are white with noncoordinating and hard anions such as hexafluorophosphate and 1,1,1-trifluoro-N-((trifluoromethyl) sulfonyl) methanesulfonamide (also known as bistriflimide) [10][11][12][13]; bromide salts are yellow [11], iodide salts are red [14] while iodide salts of N-N'-dialkyl-2-2 -6-6 -tetrametyl-4,4 -bipyridinum are orange [10], to mention but a few examples which nevertheless highlight how relatively minor chemical differences can affect the UV-Vis absorption spectrum. When dissolved in low-or nonpolar solvents the corresponding solutions have a similar color to the solid salt [9].
The prediction of the absorption spectrum of organic dyes, by quantum chemical protocols, is nowadays quite advanced: covalent systems do not pose significant problems and it is possible to successfully predict the UV-Vis spectrum of organic compounds, using Time-Dependent Density Functional Theory (TD-DFT), with very high accuracy [15][16][17]. On the other hand, the application of the same protocols to noncovalent systems is much more challenging [18]: Which is the correct geometry to be considered for a noncovalent compound? Is DFT capable to correctly account for long-range effects and electron density in noncovalently bound systems? Moreover, especially regarding viologen salts in solution, can intermolecular charge transfer (CT) bands be accurately modelled by DFT methods?
In this work we will test several TD-DFT protocols for the calculation of the absorption spectra of viologen salts in medium and low polar solvents in order to provide a recipe for the prediction of the maximum of the CT absorption. To this end we select here few experimental data, mostly from the author, concerning different viologen salts in different solvents.
The experimental data of the maximum of the absorption of the lowest energy band in the UV-Vis spectrum in dichloromethane (CH 2 Cl 2 ) for four representative examples of viologen salts, have been taken from the literature and listed in Table 1. For octylviolgen iodide, we also have experimental data in acetone. As discussed in Ref. [9] and also consistent with the observations reported by Arduini et al. [19,20], there is no evidence of significant aggregation of bromides and iodides viologens salts in solvents of low and medium polarity, except for the formation of neutral clusters. This is in striking contrast with the aggregation behavior of viologens salts of "hydrophobic" anions, like bistriflimide, which show a strong dependence of the spectroscopic properties on the concentration [21,22]. However, for the systems investigated here it is reasonable to assume the formation of just neutral clusters of 1:2 stoichiometry. OV = octylviologen, 1,1 -dioctyl-4,4 -bipyridinium; DTMV = decyltetramethylviologen, 1,1 -didecyl-2,2 ,6,6tetramethyl-1,1 -bipyridinium.

Materials and Methods
The structures of the model systems were built with the Molden software package [23]. Neutral clusters made of a dimethylviologen cation and two anions have been considered since the alkyl chains, although important in order to guarantee some solubility in low-polar solvent, does not give any significant contribution to the charge-transfer band between the anions and the bipyridinium core. Energy minimization was run at the B3LYP-D3(PCM)/6-31+G(d,p), that is including the dispersion correction of Grimme [24] and the long-range solvent dielectric response [25]. n-Hexane (ε = 1.8819), dichloromethane (ε = 8.93) and acetone (ε = 20.493) were selected as solvents. All structures were checked, by frequency calculations, to be true minima with no imaginary frequencies. Xyz coordinates can be found in Supporting Information. Subsequent TD-DFT calculations were run with a variety of functionals and two basis sets of double and triple-ζ quality: 6-31+G(d,p) and aug-cc-pVTZ. Functionals used for the TD-DFT, again always including the Polarizable Continuum Model (PCM model) of solvation, were: PBE0 [26], M06-2X [27], CAM-B3LYP [28], ωB97-XD [29] and LC-ωHPBE [30,31]. The number of states used for the calculations was set to 10. All calculations have been run with the software package Gaussian16 [32]. GaussView [33] software was used to simulate the UV-Vis spectrum applying a line broadening of 0.4 eV.

Geometry
As mentioned in the Materials and Methods Section, while experimental systems consist of dialkylviologens with relatively long alkyl chains (C 8 and C 10 ), the model systems used in the calculations are the simpler methylviologens (C 1 ). This structural difference is expected not to affect the comparison between calculated and experimental results since the optical properties of the viologen salts are not influenced by the length of the alkyl chain. In fact, recent experiments concerning the electrochemistry and UV-Vis properties of several alkyl and polyfluoroalkyl viologens revealed that the CV (Cyclic Voltammetry) potentials and spectral properties are not affected by the chain length (though they are affected by the type of chain, alkyl vs. polyfluoroalkyl) [3]. Moreover, as we will see below, none of the relevant orbitals involved in the electronic transitions are localized on the methyl group of the model systems, but only on the bipyridinium core and anions, see also [34]. On the other hand, the weak electron-donor effect of the alkyl chain can be well-modelled by a simple methyl group.
We first discuss the geometry of the model systems and its dependence on the solvent reaction field. In Table 2 we report some relevant structural parameters, that is the average X-N distance (with respect to the closest nitrogen) in the neutral cluster MVX 2 , where MV represents a model viologen, and MTMV is a model 2,2 ,6,6 -tetramethyl methylviologen, and the anion X (Br − , I − or P for PF 6 − ), and the dihedral angle between the two aromatic rings of the cation. Table 2. Some relevant geometrical parameters of the systems investigated. The optimized structures in CH 2 Cl 2 for the four model systems are shown in Figure 1. Finally, in Figure 2 we graphically show the trend of the geometrical parameters of Table 2 as a function of the dielectric constant of the solvent.

R(X-N)/Å ϕ (C-C-C-C)/
First we note that, qualitatively, the structure in CH 2 Cl 2 of the neutral cluster MVBr 2 , in Figure 1b, is quite different from the other ones. For the bromide case, each anion is approximately sitting on top of the pyridinium ring and the two pyridinium rings are strongly twisted (ϕ~70 • ). The other systems, instead, feature the anion closer to the nitrogen and the twist of the two pyridinium units is smaller, see Table 2. The structure of the iodide viologen salt is qualitatively rather similar to the structure of iodide methylviologen in the solid state [36] since the anion is close to the positive nitrogen and to the ortho carbons of the pyridinium ring. However, the bipyriinium is planar in the solid state. The similarity with the solid state structure is much less for the bromide salt, both concerning the position of the anion as well as for the twist angle between the pyridinium units [36].
The geometry also appears to be influenced by the inclusion of a solvent reaction field in the optimization run. As expected for charged systems, the presence of a solvent decreases the strength of the electrostatic interaction compared to the gas phase. The effect is the largest for the bromide salt, with bromide being the hardest anion considered, and almost negligible for the PF 6 case, since the charge density is largely reduced in the hexafluorophosphate anion. Therefore the inclusion of the solvent in the optimization of the geometry is strongly recommended, and, especially for relatively low polar solvents, it makes a significant effect on the geometry, compared to vacuum. The effect, being a nonspecific reaction field only dependent on the dielectric constant (and geometric details of the cavity) levels off for dielectric constants larger than about 10 a result which is consistent with the dielectric response in the simpler Onsager model [37] of a dipole in a spherical cavity, which is 2(ε − 1)/(2ε + 1). The optimized structures in CH2Cl2 for the four model systems are shown in Figure 1. Finally, in Figure 2 we graphically show the trend of the geometrical parameters of Table 2 as a function of the dielectric constant of the solvent. First we note that, qualitatively, the structure in CH2Cl2 of the neutral cluster MVBr2, in Figure  1b, is quite different from the other ones. For the bromide case, each anion is approximately sitting on top of the pyridinium ring and the two pyridinium rings are strongly twisted (φ ~ 70°). The other systems, instead, feature the anion closer to the nitrogen and the twist of the two pyridinium units is smaller, see Table 2. The structure of the iodide viologen salt is qualitatively rather similar to the structure of iodide methylviologen in the solid state [36] since the anion is close to the positive nitrogen and to the ortho carbons of the pyridinium ring. However, the bipyriinium is planar in the solid state. The similarity with the solid state structure is much less for the bromide salt, both concerning the position of the anion as well as for the twist angle between the pyridinium units [36]. First we note that, qualitatively, the structure in CH2Cl2 of the neutral cluster MVBr2, in Figure  1b, is quite different from the other ones. For the bromide case, each anion is approximately sitting on top of the pyridinium ring and the two pyridinium rings are strongly twisted (φ ~ 70°). The other systems, instead, feature the anion closer to the nitrogen and the twist of the two pyridinium units is smaller, see Table 2. The structure of the iodide viologen salt is qualitatively rather similar to the structure of iodide methylviologen in the solid state [36] since the anion is close to the positive nitrogen and to the ortho carbons of the pyridinium ring. However, the bipyriinium is planar in the solid state. The similarity with the solid state structure is much less for the bromide salt, both concerning the position of the anion as well as for the twist angle between the pyridinium units [36].
The geometry also appears to be influenced by the inclusion of a solvent reaction field in the

TDDFT Results
In this section, we will compare the results of the TD-DFT calculations, as described in the computational section. In Figure 3, we show the correlation between the calculated wavelength of the maximum of the absorption CT band of the model cluster and the experimental values reported in Tables 1 and 3. First we note that by using a double-ζ and or a triple-ζ basis set, namely the 6-31+G(d,p) and the aug-cc-pVTZ with the CAM-B3LYP functional, the calculated absorption spectra change very little, as we can appreciate by inspecting the corresponding columns in Table 3. Therefore, all the remaining calculations were run using the smaller basis set of double-ζ quality which strongly reduces the computational burden.

TDDFT Results
In this section, we will compare the results of the TD-DFT calculations, as described in the computational section. In Figure 3, we show the correlation between the calculated wavelength of the maximum of the absorption CT band of the model cluster and the experimental values reported in Tables 1 and 3. First we note that by using a double-ζ and or a triple-ζ basis set, namely the 6-31+G(d,p) and the aug-cc-pVTZ with the CAM-B3LYP functional, the calculated absorption spectra change very little, as we can appreciate by inspecting the corresponding columns in Table 3. Therefore, all the remaining calculations were run using the smaller basis set of double-ζ quality which strongly reduces the computational burden.  Table 1 and Table 3 at three different levels of theory, TD-DFT(PCM)/6-31+G(d,p). See Table  3 for complete data. Solvent reaction field is included in the TD-DFT calculation (either CH2Cl2 or acetone). Table 3. Experimental (see Table 1) and calculated λmax data (nm) and statistical parameters of the linear fit y = a·x + b. MAE is the mean absolute error; R 2 the correlation coefficient. In the last two columns we report the data of the absorption maxima and statistical parameters in eV. Concerning the functionals used, we note that the correlation itself is always rather good, with relatively high correlation coefficients R 2 for all levels of theory used. Nevertheless, some of the functionals strongly overestimate the slope of the correlation, namely PBE0, while others  Table 3 for complete data. Solvent reaction field is included in the TD-DFT calculation (either CH 2 Cl 2 or acetone). Table 3. Experimental (see Table 1) and calculated λ max data (nm) and statistical parameters of the linear fit y = a·x + b. MAE is the mean absolute error; R 2 the correlation coefficient. In the last two columns we report the data of the absorption maxima and statistical parameters in eV.

6-31+G(d,p)
aug-cc-pVTZ 6 Concerning the functionals used, we note that the correlation itself is always rather good, with relatively high correlation coefficients R 2 for all levels of theory used. Nevertheless, some of the functionals strongly overestimate the slope of the correlation, namely PBE0, while others underestimate it, for example, LC-ωHPBE. Also CAM-B3LYP and ωB97-XD do not perform particularly well. In contrast, M06-2X gives calculated values of λ max in very good agreement with the experimental values. Not only the correlation is good, but also the slope of the fitting line is close to unity and the mean absolute error (MAE) is significantly smaller than for the other cases. The calculated results, at the M06-2X level, appear to be only systematically shifted to lower values, compared to experiments, by a roughly constant value.
The result obtained at the M06-2X(PCM)/6-31+G(d,p)//B3LYP-D3(PCM)/6-31G+(d,p), represent, therefore, the ones in better agreement with experiments among the protocols tested. At the same time, this is a relatively cheap protocol from a computational point of view. The predicting power of the protocol is extremely high if we consider that we are concerned with intermolecular CT bands between anions and cations, therefore noncovalently bound donors and acceptors. Moreover, the experimental values are obtained from systems differing because of the anion (Br − , I − and PF 6 − ), the cation (viologen and tetramethylviologen) and the solvent (dicholoromethane and acetone). Therefore the protocol appears robust. From a practical point of view, in the visible range of wavelengths, once the absorption maximum has been calculated with the above TD-DFT computational protocol, it is sufficient to rescale the calculated value by adding about 35 nm in order to get the experimental result: λ max = λ max (TD-DFT) + 35.2 nm. The additive constant proposed, 35.2 nm, is obtained directly from the MAE reported in Table 3: since all calculated data are underestimated with respect to the experimental wavelength, but nevertheless well linearly correlated, adding the MAE brings the DFT predictions in closest match with experiments.
In order to get more insights on the electronic transitions involved in the calculated absorption bands, we present an analysis of the orbitals involved.
In Figure 4, we show the simulated lowest energy absorption band of four salts in dicholomethane together with the two electronic orbitals giving the largest contribution to the electronic transition(s) which are responsible for the absorption band itself. We limit our analysis to the four cases in CH 2 Cl 2 since the results for MVI 2 in acetone are qualitatively similar. In all cases, the largest contribution comes from a HOMO → LUMO transition, the two orbitals also represented in Figure 4, these two being the Highest Occupied Molecular Orbital (HOMO) and the Lowest Unoccupied Molecular Orbital (LUMO). However, while for the iodide and bromide cases, the HOMO is represented by the p orbitals of the anion and the LUMO by a delocalized π system of the cation, for the hexafluorophosphate salt both HOMO and LUMO are localized on the cation, therefore in this case no CT band is predicted. This is, in fact, in agreement with the absence of color of the corresponding solutions and solid salts.
In Figure 5, we show the orbital energies of the HOMO orbitals of the four representative samples in CH 2 Cl 2 (between −0.40 and −0.25 a.u.), and the corresponding LUMO's (between −0.15 and −0.05 a.u.). It is clear from the figure that the tetramethyl substitution of the bipyridinium core destabilizes the LUMO level, while the HOMO level (being that of the iodide anions) is almost unaffected. This leads to a blue shift of the band. On the other hand, replacing the iodide anion with bromide, induces both a stabilization of the HOMO level (now localized on the bromide anion) as well as a destabilization of the LUMO level of the bipyridinium core leading to an even larger blue shift of the CT band. Finally, the HOMO level in the case of the hexafluorophosphate salt, is much lower in energy being also localized on the bipyridinium ring, thus producing an absorption band at much higher energy without any significant CT contribution.
Unoccupied Molecular Orbital (LUMO). However, while for the iodide and bromide cases, the HOMO is represented by the p orbitals of the anion and the LUMO by a delocalized π system of the cation, for the hexafluorophosphate salt both HOMO and LUMO are localized on the cation, therefore in this case no CT band is predicted. This is, in fact, in agreement with the absence of color of the corresponding solutions and solid salts.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 7 of 11 In Figure 5, we show the orbital energies of the HOMO orbitals of the four representative samples in CH2Cl2 (between −0.40 and −0.25 a.u.), and the corresponding LUMO's (between −0.15 and −0.05 a.u.). It is clear from the figure that the tetramethyl substitution of the bipyridinium core destabilizes the LUMO level, while the HOMO level (being that of the iodide anions) is almost unaffected. This leads to a blue shift of the band. On the other hand, replacing the iodide anion with bromide, induces both a stabilization of the HOMO level (now localized on the bromide anion) as well as a destabilization of the LUMO level of the bipyridinium core leading to an even larger blue shift of the CT band. Finally, the HOMO level in the case of the hexafluorophosphate salt, is much lower in energy being also localized on the bipyridinium ring, thus producing an absorption band at much higher energy without any significant CT contribution.

Direct and Indirect Solvent Effects
Having obtained a rather good correlation between calculated and experimental CT bands of viologens in low and medium polar solvents, we now proceed to analyze the direct and indirect solvent effects. Direct solvent effects are the ones directly affecting the TD-DFT calculation of the absorption band, for a fixed geometry, while indirect solvent effects are those due to a different

Conclusions
We have analyzed the results of TD-DFT calculations of the absorption charge transfer band of viologen salts in solution of medium to low polar solvents. Comparison of the calculated results with experimental data, using several computational protocols, allowed us to select the TD-DFT = M06-2X(PCM)/6-31+G(d,p)//B3LYP-D3(PCM)/6-31G+(d,p), as the protocol with the best accuracy to computational cost ratio among the ones investigated. A very good linear dependence of the calculated results in the UV-Vis region allowed us to propose the following empirical equation, λ max = λ max (TD-DFT) + 35.2 nm, as a tool to predict the absorption of the viologen salts in a given solvent.
In this work, we have not included highly polar and protic solvents, such as methanol or water, where we expect our computational protocol to fail because the continuum model of solvation cannot capture the additional effect due to hydrogen bonds between water and the anions. Therefore, too tightly bound species would be predicted in this case not corresponding to the real structure in solution. On the other hand, considering explicit solvent molecules for these cases, is not expected to allow quantitative prediction of the kind discussed above. The reason, as shown in Ref. [9], is also related to the fact that when the ions are well solvated and separated, the CT band is hardly detectable and, even if a weak absorption can be observed at some relatively low wavelengths (high energy due to the large separation of donor and acceptor), its position strongly depends on the solvent, as well as ion cluster, dynamics. Such dynamic behavior, however, cannot be properly described by DFT static calculations.