A combined molecular dynamics and computational spectroscopy study of a dye-sensitized solar cell

An organic dye-sensitized solar cell consisting of a squaraine molecule attached to a TiO2 surface is modeled using first-principles molecular dynamics and time-dependent density functional theory. The system is surrounded by solvent molecules that are treated at the same level of theory as the dye molecule and the surface. The effect of the solvent on optical properties is investigated by computing many absorption spectra for various configurations along a molecular dynamics trajectory. It is shown that the dynamical effects induced by thermal fluctuations have a strong effect on the optical properties and that satisfactory agreement with experiments is achieved only when those thermal effects are accounted for explicitly.


Introduction
In recent years, shipments of photovoltaic panels have grown drastically and reached an estimated 17.5 GW in 2010 [1]. At the same time, the average price per watt has fallen sharply. The vast majority of photovoltaic devices use semiconductor p-n juctions, and most devices are based on either silicon or CdTe thin films. Although these technologies are well established on an industrial scale, there are also alternative types of solar cells, which are currently attracting considerable attention. Such alternative technologies cannot yet compete with the more traditional devices in terms of manufacturing cost, stability, efficiency and mass production, but they hold, at least in principle, the potential to compete in future in this rapidly expanding market.
One prominent example of such alternative systems is dye-sensitized solar cells (DSSCs) [2][3][4], where light is absorbed by dye molecules that are adsorbed on the surface of a nano-porous wide-gap semiconductor (usually TiO 2 in the anatase phase). Following light absorption by the dye molecules, the excited electron is transferred into the conduction band of the semiconductor, while the remaining hole is filled by electron transfer, typically from an electrolyte or a hole-conducting polymer. Such DSSCs have achieved efficiencies of up to 11% [5,6]. In those high-efficiency DSSCs, dyes containing transition metal complexes, often based on Ru, are employed. The choice of the dye complex is obviously fundamental for the functioning, efficiency and lifetime of the cell, and most studies of DSSCs are centered on this topic.
The relative alignment of the dye molecular orbitals and the substrate energy levels determine the cell's voltage V oc , while the geometric overlap of the orbitals determines the lifetime of the excited state and the charge injection efficiency. Those are the key properties determining the overall efficiency of the device, and a detailed understanding of those properties at the atomistic level might guide the development of better dye-substrate combinations. But despite the potential benefits of precise and systematic computational investigations on DSSC devices, such studies are very often limited to dye complexes alone or to minimal models of the substrate in terms of very small TiO 2 clusters.
More comprehensive computational studies face the complex nature of DSSC systems: a realistic model of a typical dye together with the semiconductor substrate contains at least 100-200 atoms. Moreover, the substrate is often best described using a two-dimensional (2D) periodic model for the surface. Such system sizes as well as the periodic boundary conditions are difficult to handle with sophisticated wavefunction-based quantum chemistry methods, which otherwise could provide precise information about the electronic structure and the electronic excitations in the systems. While being less precise than wavefunction-based methods, density functional theory (DFT)-based computations can routinely handle systems of that size, but they are limited to the electronic ground state. Time-dependent DFT (TDDFT), an extension of standard DFT that is capable of describing excited states, is computationally much more demanding than basic ground state DFT. Still, TDDFT studies of DSSC systems are possible, but the large system sizes and the broad spectral region of interest (from infrared to ultraviolet over the whole visible range) bring TDDFT computations quickly to its limits.
With few exceptions [7][8][9], computational models of solar cell devices do not normally include the electrolyte surrounding the sensitized surfaces or include it only in an approximate way using continuum solvation models [10]. While such effective models can capture many essential features of the solvent, they are nevertheless limited because they only model the solvent's dielectric constant, but do not take into account the precise interactions between the solvent and the solute at the atomic level.
In this paper, we present a computational study of a dye-sensitized semiconductor surface that explicitly includes water molecules [11][12][13][14][15] surrounding the dye. The water is treated at the same level of theory as the surface and the dye, and its complex interaction with the solute is unraveled using Car-Parrinello [16] (CP) molecular dynamics (MD). From the MD trajectory, snapshots of the system at different times are taken, and the optical properties of each snapshot configuration is analyzed using TDDFT. Such a systematic study including dynamical effects and solvent molecules allows us to unravel the relative influence of these factors on the experimentally observed spectra. The need to compute many optical spectra over the whole visible frequency range in a system containing more than 400 atoms has pushed this study to the limits of present-day computers and algorithms.

Model system
Our computational model consists of an organic dye molecule, squaraine, adsorbed on a two-layer periodically repeated TiO 2 slab. We have chosen squaraine as a dye in our model computations, because it is representative of all organic dyes that were studied recently as alternatives to dyes containing transition metal complexes. Another reason to choose this molecule is that the relatively simple adiabatic approximations to TDDFT are able to describe its absorption properties very well, which is not the case for other molecules, including Ru-based complexes.
The exposed surface of the slab is the (101) surface of the anatase phase, using a (1 × 4) primitive cell. We have chosen this surface because it is the most stable anatase surface [17,18], and is commonly investigated in the context of DSSC models. The squaraine molecule as described in [19] was simplified by replacing the octyl substituents with methyls. Such a simplification does no lead to significant changes in the electronic structure of the molecule [19]. The dimensions of the orthorhombic simulation cell are 19.35 × 28.61 × 54.28 a 3 0 , where a 0 is the Bohr radius. Periodic boundary conditions are applied to this unit cell. We have checked that these dimensions are sufficient to guarantee a good compromise between numerical cost and minimizing spurious interactions between the repeated images of the dye molecules in the directions parallel to the surface. The minimal distance between repeated images of the dye is about 5 Å. In the direction perpendicular to the surface, the distance between the lowermost TiO 2 layer of the repeated slab and the top of the dye molecule is approximately 6 Å. The space surrounding the dye and the surface slab has been filled with a total of 90 water molecules, representing liquid water in ambient conditions. In figure 1, this system is depicted. The surface slab plus dye system consists of 159 atoms, and together with the 90 water molecules, this computational model consists therefore of 429 atoms.
Please note that the squaraine model system described above, but without the surrounding water, is exactly equivalent to the model used previously in a study of electronic transitions using TDDFT [20]. Adding the surrounding water does not simply mean to repeat the vacuumphase computations with more atoms present. Instead, the liquid solvent introduces complex interactions with the solute, and its global effect should be accounted for by a statistical sampling over a sufficiently long time span. While it was enough in [20] to find the minimum energy configuration of the surface-dye system and to calculate one optical spectrum of this configuration, we aim here at obtaining the statistically averaged spectrum by performing many TDDFT computations along an MD trajectory.

Computational method
The computations presented here use a plane-wave basis set. The interactions between the ions and electrons are described using ultrasoft pseudopotentials [21] 5 . We use a kinetic energy cutoff of 25 Ry for the wavefunctions and 200 Ry for the charge density. In the given model system, these cutoff values lead to a basis set of 181 600 plane waves for the orbitals and 717 700 plane waves for the charge density.
The Brillouin zone is sampled using the -point only. The Perdew-Burke-Ernzerhof (PBE) functional [24] is employed for the DFT ground state computations, the CP MD and the TDDFT spectral calculations, where in the last case the adiabatic approximation [25] is used. All calculations have been performed using the Quantum ESPRESSO suite of programs [22,23].
The CP molecular dynamics was carried out at a fixed unit cell volume and an ionic temperature of 300 K. The fictitious electron mass in the CP Langrangian was 700 a.u. and a time step of 0.121 fs was used. After an initial equilibration run of 1 ps, a total of 120 000 time steps were performed, leading to an ionic trajectory of 14.5 ps. During this trajectory, 78 electronic absorption spectra were computed, one approximately every 181 fs.

5
TDDFT computation in systems comprised of 429 atoms using a plane-wave basis set containing more than 180 000 basis functions over the whole optical range is a daunting task. It has become possible thanks to a recently developed recursive scheme [23], [26][27][28][29], where instead of individual excitation energies the complete dipole response function is computed using a Lanczos algorithm. Following other applications of this method [20,30], 2500 Lanczos steps were performed for each of the three Cartesian directions of the polarization, leading to the frequency-dependent dipole polarizabilities α x x (ω), α yy (ω) and α zz (ω), which in turn determine the photoabsorption spectra. For the spectra plotted in the following, a small imaginary component of 0.027 eV is added to the frequency, leading to a broadening of the spectral peaks. The depicted absorption spectra are the average of the three directions of possible light polarization.
The employed Liouville-Lanczos algorithm, although computationally more efficient, is equivalent to other linear-response implementations of TDDFT. It has also been shown that the resulting spectra are essentially identical to absorption spectra obtained using real-time implementations of TDDFT [23,27,30].

Results and discussion
An analysis of the electronic structure projected onto atomic orbitals allows one to examine the rearrangement of electronic charges in the presence of water. In fact, the addition of the solvent molecules leads to a slight charge transfer to the TiO 2 slab: about 0.3 e are transferred to the slab when water is added. Given that the slab consists of a total of 32 TiO 2 formula units, this is a very small quantity. At the same time, about 0.15 e are transferred to the squaraine molecule. These charge transfers lead to a shift in the relative level alignment: while in the dry system the dye highest occupied molecular orbital (HOMO) lies about 0.55 eV beneath the conduction band edge of the semiconductor, this energy difference is strongly increased to 1.4 eV in the solvated case. This shift has, however, little consequence for the absorption spectrum: the strongest absorption lines determining the shape of the spectrum are intra-dye transitions. Direct transitions from the dye HOMO to the conduction band edge cannot be discerned in the spectra, as shown below.
In figure 2, we show the computed photoabsorption spectrum of the squaraine dye in the presence of the solvent. The spectrum is computed from a configuration of the system after an initial equilibration dynamics of about 1 ps. For comparison, the figure includes the experimental spectrum from [19] and the spectrum of the adsorbed dye in vacuum, in its minimum energy configuration. The latter two spectra have already been discussed in [20]. Upon contact with the solvent and subsequent equilibration, the main absorption peak of the dye is red-shifted by about 0.12 eV and now falls almost exactly on top of the experimental peak at 1.92 eV. Such precise agreement between experiment and theory is certainly fortuitous, since it is well beyond the expected precision of the employed adiabatic functionals.
While the experimental spectrum shows a pronounced shoulder at about 2.03 eV at a higher frequency than the main absorption peak, the computed spectrum exhibits a shoulder at lower frequencies, around 1.76 eV. In the spectrum computed in vacuum, no shoulder is present. One important point to note, however, is that the spectrum obtained in vacuum shows three unphysical peaks at low frequency (not shown in figure 2), at 0.72, 1.05 and 1.40 eV. Those peaks have been discussed in [20], where it was shown that they correspond to direct charge-transfer excitations from the dye to the semiconductor. It is well known that such 6  Absorption spectra for the squaraine dye adsorbed on the TiO 2 surface slab. Black curve: absorption in the dry system, where the geometry has been determined by energy minimization, as discussed in [20]. Red: absorption in the water solution after initial equilibration at 300 K. Green: the experimental spectrum from [19].
charge-transfer excitations are wrongly described in adiabatic GGA [31][32][33] which systematically underestimates the excitation energies in such cases, sometimes by as much as 1 eV or more. In our computations for the solvated complex, no such direct charge-transfer excitations are present, and the low-energy shoulder at 1.76 eV is the spectral feature with the lowest energy in that particular configuration.
At the level of a single equilibrated snapshot of the solvated system, we can thus conclude that the main effect introduced by the solvent is a red shift of the main absorption peak, a broadening of the spectrum and a suppression of direct charge-transfer excitations.
More generally, the solvent influences the solute and its optical properties in various ways. On the one hand, the presence of the solvent determines the effective dielectric constant around the solute. This dielectric constant is crucial for the electrostatics and influences both the ground state and excited electronic states. The solvent-determined dielectric constant can be easily accounted for using standard continuum solvent models [10]. A typical effect of this dielectric screening is the common 'bathochromic shift', where spectra undergo a red shift upon solvation. On the other hand, due to hybridization, the solvent's orbitals may participate actively in the photoabsorption, even at frequencies where the pure solvent is not absorbing. In such cases, the presence of the solvent may lead to enhanced absorption strength at given frequencies.
A third important influence of the solvent is through the nuclear forces exerted by the solvent molecules on the solute and affecting its structure. For this last effect, it is not sufficient to study a single configuration of the solvated system like above, but it should be observed along an MD trajectory.
During the MD trajectory for our system, the absorption spectra maintain their main features, i.e. most configurational snapshots lead to one main peak in the optical range, but the frequencies of the principal absorption peak are highly dependent on the precise molecular  configuration and vary between 1.84 and 2.07 eV. Obviously, it is tempting to look for a correlation between these frequencies and geometric properties of the corresponding snapshots. In this respect, one of the most striking features during the dynamics is that the squaraine's geometry varies considerably over time. While the molecule in its equilibrium configuration without solvent is essentially planar, it undergoes slow but strong deviations from the planar configuration during the dynamics. In the middle and right panels of figure 1, we present a side view of the dye in its planar and bended configurations. The time evolution of the angle is depicted in figure 3, where it can be seen that during the dynamics, the dye deviates by up to 47 • from a planar configuration. Examining the dependence of the main absorption peak on the bending angle, however, no correlation could be found between these features. It is therefore interesting to note that the most striking characteristic of the geometric evolution of the solute is not directly related to variations of the absorption frequencies.
In figure 3, it can also be seen that typical time scales of fluctuations of the solute are in the range between 0.5 and 5 ps. Our sampling of system geometries every 181 fs leads therefore to a faithful representation of the influence of dynamics on the optical properties.
A second possible influence on the absorption frequency is the surrounding water. Are fluctuations in the first solvation shell responsible for the observed shifts in the absorption? In order to examine this point, we have chosen a subset of 37 snapshots from the MD trajectory, spaced equally in time. For those snapshots, we have removed all the surrounding 90 water molecules and thus computed the optical spectrum in vacuum, but in the instantaneous molecular geometry without any further molecular relaxation. In this way, we could compare the frequencies in identical geometries, with and without the surrounding water. The result is depicted in figure 4, where we plotted the frequency of the main absorption peak in the dry system plotted against the frequency in the solvated system for each of the 37 snapshots. It is evident from that figure that the main effect of water is to red shift the frequencies. But apart from this global shift, the frequencies with and without water are very well correlated, showing that the variability of the absorption frequency is due to intrinsic variations in the geometry of the solute and not due to changes in the instantaneous solvation geometry.  We therefore infer that the global effect of the solvent is rather subtle: on the one hand, it is responsible for a general red shift of the absorption frequencies, a fact that can probably be accounted for by using standard continuum solvent models. On the other hand, the solvent is inducing continuous geometrical changes of the solute, which in turn influence the absorption frequency. The variations of the absorption frequency thus induced are larger than the overall red shift due to the presence of water molecules. There is no obvious dependence of the frequency on any one structural characteristic (angles, bond distances, etc). This seems to indicate that actual computation of the spectra for many configurations is the only way of systematically accounting for the observed variance in the spectra.
Finally, in figure 5, we show the resulting absorption spectrum obtained by averaging the individual spectra of the snapshots over all configurations. We also show the averaged spectra of the dry system as described above and the experimental absorption from [19]. Recall that the solute geometries in the dry system are exactly the same as those in the solvated system, being obtained from the same MD trajectory.
It is remarkable how well the averaged spectra reproduce the experimental curve: both computed spectra present one main absorption peak and one shoulder at higher energies. This is at variance with the spectra of the single snapshots shown in figure 2, where no shoulder at higher frequencies is present. The main absorption peak is found at 1.97 eV in the averaged solvated system and at 2.0 eV in the averaged dry system, in comparison with the experimental value of 1.92 eV. We ascribe this deviation of 0.05 eV at least partly to the semi-local adiabatic functional employed. The averaged spectrum also reveals that the precise coincidence of the main peaks in figure 2 is true only for the shown snapshot. The high-energy shoulder of the averaged spectrum is found at 2.06 eV in the solvated system and at 2.07 eV in the dry system. The experimental shoulder is found at 2.03 eV.
In figure 5, the effect of water molecules on the absorption spectra is also shown. Since all the snapshots that underlie these spectra have the same solute geometry, the difference between the solvated and dry spectra can be clearly attributed to the electronic effects of water molecules. On the one hand, the strong red shift upon solvation is a typical bathochromic shift, as discussed above. But the overall reduction of absorption intensity when water is removed suggests that the water molecules-interacting with the squaraine molecule-contribute to the optical response of the system also in the visible region, where pure liquid water is transparent.

Conclusions
We have performed CP MD simulations of an organic dye attached to a TiO 2 surface and surrounded by water in ambient conditions. The optical properties of this system are investigated by TDDFT computations at geometrical snapshots along the MD trajectory.
We find that the inclusion of solvent molecules, which are treated at the same level of theory as the solute and TiO 2 , considerably changes the optical absorption spectra and that a meaningful comparison with experiment is possible only in the solvated case. Moreover, the spectra of single snapshots of this system do not show the characteristic structure of the measured spectrum, namely one main absorption peak and a shoulder at higher frequencies.
These features appear only in the averaged spectra.
These results show that the optical activity of such a system may dramatically depend on the dynamical modifications of the molecular configuration induced by thermal fluctuations. No static model can properly describe these effects, which are much stronger than typical errors induced by the approximate functionals used in practical TDDFT computations.
The computation of the optical spectra by time-averaging individual spectra calculated for selected configurations is a demanding task, which can however be addressed using recent developments in the implementation of TDDFT. A major issue to be addressed and understood is the extent to which the optical properties calculated for static, minimum-energy configurations of isolated molecules are robust with respect to the combined effects of thermal fluctuations and the interaction with the solvent molecule. This will certainly be the theme for further investigations.