Theoretical insights into the photo-protective mechanisms of natural biological sunscreens: building blocks of eumelanin and pheomelanin

Abstract

Eumelanin (EM) and pheomelanin (PM) are ubiquitous in mammalian skin and hair – protecting against harmful radiation from the sun. Their primary roles are to absorb solar radiation and efficiently dissipate the excess excited state energy in the form of heat without detriment to the polymeric structure. EU and PM exist as polymeric chains consisting of exotic arrangements of functionalised heteroaromatic molecules. Here we have used state-of-the-art electronic structure calculations and on-the-fly surface hopping molecular dynamics simulations to study the intrinsic deactivation paths of various building blocks of EU and PM. Ultrafast excited state decay, via electron-driven proton transfer (in EU and PM) and proton-transfer coupled ring-opening (in PM) reactions, have been identified to proceed along hitherto unknown charge-separated states in EU and PM oligomers. These results shed light on the possible relaxation pathways that dominate the photochemistry of natural skin melanins. Extrapolation of such findings could provide a gateway into engineering more effective molecular constituents in commercial sunscreens – with reduced phototoxicity.

---

1. Introduction

Ultrafast internal conversion (IC) in photoexcited molecules has received much attention in the past two decades.^1–12 Such phenomena are spin-allowed processes and constitute efficient radiationless decay pathways in many photoactive biological molecules.^10,13–16 Degeneracies between potential energy surfaces (PESs), in particular conical intersections (CIs), are crucial to enabling population transfer between different electronic states – thus mediating efficient IC between PESs of the same spin multiplicity.^17,18 To date, diverse ranges of nuclear geometries at CIs have been identified in cyclic and aromatic molecules – all of which involve ring-deformations, isomerisms and/or bond extensions.^2,6–8,15 In this article, the topographic details of the PES in the vicinity of curve crossings along bond extension coordinates are reported in the context of photoprotection mechanisms occurring in natural biological sunscreens present in mammals.

Natural biological sunscreens are present on parts of the human body (e.g. the skin or the surface of the eye) that are susceptible to solar irradiation and, as such, attenuate the near-IR, visible and most importantly, the near-UV components within the solar spectrum. Their presence helps protect against skin cancers like melanoma or squamous cell carcinoma. Natural sunscreens generally comprise long chain helices of coiled conjugated organic molecules that offer photoprotection without degradation. These molecules offer photoprotection by virtue of their high absorption cross-sections and high IC efficiencies.¹⁹ Upon near-UV, visible or near-IR absorption by a biological sunscreen, rapid dissipation of the excess excited state energy is required in order to ensure photostability. This energy dissipation is usually in the form of IC to the ground electronic state, followed by cooling of vibrational energy – which dissipates as heat.

The present study explores the anatomy of possible IC routes in selected chromophoric sections of two of the most abundant biological sunscreens present on mammalian skin: eumelanin and pheomelanin (henceforth referred to as EU and PM, respectively). EU is a dark pigment, the most abundant constituent of melanin present in human skin and hair and comprises of polymeric chains of straight and cross-linked 5,6-dihydroxyindole, indole-5,6-dione, hydroxyindolone and 5,6-dihydroxyindole-2-carboxlic acid motifs – some examples of which are displayed in Fig. 1(a). It is responsible for attenuating harmful UV radiation from the sun. The precise polymeric sequence of EU is still under debate but an example of a widely accepted motif is E/F in Fig. 1(a), which is the linear trimer present at the end chain of the helical structure of EU. Various tetrameric (e.g.G), pentameric and octameric cyclic structures have also been previously proposed as likely sub-units within the three-dimensional crystal structure of EU.^20–23 Oligomer G, for example, represents a key geometric form by virtue of previous theoretical investigations suggesting that the possible presence of such oligomers allow for the efficient uptake and delivery of metal cations within the polymeric chain.²¹ Selected oligomers of EU are also known to π-stack – giving rise to low-lying exciton states following photo-irradiation.²² Such π-stacked melanins are beyond the scope of the present study but will form the basis of an in-depth future study.

Fig. 1 Depictions of exemplar sections through the complex polymeric structures of (a) eumelanin and (b) pheomelanin the ground state minima of which were optimized at the MP2/cc-pVDZ (for A, B, C, D, E, F and PM) and DFT/B3LYP/6-31G (for G) levels of theory. The atom colourings indicate carbon (grey), oxygen (red), hydrogen (white), nitrogen (blue) and sulphur (yellow). Structures within (a) represent possible arrangements (i.e. linear trimer and cyclic tetramer, respectively) of EU – on which ab initio calculations were undertaken. The section through PM was scaled down to save computational expense. The exact structure of the PM section on which in-depth theoretical analysis was undertaken is shown in (b).

In contrast, PM is much less abundant in natural melanin and comprises a pink pigment that is found on the epidermal layer of the mammalian body. It consists of oligomer chains comprising benzothiazone and benzothiazole units – a common repeating motif of which is displayed in Fig. 1(b).

EU and PM both have extensive π-conjugation, which has the effect of red-shifting the absorption so as to capture the near-IR, visible and near-UV portions of the solar spectrum (cf. monomeric aromatic systems). Computational studies on the 5,6-dihydroxyindole monomer, as well as their associated dimers and oligomers have attracted much attention – in each case contributing towards revealing the dominant short wavelength (mid to near-UV) photophysics of EU.^21,24–26 Experimentally, EU has attracted vast attention, the most notable of which includes the ultrafast time-resolved experiments of Nofsinger et al.^27–29 and Sundström and co-workers^30–34 in solution. In contrast, the possible fates of the excited states of building blocks of PM have received comparatively less photochemical attention from the ultrafast community. That said, time-resolved fluorescence up-conversion experiments by Sundström and co-workers on derivatives of benzothiazole show rapid deactivation of the excited state on an ultrafast time-scale.³⁵ Simon and co-workers have also reported similar short lifetimes for the excited states of polymeric PM by time-resolved transient-absorption spectroscopy measurements.^36,37

We note that all melanin constituents in the present study contain proton donors (a hydroxyl group) and acceptors (the carbonyl oxygen atom) aligned in such a way as to encourage intramolecular proton transfer – which is the primary reaction path of focus in the present study. Such excited-state proton transfer reactions have already been suggested in monomeric building blocks of EU^25,38 and have been identified as key relaxation pathways underpinning the photostability of many biomolecules.^{1,2,5,13,14,39}

Though we recognize that the presently studied molecular systems are small sections through the polymeric structures of EU and PM, the present study is nevertheless informative as we have considered all key molecular constituents associated with possible excited-state proton-transfer reactions. In PM, we have deliberately ignored the benzothiazole unit in an attempt to reduce computational expense. Whilst we recognize that the benzothiazole moiety may encourage alternative intrinsic decay paths, its local position is too remote from the reaction center associated with proton-transfer (see the full molecular structure of PM in Fig. S1.1 of the ESI†).

The present work reports a computational study of the energetics associated with the dominant decay pathways in possible oligomeric building blocks of isolated EU and PM sub-units, the results of which enable predictions of the most probable photoprotection mechanisms in each case. Static electronic structure calculations, using the second-order algebraic diagrammatic construction (ADC(2)) method were employed in order to ascertain the minimum energy paths (MEP) towards the intrinsic decay mechanisms. In addition, nonadiabatic surface hopping molecular dynamics simulations were employed in order to determine both the relevance of the paths identified by electronic structure theory and to provide timescales for these possible IC routes.

2. Results and discussion

2.1 Ground state geometries

Eumelanin. Fig. 1 presents the possible ground state minimum energy conformers of selected dimeric (A–D), trimeric (E and F) and cyclic tetramer (G) forms of EU. These assigned letters will be used henceforth to identify the various structures where relevant.

Starting from the cyclic tetramer (G), such porphyrin-like arrangements within the polymeric sequence have been suggested to account for the efficient ion uptake and delivery measured in EU.^20,21 The minimum energy geometry of this cyclic tetramer is slightly puckered in order to accommodate the strain caused by cyclisation. Smaller subunits of G are also present in EU, an example of which is dimer C. Structure C can be viewed as the top half of G but, in contrast, retains a planar geometry and is predicted to exist in a linear or branched configuration within the EU polymer. NH–N hydrogen bonding also helps maintain the planarity of C, keeping this molecular arrangement rigid. Upon double-tautomerisation and subsequent E–Z isomerism, structure C can form D, which is calculated to be 830 cm⁻¹ more stable than C.

Structures A and B are syn–anti rotational isomers of each other and are positional isomers of C and D. A is 2570 cm⁻¹ more stable than B which is attributable to the loss of hydrogen bonding upon geometric isomerism, as has been observed before in, for example, catechol and caffeic/ferulic acid.^2,40 In contrast to C and D, the minimum energy geometries of A and B are non-planar to amend a steric clash between the bi-indole moieties. Such rotational isomerism is also present in structures E and F, which are analogues of A and B, with an addition of a pyrrole-2-carboxylic acid functionality – representing trimeric units. E and F represent the two lowest energy rotamers in an ensemble of isomers studied (see Fig. S1.2 of the ESI†). Within this ensemble, structure E is the most stable and 660 cm⁻¹ more stable than F.

Pheomelanin. Fig. 1(b) depicts the ground state minimum energy configuration of a simple repeating unit present within PM. As shown in Fig. 1(b), the chosen section consists of a thiazinoquinolinol (TQ) and benzothiazinone (BT) moiety both of which contain a single sulphur atom. The TQ and BT units contain, respectively, a hydroxyl (OH) group and a carbonyl group in close proximity. As a result, the OH and CO groups act, respectively, as a hydrogen bond donor and acceptor – thereby instigating an intramolecular hydrogen bond. In this configuration, the TQ and BT moieties contain non-planar geometries in which the cyclic ring is deformed at the sulphur atom due to monosaturation. In addition, both TQ and BT are tilted with respect to each other in order to reduce the steric hindrance that would otherwise arise at a planar configuration. As a result, in order to maintain the strong intramolecular OH–OC hydrogen bond, the OH moiety is orientated out of the plane of the TQ unit.

2.2 Vertical excitation energies and oscillator strengths

Eumelanin. Table 1 lists calculated vertical excitation energies and accompanying oscillator strengths of transitions to the lowest lying singlet excited electronic states of structural motifs A–G. A select sub-set of associated orbitals and low energy orbital promotions for different polymeric structures are shown in Fig. 2 for A, F and G. Orbitals and orbital promotions associated with electronic excitations of the other structural forms of EU are displayed in Fig. S2.1 of the ESI.†

Table 1 Vertical excitation energies (in eV) and accompanying oscillator strengths (in parentheses) to the lowest lying singlet excited electronic states of dimeric, trimeric and cyclic tetramer forms of EU. A, B, C, D, E, F and PM were computed at the ADC(2)/cc-pVDZ level of theory whereas G was computed at the TD-DFT/CAM-B3LYP/cc-pVDZ level of theory

Electronic transition	Vertical excitation energy (oscillator strengths)
	A	B	C	D
Eumelanin dimer
S1–S0	1.78 (0.0537)	1.78 (0.0436)	1.34 (0.0163)	1.71 (0.0819)
S2–S0	2.03 (0.0405)	2.05 (0.0477)	1.77 (0.0000)	1.78 (0.0005)
S3–S0	2.29 (0.0062)	2.36 (0.0089)	2.05 (0.0592)	1.96 (0.0175)
S4–S0	2.67 (0.0018)	2.78 (0.0040)	2.23 (0.0167)	2.23 (0.0759)

---

Electronic transition	Vertical excitation energy (oscillator strengths)
	E	F
Eumelanin trimer
S1–S0	1.79 (0.0620)	1.78 (0.0448)
S2–S0	2.04 (0.0470)	2.04 (0.0505)
S3–S0	2.21 (0.0089)	2.19 (0.0240)
S4–S0	2.60 (0.0039)	2.74 (0.0092)

---

Electronic transition	Vertical excitation energy (oscillator strengths)
	G
Eumelanin cyclic tetramer
S1–S0	1.22 (0.0056)
S2–S0	1.31 (0.0611)
S3–S0	1.78 (0.0413)
S4–S0	1.90 (0.0079)
S5–S0	2.09 (0.0001)
S6–S0	2.13 (0.1195)
S7–S0	2.26 (0.1067)
S8–S0	2.63 (0.0758)
Pheomelanin
S1–S0	3.57 (0.0799)
S2–S0	3.86 (0.1181)
S3–S0	4.17 (0.0075)
S4–S0	4.34 (0.0515)

---

Fig. 2 Orbitals and orbital promotions involved in forming the lowest four singlet excited singlet states of the (a) dimeric, (b) trimeric and (c) cyclic tetramer forms of EU computed at the (a and b) ADC(2)/cc-pVDZ and (c) TDDFT/CAM-B3LYP/cc-pVDZ levels of theory. For compactness, only structures A, F and G are displayed in this figure. The orbitals and orbital promotions of the other structures are displayed in Fig. S2.1 of the ESI.†

In all studied polymeric structures, the lowest lying excited states are of ¹ππ* character – involving electron promotion from ring-centred bonding π to anti-bonding π* orbitals. In dimeric, trimeric and tetrameric units, the vertical excitations involve a mixture of states that are of either (i) locally excited character – i.e. there is no net change in orbital localisation or (ii) charge-transfer (CT) character – i.e. that the π and π* orbitals are localised on differing indole rings, leading to charge separation (vide infra).

In the presently studied dimers and trimers, the first four singlet excited states absorb in the near-IR and visible – the lowest of which involves a π* ← π orbital promotion across two different indole moieties – leading to the formation of a charge-separated state. As a result, the S₁ state contains charge transfer (CT) character. These excitation energies agree well with those computed previously.^41,42

The lowest three singlet excited states of the cyclic tetramer G are also of ¹ππ* character. Due to the highly extended conjugation of G, the lowest three states absorb in the near-infrared, but higher states of ¹ππ* character (i.e. ≥S₄), absorb in the visible.

We note however that oligomer chains that solely contain fully hydrogenated dihydroxyindoles absorb in the near-UV whereas partially (i.e. hydroxylindolone) or fully dehydrated (i.e. indole-5,6-dione) oligomers are known to contribute the near-IR and visible – as outlined recently by Prampolini et al.⁴² In the present study we focus on oligomers that contain mixtures of dihydroxyindole, hydroxyindolone and indole-5,6-dione moieties – but note that the presence of even a single dehydrated form (i.e. hydroxylindolone or indol-5,6-dione), alongside 5,6-dihydroxyindole in the oligomer chain, leads to a large bathochromic shift in the vertical excitation (cf. bare 5,6-dihydroxyindole^25,42).

Pheomelanin. Table 1 also lists calculated vertical excitation energies and accompanying oscillator strengths of transitions to the lowest lying singlet excited electronic states of the presently studied PM sub-unit (henceforth simply referred to as PM). The associated orbitals and orbital promotions are shown in Fig. 3. As depicted, the lowest two excited electronic states of PM absorb in the near-UV and involve π* ← π electron promotions, whilst the third and fourth excited states involve a mixture of π* ← n and π* ← π electron promotions. All transitions exhibit appreciable oscillator strengths. As with the polymeric structures of EU, the S₁ state develops increasing CT character upon R_O–H bond extension, thus changing electronic character upon excited state intramolecular hydrogen transfer (vide infra).

Fig. 3 Hartree–Fock orbitals and orbital promotions involved in forming the lowest four singlet excited states of PM computed at the ADC(2)/cc-pVDZ level of theory.

As Fig. 4 indicates, the calculated vertical excitation energies for structures of EU and PM and the associated profile of the continuous spectrum for EU agree moderately well with the experimentally observed absorption spectra of synthetic EU and PM³⁷ – potentially protecting against harmful UV-A and UV-B radiation from the sun that can damage DNA.

Fig. 4 Experimentally measured absorption spectra of eumelanin (solid black curve) and pheomelanin (dashed black curve) – reproduced from ref. 39. The calculated vertical excitation energies of the various structures are given as stick spectra. The vertical left axis accompanies the experimentally measured spectrum whereas the vertical right axis accompanies the calculated vertical excitation energies. Though absorptions of the EU trimers exist at wavelengths shorter than 350 nm, for simplicity we have not included them since they appear at very similar peak position to those of the EU dimers. The grey profile shows the continuous absorption spectrum derived from the corresponding vertical excitation energies of the dimeric and trimeric forms of EU (see Computational methodology for details and Fig. S2.2 of the ESI†).

2.3 A static ab initio picture of the potential ultrafast non-radiative decay paths

Eumelamin. The intrinsic excited state relaxation process, in selected polymeric forms of EU, was studied in the context of intramolecular proton transfer – considering the dimeric (A, B and D) and trimeric (E and F) oligomers. Such a proton transfer photo-reaction (also referred to as electron-driven proton transfer) could be viewed as an enol(imino) → keto tautomerism – driven by electronic excitation.⁴³Fig. 5 presents the potential energy profiles along the excited state intramolecular proton transfer coordinate (i.e. the R_O–H or R_N–H driving coordinate) for structures A, B and D. In each case, the filled black circles represent the relaxed profile of the ground electronic state along R_O–H or R_N–H. The open red, blue, green and pink circles represent, respectively, the energies of the S₁, S₂, S₃ and S₄ states, calculated along the corresponding ground state MEP along R_O–H or R_N–H and thus represent the unrelaxed vertical energies of the excited states along R_O–H or R_N–H. The filled red circles represent the relaxed MEP in the S₁ state along the R_O–H or R_N–H driving coordinate, whilst the open black circles represent the corresponding ground state energies computed using the returned geometries of the S₁ MEP along R_O–H or R_N–H.

Fig. 5 ADC(2)/cc-pVDZ potential energy profiles of the S₀ (black), S₁ (red), S₂ (blue), S₃ (pink) and S₄ (green) states along the R_O–H (or R_N–H) coordinate – representing intramolecular H transfer in structures (a) A, (b) B, and (c) D.

To test the validity of the ADC(2) PE profiles in Fig. 5, benchmark CASPT2 calculations were performed for A at selected R_O–H values along the relaxed ground state profile (see Fig. S3.1 of the ESI†). These benchmark CASPT2 calculations show (i) that the CASPT2 S₁ vertical excitation energy at equilibrium R_O–H is in excellent agreement with that derived using ADC(2) (ADC(2) = 1.78 eV, CASPT2 = 1.73 eV) and (ii) that the unrelaxed CASPT2 PE profile of S₁ mirrors that of ADC(2). We are therefore confident that the present ADC(2) calculations appropriately describe the PE profiles of the present systems.

Focusing on Fig. 5(a) and (b) first, structures A and B return essentially identical profiles along R_O–H – since the only molecular distinction between A and B is OH rotational isomerism. The relaxed potential energy profile of the ground state is endoergic with respect to electron-driven proton transfer – thus increases in potential energy upon increasing R_O–H. This indicates a lack of driving force for intramolecular proton transfer in the stable ground state electronic configuration. In contrast, the unrelaxed and relaxed S₁ paths are exoergic with respect to proton transfer and decrease in energy upon increasing R_O–H. We note that the stabilization of the S₁ state, along the unrelaxed R_O–H reaction path, is unique to the present systems but can be understood by inspecting the orbitals (in Fig. 2(a)) that constituent the electronic configuration of the S₁ state along the unrelaxed and relaxed paths. In both cases, the electronic character of the S₁ state at R_O–H = 1 Å is described as a ππ* electronic configuration, in which the participating π and π* orbitals are localized on different indole rings. This creates a charge separation in the S₁ state – thus constituting an electronic configuration that contains CT character. Proton transfer, in an electron/proton coupled reaction, counters the charge separation caused by the initial excitation step – thus imparting a favourable reaction along R_O–H elongation. The respective decrease and increase in energy of the S₁ and S₀ states leads to an inevitable curve crossing which may constituent a low energy CI – motion through which could lead to ultrafast internal conversion of the excited state population. Along the unrelaxed S₁ path, this crossing occurs at R_O–H = 1.5 Å, the position of which shifts to shorter R_O–H (R_O–H = 1.1 Å) along the relaxed S₁ path – thus encouraging internal conversion closer to the FC region.

Double enol → keto tautomerisation and subsequent E → Z isomerization of C, forms dimer D – which, in contrast to A, B and C, comprises an intramolecular hydrogen bond with an acidic N–H donor rather than an O–H donor. The potential energy profile, along the R_N–H driving coordinate of D is depicted in Fig. 5(c). The relaxed profile of the ground electronic state is, as before, endoergic with respect to intramolecular proton transfer – thus increasing in potential energy upon R_N–H bond elongation. The profiles of the unrelaxed S₁, S₂, S₃ and S₄ states, mirror that of the ground state, indicating no driving force for proton transfer in the vertical electronic configuration. This is understood by considering the orbitals involved in the vertical excitation – i.e. that they involve local π and π* excitations that are localized on the same indole moiety and do not lead to a net charge separation (see Fig. S2.1 of the ESI†). In contrast, the relaxed S₁ state at R_O–H ∼ 1.0 Å contains CT character, which arises from an electron promotion between sets of π and π* orbitals localized on different indole moieties – thus leading to a charge separated state (see orbitals in Fig. S2.3, ESI†). As a result the S₁ relaxed profile is exoergic with respect to intramolecular proton transfer and, as such, decreases in potential energy upon extending R_N–H. Proton transfer, via R_N–H bond elongation leads to conservation of the neutral configuration of the molecule. A curve crossing at R_O–H ∼ 1.06 Å is returned along the relaxed S₁ path en route to H transfer – indicating, as before, a possible ultrafast internal conversion route close to the FC region.

Possible trimeric units within EU are exemplified by trimers E and F – which are distinguishable only by syn–anti rotational isomerism of O–H. As such, the potential energy profiles along the equivalent R_O–H driving coordinate for E and F (Fig. 6(a) and (b), respectively) are very similar to each other and to those obtained for dimers A and B. One caveat with dimer F is the possibility of double intramolecular proton transfer – due to the availability of two labile hydrogen bonds. The two distinct R_O–H driving coordinates in trimer F are henceforth referred to as R_O(1)–H and R_O(2)–H – as indicated by the molecular depictions in Fig. 6(b) and (c). As shown, neither the unrelaxed nor relaxed S₁ profiles cross the S₀ state along the R_O(2)–H coordinate. This can be understood by considering the orbitals involved in preparing the S₁ state (i.e. π(1) → π*(1) in Fig. 2), which indicates an electron promotion from the central indole moiety to the neighbouring indole moiety, i.e. that the net movement of charge is in the opposite direction to that of proton-transfer along R_O(2)–H. Proton transfer, along R_O(2)–H would seek to further destabilize the central indole moiety and leave a net negative charge on the neighbouring indole ring on which the π*(1) acceptor orbital is localised. This lack of driving force for proton-transfer along R_O(2)–H can be seen more clearly in the two-dimension potential energy surface in Fig. S3.2 of the ESI.†

Fig. 6 Potential energy profiles of the S₀ (black), S₁ (red), S₂ (blue), S₃ (pink) and S₄ (green) states along the R_O–H coordinate – representing intramolecular proton transfer in structures E (a) and F ((b) & (c)).

Pheomelanin. PM also has the possibility to undergo photoinduced proton transfer. In much the same way as EU, this mechanism requires that, upon photoexcitation of the enol tautomer, proton migration from the donor to acceptor O atoms, leads to the formation of a keto-tautomer. Fig. 7(a) depicts the potential energy profiles of the lowest five singlet states of PM along the R_O–H coordinate (indicated by the purple arrow above the molecular structure in Fig. 7(a)). The relaxed potential energy profile of the ground electronic state was obtained, as before, by holding R_O–H fixed at various values along the path, and allowing all other internal degrees of freedom to relax to their respective minima. This procedure was then repeated on S₁ in order to obtain the S₁ relaxed potential energy profile across the range 1.0 Å ≥ R_O–H ≥ 1.2 Å. At R_O–H ≥ 1.3 Å, fixing only the R_O–H internal coordinate on S₁ returned a negative gradient along an orthogonal ring-opening C–S bond stretch coordinate – representing an alternative region of configuration space with a lower energy reaction path (vide infra). Therefore, in order to fully describe motion along the intrinsic proton transfer coordinate, the remainder of the relaxed S₁ relaxed profile (i.e. R_O–H ≥ 1.3 Å) along R_O–H required the ring centred C–S bond stretch (R_C–S) coordinate to be held fixed at its equilibrium value.

Fig. 7 (a) Potential energy profiles of the S₀, S₁, S₂, S₃ and S₄ electronic states of PM along the R_O–H driving coordinate – which is indicated by the purple arrow above the molecular structure. The filled black and red and open black, red, blue, pink and green circles are as described before. (b) Potential energy profiles of the S₀, S₁, S₂, S₃ and S₄ electronic states along the coupled R_O–H and R_C–S driving coordinates (as indicated by the purple arrows above the molecular structure) – separated by the left and right panels. The left panel displays a reprise of the potential energy profiles between 1.0 Å ≥ R_O–H ≥ 1.3 Å, displayed in Fig. 7(a). The right panel indicated the potential energy profile along the R_C–S coordinate, whilst holding R_O–H fixed at 1.3 Å.

Initially focussing on the intramolecular proton transfer depicted in Fig. 7(a), the ground state relaxed profile (filled black circles), and the corresponding vertical excited states energies along this path (open red, blue, pink and green circles), all show no significant driving force for intramolecular proton transfer. In contrast the relaxed profile of the S₁ state (filled red circles) shows a barrier of ∼0.3 eV at R_O–H = 1.2 Å, after which the S₁ relaxed profile is reactive with respect to proton transfer – returning an S₁/S₀ crossing at R_O–H = 1.8 Å.

Following photo-excitation, an additional path by which internal conversion could occur is via motion along a ring-opening coordinate. The right hand panel of Fig. 7(b) presents the potential energy profiles of the ground and lowest four excited singlet states of PM along the intramolecular ring-centred C–S bond extension coordinate (R_C–S, as indicated in the molecular inset in Fig. 7(b)). Again, the ground electronic state configuration contains no driving force for C–S bond extension and thus rises in potential energy, whilst the relaxed and unrelaxed potential energy profiles of the S₁ state are labile with respect to C–S bond extension of ring-opening. As a result, a long range S₁/S₀ crossing occurs at R_C–S = 2.8 Å and R_C–S = 3.2 Å, along the S₁ relaxed and unrelaxed paths, respectively. This crossing could represent another example of conical intersection, through which ultrafast internal conversion could occur. As evident in Fig. 7(b), the MEP by which this long range CI is reached occurs via coupled O–H and C–S bond extension motions – the former of which returns a small activation barrier as described above. To understand the topography of the PES associated with the O–H and C–S coupled stretching motions, we have scanned the two-dimensional unrelaxed PES along these two driving coordinates (depicted in Fig. 8; see methodology for details). This PES shows a local minimum at the FC region (at R_O–H = 1.0 Å and R_C–S = 1.88 Å) and an unrelaxed MEP that initially involves O–H bond elongation to ∼1.3 Å – along which the energy rises by ∼0.3 eV – followed by favourable motion along the C–S bond extension coordinate, which is exoergic by ∼1.5 eV (purple arrow indicates the MEP in Fig. 8). Since the excited state in Fig. 8 is unrelaxed, we note that the ∼0.3 eV barrier represents an upper limit to the true MEP. At the S₁/S₀ crossing at R_C–S = 3.2 Å, a seam of intersection is created along 1.0 Å ≥ R_O–H ≥ 1.35 Å. To the best of our knowledge, this so-called proton-transfer coupled ring-opening reaction has not been identified hitherto and could be an interesting mechanism that drives the photochemistry of related thiol oligomers.

Fig. 8 Two-dimensional potential energy surface of the ground and first excited electronic states of PM along the R_O–H and R_C–S driving coordinates.

2.4 Dynamics of melanin oligomers

Non-adiabatic surface-hopping molecular dynamics (MD) simulations were computed ‘on-the-fly’ on structure A of EU and for PM. Attempts at undertaking the MD simulations using the ADC(2)/cc-pVDZ level of theory were significantly limited by the large molecular sizes and number of basis functions of the present systems (especially for PM). As such the B3LYP functional of TD-DFT (see Computational methodology) was used in order to compute the energies and gradients required for the ensuing MD simulations. Though we recognize that TDDFT and Kohn–Sham DFT are limited near intersections with the ground state, nonadiabatic couplings between S₁/S₀ and thus the hopping probability were nevertheless computed at each time step. In the present case it was appropriate to do so as our primary aim is to assess the degree to which a particular melanin chromophore affords photostability following IC to S₀. As such the present approach is at least appropriate for assessing the approximate timescales for IC as well as elucidating the fate of the molecule on S₀ following IC.

For benchmarking purposes, the vertical excitation energies derived by TDDFT were compared with those obtained from ADC(2)/cc-pVDZ (see Table S1 of the ESI†). These benchmark calculations show that whilst the ADC(2) and TD-DFT/B3LYP/6-31G vertical excitation energies of PM are broadly in agreement, those for structure A derived by TD-DFT/B3LYP/6-31G underestimate the ADC(2) value. For comparative purposes the ab initio MD simulations for A were also repeated using the TDDFT/B3LYP/cc-pVDZ level of theory, the results of which show only modest differences from those obtained from TDDFT/B3LYP/6-31G (vide infra). We also note that at all levels of theory, the S₁ state contains charge-transfer character with comparable oscillator strengths; therefore the TDDFT/B3LYP/6-31G level of theory is a good compromise between accuracy and computational expense – which is essential for systems of the present molecular size. The TDDFT/B3LYP/6-31G level of theory was thus used in order to obtain the energies and gradients for the subsequent molecular dynamics (MD) simulations which treat the nuclei within the framework of Newton's classical equations of motion whilst treating the electrons quantum mechanically. Given the classical treatment of the nuclei, they cannot return quantum behavior such as proton tunneling – which may be important when describing proton transfer reactions. That said, the present MD simulations are powerful since they include the effects of non-adiabaticity between PESs and describe the dynamics in full dimensionality – a property that is clearly desirable in such large systems where full quantum dynamics simulations would fail. In order to obtain the relevance of the photochemical mechanism described in Section 2.3 and the associated timescales on which they occur, an ensemble of 100 trajectories for both dimer A and PM were computed – each starting from the S₁ state (see Computational methodology for more details).

Fig. 9 presents the average population in the S₁ state for (a) dimer A of EU and (b) PM as a function of propagation time – representing the decay of the S₁ state. For both A and PM, the lifetime of the S₁ state was obtained by fitting the decay to the function:

which was shown to work well by Ruckenbauer et al. and to perform well for related systems.^13,44 In eqn (1), t₁ is a time offset that represents the time taken for the first trajectory to decay to S₀ and t₂ is the decay constant. The overall lifetime (τ) of the S₁ state is given by the sum t₁ + t₂ which returns 15 fs for A and 59 fs for PM. As an additional benchmark, the ab initio MD simulations were repeated for A using the TD-DFT/B3LYP/cc-pVDZ level of theory. The subsequent decay profile is shown in Fig. S4 of the ESI,† which when fitted to eqn (1), returns a lifetime of 16 fs – revealing that the 6-31G and cc-pVDZ basis sets show no significant discrepancies to the obtained dynamics of A. Whilst we note that 100 trajectories may not yield quantitatively converged statistics, the presently selected ab initio methods and molecular structures limit the feasibility to compute many more. Nevertheless, the smooth variation of the S₁ decay profile as a function of time indicates that the statistics are at least qualitatively converged.

Fig. 9 Fraction of trajectories in the S₁ electronic state as a function of time. The continuous curves represent the decay profiles of (a) dimers A of EU and (b) PM. The corresponding dashed curves through the continuous profiles represent the fits to the decay profile using eqn (1).

For dimer A of EU, all trajectories initiated on the S₁ state maintained the ππ*(CT) electronic configuration shown in Fig. 2 and decayed to S₀ within 50 fs – all of which deactivated via the intramolecular proton transfer path detailed in Section 2.3. Such a short lifetime for A is consistent with the calculated potential energy profiles displayed in Fig. 5(a) – i.e. that an S₁/S₀ CI occurs close to the FC region, supporting previous conclusions that excited state proton/hydrogen detachment and transfer reactions generally occur faster than the pump–probe time-response of many time-resolved experiments.^35,45,46 This fast decay can be observed more clearly in Fig. 10 which displays the variation in the energy profile as a function of time of a representative trajectory of dimer A. The profile shows an initial S₁ energy of ∼2.0 eV which displays a net decline at t > 0 fs, whilst that of the S₀ state steadily rises. At short propagation times (0 fs > t > 20 fs), large amplitude motions of O–H stretch are observed, indicating a strong driving force for R_O–H extension. This first attempt at proton transfer occurs at t = 8 fs where ΔE(S₁ − S₀) ∼ 0.5 eV – representing an energy difference that is too large for a sufficient transmission probability to encourage surface-hopping. Upon a second attempt and at R_O–H ∼ 1.15 Å, the population in the S₁ state hops to the S₀ state at t = 20 fs – signifying internal conversion to the S₀ state. This agrees very well with the returned geometry at the approximate S₁/S₀ crossing point along the calculated ADC(2) profile in Fig. 5(a) –highlighting that the energies derived by TDDFT are appropriate for describing the photodynamics of A.

Fig. 10 Energy profiles of a representative trajectory of dimer A as a function of time. The black and red curves represent the profiles of the S₀ and S₁ states, respectively. From left to right, the molecular depictions represent the geometry at t = 0 fs, t = 20 fs (time at which S₁ → S₀ hopping occurs), t = 40 fs and t = 90 fs. The total energy was conserved at ∼5.02 throughout the entire propagation time.

Following S₁ → S₀ surface-hopping at t = 20 fs, the molecular structure is restored to that of the ground state geometry within 20 fs of propagation on the adiabatic S₀ potential energy surface – with no further attempt of proton transfer. This representative trajectory is a strong indication that dimer A is photostable – driving the thesis that following photoexcitation, rapid internal conversion leads to swift reformation of the starting ground state molecule, without detriment to the molecular structure.

For PM, an ensemble of 100 trajectories was initiated on the S₁ state. The lifetime of the S₁ state, obtained by fitting the decay curve displayed in Fig. 9(b) to eqn (1), returns τ = t₁ + t₂ = ∼59 fs – which is about half the value derived experimentally and theoretically for C–S ring opening in bare thiophene in the gas phase.^47,48 On this timescale, the majority of the evolving trajectories reached the S₁/S₀ crossing via a variety of R_O–H and R_C–S stretch motions. Out of the 100 computed trajectories, 42 deactivated via the mechanism outlined in Fig. 7(b), 16 via the mechanism outlined in Fig. 7(a), whilst 25 involved the C–S ring opening of the other heterocyclic ring (i.e. the terminal TQ moiety – henceforth referred to as R_C–S′). These are the three dominant relaxation pathways identified and will henceforth be referred to as pathway 1, pathway 2 and pathway 3, respectively. The remaining 17 trajectories continued on the excited state potential energy surface up to the maximum propagation time of 100 fs. Example trajectories associated with pathways 1, 2 and 3 are depicted in Fig. 11(a)–(c).

Fig. 11 Energy profiles of representative trajectories of PM as a function of time for the three major S₁ state deactivation pathways. The black and red curves represent the profiles of the S₀ and S₁ states, respectively. The open red squares represent the population at a given time. The total energy in each case was conserved at ∼8.80 eV.

Fig. 11(a) shows that the dominant deactivation via pathway 1 involves coupled R_O–H and R_C–S stretch motions – in a proton-transfer coupled ring-opening reaction. Within the initial 20 fs, trajectories evolving along pathway 1 show large amplitude motions in R_O–H stretch, oscillating in the range 1.0 Å > R_O–H > 1.3 Å. The upper limit of the R_O–H range is consistent with the R_O–H value required to circumvent the barrier to R_C–S ring-opening along the MEP reported in Fig. 7(b). This R_O–H oscillatory motion triggers rapid R_C–S elongation (i.e. ring-opening) in the time range ∼20 fs > t > ∼100 fs. In the case of the exemplar trajectory in Fig. 11(a), the population hops onto the ground electronic state PES at t = 36.5 fs and at R_C–S ∼ 3.0 Å, which is in excellent agreement with the calculated ADC(2) potential energy profile in Fig. 7(b). Following S₁ → S₀ hopping, the molecule continues on S₀ eventually reforming the starting PM molecule by 200 fs – via C–S ring-closure. Ring closure is driven by the steep gradient associated with shortening of R_C–S on the ground state PES – since the ring-opened bi-radical introduces large steric hindrances and loss of conjugation. On entropic grounds, ring closure is likely to be driven by the dampened local modes of C–S stretch due to the large π-conjugation. Trajectories show that the regaining of conjugation and restoration of PM leads to large vibrational excitation predominantly in O–H, C–C and C–H stretching modes. Such large vibrational excitations are expected to cool in a complex solvent environment – which is of great interest but beyond the scope of the present study. Given these findings, we therefore expect that deactivation of photoexcited PM via the coupled R_O–H and R_C–S path will likely result in photostability.

Fig. 11(b) depicts the energy as a function of time for a representative trajectory evolving along the alternative decay path summarized by the potential energy profiles in Fig. 7(a) (i.e. predominant excited state proton transfer or pathway 2). As with dimer A of EU, trajectories evolving along pathway 2 deactivate to S₀ within the first 30 fs, representing the light H-atom motion. In the specific trajectory depicted in Fig. 11(b), S₁ → S₀ surface hopping occurs at 20 fs – after which the trajectory continues on the ground electronic state potential energy surface until the original molecular structure is restored. As such, deactivation via pathway 2 also shows photostability.

The third identified deactivation route of photoexcited PM is via R_C–S′ bond extension along pathway 3 – i.e. ring-opening of the TQ moiety. Trajectories evolving along this path undergo S₀ → S₁ surface hopping in the timescale range ∼30 fs > t > ∼100 fs – reflecting the heavy atom C and S motions required to access the S₁/S₀ crossing. The particular trajectory presented in Fig. 11(c) shows that S₁ → S₀ surface hopping occurs at t = 37 fs and at R_C–S ∼ 3.5 Å. In contrast to pathways 1 and 2, S₁ → S₀ surface hopping via pathway 3 leads to continued R_C–S′ bond extension and eventual ring-opening (see right hand molecular depiction at t = 200 fs in Fig. 11(c)) without reformation of the starting PM molecule in 200 fs of propagation time. As a result, internal conversion via R_C–S′ is predicted to be detrimental to the PM chromophore – showing no photostability. This is in stark contrast to pathway 1 in which S₁ → S₀ surface hopping, mediated via R_C–S bond extension, leads to prompt ring closure within 200 fs. There are many factors that provide plausible explanations for the differences in the fates of PM following deactivation via pathways 1 and 3. One such factor is that evolution along pathway 3 leads to a continuation in ring-opening on the ground state PES on entropic grounds, i.e. ballistic ring-opening will result in a highly vibrationally excited biradical disfavouring subsequent ring-closure. Unlike pathway 1, ring-opening via R_C–S′ does not lead to substantial molecular strain – thus resulting in a lack of entropic driving force for ring-closure.

Motivated by the outcomes of the MD simulations, additional potential energy profiles were constructed along R_C–S′ in order to assess the local topography along this coordinate (see Fig. 12). Upon vertical excitation to S₁, Fig. 12 shows an S₁ unrelaxed profile (i.e. the profile containing the open red circles) that contains a net decrease in potential energy upon R_C–S′ ring-opening – en route to which a small barrier of ∼0.03 eV exists. The S₁ relaxed profile (i.e. the profile containing the filled red circles) shows an S₁/S₀ crossing at R_C–S′ = 2.6 Å, motion through which leads to progressive ring-opening on S₀ – as indicated by the acquired MD simulation. As such we expect PM to be less photostable than EU – which is in good agreement with previous dermatological studies of melanins.³⁶

Fig. 12 Potential energy profiles of the lowest five singlet states of PM along the R_C–S′ ring-opening coordinate. The filled black and red and open black, red, blue, pink and green circles are as described before.

3. General discussion and conclusions

In the present study we have characterised the most probable nuclear motions that drive IC in two classes of molecule that contribute to natural sunscreens found on mammalian skin. We have employed ab initio electronic structure calculations and molecular dynamics simulations which are arguably the method of choice for providing detailed insight into such radiationless transitions in such systems.

We acknowledge that the present calculations do not account for possible intersystem crossings and relate solely to the isolated molecules. We are thus blind to additional complexities (e.g. interactions with proximal solvent molecules and epidermal ions or bimolecular reactions with molecules such as O₂) that would form the detailed parts of a biochemical study.

It is also noteworthy that weak hydrogen-bonds can be heavily affected by polar-protic environments. A notable example is that of catechol in methanol solution in which the classic OH–OMe intramolecular hydrogen-bond observed in the gas phase, is broken – resulting in an altered photochemistry in methanol solution (cf. gas phase). In contrast, strong intramolecular hydrogen-bonds (such as those of the present systems) are largely retained in polar-protic environments – often contributing the ensuing photochemistry. Notable examples include oxybenzone (an active ingredient in commercial sunscreens) in methanol,⁴⁹ the guanine–cytosine base-pairs in D₂O and chloroform,⁵⁰ and d(adenine–thymine)₈ and d(guanine–cytosine)₈ duplexes in D₂O.⁴⁶

Therefore, the current findings are nevertheless revealing and a good starting premise for more complex computational studies. As noted in the Introduction, the primary role of the conjugated organic molecules present in melanin constituents is to absorb solar radiation and to convert the excess electronic energy into vibrational energy (which is then dissipated as heat) without detriment to starting molecular structure. Thus the molecular sunscreen is required to provide photoprotection via ensuring structural photostability.

Focussing primarily on EU, the calculated vertical excitation energies, of the dimers (A, B, C and D), trimers (E and F) and tetramer (G), to the various ¹ππ* states accord well with the experimental absorption spectrum – placing their absorption maxima well within the solar spectrum. We also note that Experimental studies have confirmed the impressive photostability of synthetic eumelanin and those extracted from human skin cells.^31,51

The present electronic structure calculations reveal a barrierless and potentially ultrafast electron driven proton transfer pathway from the S₁(1¹ππ*) state to an S₁/S₀ crossing – providing a static but clear explanation for the impressive photostability of the polymeric forms of EU. Though the energies of these reactive CT states are in the near-IR, the large density of electronic states absorbing at somewhat higher energies (i.e. those depicted in Fig. 4) could principally couple to (the presently identified) low lying CT states that show ultrafast proton-transfer. As outlined in the Introduction, previous experimental studies in aqueous solution and electronic structure calculations have identified possible proton-transfer driven relaxation paths in monomeric structures of photoexcited eumelanin building blocks. Sobolewski et al.²⁵ and Huijser et al.³⁸ both identified intrinsic proton-transfer reaction paths in 5,6-dihydroxyindole and indole-2-carboxylic acid, respectively. The present work on EU builds on previous high profile investigations of EU by expanding a given EU building block into a series of known donor–acceptor oligomers. In so doing, we have shown the details of the excited state topography along possible proton-transfer coordinates which could potentially promote ultrafast relaxation in the oligomers of EU.

The accompanying MD simulations confirm the ultrafast decay of dimer A of EU – showing IC to occur in 20 fs in both cases with no detriment to the original molecule. Within the cellular medium, the demonstrable photostability of EU leads us to conclude that almost all of the evolving population on S₁ must be funnelled towards the parent S₀ structure, which is expected to vibrationally cool by energy transfer to the bulk cellular environment. The present electronic structure calculations and associated MD simulations reveal that the dominant molecular distortions that drive IC in structures A, B and D is via electron driven proton transfer – thus adding these molecules to the multitude of systems whose non-radiative decay is dominated by intramolecular electron driven proton transfer.^{4,5,13,14,52,53} Other than intramolecular proton transfer, no further low-energy and potentially competing decay paths were identified in the presently studied structures of EU.

The mechanisms by which the presently studied section through PM ensures photostability is less clear cut than EU. The experimentally reported absorption spectrum of synthetic PM exists well within the wavelength range required for solar light protection – as do the presently calculated excitation energies of the sensibly selected section through PM (Table 1 and Fig. 4).²⁴ As with the EU oligomer, static ab initio calculations suggest that excitation to the first excited (1¹ππ*) state of PM shows a favourable driving force with respect to intramolecular proton transfer. Dynamic simulations reinforce the earlier conclusion that the excited state decay of PM is ultrafast – suggesting that intramolecular proton transfer is a contributing decay path (i.e. pathway 2). In addition, the dynamics and electronic structure calculations show that deactivation via R_C–S stretch motions is also a strong contributor for the ultrafast deactivation of photoexcited PM. Pathways 1 and 3 involve ring-opening nuclear motions that mediate IC of the evolving S₁ population back to S₀. In the former case, ring-opening is coupled to large amplitude motions in O–H stretch, which ultimately mediates S₁ → S₀ IC. Following S₁ → S₀ IC via pathway 1, the resulting S₀ bi-radical is rapidly quenched and thus leads to structural reformation of the starting parent molecule via ring-closure within 200 fs. In contrast, pathway 3 involves isolated ring-opening of the TQ moiety which also mediates S₁ → S₀ IC. The subsequent S₀ bi-radical remains structurally active on the ground state and continues to ring-open – leading to twisting of the terminal TQ moiety and thus molecular denaturation. We therefore conclude that PM is less photostable than EU, reinforcing earlier conclusions that PM is a less effective biological sunscreen (cf. EU).

The present study is an addition to the small but fast growing literature detailing the intrinsic excited state photochemistry and photophysics of molecular systems present in biological and commercial sunscreens. Such studies can help to determine the extent to which a given molecular building block affords photostability and photoprotection in response to UV excitation (if at all). We can anticipate many more such studies in the next few years. Since most aromatic chromophores absorb strongly in the mid-UV (200–300 nm), an important challenge is to design molecular constituents in commercial sunscreens with absorption maxima in the range of 280–400 nm – so as to capture the UV-A and UV-B components of the solar spectrum. We and others⁴² have shown that dehydrogenation of 5,6-dihydroxyindole (to form hydroxylindolone and indol-5,6-dione) leads to a bathochromic shift in absorption maximum (cf. 5,6-dihydroxindole) – the onsets of which ranges from the near-IR to the visible. Polymers of EU containing hydroxylindolone, indole-5,6-dione and 5,6-dihydroxyindole retain this long wavelength absorption and contain low lying reactive states that show molecular photostability. We have also shown that ultrafast deactivation of PM could be driven by a proton-transfer coupled ring-opening reaction which has (to the best of our knowledge) been unexplored hitherto. Such a process could provide a plausible route to lowering the excited state reaction barriers associated with C–S ring-opening in aromatic thiols (e.g. in oligomeric thiophenes that act as molecular photoswitches).

The present findings help to better understand the intrinsic deactivation mechanisms of natural biological sunscreen. Such outcomes could draw synergies with future experimental studies on commercial sunscreens – leading to the design of better and more efficient alternatives in commercial sunscreen products.

4. Computational methodology

The ground state minimum energy geometries of selected low energy dimeric and trimeric forms of EU (depicted in Fig. 1(a)) were optimized at the MP2/cc-pVDZ level of theory.^54,55 At this same level of theory, only the simplest repeating unit of PM was optimized (see Fig. 1(b)). Using these ground state minimum energy geometries, the corresponding vertical excitation energies and oscillator strengths of the lowest four singlet excited states were calculated at the ADC(2)/cc-pVDZ level of theory.^56,57 A continuous absorption spectrum was constructed for EU by broadening the calculated vertical excitation energies of A–F with Gaussian functions at a fixed FWHM of 100 nm. Owing to the large molecular size, attempts at computing the ADC(2) vertical excitation energies of G proved challenging. As such, the ground state equilibrium geometry of the cyclic form (G) of EU was optimized at the DFT/-B3LYP/6-31G level of theory and the corresponding vertical excitation energy calculated at the TD-DFT/CAM-B3LYP/cc-pVDZ level of theory. Since the levels of theory at which the vertical excitation energies of structures A–F and G were calculated were different – G was excluded from the construction of the continuous absorption spectrum.

For selected structures, one-dimensional relaxed potential energy profiles along the S₀ and S₁ states were computed using, respectively, the MP2/cc-pVDZ and ADC(2)/cc-pVDZ levels of theory – using R_O–H, R_N–H or R_C–S (i.e. ring opening coordinate) as the driving coordinate. The corresponding ground state energy, along the relaxed S₁ path was computed at the MP2/cc-pVDZ level of theory.

Where appropriate, two-dimensional (2D) relaxed PES of the S₀ electronic state was computed using the MP2/cc-pVDZ level of theory and two independent driving coordinates (i.e. R_O–H and R_C–S for PM). The corresponding S₁ state vertical excitation energies, along the relaxed 2D path, were computed using the ADC(2)/cc-pVDZ; thus generating the corresponding S₁ surface at the 2D S₀ relaxed PES. The ADC(2) and (TD)DFT calculations were computed using the Turbomole⁵⁸ and Gaussian 09⁵⁹ computational packages, respectively.

For benchmarking purpose, the S₀ and S₁ energies of A, along the S₀ relaxed path, were also calculated using the complete active space with second-order perturbation theory (CASPT2),⁶⁰ based on a state-averaged complete active space self-consistent field reference wavefunction (SA2-CASSCF) – alongside a cc-pVDZ basis set. The active space comprised of 6 electrons distributed in 6 orbitals – consisting of the highest three occupied π and the lowest unoccupied π* orbitals. An imaginary level shift of 0.5 E_H was used in order to aid convergence and avoid intruder states. The CASPT2/CASSCF calculations were performed in Molpro 2010.1.⁶¹

Using Newton-X,^62,63 non-adiabatic surface hopping molecular dynamics simulations were performed on dimer A of EU and on PM using the molecular dynamics with quantum transitions surface hopping method devised by Hammes-Schiffer and Tully.⁶⁴ An ensemble of 100 trajectories for both A and PM were created in order to obtain a statistical average of the excited state fate of each system. The initial conditions were obtained by sampling a Wigner distribution based on ground state harmonic frequencies calculated at the DFT/B3LYP/cc-pVDZ level of theory. At each Wigner geometry the vertical excitation energies and oscillator strength were calculated at the TD-DFT/B3LYP/6-31G level of theory. For the nuclear motions, Newton's classical equations were integrated in steps of 0.5 fs using the velocity Verlet algorithm. The semi-classical electronic time-dependent Schrödinger equation was integrated using Butcher's fifth-order Runge–Kutta method in steps of 0.025 fs.⁶⁵ A standard decoherence correction of 0.1 E_H was used as proposed by Granucci and Persico.⁶⁶ In order to evaluate the hopping probability, nonadiabatic couplings were computed by finite differences as proposed in ref. 64 and implemented in ref. 67. All trajectories were initiated on the S₁ states upon which their associated gradients were computed ‘on-the-fly’ using the TD-DFT/B3LYP/6-31G level of theory in the Gaussian 09 computational package. Test ab initio MD simulations for A were also conducted at the TD-DFT/B3LYP/cc-pVDZ level of theory – but contained only 50 initial conditions, which were also sampled using a Wigner distribution of the ground state harmonic wavenumbers calculated at the DFT/B3LYP/cc-pVDZ level of theory. Trajectories were propagated in time with a step size of 0.5 fs until the S₁–S₀ energy gap (ΔE(S₁ − S₀)) fell below 0.10 eV or when the S₁ state population hopped onto S₀ (whichever occurred first). The decay profiles of the S₁ state in all cases therefore indicate a specific trajectory fulfilling the ΔE(S₁ − S₀) < 0.1 eV threshold. Trajectories that did not fulfil this termination criterion were propagated on S₁ for a maximum time of 100 fs. Trajectories on S₀ were propagated for a maximum time of 100 fs for A and 200 fs for PM. The lifetimes of A and PM were carefully benchmarked using different values of ΔE(S₁ − S₀) of 0.10 eV, 0.13 eV and 0.15 eV, the results of which showed only modest changes to the S₁ lifetime (see Table S2 of the ESI†). Benchmark calculations of the vertical excitation energies using (TD)DFT/B3LYP/6-31G were performed. These are displayed in Table S1 of the ESI† and show no major discrepancy to those calculated using ADC(2) theory – justifying the using of DFT with the B3LYP functional.

Acknowledgements

The authors would like to thank Michael Ashfold (Bristol), Wolfgang Domcke (TUM), Rachel Crespo-Otero (Queen Mary University of London) and Deniz Tuna (Max-Planck-Institut für Kohlenforschung) for fruitful discussions. BM acknowledges the University of Bristol for the award of a postgraduate scholarship. TNVK is grateful to TUM for the award of a postdoctoral fellowship and to the EPSRC (Programme Grant EP/L005913) for funding.

References

1. A. L. Sobolewski, W. Domcke and C. Hättig, Proc. Natl. Acad. Sci. U. S. A., 2005, 102, 17903–17906.
2. T. N. V. Karsili, B. Marchetti, M. N. R. Ashfold and W. Domcke, J. Phys. Chem. A, 2014, 118, 11999–12010.
3. S. Perun, A. L. Sobolewski and W. Domcke, Chem. Phys., 2005, 313, 107–112.
4. A. L. Sobolewski and W. Domcke, J. Phys. Chem. A, 2007, 111, 11725–11735.
5. A. L. Sobolewski and W. Domcke, Phys. Chem. Chem. Phys., 2006, 8, 3410–3417.
6. I. J. Palmer, I. N. Ragazos, F. Bernardi, M. Olivucci and M. A. Robb, J. Am. Chem. Soc., 1993, 115, 673–682.
7. S. Clifford, M. J. Bearpark, F. Bernardi, M. Olivucci, M. A. Robb and B. R. Smith, J. Am. Chem. Soc., 1996, 118, 7353–7360.
8. S. Yamazaki, A. L. Sobolewski and W. Domcke, Phys. Chem. Chem. Phys., 2011, 13, 1618–1628.
9. R. S. Minns, D. S. N. Parker, T. J. Penfold, G. A. Worth and H. H. Fielding, Phys. Chem. Chem. Phys., 2010, 12, 15607–15615.
10. M. Barbatti, A. J. A. Aquino, J. J. Szymczak, D. Nachtigallová, P. Hobza and H. Lischka, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 21453–21458.
11. C. T. Middleton, K. d. L. Harpe, C. Su, Y. K. Law, C. E. Crespo-Hernández and B. Kohler, Annu. Rev. Phys. Chem., 2009, 60, 217–239.
12. T. Gustavsson, R. Improta and D. Markovitsi, J. Phys. Chem. Lett., 2010, 1, 2025–2030.
13. D. Tuna, N. Došlić, M. Mališ, A. L. Sobolewski and W. Domcke, J. Phys. Chem. B, 2015, 119, 2112–2124.
14. D. Tuna, A. L. Sobolewski and W. Domcke, J. Phys. Chem. A, 2014, 118, 122–127.
15. D. Tuna, A. L. Sobolewski and W. Domcke, J. Phys. Chem. B, 2014, 118, 976–985.
16. M. Barbatti, B. Sellner, A. J. A. Aquino and H. Lischka, Handbook of Computational Chemistry, 2008, pp. 209–235.
17. W. Domcke, D. R. Yarkony and H. Koppel, Conical Intersections: Electronic Structure, Dynamics and Spectroscopy, World Scientific, Singapore, 2004.
18. W. Domcke, D. R. Yarkony and H. Köppel, Conical Intersections: Theory, Computation and Experiment, Wiley, 2011, vol. 17.
19. N. Tarras-Wahlberg, G. Stenhagen, O. Larko, A. Rosen, A. M. Wennberg and O. Wennerstrom, J. Invest. Dermatol., 1999, 113, 547–553.
20. S. Meng and E. Kaxiras, Biophys. J., 2008, 94, 2095–2105.
21. E. Kaxiras, A. Tsolakidis, G. Zonios and S. Meng, Phys. Rev. Lett., 2006, 97, 218102.
22. C.-T. Chen, C. Chuang, J. Cao, V. Ball, D. Ruch and M. J. Buehler, Nat. Commun., 2014, 5, 3859.
23. J. I. N. Cheng, S. C. Moss, M. Eisner and P. Zschack, Pigm. Cell Res., 1994, 7, 255–262.
24. J. D. Simon, Acc. Chem. Res., 2000, 33, 307–313.
25. A. L. Sobolewski and W. Domcke, ChemPhysChem, 2007, 8, 756–762.
26. Y. V. Il'ichev and J. D. Simon, J. Phys. Chem. B, 2003, 107, 7162–7171.
27. J. B. Nofsinger, Y. Liu and J. D. Simon, Free Radicals Biol. Med., 2002, 32, 720–730.
28. J. B. Nofsinger, S. E. Forest and J. D. Simon, J. Phys. Chem. B, 1999, 103, 11428–11432.
29. J. Brian Nofsinger and J. D. Simon, Photochem. Photobiol., 2001, 74, 31–37.
30. A. Corani, A. Pezzella, T. Pascher, T. Gustavsson, D. Markovitsi, A. Huijser, M. d'Ischia and V. Sundström, J. Phys. Chem. Lett., 2013, 4, 1383–1388.
31. A. Huijser, A. Pezzella and V. Sundstrom, Phys. Chem. Chem. Phys., 2011, 13, 9119–9127.
32. M. Gauden, A. Pezzella, L. Panzella, M. T. Neves-Petersen, E. Skovsen, S. B. Petersen, K. M. Mullen, A. Napolitano, M. d'Ischia and V. Sundström, J. Am. Chem. Soc., 2008, 130, 17038–17043.
33. M. Gauden, A. Pezzella, L. Panzella, A. Napolitano, M. d'Ischia and V. Sundström, J. Phys. Chem. B, 2009, 113, 12575–12580.
34. A. Huijser, A. Pezzella, J. K. Hannestad, L. Panzella, A. Napolitano, M. d'Ischia and V. Sundström, ChemPhysChem, 2010, 11, 2424–2431.
35. A. El Nahhas, T. Pascher, L. Leone, L. Panzella, A. Napolitano and V. Sundström, J. Phys. Chem. Lett., 2014, 5, 2094–2100.
36. T. Ye and J. D. Simon, J. Phys. Chem. B, 2003, 107, 11240–11244.
37. D. N. Peles and J. D. Simon, J. Phys. Chem. B, 2010, 114, 9677–9683.
38. A. Huijser, M. F. Rode, A. Corani, A. L. Sobolewski and V. Sundstrom, Phys. Chem. Chem. Phys., 2012, 14, 2078–2086.
39. S. Perun, A. L. Sobolewski and W. Domcke, J. Phys. Chem. A, 2006, 110, 9031–9038.
40. A. S. Chatterley, J. D. Young, D. Townsend, J. M. Zurek, M. J. Paterson, G. M. Roberts and V. G. Stavros, Phys. Chem. Chem. Phys., 2013, 15, 6879–6892.
41. A. A. R. Watt, J. P. Bothma and P. Meredith, Soft Matter, 2009, 5, 3754–3760.
42. G. Prampolini, I. Cacelli and A. Ferretti, RSC Adv., 2015, 5, 38513–38526.
43. S. Lochbrunner, K. Stock and E. Riedle, J. Mol. Struct., 2004, 700, 13–18.
44. M. Ruckenbauer, M. Barbatti, T. Müller and H. Lischka, J. Phys. Chem. A, 2013, 117, 2790–2799.
45. A. Iqbal and V. G. Stavros, Abstracts of Papers of the American Chemical Society, 2009, 238.
46. Y. Zhang, K. de La Harpe, A. A. Beckstead, R. Improta and B. Kohler, J. Am. Chem. Soc., 2015, 137, 7059–7062.
47. D. Fazzi, M. Barbatti and W. Thiel, Phys. Chem. Chem. Phys., 2015, 17, 7787–7799.
48. S. Salzmann, M. Kleinschmidt, J. Tatchen, R. Weinkauf and C. M. Marian, Phys. Chem. Chem. Phys., 2008, 10, 380–392.
49. L. A. Baker, M. D. Horbury, S. E. Greenough, P. M. Coulter, T. N. V. Karsili, G. M. Roberts, A. J. Orr-Ewing, M. N. R. Ashfold and V. G. Stavros, J. Phys. Chem. Lett., 2015, 6, 1363–1368.
50. K. Röttger, H. J. B. Marroux, M. P. Grubb, P. M. Coulter, H. Böhnke, A. S. Henderson, M. C. Galan, F. Temps, A. J. Orr-Ewing and G. M. Roberts, Angew. Chem., Int. Ed., 2015, 54, 14588.
51. S. Olsen, J. Riesz, I. Mahadevan, A. Coutts, J. P. Bothma, B. J. Powell, R. H. McKenzie, S. C. Smith and P. Meredith, J. Am. Chem. Soc., 2007, 129, 6672–6673.
52. A. L. Sobolewski, W. Domcke, C. Dedonder-Lardeux and C. Jouvet, Phys. Chem. Chem. Phys., 2002, 4, 1093–1100.
53. A. L. Sobolewski, W. Domcke and C. Hättig, J. Phys. Chem. A, 2006, 110, 6301–6306.
54. C. Møller and M. S. Plesset, Phys. Rev., 1934, 46, 618–622.
55. W. J. Hehre, R. F. Stewart and J. A. Pople, J. Chem. Phys., 1969, 4, 2657–2664.
56. A. Dreuw and M. Wormit, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2015, 5, 82–95.
57. T. H. Dunning, Jr., J. Chem. Phys., 1989, 90, 1007–1023.
58. TURBOMOLE V6.4 2012, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989–2007, TURBOMOLE GmbH, since 2007, available from http://www.turbomole.com.
59. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennuchi and G. Petersson, et al., Gaussian 09, revision B.01, Gaussian Inc., Wallingford, CT, 2010.
60. B. O. Roos and K. Andersson, Chem. Phys. Lett., 1995, 245, 215–223.
61. H. J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, M. Schütz, P. Celani, T. Korona, R. Lindh, A. Mitrushenkov, G. Rauhut, K. R. Shamasundar, T. B. Adler and R. D. Amos, et al., MOLPRO, University of Cardiff: Cardiff, UK, 2010.
62. M. Barbatti, G. Granucci, M. Persico, M. Ruckenbauer, M. Vazdar, M. Eckert-Maksić and H. Lischka, J. Photochem. Photobiol., A, 2007, 190, 228–240.
63. M. Barbatti, M. Ruckenbauer, F. Plasser, J. Pittner, G. Granucci, M. Persico and H. Lischka, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2014, 4, 26–33.
64. S. Hammes-Schiffer and J. C. Tully, J. Chem. Phys., 1994, 101, 4657–4667.
65. J. Butcher, J. Assoc. Comput. Mach., 1965, 12, 124.
66. G. Granucci and M. Persico, J. Chem. Phys., 2007, 126, 134114.
67. F. Plasser, R. Crespo-Otero, M. Pederzoli, J. Pittner, H. Lischka and M. Barbatti, J. Chem. Theory Comput., 2014, 10, 1395–1405.

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c5cp06767g

‡ Both the authors contributed equally to this publication.

---