RECENT ADVANCES IN THEORETICAL INVESTIGATION OF TITANIUM DIOXIDE NANOMATERIALS. A REVIEW

© Blazhynska M.M., Kyrychenko A.V., Stepaniuk D.S., Korsun O.M., Kovalenko S.M., Ivanov V.V., Miannay F.-A., Idrissi A., Kalugin O.N., 2020 6 УДК 541.2+544.7 RECENT ADVANCES IN THEORETICAL INVESTIGATION OF TITANIUM DIOXIDE NANOMATERIALS. A REVIEW M.M. Blazhynska†,a, A.V. Kyrychenko†,b, D.S. Stepaniuk†,c, O.M. Korsun†,d, S.M. Kovalenko†,e, V.V. Ivanov†,f, F.-A. Miannay‡,g, A. Idrissi‡,h, O.N. Kalugin†,i


Introduction
Nanocrystalline titanium dioxide (TiO 2 ) is one of the most widely used nanomaterials because of its chemical stability and environmental compatibility. In addition to many traditional applications in pigment chemistry, metallurgy, and catalysis, nano-sized TiO 2 has become a very appealing material for novel emerging areas of material science, including solar energy harvesting and hydrogen generation from water through photochemical reactions [1][2]. The performance of these applications depends critically upon size, shape, crystallinity, and morphology of TiO 2 nanomaterials. TiO 2 exists in various polymorphs, such as rutile, anatase, brookite, and other rare phases ( Fig. 1) [3]. Rutile and anatase have been widely used in the industry over the last decades [2,4]. Due to its large dielectric constant and its high reflectivity across the visible spectrum, rutile is commonly used in electronic devices, such as thin-film capacitors, interference filters, and optical waveguides. Anatase https://doi.org/10.26565/2220-637X-2020-34-01 is frequently used as a photocatalyst [1] because of its high photoreactivity in homogeneous and heterogeneous catalysis. Chemical reactions on TiO 2 surfaces play decisive roles in photocatalysis and photoconduction reactions. Figure 1. Crystallographic structure of anatase, rutile, and brookite TiO 2 . Visualization and graphics were drawn with the VESTA program package [5].
One of the advantages of TiO 2 nanomaterials is their biocompatibility. Engineered TiO 2 nanostructures can exist in contact with human body tissue without causing an unacceptable degree of harm to the body [6]. Due to this unique capability, TiO 2 has emerged in many biomedical applications, such as artificial bones, joints and dental implants [7]. Therefore, significant experimental and computational efforts have been invested in studying phenomena occurring at the liquid-solid interface because their mechanisms and driving forces can point to new directions for biomaterial design and evaluation [8]. It has been shown that TiO 2 contributed to the development of a novel class of composite materials that physically integrates inorganic nanoparticles and biological molecules [9][10]. These hybrid conjugates can offer new opportunities for diverse applications in biosensing, site-selective recognition, and may potentially be incorporated into living organisms and biological cells [11]. Therefore, understanding protein adsorption on biomaterial surfaces is of great importance since, alongside water-surface interactions, it can significantly affect the performance of a biomaterial [10].
In this contribution, we discuss recent progress and achievements in the theory and computational modeling of TiO 2 nanomaterials and their applications across a wide range of research fields [12][13]. The effects of crystal morphology, crystallographic planes, surface defects, on the adsorption process are discussed. The role of the environmental effects, such as solvent, pH, and temperature, in the adsorption behavior of small organic ligands and biopolymers, is analyzed in detail. In addition, we overview the recent progress in classical MD simulations and some advanced MD sampling techniques, such as replica exchange and umbrella sampling, examining the energetics and dynamics of the interactions of single amino acids, small peptides, and proteins with TiO 2 surfaces [6]. Some quantitative insights into molecular recognition mechanisms between biopolymers and inorganic interfaces, including structure-function relationships are also considered. Particular focus is given to how MD simulations have contributed to our current understanding of peptide-surface adsorption at atomic detail. Finally, the perspectives, challenges, and promising directions in computational modeling interactions between inorganic materials and natural biomacromolecules are outlined.

Modeling of TiO 2 Surfaces and Aqueous Interfaces
The first empirical force-field (FF) for bulk TiO 2 was developed by Matsui and Akaogi (MA) [14] in 1991, in which the non-bonded interaction parameters between Ti and O were determined by reproducing the crystal structure of various polymorphic forms of TiO 2 (Fig. 1). In the MA force-field, the energetic interactions in TiO 2 were assumed to be completely non-covalent. The dispersion and repulsion interactions were represented by the Buckingham potential (1) ( Table 1) [14]. The electrostatic contributions were treated with the pairwise additive Coulomb potential.

ion-ion
A ij (kcal/mol) B ij (Å) C ij (kcalꞏÅ 6  In succeeding years, the MA force-field has been used as an initial guess for further extension and the development of more complex computational models, suitable for simulations of TiO 2 surfaces, interfaces, nanoparticles and small clusters [26][27]. The first attempt to follow this approach has been performed by Bandura et al. [28][29] by using periodic density functional theory (DFT) calculations to derive new empirical FF parameters for modeling the TiO 2 (110) surface. A primary interest in modeling and simulations of TiO 2 interfaces was focused on studying of adsorption of biological molecules, especially peptides and proteins [30][31].
Because the most experimental studies were conducted in aqueous solution, it was vital that the adsorption interactions between TiO 2 and water molecules were adequately described. Since then, numerous MD simulations have been focused on modeling TiO 2 -water interfaces [32][33][34][35].
The adsorption behavior and the interfacial structure of water molecules on the uncharged and negatively charged rutile (110) surface were studied by classical MD simulations at 298 K and 1 atm. [32][33]. The original interaction parameters of the MA force-field for TiO 2 were further adjusted by ab initio calculations of relaxed structures, and charges of the ideal rutile (110) surface. As a result, new modified interaction parameters were suggested for interfacial and deeper layers of TiO 2 atoms. Water molecules were treated with the SPC/E water model [36]. Therefore, the interactions between Ti atoms and the oxygen atoms of water molecules (Ti-O H2O ) were optimized by ab initio structure optimization and are described by a Buckingham potential (1). It was found that the positions of water molecules at the interface were found to be above each terminal titanium atom. A single additional layer of adsorbed water molecules occupies distinct sites related to the underlying crystal surface structure [32]. The water structure and mobility quickly decay to the bulk liquid properties beyond the second water layer. It was concluded that the influence of the interface for all simulated physicochemical properties, including spatial orientation and diffusivity of water molecules, is negligible beyond distances of about 15 Å from the surface [32].
Starting from the interatomic potentials proposed by Bandura and Kubicki [28], Alimohammadi and Fichthorn [34] have developed a new set of FF parameters for modeling the interaction of water molecules adsorbed either molecularly or dissociatively onto the TiO 2 surface. The proposed FF was validated by computing binding energies and conformations for molecular adsorption of water molecules on the (101) and (112) surfaces of anatase, and the (110) surface of rutile. Moreover, the dissociative adsorption on water molecules onto the bulk-terminated (001) surface of anatase and the (110) surface of rutile was considered for various water surface coverages. The MD simulation results were found to be in good agreement with those of first-principles DFT calculations and experimental observations [34].  (110) surface: DFT (black diamonds, ◊) and classical MD calculations (red circles, ○). The corresponding structure is displayed in the inset. Adapted with permission from [37]. Copyright © 2011, American Chemical Society. Schneider and Ciacchi have performed a robust computational procedure of the development and validation of empirical potentials for classical MD simulations of the Ti/TiO x /water interface employing ab initio MD simulations [37]. The procedure was based on the mapping of the adsorption energy profile of a single water molecule on the rutile TiO 2 (110) surface with Coulomb and Lennard-Jones potentials (Fig. 2). The following LJ parameters were suggested: σ(Ti)=0.07827 nm, σ(O)=0.1615 nm and ε(Ti)=1.402557 kJ/mol and ε(O)=1.91151 kJ/mol, respectively. The proposed FF was then validated over a representative series of small organic molecules, such as CH 3 OH, CH 3 NH 2 , HCOOH, as well as the tripeptide Arg-Gly-Asp (RGD) adsorbed onto the oxidized Ti surface [37].  (110) surface: a) non-hydroxylated surface, b) fully hydroxylated surface showing the bridging and terminal hydroxyl groups, c) the surface partially covered by terminal hydroxyl groups, d) the surface covered completely by terminal hydroxyl groups and partially covered of bridging hydroxyl groups. Adapted from with permission from [32]. Copyright © 2004, American Chemical Society.
Skelton and Walsh have studied the interactions of liquid water and dissolved Na + and Clions with the non-hydroxylated and hydroxylated rutile (110) surface (Fig. 3) by using a modified TIP3P water model, which is compatible with the CHARMM FF [38]. Moreover, further extension of the original MA FF has been released by MD simulations of the adsorption of water molecules onto different surfaces of TiO 2 , such rutile (110,100,101,001) and anatase (101,001), in which water molecules were represented with a flexible SPC model. A primary goal of this study was focused on studying the structure, ordering, and orientation of water molecules located in the first and second interfacial layers [39][40][41]. In 2006, Köppen and Langel developed a new set of empirical FF parameters, suitable for modeling of the interfacial phenomena occurring at the rutile (100) surface in contact with NaCl aqueous solutions by using a two-particle interaction model with the Morse potential and Coulombic interactions [42]. Their FF model was able to mimic acid-base properties, such as the protonation and hydroxylation equilibriums of the TiO 2 interface, which could be adjusted to specific pH values (4, 7.4, and 9). Interfacial pH changes were mimicked by varying a ratio of protonated and hydroxylated residues at the interface, according to Fig. 4. The primary focus of their study was on validation of the Stern model by confirming of a dense layer of counterions on the charged TiO 2 surface and a diffuse layer propagated into the bulk water [42].
Conventional MD simulations of the TiO 2 interface are able to sample statistical mechanical ensembles with a fixed surface composition, which are well described by classical non-bonded interactions or fixed valence FF models. In classical MD simulations, an initially assigned protonation state of the TiO 2 interfacial atoms can not be altered during the simulation, which imposes some limitations for modeling the acid-base interfacial phenomena. One of the promising approaches, overcoming these restrictions, has recently been proposed by the development of a new ReaxFF reactive FF [43], capable of modeling chemosorption reactions in the Ti−O−H system [44][45][46]. The ReaxFF force field parameters have been fitted to a quantum mechanical (QM) structures and energies related to bond dissociation energies, angle, and dihedral distortions. Fig. 5 shows snapshots of the rutile (110) surface with layers of water molecules [46]. Moreover, reactions between water and TiO 2 , as well as experimental crystal structures, heats of formation, and bulk modulus data were also compared. The authors demonstrated that the results of the ReaxFF calculations match reasonably well with those of DFT results for water binding energies, surface energies, and H 2 O dissociation energy barriers [46]. The agreement between the DFT/MD simulation and ReaxFF simulations of water dissociation levels was reported to be within a 10% error [44,46].  (110) surface taken at the beginning (a) and the end (b) of the sampling. Adapted with permission from [46]. Copyright © 2013, American Chemical Society.
One of fundamental simplification of classical MD simulations of the water-TiO 2 interfaces is that a site charge located at each Ti and O atom kept fixed during modeling so that the variations of the surface dipoles in TiO 2 induced by water and the dipole moments of water molecules near the TiO 2 surfaces are not taken into consideration. Recently, a novel polarizable FF for water/rutile interfaces was developed [47]. This polarizable model takes into account the effect of the surface polarity induced by liquid−solid interactions explicitly. The model represents some extension of the charge response kernel (CRK) method [48][49][50] for molecules to solid surfaces by introducing the surface CRK (SCRK) (Fig. 6a-b), in which the CRK parameters were systematically derived by the first-principles calculations in the TiO 2 slab with the dipole-correction method.
The proposed polarizable FF was validated for the water/clean rutile TiO 2 (110) interface (Fig. 6c). Good agreement was observed between structures and induced charges of a single water molecule, adsorbed onto the TiO 2 surface, obtained by the polarizable FF and those predicted by the firstprinciples calculations [50][51]. In addition, MD simulations of the liquid water/TiO 2 interface were benchmarked over stable structures of water adsorbed onto the TiO 2 surface, indicating that the dipole moments of water and TiO 2 induced by the water−TiO 2 interactions have a significant impact on the structure and dynamics of the water/TiO 2 interface [47,52]. It should be noted that the numerous studies of the absorption behavior and the structure of water layers at water/TiO 2 interfaces have many similar features with computational studies of other liquid/solid interfaces [53]. To identify interface criteria, it would be instructive to apply a concept of the "10-90 width", which is defined as the distance from the solid surface over which a specific interfacial order parameter changes from 10% to 90% of the bulk value [53]. This concept was successfully utilized in MD simulations of the basal and prism interfaces formed by ice 1h and water molecules, approximated with the rigid SPC/E model [54]. The use of this concept by analyzing of the translational and window average-density order parameters allowed the authors to show that the interfaces consist of a series of relaxed double layers, which differ in diffusion coefficients and bond arrangement. The thickness of the interfaces was found to be 0.9-1.1 nm, which was consistent with ellipsometry measurements [54]. Importance of the "10-90 width" concept was recently highlighted for the proper setup of ab initio/DFT MD simulations of water/ice interfaces, for which a relatively large size of ~1.1 nm of water interface should be utilized [55].
It has long been recognized that the structure of water molecules confined into nano-sized volumes differs significantly from that of bulk water. In ref [56], MD simulations of the SPC/E water molecules, confined in TiO 2 -rutile pores of diameters 1.3, 2.8, and 5.1 nm, were carried out at various water contents. It has been found that the proximity to the interface affects density and diffusivity of water molecules within a distance of around 10 Å from the walls, beyond which all structural and dynamics properties tend to converge [56].
To understand the difference between polar TiO 2 and hydrophobic graphite surfaces, the microstructure of water molecules absorbed at the nano-slits of the rutile (110) surfaces with the separation distance ranging from 0.8 to 2.0 nm were analyzed by using MD simulations [57]. MD simulation results show that the residence time of water at a TiO 2 surface is considerably longer than that at the graphite surface. Besides, the bound water molecules at the TiO 2 surface and the hydrogen bond network reduce the diffusivity of water through the TiO 2 slits. These MD simulations show that the surface chemistry is crucial to the diffusion of water inside nanoscale pores as compared to geometric size effects [57]. To investigate the microstructure and lubrication of water molecules confined in TiO 2 nano-slits under shearing, a series of non-equilibrium MD simulations have also been carried out [58]. To shed light on the role of the confined water molecules on lubrication, the effects of varying slit gap widths (0.8, 1.2, 1.6, and 2.0 nm) and shear velocities (200,100,50, and 10 m/s) on the friction coefficients between TiO 2 and water molecules (SPC/E) were evaluated. It has been shown that the friction coefficient decreased as the slit width increased. The detailed analysis of the microstructure of water layers revealed that water molecules confined in the slits were ordered up to two interfacial layers [58].
Using MD simulations with the LAMMPS package, the static and dynamic properties of water on TiO 2 nanotubes (TiO 2 -NTs) with a diameter of ~1.0 nm ( Fig. 7a-b) were studied [59]. The MA forcefield was used for TiO 2 , whereas the SPC/E model was utilized for water. The MD simulation results showed that the water molecules outside TiO 2 -NT conform to the two-layer model for water, which is typical for water adsorbed onto a planar surface. The difference was that the first water layer was further affected by the surface curvature, indicating the more facile water desorption from the surface. This layer disappears for water inside TiO 2 -NT, leaving a water layer with the hydrogen atoms pointing to the oxygen atom O 2c to form hydrogen bonds ( Fig. 7c-d) [59].
The mechanical behavior of anatase TiO 2 nanotubes in axial loading was examined by using MD simulations with the MA force-field [60]. The nanotubes of different radii and lengths were modeled at 300 K and subjected to both stretching and compression. The MD simulations demonstrated that, as the radius of the nanotubes increases, their surface was wrinkled axially, which in turn caused their axial behavior to change dramatically. The threshold radius was observed to be 50 Å for a length of 300 Å [60]. The phase behavior of the n-octacosane (n-C 28 )−water mixture inside TiO 2 nanopores was considered with classical MD simulations [61]. Three model systems having different TiO 2 (101) anatase pore sizes of 34, 40, and 55 Å were examined ( Fig. 8a and 8c). TiO 2 atoms were modeled as Lennard-Jones spheres with σ(Ti)=0.092 nm, σ(O)=0.303 nm and ε(Ti)=0.17138 kJ/mol and ε(O)=0.5016 kJ/mol as taken from [62]. Water molecules were treated with the SPC/E model, whereas the united-atom TraPPE (transferable potential for phase equilibria) FF was used for n-C 28 . The standard Lorentz−Berthelot combining rules were used for calculations of nonbonded Lennard-Jones interactions between sites of a different type. It was found that even though phase segregated mixtures retained their structural properties, as compared to their bulk counterparts, some significant deviations were observed in the density profiles inside the nanopore ( Fig. 8b and 8d). Water molecules were organized into two discrete layers on the TiO 2 surface so that n-C 28 was shielded from the nanopore walls (Fig. 8d). It was suggested that octacosane's self-diffusion was influenced only slightly by confinement. However, water molecules were severely hindered by the TiO 2 nanopore surface. Similarly to other studies, it was found that water diffusivity revealed a strong dependence on the distance from the nanopore center [61].
Recently, the MA force-field has been utilized as the initial guess for the development of interatomic potentials for modeling the deposition of calcium and phosphate species (calcium phosphate, Ca 3 (PO 4 ) 2 ) on the anatase TiO 2 (101) and (100) surfaces, respectively [63]. The authors have parameterized the missing two-body interactions between Ca 2+ , PO 4 3and TiO 2 , and performed classical MD simulations of the deposition of calcium phosphate on the anatase surface in aqueous solution, which 13 were further benchmarked over DFT results. The comparison between the DFT and MD simulations results have shown the outstanding performance of the new set of FF parameters in describing the adsorption of Ca 3 (PO 4 ) 2 on the most stable (101) surface of anatase [63]. New insights on the molecular processes in the initial nucleation of Ca 3 (PO 4 ) 2 on the rutile TiO 2 surfaces were obtained by using classical (MD) simulations [64]. Five kinds of TiO 2 substrates were considered: nonhydroxylated and fully hydroxylated perfect TiO 2 surfaces, and the stepped-, grooved-, and ridged-structured TiO 2 surfaces with full surface hydroxylation. The FF parameters for the surface hydroxyls were obtained from Předota's work [32]. The FF interaction parameters between Ca 3 (PO 4 ) 2 and TiO 2 were previously derived from the ab initio calculations [65]. The MD simulation results suggested that the surface hydroxylation and the shape and the structure of the TiO 2 substrate contributed significantly to the initial nucleation of calcium phosphate. Surface hydroxyls on TiO 2 provided active sites for the aggregation of calcium phosphate. Both Ca 2+ and PO 4 3ions were able to bind to the hydroxylated TiO 2 surface directly or indirectly via the first water layer ( Fig. 9a-b). Surface nanotopographies (e.g., grooves or ridges) were able to restrict the diffusion of Ca 2+ and PO 4 3ions ( Fig. 9c) [64].  (110) surface shown in the direct adsorption conformations (a) and the hydrated adsorption conformations (b). Adsorption conformation of calcium phosphate ions on the ridge-structured TiO 2 surface (c). Adapted with permission from [64]. Copyright © 2018, American Chemical Society.

Adsorption of small additives onto TiO 2
Crystal growth and morphology of crystalline TiO 2 are governed by adsorption of small organic additives. These parameters are also essential for the development of novel functional materials. Therefore, MD simulations of the binding modes and thermodynamically stable conformations of small ligands at the surfaces of TiO 2 crystals have become attractive tools for understanding the mechanisms underlying additive-surface interactions. As results, numerous computational studied have been devoted to examining the adsorption behavior of small organic compounds, such as ammonium [66], methane [67], methanol [66], formic acid [67], benzene [66][67], and dimethyl methylphosphonate [68], on TiO 2 surfaces. Adsorption of a rhodamine dye [69] and some biologically-relevant molecules, such as phospholipids (DOPC, DOPS, and DMTAP) [70], onto TiO 2 was also examined.
Recently, the use of the electronic continuum correction (ECC) theory, which suggests some adjustment of point charges at Ti and O atoms to reproduce the adsorption phenomena occurring at the TiO 2 −liquid interface, has received significant attention [32][33][71][72][73][74][75]. MD simulations, augmented by free energy calculations and supported by ab initio calculations, were used to study the adsorption of oxalic acid ions (oxalate and hydrogen oxalate) on the rutile (110) surface [76][77]. The predicted adsorption on perfect non-hydroxylated and hydroxylated surfaces (Fig. 3) was found in agreement with experimental adsorption data and predictions of the charge-distribution multisite ion complexation model [76]. The most favorable surface species were identified to be outer-sphere complexes because of the strong hydrogen bonding of oxalic acid ions with surface hydroxyls and adsorbed water molecules [76]. Moreover, MD simulations were performed to identify the thermodynamically stable conformations of glycolate (CH 2 (OH)COO − ), lactate, and 2-hydroxybutyrate ions ( Fig. 10) at the (001) and (110) planes of rutile crystals by using the MA force field [78][79]. The MD results suggested that, for both planes, the conformation and dynamics of the glycolate ion were strongly dominated by a layered structure of water on the surface. It was found that the glycolate ion was bound to the surface more stably on the (110) plane than on the (001) plane [78]. Besides, the bound conformation in which the COO − group was oriented toward the rutile surface was found to be the most stable for the lactate and 2-hydroxybutyrate ions. In contrast, a conformation, in which the COO − group was oriented toward bulk water, was found to be the most stable for the glycolate ion [79]. The adsorption process of a nanodroplet of cysteine molecules onto a perfect and the defective rutile (110) surface was simulated in the gas phase using classical MD approaches, based on reactive and nonreactive FF parameterizations [44]. Three molecular forms of cysteine species, namely, one neutral and two zwitterionic ones, were considered. It has been shown that the adsorption of cysteine molecules on perfect and defective rutile surfaces indicates that both of the computational approaches are able to recognize the importance of the sulfur atom in the binding process. However, it was also pointed out that only through a reactive FF description, the protonation state of cysteine at the interfaces produced by environmental effects could be correctly identified [44].

Interaction of biomacromolecules with TiO 2 surfaces
Understanding the interaction mechanism between metal and metal oxide surfaces and biomolecules, such as proteins/peptides/aminoacids, is crucial to the successful usage of hybrid biomaterials in complex physiological environments [13,15,[80][81]. Therefore, the development of theory and computational approaches for modeling organic/inorganic interfaces may accelerate the development of novel biomaterials [10,[81][82][83].
The MA force-field has attracted essential interest in MD simulations TiO 2 nanomaterials: however, its applicability for nano-bio MD simulations remains limited because of its discrepancy with the standard force fields for biomolecular simulations already available in the literature (e.g., AMBER, CHARMM, OPLS). So, some long-standing problems of accurate sampling of biomolecules near inorganic surfaces and calculation of binding free energies of amino acid side chain analogs and a peptide to the TiO 2 surface require the development of the compatible FF parameters. To accomplish this goal, Brandt and Lyubartsev [84] developed bonded and non-bonded interaction parameters for classical MD simulations of TiO 2 surfaces in contact with aqueous media (including biomatter) (Fig. 11a).
The proposed FF was designed in agreement with the following criteria: (1) compatibility with standard FFs for biomolecular simulations; (2) transferability of TiO 2 interactions to a wide variety of biomolecules: proteins, lipids, and nucleic acids so that Lennard-Jones nonbonded interactions with new molecules could be estimated by standard combination rules; (3) the force field should reproduce interactions at biologically relevant conditions, such as room temperature and atmospheric pressure; (4) the FF parameters were derived in a thermodynamically consistent manner from known experimental data. The modeling was performed by quantum mechanical geometry optimizations of water on TiO 2 surfaces using density functional theory using the Perdew−Wang exchange-correlation functional and ultrasoft pseudopotentials (Fig. 11b). A bonded TiO 2 −surface model was developed based on these optimized geometries [84]. The performance of the proposed bonded TiO 2 FF was validated by calculations of adsorption profiles and adsorption free energies for the side chain analogs of the 20 naturally occurring amino acids and a titanium binding peptide on the TiO 2 (100) surface [85]. Finally, adsorption simulations of a titanium binding RKLPDA-peptide were benchmarked to identify peptide binding modes on a TiO 2 (100) surface ( Fig. 11c-d) [85].
The adsorption behavior of the four DNA bases (adenine, guanine, thymine, and cytosine) on a charged partially hydroxylated rutile (110) surface was studied through classical MD simulations [86]. The FF parameters reported by Predota et al. was used for the hydroxylated rutile (110) slab [32]. The AMBER FF was used for the DNA bases and the modified TIP3P FF for water [87].
To identify the best binding arrangement, and quantify the strength of their interactions with the inorganic surface, the potential of mean constrained force (PMF) was calculated. The general tendency was observed so that each DNA base adsorbed onto structured solvent layers, which were in direct contact with the TiO 2 substrate. Upon the surface binding, the DNA bases were located to be perpendicular to the slab without adopting a preferential orientation. The DNA base binding to the TiO 2 surface was rather weak, so that desorption back into the bulk solvent was observed [86]. Single amino acids or short peptides have been the first molecules to be computationally studied [13,[80][81][88][89]. Ab initio periodic calculations and classical MD simulations were performed to investigate the adsorption mode, conformational characteristics, and inter-adsorbate interactions of alanine and several dipeptides, such as alanine-glutamic acid and alanine-lysine onto the TiO 2 (110) rutile surface [30,90].
Among the many possible peptides, the RGD (Arg-Gly-Asp) and RKLPDA sequences of amino acids are of great interest [6]. After the placement of a Ti implant inside the body, integrin receptors at the cell membrane will search for specific ligands on the surface for binding. If the ligand is present and its conformation on the surface is suitable, further interaction between the cell and the implant can occur. Protein ligands, such as fibronectin, vitronectin, and collagen, are present in the extracellular matrix (ECM). The cellular response induced by these extracellular matrix proteins, however, is mainly through the Arg-Gly-Asp (RGD) sequence (Fig. 12a) [91].
The adsorption of RGD-peptide onto the rutile TiO 2 (110) surface was modeled by using the MA force-field for TiO 2 and TIP3P for waters [92][93]. It was found that RGD-peptide was able to adsorb rapidly on the surfaces by amino groups NH 2 and NH 3 + , and the carboxyl group COOcould edge out the adsorbed water molecules and bond to the surface titanium atom (Fig. 12b). Since earlier computational studies [93], RGD-peptide has become a benchmark system for the understanding of the origin and molecular mechanism of the molecular recognition between proteins and peptides, and the TiO 2 -aqueous interface [37,65,[94][95][96][97][98][99]. While the physicochemical aspects of the adsorption of RGD-peptide onto TiO 2 surfaces were studied in detail, the mechanism by which other peptides adsorb to titanium surfaces is not yet fully understood [65,[100][101][102][103].
Two experimentally identified 12-mer peptide sequences referred to as the titania-binding peptides Ti-1 (QPYLFATDSLIK) and Ti-2 (GHTHYHAVRTQT) [88,100] revealed the unusual adsorption affinity to for titania that efficiently promote titania mineralization from an aqueous titanium bisammonium lactatodihydroxide (TiBALDH) solution [104]. Influences of Ca 2+ ions at the aqueous TiO 2 interface on the binding modes of Ti-1 and Ti-2 peptides was investigated using replica exchange MD simulations combined with solute tempering MD simulations [105]. It was observed that for Ti-1, Ca 2+ ions enhanced the adsorption of the negatively charged Asp8 residue in this sequence to the negatively charged surface, via Asp-Ca 2+ -TiO 2 bridging. However, for Ti-2, fewer residues were predicted to adsorb directly to the surface in the presence of Ca 2+ , possibly due to competition between the other peptide residues and Ca 2+ ions to adsorb to the surface [105]. Recently, a new Ti-binding peptide of 13 amino acids (AMRKLPDAPGMHC) was identified by using a series of multiple 300 ns MD samplings of adsorption dynamics of the peptide at the aqueous interface of a TiO 2 anatase surface, started from different peptide orientations with respect to the surface [106]. It was suggested that the peptide anchoring was mediated by the interfacial water layers by means of the charged groups of the side chains of the peptide [106].
Bone morphogenetic protein-2 (BMP-2), a dimeric cysteine knot protein, has therapeutic effects for ectopic bone [107]. Moreover, it plays a crucial role in osteoblast differentiation and proliferation [107]. Classical MD simulations and a hybrid approach of steered molecular dynamics (SMD) combined with MD simulations were employed to investigate the initial stages of the adsorption of BMP-2 upon approaching hydrophilic surfaces of rutile TiO 2 [62]. Surface adsorption was evaluated for six different orientations of the protein, two end-on, and four side-on, in an explicit water environment. It was observed that BMP-2 adsorbed only loosely to the hydrophilic TiO 2 surface. Despite favorable interaction energy between BMP-2 and the TiO 2 surface, the rapid formation of a two-layer water structure prevented the direct interaction between protein and TiO 2 . It was suggested that the first water adlayer had a strong repulsive effect on BMP-2, while the second attracted the protein toward the surface [62].
Adsorption orientation and conformation of myoglobin on perfect rutile (110) and (001) surfaces were investigated by combining Monte Carlo and MD simulations [108]. Myoglobin was modeled with the CHARMM27 FF, whereas the modified MA force-field [103] re-fitted to the Lennard-Jones potential, was utilized for rigid TiO 2 surfaces. It was found that the adsorbed myoglobin conformations were not influenced by the physicochemical properties of the rutile (110) and (001) surfaces (Fig.  13). These findings were understood in terms of the strong hydrophilicity of both surfaces. The adsorption behavior of myoglobin was governed by two layers of adsorbed water molecules formed on both rutile surfaces. After 80 ns MD simulation, the heme of myoglobin was found to be close to the rutile (001) surface and far away from the rutile (110) surface [108]. RNA aptamers with unique tertiary structure and high affinity for binding to a specific target with high affinity and specificity targets are used as therapeutic agents [109]. The adsorption behavior of an RNA aptamer with (101), (100) and (110) surfaces of TiO 2 anatase was investigated by using MD simulations [110]. Analysis of the structural parameters and interaction energies demonstrated that the (110) surface was energetically more favorable for the adsorption of the RNA aptamer than the (100) and (101) surfaces, respectively [110].
Collagen is the main structural protein in the extracellular matrix, which plays an important and irreplaceable role in biomineralization [111]. Collagen is used to modify biomaterial surfaces so that the pre-adsorbed collagen on TiO 2 could improve adhesion. It stimulates cell behaviors mainly through the amino acid sequence Arg-Gly-Asp interacting with the integrin at the cell membrane in cell adhesion. It is believed that the carboxyl groups in collagen play a key role in hydroxyapatite nucleation [112]. The adsorption of collagen on the perfect and grooved rutile (110) surface in aqueous solution was modeled by classical MD simulation [113][114]. Preferable binding modes were examined, including direct and indirect binding modes, in which the peptide interacted with the substrate surface directly or via the first layer water molecules (Fig. 14). The MD simulations suggested that the initial poses of collagen on the rutile surface could influence the final adsorption conformation of collagen. The reduced rutile surface provided active sites for collagen adsorption. It was concluded that the direct binding mode was responsible for the stable adsorption of collagen. The indirect binding mode might also play an important part at the initial adsorption stage, but itself alone could not 'trap' the collagen on the surface stably unless the direct binding mode had already been formed [113]. The nonspecific adsorption of noncollagenous bone matrix proteins (NCPs) onto a titania interface was examined for the case of the adsorption of the α-1 domain of osteocalcin [115].
Fibronectin (FN) is a major component of the extracellular matrix in tissue, which influences cell adhesion and migration through interactions with other extracellular components [116]. FN contains the important cell-binding sequence -arginine-glycine-aspartic acid (RGD), which can serve as a primary attachment site for TiO 2 . To investigate the topographical dependency of FN adsorption, MD simulations were employed to describe its adsorption behavior on nanostructured rutile (110) surfaces ( Fig. 15) [117]. It was found that the residence time of adsorbed FN was relied on its binding mode (direct or indirect) with the substrate and the surface interior. In the direct binding mode, FN molecules were found to be 'trapped' at the anchoring sites of the rutile surface, or even penetrate deep into the interior of nanostructures, regardless of the initial configuration on the surface [117].

Interaction of polymers with TiO 2 surfaces
Passivation of metal and metal oxide nanoparticles by polymer coating is the conventional approach to improve their colloidal stability and biocompatibility [118][119][120][121][122].
MD simulations of the interface between poly(ethylene oxide) (PEO 54 ) and TiO 2 were performed at 423 K using a quantum chemistry-based FF [123]. A representative model system, used for the FF parameterization, consisted of a PEO monomer (dimethyl ether) complexed with the TiO 5 H 9 cluster. It was found that the PEO density was significantly perturbed by TiO 2 surfaces so that layers of highly dense polymer were formed that persisted up to 15 Å from the surface. Conformational and structural relaxations of the interfacial PEO were found to be dramatically slower than those of bulk PEO. It was suggested that the surface structure and electrostatic interactions between PEO and TiO 2 governed the nature of PEO relaxation at the TiO 2 interface [123]. Adsorption of poly(3-hexylthiophene) on the nanostructured titania (110) rutile surface in vacuo was studied using MD simulations as a function of local surface curvature and roughness [124].
Coating of TiO 2 with polyethylene glycol (PEG) is particularly convenient because it is inexpensive, FDA-approved prevents opsonization from the immune system and therefore increases the in vivo circulation time [125][126][127]. In addition, it has been reported that poly(vinylpyrrolidone) (PVP) and PEG oligomers and may be used as co-adsorbents, which reduce electron recombination by the effective coverage of vacant sites on TiO 2 surfaces that improve the performance of dye-sensitized solar cells [128][129]. 13 C NMR measurements combined with dispersion-corrected density functional theory (DFT-D2) calculations [130] and classical MD simulations were performed to determine the operating condition that guarantee very high grafting densities of PEG on amorphous TiO 2 -NPs, which are desirable in any biomedical application [131]. For TiO 2 , a simplified version of the MA FF model [132] was used, while as for water, the TIP3P model was utilized. Polyethylene PEG was described with Generalized Amber Force Field (GAFF). It was shown that classical and popular models of polymer conformation on surfaces failed in determining the mushroom-to-brush transition point and proved that it takes place only at a rather high grafting density value (Fig. 16) [131].

Modeling of TiO 2 in Non-Aqueous Environment
Chemisorption of solvent molecules and ligands onto the TiO 2 surface may strongly affect its electronic structure [133]. For these reasons, growing interest has been focused on modeling and simulation of TiO 2 in the non-aqueous environment [134][135][136][137][138][139].
Acetonitrile is commonly used as a polar organic solvent for Grätzel DSSCs. Schiffmann and coworkers have studied the structure and dynamics of the anatase(101)-acetonitrile interface by combined ab initio and classical MD simulations [134]. The CP2K simulation package and the Gaussian plane wave (GPW) approach were employed to derive the empirical FF for the TiO 2 /acetonitrile interactions. A rigid slab of TiO 2 was modeled in contact with neat acetonitrile. It has been found that the strong dipolar interactions between acetonitrile molecules result in pronounced layering of adsorbed molecules propagated further away from the interface to bulk solvent up to seven distinct solvent layers (~12 Å). Each successive layer was characterized by an alternate orientation of the acetonitrile dipoles. Such capability of the solvent to make a strong interaction with the TiO 2 wall rather than with bulk solvent might have several significant practical effects, which, for instance, could lead to passivation of the solar cell surface [134].
Singh et al. examined the dynamics of ionic liquids (ILs) inside the nanopores of TiO 2 and suggested that the pore size plays a crucial role in the determination of the macroscopic electrical resistance of IL-bases DSSCs [135]. The object of their study was the interplay of the rutile wall, containing nanopores of various shapes, which were confined with 1-ethyl-3-methylimidazolium bis(trifluoromethylsulfonyl)imide [EMIM + ] [TFMSI -]. The previously developed all-atom force field parameters [140] were used for modeling the cations [EMIM + ] and anions [TFMSI -], respectively. The authors were able to identify major factors, which influence on the interaction TiO 2 -IL: the IL density inside the nanopores, distances of the ions to the pore walls, as well as the effective pore size, pore shape, surface roughness, and the surface density of electrical charges within the pore walls [135].
Recently, particular interest has been paid to studying the interaction mode between TiO 2 and various ILs at the molecular level [136]. Dai et al. studied the role of interfacial hydrogen bonding between the TiO 2 slits of different widths and ethylammoinium nitrate (EAN), determining the orientation of nanoconfined IL ions on the (100) rutile surface by combining of MD simulations and vibrational spectroscopy [137]. MD simulations were carried out with the NVT ensemble at T=353K by using the all-atom FF for the EAN IL, developed by Acevedo and Tirado-Rives [141]. It was found out that NH 3 + is more likely to adsorb to the TiO 2 wall than NO 3 -. However, both NH 3 + and NO 3 tend to be closer to the surface than CH 3 + ions, which could be explained by the hydrophilicity of the TiO 2 slit. Moreover, the orientation of IL ions was depended on the width of TiO 2 , which significantly affects the microstructure of the EAN: the smaller is the TiO 2 wall, the stronger is hydrogen bonding of NH 3 + with the TiO 2 surface that results in lying down of cations along the TiO 2 wall. The increase of the width of the slits leaded to decreasing of the hydrogen bond strength and proportion between IL and TiO 2 and, consequently, to the interfacial reorientation of IL ions [137].
Malali and Foroutan have studied the ionic liquid droplets containing 1-butyl-3-methylimidazolium and fluorophosphate on the rutile (110) surface by MD simulations [142]. Using the OPLS-AA FF customized for ionic liquids [141], it was found that the closest surrounding of IL to the nanosurface could be represented by two differently oriented layers of cations. The imidazolium rings are parallel to the surface in the first layer, and they tend to orient perpendicular to the surface in the second region. Such orientation creates some vacancies and empty spaces between these two layers, making the spherical shape gaps, which favors the IL becomes non-wetting at the upper layers of the droplets [142]. Consequently, such non-wetting behavior depends on the covering with a full layer of cations, but not on the density on the TiO 2 surface and this effect of the interaction energy of the substrate on ionic liquids leads to a circular configuration of ions in the center of the ionic liquid droplet (Fig. 17).  (Fig. 18a), in the aqueous environment confined between of amorphous slit nanopores [143]. The MD simulations have been performed in an NVT ensemble at 300 K by using the GROMACS package.  6 ] ILs were taken from [141].
The authors have found that the behavior of IL ions inside the nanopores differs significantly from the adsorption behavior at the TiO 2 anatase surface so that the extent of an ion-pairing lifetime under such nanoscale confinement is substantially lower than its value found in the bulk solution. The asymmetric distribution of confined IL ions inside amorphous walls was attributed to the formation of ion-pairs ( Fig. 18b-f), whereas there is a more homogenous-like distribution of cations and anions inside anatase pores. It has been suggested that the ion-pairing decreases in dilute aqueous mixtures of these ILs [143]. Considering the fact that a key element of any DSSC devices is the TiO 2 /electrolyte interface, Kislenko et al. have studied the influence of different electrolytes (lithium ion (Li + ), 1,2dimethyl-3propylimidazolium (DMPIM + ), tetrabutylammonium (TBA + ) as a cation and iodide ion (I -) as an anion) (Fig. 19) on the formation of the anatase (101) double-layer structure in acetonitrile media [144]. Acetonitrile was described with the model proposed by Nikitin et al. [145]. The Amber94 FF was used for the TBA + cations. FF parameters for DMPIM + were taken from the work of Liu et al. [146], whereas some missing FF parameters were derived by ab initio calculations. Finally, the interaction potentials between TiO 2 , acetonitrile, and IL ions were adapted from the earlier works of Schiffmann et al. [134].
It has been found that the physicochemical properties of such systems differ significantly, because of different size, symmetry, rigidity, and adsorptive properties of each cation. In the case of Li + , which is prone to adsorption and intercalation, ions in bulk tend to aggregate in pairs, triplets, and larger clusters, representing chains of alternating Li + and I − ions, which leads further to decrease of their diffusion. In the case of DMPIM + , the formation of ion pairs and ion clusters of irregular shapes were observed, while the aggregation of TBA + ions in bulk was not frequently occurred, because of its low absorptiveness. The acetonitrile molecules form a self-assembled monolayer on the TiO 2 surface, preventing adsorption/desorption phenomena for electrolytes. However, it has been pointed out that the adsorbed I − anions produce local "defects" in the solvent monolayer, which creates adsorption "channels" for Li + and DMPIM + cations and results in the formation of the adsorbed ion pairs [144]. Figure 19. Snapshots of the MD simulation boxes with the electrolytes containing different cations: Li + (a), DMPIM + (b), and TBA + (c). Adapted with permission from [144]. Copyright © 2013 American Chemical Society.

Amorphous TiO 2 Nanoparticles
Despite many advantages and popularity of the MA force-field of TiO 2 in material science modeling, its functional form has rather limited transferability to conventional biomolecular MD simulations. The majority of popular MD packages, such as GROMACS, ABMER, NAMD, etc., have welldocumented FFs and numerous libraries of interaction parameters customized for describing the van der Waals interaction between atoms by the modified Lennard-Jones (LJ) potential (2). To overcome this difficulty, Zhou and co-workers [132] have re-parameterized the MA FF parameters for TiO 2 by fitting of the Buckingham potentials for TiO 2 with the Lennard-Jones one ( Fig. 20 and Table 2).  The customized FF parameters, adopted to the LJ form, were tested in MD simulating of the bulk structure of TiO 2 , water interaction, and the adsorption of the C-terminal (from Met 548 to Glu 556 ) fragment of the human serum albumin (HSA) at a spherical TiO 2 -NP [132].

Covalently Bound Models of Ti-O Interaction
A bonded FF, based on ab initio potentials, has been developed for interactions of water molecules with the (110) surface of rutile TiO 2 surface slabs [32]. This FF was suitable for atomistic MD simulations of the interfacial structure of the uncharged surface in contact with liquid SPC/E water and negatively charged surfaces in contact with dissolved electrolyte ions (Rb + , Sr 2+ , Zn 2+ , Na + , Ca 2+ , Cl -) [32][33]147].

TiO 2 Nanoparticles
Recent progress and achievement in MD simulations of physicochemical properties of TiO 2 nanomaterials have been reviewed elsewhere [6]. Here, we provide the most significant advances and critical examples in this field.
MD simulations of isolated TiO 2 -NPs in the three characteristic morphology (anatase, brookite, and rutile) have been reported, in which FF parameters were developed based on the MA force-field and fine-tuned with simulated X-ray diffraction data [19]. It has been suggested that the Ti-O bond length was found to be dependent on the coordination environment of Ti and independent of morphology and the size of NPs. The equilibrium Ti-O bond length for a titanium ion was assigned to be 0.186 nm for a four-coordinated ion, 0.192 nm for a five-coordinated ion, and 0.194 nm for an octahedral Ti ion, respectively [19]. Fig. 21 shows MD snapshots of 3 nm size NPs of quasi-spherical anatase and rutile NPs, respectively. To get a better understanding of the size TiO 2 nanoparticles below which anatase NPs become more stable than rutile NPs (crossover diameter), MD simulations of TiO 2 nanoparticles were performed based on the MA potential [22]. Sintering-induced phase transformations between various polymorphs, such as anatase and rutile, were considered. It was found that phase transformation between anatase and rutile phases was facilitated by enhanced ionic mobility at temperatures near 1473 K, corresponding to the melting point of the NPs. The final sintering agglomerate transformed to the rutile phase provided one of the sintering NPs was rutile, while sintering of anatase and amorphous nanoparticles resulted in a brookite agglomerate. No such phase transformations were observed at temperatures away from NPs melting temperatures [22].
Alimohammadi and Fichthorn used MD simulations to study the aggregation of various TiO 2 (anatase) nanocrystals in vacuum [24]. They found that the aggregation occurred during the picosecond time-scale (Fig. 22). It was not driven by dipole-dipole interactions so that higher-order multipole moments were suggested as the driving force for the preferential alignment of TiO 2 -NPs. These highorder multipole moments originate from under-coordinated O and Ti surface atoms on the edges between nanocrystal facets, which create localized regions of positive and negative charge [24].
MD simulations of the aggregation behavior of amorphous spherical TiO 2 of the size 6 nm were examined at 350 K and 1400 K by using the MA force-field [148]. Some localized atomic rearrangements were found to induce preferential crystallization near the surfaces of the NPs. The deposition of amorphous TiO 2 -NP of the size 3 nm onto a flat rigid TiO 2 (B) nanocrystal surface was studied by using spatially resolved pulsed laser deposition combined with MD simulations to gain insight into the binding dynamic and the structural behavior of TiO 2 -NP [149]. The simulations were based on the MA potential and carried out at discrete temperatures of 673 K, 873 K, and 1073 K, respectively. MD simulations revealed that an amorphous TiO 2 -NP was attached to the surface within a picosecond time scale, during which some interdiffusion of its atoms into the topmost layer of the TiO 2 (B) crystal occurred. Starting from a metastable amorphous nanoparticle with 8% TiO 6 octahedra, the TiO 2 -NP tended to crystallize within the next 2 ns so that the Ti atoms have formed octahedral oxygen bonds [149].

Quantum-Chemical Calculations of Dye-Sensitized TiO 2
Dye-sensitized solar cells (DSSC), which are based on dyes adsorbed (or covalently connected) on the surface of TiO 2 , have been intensively studied in recent years [150]. High efficiency of solar energy conversion, low cost of structural elements, and production of solar cells made them attractive for researchers. Many organic dyes have been studied as constituent elements of DSSC and the relationship between chemical structure and performance of corresponding solar cells has been intensively investigated [150]. In the current section, we present a brief description of the operational principle of DSSC. A more detailed overview can be found elsewhere [151][152][153][154][155].
The main constructive elements of DSSC are an organic dye and a TiO 2 semiconductor (Fig. 23). An organic dye (usually π-conjugated system) contains two structural parts connected by a bridge: a donor and an acceptor. The acceptor moiety is linked to the surface of the semiconductor.
e TiO e (TiO ) The oxidized dye ( S  ) regenerates the ground state by I  ion (mediator).

S e S
    .
where LHE (light-harvesting efficiency) can be obtained using oscillator strength (f) directly from quantum chemistry calculations of electronic spectra. The  It is clear that the light conversion efficiency of a DSSC element is heavily influenced by the chemical structure of a dye. It is why quantum chemical calculations play an important role in the investigation and development of high-performance DSSC devices.
During the last decade, the complexes "organic dye-TiO 2 " were intensively studied by using different level quantum chemistry (ab initio) methods. The main problems of such studies are connected, firstly, with the necessity to consider large and "heavy" systems, and secondly, the study of spectral characteristics usually requires rather sophisticated methods. In recent years, many works have been devoted to the study of the electronic structure, geometry, and effects of energy transfer in "organic system-TiO 2 " complexes. In the present review, we consider the essential articles only, in which, besides the actual computational data concerned with specific molecular systems, the general theoretic methodological aspects were discussed as well.
TiO 2 is one of the most investigated single-crystalline systems. The geometrical structure and properties of rutile and anatase surfaces are described in detail in many reviews (see, for instance, [156]). The results of quantum chemical calculations for TiO 2 are presented in a series of articles. Among the most interesting works in this field is the article [157], where the quantum chemical calculations for cluster models (TiO 2 ) n (rutile n=20, and anatase n=92) were performed. Different structural and electronic properties were obtained by using several quantum chemistry approaches. Among them, the semiempirical PM6 method, Hartree-Fock (HF), DFT (BLYP, B3LYP, PBE, PBE0, CAM-B3LYP, and ωB97 functionals) approaches. The LANL2DZ basis set was used. The presented data demonstrated significant sensitivity of some critical parameters of the electronic structure to the chosen functional. However, the parameters of the electronic structure demonstrated small sensitivity to the geometrical structure optimization method. For instance, the values of the HOMO-LUMO gap (HL) ob-tained using the B3LYP functional vary by only 0.08 eV for differently optimized structures (different DFT functionals were used for optimization). Surprisingly, the HL gap, obtained, for instance, with B3LYP and CAM-B3LYP, differs by as much as 2.91 eV! Also, the performance of different exchange-correlation functionals for calculations of rutile and anatase was compared in ref [158].
Another structural aspect of anatase TiO 2 NPs was discussed [159]. In this study, the nanospheres and nanocrystals of various sizes were compared (Fig. 24). The HSE06 functional, which is considered as the best suited to describe semiconductors, was used. The comparison of HSE06 with the standard B3LYP functional was also presented. Faceted and spherical nanoparticles have been cut from the bulk anatase crystal (Fig. 24a-b). The initial geometry of anatase was derived from DFT calculation. It was shown that the calculated surface energies for NPs higher than for a regular surface. In general, it was demonstrated that these NPs have essentially different structural and electronic properties. The electron and geometrical structure of TiO 2 nanocrystals with an excess of electrons introduced by defects were investigated in [160]. The spin-polarized variant of the DFT theory and the exchangecorrelation functional HSE06 were used to characterize the electronic properties of the nanocrystals. The various oxygen and titanium atom (interstitial) defects were created by removing an O atom (for oxygen vacancy), and Ti was inserted into interstitial location between oxygen atoms (Ti int defect). According to the presented data, the oxygen vacancy demonstrated similar stability at the surface and sub-surface sites, whereas Ti int was more stable at sub-surface sites. In general, the result of the study reveals the possibility of tuning the band structure through defects.
A water-TiO 2 (anatase) interface was investigated in [161][162]. The calculations were performed by using a semiempirical variant of the DFT theory. Namely, the tide binding model, as an approximate DFT method (DFTB), was used for large cluster calculations. The DFTB method is more than a thousand times faster than the conventional B3LYP approach. Thus, the DFTB method can be used to calculate large clusters composed of thousands of atoms. The test calculations that were performed in [161][162] showed a rather good correspondence between the geometric parameters (water molecules and anatase slab cluster) obtained by the DFTB method and those obtained by the DFT/PBE method. In general, it proved that the DFTB set of parameters was able to reproduce the main features of the Ti/(water multilayers) interface with the DFT accuracy. However, it was pointed out that the DFT level of theory is still required when the electron excitations and electron transfer were under investigation. The benchmark study of DFTB against DFT is given in [163].
The accuracy of quantum chemical calculations of excited states was recently discussed in [164]. The results of the time-dependent DFT theory (TD-DFT) were compared with the equation of motion coupled-cluster singles and doubles (EOM-CCSD) data for rather small (TiO 2 ) n with n varying from 1 up to 13 clusters, respectively. High-quality basis sets were used, such as def2-SV(P), def2-TZVP, and def2-TZVPP, respectively. According to the computational data of the TD-DFT (with different functionals used) and EOM-CCSD methods demonstrated qualitatively similar results for excitation energies of TiO 2 NPs. However, for systems with a significant charge-transfer character, the standard functionals (e. g. B3LYP) demonstrated a rather poor agreement with the exact EOM-CCSD approach [164]. Such a situation can be illustrated by the data presented in Table 3. Table 3. Excitation energies (in eV) for the four lowest states ((TiO 2 ) 3 cluster, C 1 symmetry) [ One can see from Table 3 that the best quantitative agreement between TD-DFT and EOM-CCSD (with the best basis set def2-TZVPP) was achieved for the CAM-B3LYP functional, while B3LYP demonstrated the rather poor accuracy [164].
Some essential structural and electronic criteria, which required (but not sufficient) for a dye to generate a photocurrent upon adsorption on a semiconductor surface, have been suggested as following [165]:  The anchoring group should be positioned on the electron-accepting group.  The HOMO and LUMO of the dye should be energetically above the valence band edge (VBE) and the conduction band edge (CBE) of the semiconductor, respectively.  The absorption spectrum should significantly overlap with the solar spectrum.  The transition involved in the photoconversion phenomenon should have a high intensity.  The fluorescence lifetime should be sufficiently high (above the nanosecond time scale) to allow an electronic injection to the semiconductor from the excited state of the dye before quenching of the excited state by radiative decay.
One of the most interesting works concerned dye-TiO 2 complex is the article [183]. Herein indicates that even for rather simple π-systems, the interaction with TiO 2 particles may lead to substantial changes in their optical properties. The catechol and dopamine molecules (Fig. 26a) in contact with TiO 2 -NPs were studied. Three modes of connection with TiO 2 were discussed. The bidentate mode corresponds to the connection of the molecule with two Ti sites via two oxygen atoms (Fig. 26b). The chelate mode is a connection of a molecule with a single Ti site via two oxygen sites. The monodentate mode corresponds to the connection of a molecule with a single Ti site, which is realized via a single oxygen site (Fig. 26b). The TD-DFT/CAM-B3LYP/6-311+G* approach was used for an adequate description of charge-transfer excitations. It was demonstrated that the interacting systems show the redshift of the absorption compared to isolated molecules and clusters. The calculation of natural orbits indicates the significant charge-transfer character of the excitations [183]. Based on synthesized polycyclic 9,14-dihydro-4H-dithieno[2′,3′:2,3;3′,2′:10,11]-piceno [1,14,13,12bcdefgh]carbazole (DTPC) the in silico analysis and design for suitable electron withdrawing structures was performed [184]. The electronic structure calculatіоns have been performed by using the B3LYP and MPW1K functionals with the 6-311G(d,p) basis set. The molecular structure of investigated systems is typical for a dye used in DSSC elements (Fig. 27). Here the main part of a molecule is electron donor fragment, while X is an electron acceptor. The carboxylic acid that was attached to acceptor group X is an anchoring group since this group can provide strong dye-TiO 2 adsorption [184]. Geometry and electronic structure of the ground and excited states for systems with different X (Fig. 27) was predicted by means of the DFT method. The calculated by using functional B3LYP maximum absorption wavelength for parent dye (λ max =596 nm) agreed well with known experimental data (λ max =609 nm) [184]. Figure 27. The molecular structure of the parent dye studied in [184].
The most critical parameters that influence the DSSC device performance were calculated. Among them, the LHE (light-harvesting efficiency), the number of injected electrons, the electron-injection lifetime, and the quantum-energy alignment of the adsorbed dye molecule to TiO 2 components. The results of work demonstrated that by modifying the acceptor group (X) the optical properties of dyes and the energy alignment of dye-TiO 2 structure could be tuned [184].
The effects of electron donor group on electronic and optical properties of arylamine-based dyes for DSSC were investigated in [185] (Fig. 28). DFTB optimized geometries of a dye and dye-TiO 2 conjugates were used for TD-DFT excited states calculations. The following functionals were used: CAM-B3LYP/3-21G(d), LC-ωPBE/6-31+G(d), and mPWDandHPW91/6-31G(d), respectively. The results show that the functional ωPBE/6-31+G(d) with the damping parameter of 0.175 au -1 provides the best agreement with experimental values. New rigid triarylamine based dyes were designed. Some critical parameters, such as η(%) and LHE, as well as other characteristics, were calculated. It was demonstrated that the investigated optical properties of new systems match with the requirement of an effective sensitizer [185]. Figure 28. Schematic illustration demonstrating that new dye Dhk4 was designed by modifying the rigidity of its donor moiety via altering the type and position of binding in the donor group. This structure modification leads to distinct planarity of the dyes, which significantly affects its properties so that Dhk4 dye showed more redshifted and broadened absorption spectra owing to the enhanced coplanarity between its donor and π-bridge moiety, as compared with a parent dye DIA3. Adapted with permission from [185]. Copyright © 2016 American Chemical Society.
In the ref. [186], the effects of various anchoring groups on the electronic coupling of dye-to-TiO 2 were investigated by means of ab initio and molecular dynamics simulations. The set of anchoring groups includes the following structures: carboxylic acid (-COOH), biscarbodithiolic acid (-CSSH), phosphonic acid (-PO(OH) 2 ), sulfonic acid (-SO 3 H) and hydrixamic acid (-CONHOH), respectively (Fig. 29a). The TD-DFT/PBE0/6-31G(d,p) with the PCM model for tetrahydrofurane solvent was used for excited-state calculations. It was shown that different anchoring groups give noticeably different contributions to LUMO. For instance, the contribution of a CSSH group is more significant than that of the COOH moiety. Hence, the CSSH group is favorable for electron injection from a dye to TiO 2 . Also, the effect of anchoring groups on an aggregation of a dye was investigated by MD simulations (Fig. 29b). In addition, the role of anchoring modes in adsorption of retinoic acid (RA) Fig. 30a on TiO 2 was studied by DFT calculations [187]. Three anchoring modes of RA, which contains a carboxylic acid group, were considered, as shown in Fig. 30b. The analysis of the B3LYP/6-31G(d) optimized structures revealed (Fig. 30b) that among the three anchoring modes, the mode c involving surface Ti−OH groups was suggested to have the highest driving force for photo-induced electron transfer. In this bound geometry, the mixing between the RA's LUMO and the conduction band states of TiO 2 was also the highest. The calculated results were supported by EPR measurements, which showed enhanced light-driven charge separation upon increasing the density of the Ti−OH groups on the surface of TiO 2 . It was found that the anchoring modes affected the bond length, degree of conjugation, energies of the HOMO and LUMO, and the mixing of the dye's excited states with the conduction band states of a semiconductor (Fig. 30e-g).
The influence of the chemical environment at the electrolyte/dye/TiO 2 interface on the adsorption geometries and electronic properties of a dye adsorbed on TiO 2 was considered by first-principles calculations [188]. Various adsorption modes of a series of organic dyes (Fig. 31a-b) on the (101) surface of anatase TiO 2 were examined in the absence and the presence of explicitly coadsorbed solvent (water or acetonitrile) molecules. Periodic DFT calculations were carried out within the generalized gradient approximation (GGA) using the PBE exchange-correlation functional [188]. It was found that coadsorbed solvent molecules introduced important modifications on the dye adsorption geometry with respect to the geometry calculated in vacuo (Fig. 31c-f). In particular, the bonding distance of the dye from the Ti anchoring atoms increased, whereas the adsorption energy decreased. Moreover, the adsorbed solvent induces the deprotonation of the dye due to the changing the acid/base properties of the system. Analysis of the electronic properties for the dye-sensitized TiO 2 structures in the presence of coadsorbed solvent molecules revealed an upward shift in the TiO 2 conduction band from 0.5 up to 0.8 eV in acetonitrile [188]. It should also be noticed that a similar shift of the TiO 2 conduction band was estimated for a solvent monolayer on unsensitized TiO 2 .  A joint experimental and computational study of the optical and charge transfer properties of D102 dye in the presence of Li-TFSI and EMIM-TFSI additives, as well as at the dye-sensitized TiO2 interface, was reported [181]. Experimental UV adsorption spectra showed that the spectral shifts of D102 are negligible for the dye in solution. However, upon moving to the TiO 2 film these additives resulted in some slight red-shifts in spectra, accompanied with the appearance of a lower-energy shoulder in the absorption spectrum. B3LYP/6-311G** and PBE0/6-311G** calculations showed the formation of weak dye/additive complexes in solution and also predicted spectral red-shifts. Moreover, the strong effect of these additives on the electronic coupling between the dye's LUMO and the TiO 2 conduction band states, approximated by (TiO 2 ) 82 (101) two-layer anatase slab [181].
The systematical theoretical study of seven D−A−π−A-based indoline (IND) dyes (Fig. 33a) was comprehensively investigated via quantitative-structure−property relationships to evaluate their prospect of application in DSSCs [195]. Adapted with permission from [195]. Copyright © 2019 American Chemical Society. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) Optoelectronic properties of the isolated dye and dyes adsorbed on a TiO 2 cluster were explored by DFT (B3LYP and CAM-B3LYP with a combined 6-31G(d,p)/LANL2DZ basis set) and TD-DFT methods (Fig. 33b-c). Light absorption spectra, vertical dipole moment, shift of the conduction band of the semiconductor, excited-state lifetime, a driving force of electron injection, the photostability of the excited state, and exciton binding energy were computed. The authors showed that the presence of an internal acceptor such as pyrido [3,4-b]pyrazine would influence greater the open-circuit voltage (V OC ), compared to the benzothiadiazole moiety (Fig. 33d). Considering the balance between the V OC and J SC (short circuit current), the IND3, IND5, and IND10 dyes were suggested as the most suited to be utilized as potential candidates for DSSC [195].
Different mechanisms for electron injection in DSSC are discussed in ref. [196]. The aim of work was to investigate structural aspects of a dye, which lead to specific electron injection to TiO 2 mechanism. The classification of sensitizers was proposed. Type I of sensitizers corresponds to the indirect mechanism. For these systems, the orbitals involved in excitations mainly localized in the molecule. Type II is characterized by virtual orbitals mostly localized on a dye that appears well into the conduc-tion band of the semiconductor. These orbitals are formed via the interaction of d z2 orbitals of Ti atom and p-orbitals of the oxygen atom of the dye [196].
A series of Ru(II)-based dyes, labeled as N dyes, such as N3 (Fig. 25), N719, N749, etc. have become prototypical benchmark systems for ab initio calculations sensitizing effect to TiO 2 [197][198][199]. Earlier computational investigations of the absorption spectrum and the alignment of the ground and excited state energies for an N719 dye, adsorbed on an extended TiO 2 anatase slab exposing the majority (101) surface, were carried out by DFT/TDDFT calculations (Fig. 34a) [198]. The geometry of the dye-TiO 2 conjugate was optimized by using the PBE exchange-correlation functional together with a plane wave basis set and ultrasoft pseudopotentials [199]. TDDFT calculations were performed on the optimized geometries, employing the hybrid B3LYP functional with a 3-21G* basis set. A polarizable continuum model of solvation (C-PCM) was employed to treat the effect of the surrounding water solvent. It was found that the calculated and experimental absorption spectra for the dye-TiO 2 conjugate were in excellent agreement with an absorption maximum deviation below 0.1 eV (Fig. 34b). The lowest optically active excited state was found to lie ca. 0.3 eV above the lowest TiO 2 state. The latter was characterized with a sizable contribution from the dye π* orbitals, strongly mixed with unoccupied TiO 2 states (Fig. 34c). Therefore, an ultrafast electron injection component was predicted on the basis of the strong coupling and of the matching of the visible absorption spectrum and density of TiO 2 unoccupied states (Fig. 34d). Comparison between the experimental (red) and calculated (blue) absorption spectra of N719 on TiO 2 . For clarity, the experimental spectrum intensity has been normalized to the calculated absorption maxima. Blue vertical lines represent unbroadened transition wavelengths and oscillator strengths. (c) Isodensity plots of the HOMO and LUMO+11 of the N719-TiO 2 conjugate. Bottom: Charge density difference between the ground state (S 0 ) and the S 18 excited state. A blue (yellow) color signifies and increases (decrease) of charge density upon electron excitation. (d) Left: alignment of the ground and excited-state energy levels for the N719-TiO 2 conjugate. Right: the calculated density of unoccupied TiO 2 states, aligned to the energy of the lowest TiO 2 state in the conjugate. Adapted with permission from [198]. Copyright © 2011 American Chemical Society. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) Successes, failures, and challenges of quantum chemical calculations of excited states and alignment of energy levels in dye-sensitized solar cells were critically analyzed in the context of the potential and limitations of current DFT and TDDFT computational methods to model DSCs [197]. Pastore, Fantacci, and De Angelis considered the ground and excited-state properties of both isolated and TiO 2adsorbed metallorganic and fully organic dyes, relevant to modeling the dye→semiconductor electron injection process. By reviewing previous data from their laboratory, combined with new calculations (Fig. 35a-b), they suggested that ruthenium dyes are typically accurately described by standard DFT approaches [165,[197][198]. However, in the case of highly conjugated organic dyes, which are often characterized by strong charge-transfer excited states, some specifically tailored exchange-correlation functionals, such as CAM-B3LYP, are required. It was shown that for donor−acceptor organic dyes, a simple picture based on the dye lowest unoccupied molecular orbital (LUMO) broadening accounts for the different interfacial electronic coupling exhibited by dyes with different anchoring groups. It was also concluded that the explored DFT/TDDFT functionals, such as B3LYP, MPW1K, and CAM-B3LYP, estimate a balanced description of the dye-TiO 2 excited states and of the alignment of the dye excited states with the semiconductor manifold of unoccupied states with different accuracy (Fig. 35c) [197].  [197]. Copyright © 2013 American Chemical Society.
A fully atomistic computational approach was suggested to capture the differences between direct (type-I) and indirect (type-II) photoinjection mechanisms by a time-dependent density functional tightbinding (TD-DFTB) approach [200][201][202][203]. The type-I photoinjection involves two steps (Fig. 36a): (i) an excitation from the ground state to the excited state of the dye, which is produced by the absorption of a photon, (ii) an electron transfer to the conduction band of the semiconductor nanoparticle (NP). Whereas the type-II photoinjection is referred to as the one-step electron injection from the ground state of the dye to the conduction band of the semiconductor upon photoabsorption. It should be noted that, in the type-II mechanism, a new charge-transfer band appears, which characterizes the direct electron excitation into the conduction band and can be observed in the UV−vis absorption spectrum. On the opposite, in the type-I mechanism, no new bands are expected in the UV−vis spectra. A model anatase TiO 2 NP (a 270 atom (90 TiO 2 units)) functionalized with different dyes, such as alizarin, coumarin C343 (C343), catechol, cresol, [Ti(Ph-tetra-t-Bu)(catechol-CO 2 H)], aniline derivative, and naphthalenediol, was considered. The equilibrium structure of the anatase TiO 2 NP was taken from an equilibrated MD simulation at 300 K by using the MA FF parameters [200][201]203]. It was shown that the nature of the photoabsorption process in the dye-TiO 2 NP can be understood in terms of orbital population dynamics occurring during photoabsorption (Fig. 36b) [200]. Fig. 36c-d shows that alizarin and coumarin C343 exhibited an indirect (type-I) mechanism for electron injection because no new absorption band was observed in the visible region upon adsorbing the dye onto the TiO 2 NP. The lowest-energy band suffers a redshift for the alizarin−NP, while coumarin C343 shows a slight shift to higher energies. The lack of energy shift in the case of C343 might be understood in terms of the weak electronic coupling between the excited states of the dye and the anatase NP [200]. The proton transfer at a dopamine-functionalized TiO 2 interface was investigated by using dispersion-corrected hybrid DFT calculations and DFT tight-binding (DFTB) MD simulations [204]. The adsorption modes, patterns of growth, and configurations of dopamine on the anatase (101) TiO 2 surface were considered with reference to the binding archetype of catechol. It was found that at low cov-erage, the isolated dopamine molecule prefers to bend toward the surface, coordinating the NH 2 group to a Ti 5c ion. Whereas, at high coverage, the packed molecules were succeeded in bending toward the surface only in some monolayer configurations so that a proton transfer from the surface to the ethylamino group occurred. It was suggested that improving the probability of dopamine molecules being free to bend toward the surface through thermodynamic versus kinetic growth conditions should lead to the formation of a monolayer of fully protonated dopamine molecules [204].
Self-assembly of organic molecules on semiconductor surfaces is of crucial importance in the performance of DSSC devices. Several computational works were focused on a strong tendency of organic dyes to form molecular aggregates on the TiO 2 surface [204][205][206].
The optical and electronic properties of dye aggregates of p-methyl red on a TiO 2 anatase (101) surface were considered as a function of aggregation order varying from a monomer up to pentameric dyes (Fig. 37) by using dispersion-corrected DFT calculations (PBE functional and ultrasoft pseudopotentials) [205]. A progressive red-shifting accompanied by the intensity increase toward the visible region in UV−vis absorption spectra were observed from monomeric-to-tetrameric dyes. The calculated density of states (DOS) and partial DOS spectra reveal similar dyeꞏꞏꞏTiO 2 nanocomposite conduction band characteristics; however, different valence band features were revealed. Associated molecular orbital distributions demonstrated dye-to-TiO 2 interfacial charge transfer in all studied aggregates.

Figure 37.
Optimized structural geometries of different adsorbed dye aggregates. Adapted with permission from [205]. Copyright © 2014 American Chemical Society. Adsorption modes, dynamics, self-aggregation, and optical properties of an indoline D102 dye adsorbed on TiO 2 anatase (101) substrates were studied by reactive MD (ReaxFF) simulations combined with TD-DFT calculations [206]. The morphology and electronic properties of dye-sensitized heterogeneous interfaces were examined at an atomistic level. Extensive MD simulations of the adsorption of D102 dyes on medium/large size TiO 2 anatase slabs allowed to monitor a variety of dye binding modes. It was suggested that the main driving forces toward the formation of ordered surface aggregates of D102 dyes are π stacking and T-shaped interactions between the aromatic rings of the donor moiety of the dyes, as well as the tendency to maximize the anchoring points with the surface. The dye aggregates were found to be organized in domains, which are characterized by a different orientation of the packing units with a certain degree of short-to-medium range order [206].
Finally, at the end of this section, we suggest a few methodological conclusions from the results described above: 1) The results of the electronic structure calculations are insensitive to the geometry optimization method. The tight-binding DFT approach (DFTB) can be a realistic low-cost alternative to regular DFT geometry optimization.

Ab initio MD Investigations
In the current section, we are briefly discussing the most interesting papers where ab initio MD simulation methods were predominately used for the theoretical investigation of the pristine TiO 2 fragments or TiO 2 -based functional materials in the last two decades.
Ab initio (quantum or first principles) MD is the most powerful molecular modeling method for investigating the structure and dynamics of many-body systems in the framework of the condensed matter physics and computation quantum chemistry. The main idea underlying in any first principle MD method is to compute the true forces acting on the nuclei from proper electronic structure calculations that are performed "on the fly" as the MD trajectory is generated [207]. The three basic techniques are traditionally distinguished among an ab initio MD simulations: The first one is the conventional MD of nuclei using the quantum chemistry predicted Born-Oppenheimer (BO) potential energy surface (BOMD method). In the second MD technique, that was proposed by Car and Parrinello (CP), instead of solving the electronic structure problem at each new nuclear configuration, the so-called fictitious dynamics for the electronic orbitals (wave functions) is introduced that, as before, allows adiabatic motion of the nuclei (CPMD method) [208][209]. Both mentioned approaches deal with the ground state of an electronic sub-system. The third specific ab initio MD technique was designed for the solution of the more complicated non-adiabatic (NA) problems, such as electronic excitation and relaxation or particle collision and scattering (NA or Ehrenfest MD method) [210].
Earlier investigations of the photoinduced electron transfer from a molecular donor to the TiO 2 semiconductor acceptor, using the NAMD simulation method, were carried out in a series of works by Stier and Prezhdo [211][212][213]. The small systems under study are typically composed of dye-sensitized semiconductor nanomaterials, used in solar cells, photocatalysis, and photoelectrolysis. The electronic structure of a dye-semiconductor system and the adiabatic dynamics were simulated using BOMD, while the NA effects were incorporated by the quantum-classical mean-field approach that neglects the quantum effects of the ionic motion, in particular, the zero-point energy and decoherence effects. A novel procedure separating the NA and adiabatic electron transfer pathways (Fig. 38) was also developed. The simulation provided a detailed picture of the ET process. It was initially assumed that ground and excited state ion dynamics are similar. So, the change in the quantum force due to the electronic photoexcitation could be neglected and the following analysis might be greatly simplified. An important result was that thermally driven adiabatic transfer is the dominant electron transfer mechanism. Compared to the earlier simulation at 50 K, the rate of NA transfer at 350 K remained almost unchanged, whereas the rate of adiabatic electron transfer increased substantially [211].
In the next two publications, the ultrafast dynamics of the photoinduced electron transfer process occurring between a model organic sensitizer molecule (isonicotinic acid) and a TiO 2 semiconductor surface was simulated and discussed in details (Fig. 39). The model used in simulation exhibited the key features of the larger real systems studied before experimentally while providing a detailed simulation picture of the electron transfer mechanism [212][213].
Finally, the authors indicated that the NA electron transfer entirely dominated at short times and could occur due to the presence of the strong localized avoided crossing region, as well as the extended region of weaker NA coupling. Although the adiabatic electron transfer contribution accumulated more slowly, it approached that of the NA transfer pathway asymptotically. As a result, the simulation predicted a complex non-single-exponential time dependence of the electron transfer process [213]. An interesting theoretical investigation in the field of phase transformation mechanisms from anatase to the TiO 2 nanosheets was carried out by Alvarez-Ramirez and Ruiz-Morales [214]. Several ab initio MD simulations with the non-self-consistent Harris functional (so-called non-self-consistent BOMD) including (i) the Γ-point approximation and (ii) the k-points were performed [214]. The three different size wire models of anatase with the periodic conditions along the (0, 1, 0) direction were used and also several different environmental conditions were investigated.
It was found the titanate-type structures are energetically favorable. The transformation (Fig. 40) was followed by monitoring the formation of four-coordinated oxygen atoms, which constituted the characteristic sequence of edge-sharing TiO 6 octahedrons of titanates. In addition, it was shown that the titanate-type structure could be transformed to anatase in a reversed mechanism. An essential factor in the transformation of the anatase-type sections to the TiO 6 octahedrons was the interaction of the (0, 0, ±1) surfaces of the anatase with positive counter ions, which stressed the surface, deforming it and promoting the transformation [214].
Outstanding NAMD investigation of the so-called "wet-electrons" on the TiO 2 surface was carried out by Prezhdo and coworkers [215]. Simulations directly mimicked real time-resolved experiments and revealed the nature of electron transfer in the "wet-electron" system. Earlier, the experimental and theoretical studies were done on the TiO 2 (110)/H 2 O interface have shown a short-lived state around 2.4 eV above the Fermi level. The electron in this state was localized on the surface and was partially hydrated (Fig. 41). Therefore, it was named the "wet-electron". The time-resolved studies of the "wetelectrons" provide a fundamental understanding of the solvent role in the DSSC processes. The "wetelectron" state is higher in energy than the excited states of the majority of chromophores used in DSSC devices. Therefore, pristine "wet-electron" states are not directly involved in the photovoltaic processes [215]. It was shown that thermal phonon motions induced a large fluctuation of the "wet-electron" state energy, generated frequent crossings of the donor and acceptor states, and drove the adiabatic mechanism. The rapid phonon-assisted NA tunneling from the "wet-electron" state to the TiO 2 surface was facilitated by the strong water-TiO 2 electronic interaction. Besides, it was demonstrated that the motions of H 2 O molecules had a greater effect on the electron transfer dynamics than the hydroxyl vibrations [215].

Figure 41.
Charge densities for the "wet-electron" state (left), the surface state (middle), and the bulk state (right) of the "wet-electron" system. The charge is transferred sequentially from the "wet-electron" state through the surface state to the bulk state. The large delocalization of charge around the water molecules in the "wetelectron" state is shown, which is due to the dangling hydrogen atoms that help stabilize the electron. Adapted with permission from [215]. Copyright © 2009, American Chemical Society In order to better understand the electron-hole relaxation of nano-scale TiO 2 structures, it was important to start with an understanding of TiO 2 synthesis building blocks. Such investigation was carried out by Vogel and Kilin [216], in which the Ti(OH) 4 complex with tetrahedral coordination was examined. To simulate the electronic properties of the Ti IV ion in solution, the model was surrounded with explicit 27 H 2 O molecules (Fig. 42). The model was explored by means of standard DFT BOMD followed by NA electron dynamics computed with a reduced density matrix approach combined with "on-the-fly coupling" [216]. According to the MD simulation results [216], the titanium ion is affected by coordination to hydroxyl groups and an explicitly modeled water environment in three ways: First, the crystal field environment stabilizes the oxidation state. Second, the crystal field provides splitting of the titanium ion dorbitals, which is responsible for the density of states in the conduction band range of energies. Third, thermal motion of near hydroxyl ligands and water coordination sphere are responsible for the nonradiative relaxation mechanisms: a photoexcitation in the solvated titanium ion is followed by the energy dissipation from electronic degrees of freedom to thermal motion. Thus, it was found that electrons dissipate energy slower than holes [216].
Chemical and physical adsorption of water on small-sized wet TiO 2 nanoparticles was investigated by ab initio MD simulations [217]. The wetting behavior of (TiO 2 ) 24 ꞏnH 2 O (n=0, 1, 3, 8, 15, 30) was studied as a function of water content, ranging from dry nanoparticles to wet nanoparticles with monolayer coverage of water. The simulations were performed with the mixed Gaussian and planewave (GPW) method, as implemented in the CP2K code. It was found that the surface reactivity was shown to be a concave function of water content and driven by surface defects [217]. It was suggested that the local coordination number at the defect of TiO 2 NPs was identified as the key factor in deciding whether water adsorption proceeds through dissociation or physisorption on the surface.

Hybrid QM-MD Simulations
In the last decade, hydride QM-MM calculations utilizing the DFTB approximation have become an attractive alternative to conventional non-empirical first principle ab initio MD simulations. Several examples of such studied are summarized below.
Small TiO 2 NPs (anatase, rutile, brookite) have been studied by ab initio MD simulations to elucidate chemical and physical adsorption of water molecules on the TiO 2 -water interface, ranging from dry NPs up to monolayer solvated NPs [217]. QM-MM simulations have been used to model the adsorption of insulin on TiO 2 (100) surfaces [218]. Recently, a multistep and multiscale investigation based on hybrid DFT, DFT tight-binding (DFTB), and QM-MM calculations have been applied to consider solvation behavior of curved TiO 2 -NPs, ranging from vacuum, multilayered solvation and up to the bulk water environment, resulting in an overall stoichiometry of (TiO 2 ) 223 ꞏ8958H 2 O (Fig. 43) [219].
Recently, DFT-based BOMD was carried out by Byrne and English to investigate the effect with which the choice of exchange-correlational functional has on the structural properties of a [BMIM] + [NTF 2 ] − RTIL solvating an N719 sensitizing dye, adsorbed onto the anatase TiO 2 (101) surface [220]. The systematic study of an N719 dye, non-solvated and solvated by RTILs, was performed by comparing the widely used BLYP and PBE functionals, which were utilized with and without D3 dispersion corrections by Grimme [220]. The above results show that neither the addition of dispersion corrections nor solvation by the strongly ionic solvent has any noticeable effect on the anatase surface structure. The structure of the surface was stable throughout the MD simulations, as is the carbon backbone of an N719 dye with neither group showing any real differences in their structure upon change of the DFT functional. However, the pair distribution functions (PDFs) for some of the carboxylate oxygens showed that the results obtained with the PBE functional led to some structural peaks shifts by 0.01 nm closer than their position derived by the BLYP functional [220]. The metal complex at the center of the dye was also largely unaffected by utilization of the dispersion correction, variations in the functionals or the presence of the solvent. The structure of the dye was found to be essentially unaffected by solvation with the bipyridine rings and anchoring groups retained their rigid geometry. A comparison of the PDFs, derived from the BLYP and PBE functionals, showed the structure of the anions and cations to be quite similar. On the other hand, some deviation was found between the results of the BLYP-D3 and the PBE-D3 so that the cations and anions are seen to be slightly closer to each other in the case of the PBE-D3 calculation by about 0.02 nm. Moreover, it was suggested that the dispersion corrections also introduced the pronounced effect of an orientational alignment in the liquid, the phenomenon which would be expected for a real liquid [220].
Stable structures and photoexcitation character of adsorption of an N749 dye to TiO 2 anatase (101) interface, solvated by explicit acetonitrile (ACN) molecules were considered by DFT-MD simulations [221]. It was found that the adsorption of an N749 dye via its deprotonated carboxylate two anchors (type d2) was the most stable at the interface (Fig. 44a). However, the adsorption mode with the one protonated carboxyl anchor (type p1) was observed to be at the average energy only slightly higher than that of type d2 (Fig. 44b).
Other possible anchoring modes, shown schematically in Fig. 44c-e were all higher in energy as compared to those of p1 and d2, respectively. It was suggested that, under equilibrium, the mode p1 might still coexist with the mode d2 [221]. These findings were in contrast with the high stability calculated for the mode p1 in vacuo. It was proposed that some stabilization of the mode d2 could be caused by inhomogeneous charge distribution and anchor fluctuation enhanced by ACN solution (Fig.  44f-g). The calculated projected densities of states and the photoabsorption spectra supported the conclusions that the mode d2 has the larger driving force of the electron injection from an N749 dye into the TiO 2 , whereas the characteristic sensitizer photoabsorption in the wavelength region over 800 nm, was mainly attributed to the mode p1 even in the bulk ACN solution. It was concluded that the better performance of a N749 dye in DSSC devices can be regionalized in terms of the co-sensitizer framework of the adsorption modes d2 and p1, respectively [197][198]221].

Figure 44.
Representative snapshots of the equilibrium trajectories of adsorption of a N749 dye to TiO 2 anatase (101)/acetonitrile (ACN) interfaces: (a) adsorption with one protonated carboxyl anchor (−COOH), labeled as p1, (b) adsorption with two deprotonated (−COO -) anchors, labeled as d2, (c) adsorption with two protonated anchors (label p2), (d) adsorption with one protonated and one deprotonated anchors (label p1d1), and (e) the desorbed state, assigned as D. Each trajectory was obtained with DFT-MD simulation with different initial configuration of the dye, which was optimized on the vacuum surface, in advance. (f-g) Average structures and distances (in angstroms) of the characteristic atoms near the interface in (f) the p1 equilibrium state and (g) the d2 state, respectively. Adapted with permission from [221]. Copyright © 2014, American Chemical Society.
Small TiO 2 NPs, functionalized with a supramolecular complex consisting of an entirely organic naphthalene-diimide (NDI) dye covalently bound to a mononuclear Ru-based water oxidation catalyst (Fig. 45a) were studied by ab initio MD simulations to elucidate the mechanistic insight and microscopic details of water oxidation at the photoanode induced by visible light absorption [222]. The MD nuclear trajectory was calculated through ab initio MD by using the Car−Parrinello MD (CPMD) code. The orthorhombic cell, containing a two-layer anatase slab (TiO 2 ) 32 functionalized with the NDI1 chromophore, was used together with periodic boundary conditions to model the TiO 2 surface. The fast photoinduced electron injection from the NDI into the semiconductor was observed, which provided the driving force for the activation of the Ru catalyst (Fig. 45b). The explicit description of the water environment was utilized to elucidate the sub-picosecond time-scale and the nature of the proton-coupled electron transfer of the catalytic reaction path. It was suggested that the explicit water treatment was essential to determine the proton diffusion channel and the free energy change along with the reaction [222].
The above photochemical reaction simulations are a promising example for the computational design and DSSC device optimization with an aim to explore the effect of different anchoring groups, which allow tuning the driving force for the subsequent steps of the catalytic water photo-oxidation cycle and reducing the charge recombination rate [223][224][225]. , the molecular chromophore NDI, and the Ru-based water oxidation catalyst; (NDI 1 ) chromophoresemiconductor subsystem employed for the photoinduced electron injection MD simulation. (b) Electron injection profile obtained through the time-dependent population analysis of the wavepacket projected only over the dye-TiO 2 conjugate (solid line). The insets show the distribution of the total wavepacket after initialization (i), and along with the time evolution of the MD trajectory (ii and iii). Adapted with permission from [222]. Copyright © 2016, American Chemical Society.

Conclusions and Outlook
The review describes recent advances in theoretical studies of TiO 2 nanomaterials. Particular attention is focused on the MD simulation approaches that are essential for the modeling of bulk TiO 2 , its interfaces with aqueous and non-aqueous environments, as well as large-scale nanoparticles and their aggregates. Starting from the earlier FF model by Matsui and Akaogi (MA) [14], rapid progress in computational modeling of TiO 2 nanomaterials resulted in the development of well-validated FFs for a broad range of TiO 2 nanomaterials, ranging from small perfect clusters and up to water/TiO 2 interfaces. The applicability of MD simulation tools has currently extended up to about 100 nm size and microseconds, and can further be scaled up by coarse-grain and finite element models to the range of micrometers. However, the critical overview of recent computational models suggested for modeling partially hydroxylated TiO 2 surfaces reveals that the existing MA approach still requires further improvements and validation. The main limitation of the MA force-field, developed originally for bulk TiO 2 , is that it utilizes only two types of atoms, such as Ti and O. However, in realistic TiO 2 nanostructures, these atoms are characterized by different coordination numbers so that they can be classified into inner-bulk and outer-surface-exposed atoms. The latter can be further divided into subgroups of corner, edge, or side atoms. Recent ab initio MD studies demonstrated that indeed these atoms bear different partial charges and, hence, they should be represented by different atom types in an empirical force field [217]. The different atoms types for Ti and O have already been implemented in some MD simulation studies of TiO 2 nanostructures [19,32,84]. Moreover, the development of computational models for some realistic systems with the simultaneous presence of several crystalline facets (100,110,111) or different phases, such as rutile and anatase, will enable modeling the adsorption phenomena that may take place at the grain boundary and at their common interface with the outer environment [226].
The understanding of the mechanisms that govern the interaction of biological systems with inorganic materials remains a challenge for both fundamental and applied material science [227][228]. The protein adsorption modulates the formation of biofilms onto surfaces, which is essential for medical implants and dental technologies. The interaction of TiO 2 nanomaterials with biomacromolecules is crucial to further response of cells to foreign inorganic materials, such as engineered implants or accidentally entrapped inorganic nanoparticles. Therefore, detailed knowledge of the TiO 2 /biological interfaces is required for the design of new biocompatible materials. MD simulations of conformation equilibria of biomolecules at perfect, defect-free, and curved interfaces have provided insight into the mechanism, specific recognition, and sequence-relation suggestions to understand their adsorption, binding, and assembly. However, due to the complexity of the configuration space, it remains challenging to explore because of the large number of possible conformations of biomacromolecules. Therefore, the critical review of MD studies in this field reveals that free, unbiased MD sampling of protein binding and adsorption onto TiO 2 may not guarantee finding of proper binding modes so that the development of some innovative sampling techniques is required. Some recent studies report new promising approaches in this direction, such as the potential mean-force calculations, free-energy umbrella sampling, and replica-exchange MD simulations [229].
In the last decade, significant progress has been achieved in quantum chemical calculations of the structural, electronic, and optical properties of TiO 2 materials in the context of their use as components of dye-sensitized solar cells (DSSC). Successes, failures, and challenges of quantum chemical, DFT, and TDDFT computational methods of dye-TiO 2 conjugates have been critically analyzed elsewhere [1,197,230]. To summarize, we found that the electronic structure calculations are little sensitive to the basis set and DFT functionals used for geometry optimization. Moreover, the tight-binding DFT approach (DFTB) can provide a low-cost alternative to regular CPU-expensive DFT geometry optimization. However, it was found that the results of electronic property calculations are crucially sensitive to the choice of the DFT functional [165,197]. The presented DFT/TDDFT computational models, employing conventional hybrid (B3LYP) and range-separated functional and based on the use of large TiO 2 clusters, PCM-continuum solvation models, allows one to reproduce to a reasonable degree of accuracy the electronic and optical properties of dye-TiO 2 systems. However, for highly conjugated organic dyes, which are often characterized by a high degree of charge-transfer excited states, specifically tailored exchange-correlation functionals, such as CAM-B3LYP or M06-2X, are required [197]. Some efforts have already devoted to solving the shortcomings of current TDDFT methodologies for dye-TiO 2 systems [231].
Finally, our review sheds light on the progress and development of modern computational tools, combining classical FF modeling and high-level ab initio MD simulations, which may assist the molecular engineering and rational design of more efficient dye sensitizers and deeper understanding of the electron transfer dynamics in dye-sensitized solar cells.