Vibrational Spectra of a N719-Chromophore/Titania Interface from Empirical-Potential Molecular-Dynamics Simulation, Solvated by a Room Temperature Ionic Liquid

Summary

A dye-sensitized solar cell was solvated by RTILs; using optimized empirical potentials, a molecular dynamics simulation was applied to compute vibrational properties. The obtained vibrational spectra were compared with experiment and ab initio molecular dynamics; various empirical potential spectra show how partial-charge charge parameterization of the ionic liquid affects vibrational spectra prediction.

Abstract

The accurate molecular simulation prediction of vibrational spectra, and other structural, energetic and spectral characteristics, of photo-active metal-oxide surfaces in contact with light-absorbing dyes is an ongoing thorny and elusive challenge in physical chemistry. With this in mind, a molecular-dynamics (MD) simulation were performed by using optimized empirical potentials for a well-representative and prototypical dye-sensitized solar cell (DSC) solvated by a widely-studied room temperature ionic liquid (RTIL), in the guise of a [bmim]⁺[NTf2]^- RTIL solvating an N719-sensitizing dye adsorbed onto 101 anatase-titania. In doing so, important insights were gleaned into how using a RTIL as the electrolytic hole acceptor modulates the dynamical and vibrational properties of a N719 dye, estimating the spectra for the DSC photo-active interface via Fourier transformation of mass-weighted velocity autocorrelation functions from MD. The acquired vibrational spectra were compared with the experiment spectra and those sampled from ab initio molecular dynamics (AIMD); in particular, various empirical-potential spectra generated from MD provide insight into how partial-charge charge parameterization of the ionic liquid affects vibrational spectra prediction. In any event, careful fitting of empirical force-field models has been shown to be an effective tool in handling DSC vibrational properties, when validated by AIMD and an experiment.

Introduction

In dye-sensitized solar cells (DSCs), the optical band gap of semiconductors is bridged by a light-absorbing, or sensitizing, dye. DSCs require continual recharging: therefore, a redox electrolyte is essential to foster this constant supply of charge (typically in the form of an I^-/I^3-, in an organic solvent). This facilitates the passage of holes from the sensitizing dye to the electrolyte, with injected photo-excited electrons into the metal-oxide substrate passing through to an external circuit, with eventual recombination taking place at the cathode¹. A crucial aspect underpinning DSCs' positive outlook for a wide variety of real world applications originates in their straightforward manufacture, without the need for raw materials high in purity; this is in stark contrast with the high capital cost and ultra-purity required for silicon-based photovoltaics. In any event, the prospect of significantly improving the working-life timescales of DSCs by swapping less stable electrolytes with room-temperature ionic liquids (RTILs) having low volatility shows significant promise. The solid-like physical properties of RTILs combined with their liquid-like electrical properties (such as low toxicity, flammability, and volatility)¹ lead these to constitute rather excellent candidate electrolytes for usage in DSC applications.

Given such prospects for RTILs in DSCs, it is hardly surprising that, in recent years, there has been a substantial boost of activity in studying DSC-prototype N719-chromophore/titania interfaces with RTILs. In particular, important work on such systems has been performed^2,3,4,5, which consider a broad suite of physico-chemical processes, including the charge-replenishment kinetics in dyes^2,5, the mechanistic steps of electron-hole dynamics and transfer³, and, of course, the effects of titania substrates' nanoscale nature upon these, and other, processes⁴.

Now, bearing in mind impressive advances in DFT-based molecular simulation, particularly AIMD⁶, as a highly useful prototypical design tool in materials science and particularly for DSCs^7,8,9,10,11, with critical assessment of optimal functional selection being vital^8,9, AIMD techniques have proven very useful previously in scrutinizing rather significant dispersion and explicit-RTIL solvation effects on dye structure, adsorption modes and vibrational properties at DSC-semiconductor surfaces. In particular, the adoption of AIMD had led to some success in attaining reasonable, semi-quantitative capture and prediction of important electronic properties, such as band gap, as well as structural binding¹³ and vibrational spectra¹⁴ In refs. 12-14, AIMD simulations were performed extensively on the photo-active N719-chromophore dye bound to (101) anatase-titania surface, assessing both electronic properties and structural properties in the presence of both [bmim]⁺[NTf2]^{- 12,13} and [bmim]⁺[I]^{- 14} RTILs, in addition to vibrational spectra for the case of [bmim]⁺[I]^{- 14}. In particular, the rigidity of the semiconductor's surface¹⁵, apart from its inherent comparative photo-activity, led the surface to slightly alter within the AIMD simulation, which makes (101) anatase interfaces^12,13,14 a suitable choice. As ref. 12 shows, the mean distance between the cations and the surface dropped by about 0.5 Å, the average separation between the cations and anions decreased by 0.6 Å, and the noticeable altering of the RTILs in the first layer around the dye, where the cation was on average 1.5 Å further from the center of the dye, were directly caused by explicit dispersion interactions in RTIL-solvated systems. Unphysical kinking of the adsorbed N719 dye's configuration was also a result of the introduction of explicit dispersion effects in vacuo. In ref. 13, analysis was conducted on whether these structural effects of explicit RTIL solvation and functional selection affected the behavior of the DSSCs, concluding that both explicit solvation and treatment of dispersion is very important. In ref. 14, with high-quality experimental vibrational-spectral data of other groups on hand, the particular effects were benchmarked systematically on both explicit [bmim]⁺[I]^- solvation and accurate handling of dispersion established in refs. 12 & 13 on the reproduction of salient spectral-mode features; this led to the conclusion that explicit solvation is important, alongside accurate treatment of dispersion interactions, echoing earlier findings for both structural and dynamical properties in the case of AIMD modeling of catalysts in explicit solvent¹⁶. Indeed, Mosconi et al. have also performed an impressive assessment of explicit-solvation effects on DFT treatment of DSC simulation¹⁷. Bahers et al.¹⁸ studied experimental absorption spectra for dyes along with the related spectra at the TD-DFT level; these TD-DFT spectra agreed very well in terms of their computed transitions with their experimental counterparts. In addition, absorption spectra of pyrrolidine (PYR) derivatives were studied by Preat et al. in several solvents¹⁹, providing significant insights into the dyes' geometrical and electronic structures, and evincing adequate structural modifications that serve to optimize the properties of the PYR-based DSSCs - a spirit of simulation-led/rationalized 'molecular design', indeed.

Having clearly established the important contribution of both DFT and AIMD towards accurate modelling of DSCs' properties and function, including such important technical matters like explicit solvation and appropriate treatment of dispersion interactions from structural, electronic and vibrational standpoints^{7,8,9,10,11,12,13,14}, now - in the present work - the focus turns towards the pragmatic question of how well empirical-potential approaches can be tailored to address the apposite and reasonable prediction of structural and vibrational properties of such prototypical DSC systems, taking the N719 dye adsorbed on anatase (101) in the [bmim]⁺[NTf2]^- RTIL as a case in point. This is important, not only because of the large corpus of forcefield-based molecular-simulation activities and methodological machinery available to tackle DSC simulation⁷, and metal-oxide surfaces more widely, but also because of their staggeringly reduced computational cost vis-à-vis DFT-based approaches, together with the possibility of very efficient coupling to biased-sampling approaches to capture more efficiently phase space and structural evolution in highly viscous RTIL solvents, dominated by solid-like physical properties at ambient temperatures. Therefore, motivated by this open question of gauging and optimizing forcefield approaches, informed by both DFT and AIMD as well as experimental data for vibrational spectra¹⁴, we turn to the pressing task of assessing empirical-potential performance at vibrational-spectra prediction from MD, using mass-weighted Fourier transforms of the N719 dye's atomic velocity autocorrelation function (VACF). One key concern is how different partial-charge parameterizations of the RTIL may affect vibrational-spectra prediction, and particular attention was given to this point, as well as the wider task of tailoring forcefields for optimal spectral-mode prediction relative to experiment and AIMD²⁰.

Protocol

1. Performing MD Simulation using DL_POLY

1. Construct the DSC-systems initial structure of the N719-dye adsorbed to an anatase-titania (101) surface solvated by [bmim]+[NTf2]^- taken from previous work^12,13. Draw the required structure using VESTA software.
2. Choose the N719 cis-di(thiocyanato)-bis (2,2'-bipyridl-4-carboxylate-4'-carboxylic acid)-ruthenium (II)-sensitizing dye with no counterions and ensure the presence of two surface-bound protons to provide for overall-system charge neutrality.
   NOTE: Indeed, in depth studies by De Angelis et al. have established this to constitute a realistic representation of N719 adsorbed to anatase-titania²¹. This is because in ref. 21, there was the most convincing level of agreement with experimental results for a number of properties; in experimental systems, the de-facto surface protonation is thought to emerge from ILs' cations and anions leading to some degree of charge transfer with the surface²¹.
3. Ensure that the dye is adsorbed chemically to the TiO₂ surface through two carboxylate groups (i.e., bidentate). This initial adsorbed dye configuration is similar to that denoted as I1 and found by Schiffman et al.²², which was determined to be the most stable with surface protonation taken into account. Consult refs. 12 & 13 for a detailed account of how this is done, including chemical-adsorption coordinates.
4. Ensure that the sensitizing dye N719 (cis-di(thiocyanato)-bis(2,20-bipyridl-4-carboxylate-40-carboxylic acid)-ruthenium(II)) has no counterions. Add two surface-bound protons for charge neutrality, as in ref. 12 & 13.
5. Select 12 cation-anion pairs of 1-butyl-3 methylimidazolium bis(trifluoromethyl sulfonyl)imide, comprising 480 atoms^12,13. These were taken from refs. 12 & 13.
6. Relax the RTIL configuration via empirical-potentials, using the well-validated forcefield of Lopes et al.²³. Model anatase using the Matsui-Akaogi (MA) force-field and include the mobility of the titania in the relaxation process. Using the DL-POLY details in step 2.1 below, perform a geometry-optimization in DL-POLY, rather than MD, with a conjugate-gradient-minimization relative termination gradient of 0.001. Here, specify optimization in the FIELD file, rather than dynamics.
7. For the anatase surface, (TiO₂)₉₆, consisting of 288 atoms, ensure it is periodic along x- and y- laboratory axes, projecting to the RTIL a pair of parallel (101) surfaces; the dimensions in x-axis 23 Å and y-axis to 21 Å. This was taken from refs. 12 & 13.
8. Ensure that the entire DSC system with an explicit solvent is composed of 827 atoms^12,13; For the 'in-vacuo' case, bereft of RTIL solvation, there should be 347 atoms in the system.

2. Performing forcefield-based MD Simulation using DL_POLY

1. Perform MD using DL-POLY with various different partial-charge sets (vide infra) for 15 ps with a 1 fs time-step and at 300 K in an NVT ensemble^24,25, using the Lopes et al.²³ forcefield parameters for the RTIL and general-purpose OPLS model for the dye²⁶, with the well-studied and reliable Matsui-Akaogi potential acting for titania²⁷ force-field. To run DL-POLY on the terminal, type DLPOLY.X & where the input files are located.
2. Perform classical MD via these above-specified empirical force-fields, as implemented in DL_POLY²⁸. Here, there is no need to use a graphical user interface (GUI) in the software, so it is recommended to input the details using the comprehensive, and easy-to-follow software manual²⁹. Here, in the CONTROL file (check Supplemental Information for the input files), specify 'Nose-Hoover' for NVT, and opt for position-velocity trajectory printing every 1 fs.
3. In the FIELD file, for Lennard-Jones parameters, apply Lorentz-Berthelot combining rules²⁵. Taking the arithmetic mean of the Lennard-Jones (LJ) radii and the geometric mean of the LJ well-depths, for the empirical force-fields, as detailed in ref. 25, and enter this in the bottom section of the FIELD file under the non-bonded-interactions tab.
4. To handle long-range electrostatics, apply the Ewald method²⁵; use a non-bonded cut-off length of r_cut = 10 Å. Consult refs. 25 & 30 for precise details as to how to optimize electrostatic parameters. Set the real-space decay parameter for the Ewald method in the CONTROL file to be ~3.14/r_cut, and choose the number of Ewald wave-vectors to ensure a relative tolerance in the Ewald evaluation of 1E-5; specify that in the CONTROL file.
5. Ensure, in the CONTROL file, state that r_cut = 10 Å; carry out a series of potential-energy evaluations with a REVCON file (renamed as CONFIG) until the system pressure in OUTPUT converges to within a few percent to choose r_cut, but avoid any r_cut below ~2.5 times the largest LJ distance^25,30.
   NOTE: This short 15 ps MD-propagation timescale is chosen to be similar to that of ~8.5 ps Born-Oppenheimer-MD (BOMD) simulations with an identical starting configuration of refs. 9, 10 & 17, so as to allow for direct comparison vibrational-spectra prediction afforded by both AI-²⁰ and forcefield-based MD (with comparison and validation against experiment also).
6. From the HISTORY file (into which both velocities and positions have been printed at each time step, as directed from the CONTROL file), extract the x-, y-, z-velocities by using python dye_atom_velocity_seperate.py (see Supplemental Information) in the terminal. It will separate the velocities at each step.
7. Compute the VACF by using vacf151005.py (see Supplemental Information). In the terminal, type ./classical_dye_autocorr.sh; it will compute the VACF of all dye atoms. Compute spectra from MD (whether AIMD^14,20 or forcefield-based) using mass-weighted Fourier transforms of the dye's atomic velocity autocorrelation function (VACF)^31,32,33 by using python MWPS.py (see Supplemental Information). In the terminal, type ./run_all_4.sh; it will compute the mass weighted power spectra.
8. Perform a Fourier transform on these VACFs using commonly available software.
9. Bear in mind that the best-quality DFT treatment in AIMD (e.g., use of explicit solvation and accurate treatment of dispersion, alongside a canny choice of functional) is important to use for benchmarking against experimental data^12,13,20 and gauging/tailoring the comparative performance of empirical-potential MD, in particular, the large influence of electrostatics and choices for partial charges¹⁶. See refs. 12 and 13 and study these to gain an in-depth appreciation, and, if intending to do AIMD (which is not the case in the present study), act accordingly in such a situation in the future, should the need to conduct AIMD arise.

3. Comparing results of each of the force-fields

NOTE: It is important to assess partial-charge sets for the RTIL for empirical-potential-based MD simulation in step 2, for ready comparison against each other, experiment and ab initio-MD results in explicit RTIL solvent (using the PBE functional with Grimme-D3 dispersion, given its superior performance for prediction of vibrational spectra)²⁰; these were as follows:

1. Note that in the case of literature-derived RTIL charges, the anion charges are to be found from Extended Hückel Theory^34,35, based on AIMD trajectories²³, owing to the absence of anion-charge parameterisation in ref. 20, with cation charges to be taken from Lopes et al.²³ Prepare a table of the literature charges, and put into FIELD-file format for DL-POLY.
2. Note that Mulliken RTIL charges are to be calculated via Mulliken population analysis. Perform the Mulliken analysis by averaging over four points of the ab initio MD trajectory²⁰, renormalize and prepare a table of the literature charges, and put into FIELD-file format for DL-POLY.
3. Note that Extended Hückel Theory (EHT) charges are to be fitted from the final configuration of AIMD trajectories²⁰, using EHT, applied to both RTIL anions and cations. Perform the EHT analysis by averaging over four points of the ab initio MD trajectory²⁰, as implemented in the MOE software package (by selecting 'Charge Analysis' menu after reading in the configuration file)³⁵, renormalize and prepare a table of the literature charges, and put into FIELD-file format for DL-POLY.
4. Note that Hirshfeld RTIL charges are to be calculated from Hirshfeld-charge analysis by averaging over four points of the ab initio MD trajectory²⁰, for anions and cations²⁰, as implemented in the MOE software package (by selecting 'Charge Analysis' menu after reading in the configuration file)³⁵. From renormalization of these so-obtained charges, tabulate these in the appropriate format in the DL-POLY FIELD file.
5. Bear in mind that the different charge sets for the [bmim]⁺[NTf2]^- atoms are presented in Table 1 & Table 2, which also show the mean amount by which some of the atomic charges needed to be modified, so as to take symmetry and overall charge conservation into account.
6. Note that the final charge sets are to have the sum of charges adding up to +1 on the cation and -1 on the anion. The cation and anion are shown in Figure 1a and Figure 1b, respectively. The underlying DFT-sampled spectra, from which these charge sets were inspired, by and large, have implicit charge transfer and tend to lead to charges closer to ±1.

Results

Structural Properties of Binding Motifs
Representative binding motifs for the four different partial-charge sets are depicted in Figure 2, after 15 ps of MD. In Figure 2a, for the (above-described) literature-derived charges, it can be seen that there is a prominent hydrogen-bonding interaction with a surface proton. From careful analyses of the trajectory, the hydrogen bonds are mostly sur...

Discussion

Ab initio simulation techniques are expensive to perform and hence to perform simulation on much longer timescales would require the use of empirical forcefields for at least some of the DSC system. Towards this end, an equivalent atomistic model was created of the [bmim]⁺[NTf2]^- solvated interface, using an empirical, classical-simulation forcefield for MD. The anatase was modeled using the Matsui-Akaogi (MA) forcefield, whilst the dye structure was handled using OPLS parameters. For the RTIL, four...

Materials

Name	Company	Catalog Number	Comments
This was a molecular simulation, so no experimental equipment was used.
The name of the software was DL-POLY (the 'Classic' version of which is available under GnuPublic Licence, via sourceforge)