Investigating photodynamics of nucleobase–water systems

Abstract

We examine the photodynamics of weakly bound hydrogen-bonded molecular systems, with a primary focus on developing a pragmatic computational protocol that incorporates flexibility and quantum distributions of initial conditions for nonadiabatic dynamics. As a case study, the photodynamics of 2,6-diaminopurine nucleobase-water clusters has been investigated, highlighting the active involvement of water in their photophysics and photochemistry.

---

Excited-state dynamics of nucleotide bases has remained a prominent research topic for several decades,^1,2 driven by both the biochemical relevance of these molecules and their implications on the origins of life.³ Although progress has been made in studying more complex systems – such as molecular aggregates and condensed-phase environments^4,5 – a great deal of theoretical research still focuses on isolated systems in vacuum. It is widely recognized that water can significantly influence the photodynamics of nucleobases^6–11 and their structurally similar analogues; however, in silico approaches for modelling nucleobase–water systems remain challenging, as reliable and widely accepted computational protocols beyond gas phase are generally lacking. The complexity arises already at the simulations outset, namely the preparation of initial conditions (ICs) for nonadiabatic dynamics.^12–14

In the context of mixed quantum–classical dynamics, ICs denote nuclear position-momentum pairs that (ideally) map the nuclear phase-space of the initially populated electronic state – typically the ground state. Within the sudden excitation approximation ICs are subsequently promoted vertically to the targeted excited state(s) where nonadiabatic trajectories are initiated. Ground-state ICs sampling is conveniently performed using the Wigner quasi-probability distribution¹⁵ in combination with harmonic frequencies at the ground-state minimum. Unlike thermal (Boltzmann) sampling, Wigner sampling mimics quantum phase-space distributions, thereby accounting for important effects such as zero-point energy (ZPE). It is well understood that, due to notable implications on calculated observables, the role of ZPE should not be neglected in nonadiabatic dynamics.¹⁶ This is because ZPE determines the minimal vibrational energy deposited in each vibrational mode (considering the molecule as a collection of harmonic oscillators), including the photoactive modes. However, the widespread use of harmonic Wigner sampling is difficult to justify for systems that are anharmonic, flexible, and characterized by numerous interconverting local minima. Nucleobase–water clusters are examples of such systems.

With the goal of studying weakly bound hydrogen-bonded molecular clusters, we pose the following question: how can the photodynamics of these flexible systems be effectively investigated? To address this challenge, here we advocate for a simple, application-oriented computational protocol consisting of four steps.

(A) Thermal sampling via molecular dynamics (MD) facilitates the exploration of the ground-state thermodynamics of flexible systems. Nucleobase–water clusters form complex hydrogen-bonding networks, rendering conventional geometry optimizations intractable even for slightly larger clusters. Moreover, calculating Boltzmann populations rely on the harmonic approximation, which is not justified for such systems. Therefore, we propose to explore conformational space by performing long-timescale MD simulations using molecular mechanics or a low-cost electronic structure, spanning time scales up to hundreds of nanoseconds. We note that exploration of the ground-state free-energy surface of nucleobases in water through classical potentials has shown reliable accuracy.¹⁷ Furthermore, data clustering of the most common hydrogen-bonding patterns allows us to identify those structures that occur the most frequently in thermal MD runs (i.e., extracting structures with thermodynamic relevance), thus, avoiding the bias introduced by arbitrary cluster design.

(B) Quantum thermostat (QT)^18,19 offers an MD-based alternative to Wigner sampling by enabling molecular modes to have different effective (i.e. frequency-dependent) temperatures, thereby recovering ZPE vibrational effects. Unlike Wigner sampling, QT is suited for anharmonic and flexible systems.²⁰ It typically produces nuclear distributions comparable to those obtained from more rigorous path-integral simulations,^21,22 while simultaneously providing both positions and momenta – which makes it suitable for sampling ICs.²² Using the most frequently occurring clusters from A, we propose conducting shorter QT runs (spanning time scales of tens of picoseconds) in combination with a more accurate electronic structure method, such as DFT. Ab initio MD simulations coupled with QT serve to “heal” the pre-sampled cluster structures by readjusting them to the higher-level electronic structure and simultaneously allowing vibrational modes to thermalize properly to mimic quantum distributions.

(C) Nuclear ensemble approach (NEA)²³ provides means to estimate a linear absorption spectrum by calculating electronic excitations from an ensemble of geometries representing the initial state nuclear distribution. We propose sampling this ensemble from QT dynamics (B) of the most common chromophore-water clusters selected in A, while preserving the relative weights of structures according to their occurrence in the data clustering (this can be seen as a substitute for Boltzmann weights). Furthermore, the NEA absorption spectrum serves to define an excitation energy window¹² from which nonadiabatic trajectories are launched.

(D) Trajectory surface hopping (TSH)²⁴ is finally used to propagate swarms of classical trajectories in the manifold of nonadiabatically coupled electronic states. Trajectories are initiated from the pool of ICs sampled via QT and selected from an excitation window superimposed to the absorption spectrum (C).

To illustrate this workflow, we investigate the photodynamics of 2,6-diaminopurine (26DAP), a non-canonical nucleobase that features rich excited-state behavior and dynamic equilibrium of structural isomers in solution. Despite its relevance in biology²⁵ and applications in crystal engineering,²⁶ 26DAP has received comparatively less attention than other nucleobases. In aqueous medium, 26DAP exists as a mixture of 9H and 7H tautomers in an estimated 60 : 40 ratio.²⁷ Up to date, the photophysics of two 26DAP tautomers have been studied through potential energy surface mapping in vacuo^27,28 and polarizable continuum.²⁹ Aside from some quantitative differences arising from variations in electronic structure, C6-puckering was identified as the primary deactivation pathway in both tautomers, while C2-puckering was deemed less likely. Furthermore, 7H was found to exhibit larger excited-state stability than 9H due to a higher excited-state barrier for reaching the S₁/S₀ conical intersection associated with puckering, suggesting its higher fluorescence²⁷ and triplet quantum yield.²⁸

Using our in silico protocol, we examined the singlet state photodynamics of two tautomers within water clusters containing 10 solvent molecules (denoted as 9H-H₂O and 7H-H₂O), as well as in vacuo for comparison. For solvated systems, we performed pre-sampling through thermal MD simulations using a previously customized force field,¹⁷ followed by structural clustering to extract the predominant hydrogen bonding patterns. ICs sampling was conducted for the most frequently occurring clusters, using QT at 298 K, in combination with the DFT electronic structure. Full methodological details are available in the SI.

The calculated absorption spectra of 7H-H₂O and 9H-H₂O (using NEA + TDDFT) were compared to the 26DAP spectrum that we measured in water (experimental details in SI). Spectra are depicted in Fig. 1. As the overall absorption represents an isomer mixture, we performed a spectral fitting and found that a 64 : 36 ratio of 9H-H₂O and 7H-H₂O best reproduces the low-energy absorption band of 26DAP in water. To account for systematic errors in the electronic structure calculations and the differences between clusters and bulk, we also allowed a rigid shift of the calculated spectra (see the SI for details). Based on this analysis, a narrow energy window for ICs selection was specified, corresponding to the experimental UV excitation energy of 4.67 eV (266 nm).²⁶

Fig. 1 Experimental, calculated and fitted absorption spectra of 26DAP.

The TSH dynamics, on one hand, confirms general trends anticipated from earlier works; on the other hand, it reveals several important observations regarding the impact of water on photodynamics. Fig. 2 shows the time evolution of the C6-puckering dihedral for trajectory ensembles of 7H-H₂O and 9H-H₂O, traced from the ICs until reaching the S₁/S₀ conical intersection (see Fig. S10 for pristine 7H and 9H). The C6-distortion serves as a primary nonadiabatic decay channel for both systems, while C2-distortion was detecable only for 9H/9H-H₂O (see Fig. S3 and S4). Despite the structural similarity of the two tautomers, 7H-H₂O exhibits significantly larger excited-state stability than 9H-H₂O. Although this trend was also claimed for isolated 9H and 7H,^27,28 water has a surprising effect: it increases excited-state stability in 7H, while decreasing it in 9H. The fractions of inactive trajectories that do not reach the vicinity of S₁/S₀ conical intersection amount to 51% in 9H-H₂O, 60% in 9H, 86% in 7H-H₂O and 65% in 7H (statistical error approx. 5%). While these variations originate from the combination of several decay channels (see Table 1), it is instructive to analyse C6-puckering alone. The average deviation from planarity calculated from the C6-distortion dihedrals of all trajectories increases in 9H-H₂O (39°) with respect to 9H (36°) and decreases in 7H-H₂O (23°) with respect to 7H (30°) – which is also reflected in C6-puckering fractions (however, a relatively small number of trajectories still limits the firmer statement).

Fig. 2 Time-evolution of C6-twisting dihedral for the TSH dynamics of 7H-H₂O (A) and 9H-H₂O (B). Structures in the vicinity of S₁/S₀ conical intersections are marked by red crosses. Incept shows a C6-puckered structure with the highlighted dihedral (C2–N1–C6–N).

Table 1 Outcomes of TSH simulations over the 2.5 ps time frame, using TDDFT electronic structure (validated against SCS-ADC(2), see SI)

	9H-H2O	7H-H2O	9H	7H
C2	4	—	2	—
C6	33 (40%)	11 (14%)	29 (33%)	14 (21%)
H diss.	—	—	4	9
PCET	1	—	—	—
CT	2	—	—	—
Other	—	—	—	1
Inactive	42 (51%)	69 (86%)	52 (60%)	44 (65%)
Total	82	80	87	68

---

The second observation is that water introduces new photochemical deactivation channels while suppressing others that may be relevant in the gas phase. Specifically, in 9H-H₂O, we identified charge transfer (CT) from water to the amino group (a mechanism first described in solvated adenine⁶) and proton-coupled electron transfer (PCET) from water to the purine nitrogen as potential deactivation pathways of solvated 26DAP (see Fig. S5 and S6). This highlights that water does not merely act as passive environment but also plays an active role in photodynamics. In contrast, water clusters completely inhibit hydrogen dissociation via N–H bond fission – a significant deactivation pathway for 7H and 9H in vacuo (see Fig. S7 and S8). Although earlier studies of 26DAP^27–29 did not consider hydrogen dissociation, its occurrence is unsurprising since it has been observed in adenine – it was shown that hydrogen dissociation proceeds via dissociative πσ* states, which, in an aqueous environment, undergo a significant blue shift relative to the onset of the same process in gas phase³⁰ – an observation which also corroborates our findings for 26DAP.

Finally, our computational protocol and the results merit some critical appraisal. While simulation quality can always be improved by refining the level of electronic structure theory or addressing the well-known concern that simulations are never sampled long enough, we highlight one specific issue that needs further investigation. Fig. 3 compares structural distributions from the ground-state sampling (QT) and excited-state (TSH) dynamics of 9H-H₂O, illustrating the substantial excited-state kinetic energy flow into the water subsystem. On one hand, this is expected as in the stable ground-state cluster the slow degrees of freedom of water are fully equilibrated with the electronic density of the nucleobase, which is not the case upon the UV-excitation. The ensuing excited-state dynamics drives the system out-of-equilibrium, with water molecules naturally re-adapting to the electronic density of excited state(s). On the other hand, classical trajectories initiated from a quantum distribution may suffer from the energy leakage problem, which indicates fast kinetic energy flow from high- to low-frequency modes, violating ZPE and ultimately driving the system toward a Boltzmann-like energy equipartition.¹⁶ This problem has been studied primarily in the ground-state dynamics – paradigmatic case being the unphysical water dimer dissociation occurring within a couple of picoseconds.³¹ In the context of QT-based dynamics, energy leakage is effectively diminished within a strong system-bath coupling regime¹⁹ (see also discussion in the SI). However, the extent of excited-state energy leakage may be an important point of interest as it may have its toll on photodynamics pathways (see also Fig. S12 for the example of water dissociation in 9H-H₂O cluster). Therefore, a potential topic for future investigation is how to distinguish (and correct) the spurious solute-to-solvent energy leakage from the genuine non-equilibrium energy redistribution in the nonadiabatic dynamics of molecular clusters.

Fig. 3 QT ground-state sampling (left) vs. microcanonical excited-state TSH dynamics (right) of 9H-H₂O, over the 2.5 ps time frame.

Overall, we proposed a set of computational “good practices” to investigate the photodynamics of flexible chromophore-water systems, accounting for structural disorder due to large amplitude motions as well as “quantum” ICs. With this, we addressed the known limitations of computational protocols relying on the harmonic approximation, as well as the shortcomings of conventional Boltzmann sampling which fails to account for ZPE.¹⁶ Our approach treats all degrees of freedom on equal footing, which seems more sensible than hybrid methods that separate quantum and classical subsystems (e.g. treated via Wigner and Boltzmann sampling, respectively), popularly used along with QM/MM.^14,16 We have demonstrated that water may influence the photophysics of a solvated nucleobase in a non-trivial way, while it can also be electronically involved in its photochemistry. The active role of water supports the notion that (at least) the first solvation shell warrants treatment as part of the quantum system, not only in terms of electronic structure but also from the perspective of ICs distribution. In other words, hydrogen-bonded water molecules should not be regarded merely as system's bath, but rather as a consistent part of the system itself. Finally, while addressing the gap related to photodynamics of solute–solvent clusters, we also anticipate further generalizations of the computational protocol extending from microsolvated systems to bulk solvent environments.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the SI.

Supplementary information available: Full computational and experimental details, supplementary analyses, and method assessment. See DOI: https://doi.org/10.1039/d5cp02344k

Acknowledgements

Research was supported by the Croatian Science Foundation under the project number HRZZ-IP-2020-02-7262. We thank Dr Tomislav Stolar for supplying the material.

Notes and references

1. R. Improta, F. Santoro and L. Blancafort, Chem. Rev., 2016, 116, 3540–3593.
2. 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.
3. C. L. Kufner, D. B. Bucher and D. D. Sasselov, ChemSystemsChem, 2023, 5, e202200019.
4. R. Improta and V. Barone, Photoinduced Phenomena in Nucleic Acids I: Nucleobases in the Gas Phase and in Solvents, Springer, 2014, pp. 329–357.
5. B. Milovanović, J. Novak, M. Etinski, W. Domcke and N. Došlić, Phys. Chem. Chem. Phys., 2022, 24, 14836–14845.
6. M. Barbatti, J. Am. Chem. Soc., 2014, 136, 10246–10249.
7. R. Szabla, R. W. Góra, M. Janicki and J. Šponer, Faraday Discuss., 2016, 195, 237–251.
8. X. Wu, T. N. V. Karsili and W. Domcke, Chem. Phys. Chem., 2016, 17, 1298–1304.
9. R. Szabla, H. Kruse, J. Šponer and R. W. Góra, Phys. Chem. Chem. Phys., 2017, 19, 17531–17537.
10. P. Chakraborty, T. N. V. Karsili, B. Marchetti and S. Matsika, Faraday Discuss., 2018, 207, 329–350.
11. M. J. Janicki, R. Szabla, J. Šponer and R. W. Góra, Phys. Chem. Chem. Phys., 2022, 24, 8217–8224.
12. J. Janoš, P. Slavíček and B. F. E. Curchod, Acc. Chem. Res., 2025, 58, 261–270.
13. E. R. Curtis, C. M. Jones and T. J. Martínez, J. Phys. Chem. B, 2025, 129, 2030–2042.
14. D. Avagliano, E. Lorini and L. González, Philos. Trans. R. Soc., A, 2022, 380, 20200381.
15. E. Wigner, Phys. Rev., 1932, 40, 749.
16. M. Barbatti and K. Sen, Int. J. Quantum Chem., 2016, 116, 762–771.
17. T. Ostojić, J. Ovčar, A. Hassanali and L. Grisanti, Phys. Chem. Chem. Phys., 2025, 27, 11869–11878.
18. M. Ceriotti, G. Bussi and M. Parrinello, Phys. Rev. Lett., 2009, 103, 030603.
19. M. Ceriotti, G. Bussi and M. Parrinello, J. Chem. Theory Comput., 2010, 6, 1170–1180.
20. A. Prlj, E. Marsili, L. Hutton, D. Hollas, D. Shchepanovska, D. R. Glowacki, P. Slavíček and B. F. E. Curchod, ACS Earth Space Chem., 2022, 6, 207–217.
21. J. Suchan, D. Hollas, B. F. E. Curchod and P. Slavíček, Faraday Discuss., 2018, 212, 307–330.
22. A. Prlj, D. Hollas and B. F. E. Curchod, J. Phys. Chem. A, 2023, 127, 7400–7409.
23. R. Crespo-Otero and M. Barbatti, Theor. Chem. Acc., 2012, 131, 1237.
24. J. C. Tully, J. Chem. Phys., 1990, 93, 1061–1071.
25. R. Szabla, M. Zdrowowicz, P. Spisz, N. J. Green, P. Stadlbauer, H. Kruse, J. Šponer and J. Rak, Nat. Commun., 2021, 12, 3018.
26. T. Stolar, B. K. D. Pearce, M. Etter, K.-N. Truong, T. Ostojić, A. Krajnc, G. Mali, B. Rossi, K. Molčanov and I. Lončarić, et al. , Iscience, 2024, 27, 109894.
27. L. M. F. Oliveira, D. Valverde, G. J. Costa and A. C. Borin, Photochem. Photobiol., 2024, 100, 323–338.
28. G. Gate, A. Williams, S. Boldissar, J. Šponer, R. Szabla and M. de Vries, Photochem. Photobiol., 2024, 100, 404–418.
29. N. E. Caldero-Rodríguez, L. A. Ortiz-Rodríguez, A. A. Gonzalez and C. E. Crespo-Hernández, Phys. Chem. Chem. Phys., 2022, 24, 4204–4211.
30. G. M. Roberts, H. J. B. Marroux, M. P. Grubb, M. N. R. Ashfold and A. J. Orr-Ewing, J. Phys. Chem. A, 2014, 118, 11211–11225.
31. S. Mukherjee and M. Barbatti, J. Chem. Theory Comput., 2022, 18, 4109–4116.

---

Footnote

† These authors contributed equally to this work.

---