Article pubs.acs.org/JPCB
Cite This: J. Phys. Chem. B 2019, 123, 5555−5567
Interdependent Electronic Structure, Protonation, and Solvatization of Aqueous 2‑Thiopyridone Jesper Norell,* Anton Ljungdahl, and Michael Odelius* Department of Physics, AlbaNova University Center, Stockholm University, SE-106 91 Stockholm, Sweden
Downloaded via KEAN UNIV on July 25, 2019 at 04:54:51 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.
S Supporting Information *
ABSTRACT: 2-Thiopyridone (2-TP), a common model system for excited-state proton transfer, has been simulated in aqueous solution with ab initio molecular dynamics. The interplay of electronic structure, protonation, and solvatization is investigated by comparison of three differently protonated molecular forms and between the lowest singlet and triplet electronic states. An interdependence clearly manifests in the mixed-character T1 state for the 2-TP form, systematic structural distortions of the 2-mercaptopyridine (2-MP) form, and photobase protolysis of the 2-TP− form, in the aqueous phase. In comparison, simplified continuum models for the solvatization are found to be significantly inaccurate for several of the species. To facilitate future computational studies, we therefore present a minimal representative solvatization complex for each stable form and electronic state. Our findings demonstrate the importance of explicit solvatization of the compound and sets the stage for including it also in future studies.
■
INTRODUCTION Excited-state proton transfer (ESPT) compounds, including photoacids and photobases, are both commonly occurring and functionally essential for many biochemical processes, motivating decades of studies on their properties and reactions.1−6 In recent years, particular interest has resurged in 2-thiopyridone (2-TP) as a compact model system for ESPT7−12 to enable detailed analysis with highly selective time-resolved X-ray spectroscopy7,9−12 and state-of-the-art quantum chemical methods,10,12 but the protonation states of other small solutes have also been studied with similar methods.13−16 In addition, 2-TP belongs to the class of heteroaromatic thione compounds, which due to their respective wide range of functionalities have also become popular targets for time-resolved spectroscopy and quantum chemical simulations.17−23 In a recent study,12 we established how the excited-state dynamics of 2-TP in aqueous solution leads to N site deprotonation in a pico- to nanosecond decay pathway through the T1 state, shared for both S2 and S4 photoexcitation. Indications were also provided that the same pathway could explain previous measurements in acetonitrile.11 Yet, the time constants appear to differ in the two solutions, which oppositely indicates a key role of the solvent in mediating the excited-state dynamics. This is in accordance with the general understanding of ESPT processes, where the solvent is well known to affect, e.g., rate constants,24 but its exact mechanistic role has both been disputed5,25−27 and remained elusive even to recent measurements.28 In fact, for 2TP, we also showed that none of the so-far reported timeresolved measurements8,10−12 appear to have been directly © 2019 American Chemical Society
sensitive to whether the proton dynamics occurs intra- or intermolecularly, i.e., whether the proton is directly transferred to the S site of the compound or extracted and subsequently redonated by the solvent. It thereby remains an open question whether 2-TP actually functions as an intramolecular-ESPT compound, as implicitly assumed in the first time-resolved study of its excited-state dynamics,8 or as a photoacid. Consequently, it is crucial to establish how the compound in its different molecular forms and electronic states interacts with the solvent to fully unravel all steps of its excited-state dynamics. In particular, such results could, in turn, further the general mechanistic understanding of proton-transfer processes through detailed insights into the interplay of nonadiabatic dynamics, solvent relaxation, and proton dynamics. The protonation states of key sites in 2-TP, related to tautomerization and protolysis, are controlled by several thermodynamic parameters and light-driven processes between three hitherto identified molecular forms as shown in Figure 1. While the thiol form 2-mercaptopyridine (2-MP) is favored in the gas phase,7,29,30 the thione tautomer 2-thiopyridone (2TP) instead dominates in the solid phase31 and most solvents,8−10,12,30−32 including water. A third related form, the deprotonated anion (2-TP−), can be stabilized in basic solutions9,10 and, as described above, is possibly formed as an intermediate in the ESPT process.10,12 In previous computational studies, the solvatization of the compound has either Received: April 2, 2019 Revised: May 31, 2019 Published: June 5, 2019 5555
DOI: 10.1021/acs.jpcb.9b03084 J. Phys. Chem. B 2019, 123, 5555−5567
Article
The Journal of Physical Chemistry B
performed at the DFT level and compared to single point electronic structure calculations based on the multiconfigurational complete active space (CAS) framework. Classical Molecular Dynamics. The classical MD simulations were performed with Gromacs (version 5.1.2).33 The aqueous solvent was described with the TIP4P/200534 model and the proton-transfer compound with optimized potential for liquid simulations (OPLS)/all-atom.35,36 The optimized geometries of 2-TP, 2-MP, and 2-TP−, as well as their partial charges were obtained at the HF/6-31G* level of theory in Gaussian 09.37 In addition to the proton-transfer compound, the simulation boxes contained in all cases 100 water molecules, and 1 Na+ counterion in the case of 2-TP−. Initial conditions were generated with standard Gromacs routines to obtain a density of 997.95 kg/m3 corresponding to ambient water within the TIP4P/2005 model. The simulations were performed in the NVT ensemble with a target temperature of 300 K. A time step of 1 fs was applied for 400 ps of simulation; the first 100 ps were discarded from the analysis to ensure equilibration. The topology files for the generated force fields and a representative simulation input file are provided in the Supporting Information for further reference. Ab Initio Molecular Dynamics. The AIMD simulations were performed with the density functional theory (DFT) as implemented in the CP2K software suite (version 5.0).38−41 The gradient-corrected BLYP functional was augmented with dispersion corrections using a 12 Å radial cutoff of the van der Waals interactions.42−45 In CP2K, the Gaussian and plane wave method41 allows one to efficiently describe the aqueous solution in periodic boundary conditions. The GTH pseudopotentials46−48 were employed for implicit treatment of the (C 1s, N 1s, O 1s, S 1s2s2p, and Na 1s) core levels, molecularly optimized Gaussian basis sets (TZVP-GTH)49 to describe the Kohn−Sham orbitals, and an auxiliary plane wave basis (with a cutoff of 400 Ry) to describe the electron density. The same system and ensemble parameters (i.e., NVT with a target temperature of 300 K) were used as in the classical MD simulations, and the end points of their respective trajectories were directly imported as initial conditions for simulation of the singlet electronic ground state (S0). The singlet systems were equilibrated for 10 ps and subsequently sampled for 50 ps. To describe the lowest triplet electronic state (T1), simulations with an unrestricted Kohn−Sham framework and triplet multiplicity were directly continued from the endpoints of the singlet trajectories and subsequently sampled for at least 50 ps. A representative simulation input file is provided in the Supporting Information for further reference. Geometry Optimization. Geometry optimizations at the DFT level of the molecular forms in the gas phase were performed with the same methodology as in the AIMD simulations. To describe the isolated proton-transfer compounds, however, a cubic simulation cell with side length 30 Å (with no additional atoms) was employed in combination with a Poisson solver for isolated systems.50 Complete Active Space Calculations. Single-state complete active space self-consistent field (CASSCF)51 calculations with subsequent second-order perturbation theory (CASPT2)52,53 were performed for the S0, T1, and T2 electronic states without explicit use of symmetry on the optimized geometries obtained from DFT for 2-TP and 2-MP with the OpenMolcas54,55 program. State averaged CASSCF calculations followed by multistate CASPT256 were also tried
Figure 1. Three proton-transfer related molecular forms 2thiopyridone (2-TP), 2-mercaptopyridine (2-MP), and their deprotonated anion (2-TP−) along with some of the thermodynamic parameters and light-driven processes, which may affect their relative concentrations.
been left undescribed7,9,10 or been approximated with classical continuum8,11,12 and other simplified models,9,29 despite the fact that all time-resolved measurements8,10−12 were performed in a solvent. In this work, we simulate and compare the solvatization of the three different forms in aqueous solution with ab initio molecular dynamics (AIMD) as a first step to bridge this apparent gap. The method, as compared to previous models, provides an explicit consideration of the solvent and its dynamics. As it also explicitly describes the electronic structure of the system, it enables a comparison of the lowest singlet and lowest triplet electronic states (S0 and T1), both of which play a key role in the ESPT process,12 to connect differences in the electronic structure to solute−solvent interactions. Furthermore, the direct inclusion of the electronic structure also naturally enables the spontaneous breaking and formation of chemical bonds, which in particular allows for solventmediated protonation or deprotonation of the N and S sites of the compound. Consequently, from the obtained trajectories, we can evaluate the stability of the different molecular forms and electronic states, to explain how it is influenced by the aqueous surrounding. For reference, the singlet simulations are also compared to simulations performed with a simple classical force-field approach. Finally, from the radial and spatial distribution functions (SDFs), we present a complete three-dimensional picture of the hydration and derive minimal representative solvatization complexes, which may be considered in future works to explicitly describe the solvatization in terms of system sizes that should prove computationally tractable for more demanding electronic structure methods.
■
COMPUTATIONAL METHODS The bulk of the results in the present study were obtained from simulations of aqueous solutions with ab initio molecular dynamics (AIMD) based on density functional theory (DFT). Classical molecular dynamics (MD) simulations based on a force-field approach were additionally performed for comparison and to generate initial conditions for the AIMD simulations. To describe differences in the gas phase, geometry optimizations of the isolated proton-transfer forms were 5556
DOI: 10.1021/acs.jpcb.9b03084 J. Phys. Chem. B 2019, 123, 5555−5567
Article
The Journal of Physical Chemistry B
Figure 2. Frontier molecular orbitals of the 2-TP S0 (a), 2-MP S0 (b), and 2-TP− S0 (c) from the DFT gas-phase geometry optimizations and frontier (spin-up) molecular orbitals of 2-TP(aq) T1 (d) from the AIMD simulation at two consecutive time steps.
radial) distribution around a reference site. The SDF is calculated by centering the reference site and rotationally aligning the system with respect to two more “rotational reference” sites, subsequently binning the other particles (averaged over the trajectory) into a three-dimensional Cartesian grid. The resulting values (of the three-dimensional grid) are typically presented as isosurfaces, in our cases isosurfaces of positive deviations from the simulation box averaged density of the binned particle (H and O atoms of the water, respectively). SDF data was generated with the trajectory analyzer and visualizer (TRAVIS),60 where the simulation boxes were aligned by centering the N atom, aligning the C1 atom onto the positive x axis, and aligning the S atom into the xy-plane with a positive y-value; binning was performed on a cubic 100 × 100 × 100 grid with a bin length of 0.07 Å. SDFs and molecular structures were visualized with VMD. Line plots were obtained with the SciPy61 Python packages NumPy62 and Matplotlib.63 All figures were built together and labeled with Inkscape.64
for the two triplet states, but they showed no appreciable differences. As in our previous study,12 the active space included the four occupied π orbitals, the occupied lone pair (S-based σnS for 2-TP and N-based σnN for 2-MP), and the three unoccupied π* orbitals, for a total of 10 active electrons; an imaginary shift of 0.1 hartree was applied to the CASPT2 iterations to accelerate convergence. Calculations were performed both with and without a polarizable continuum model (PCM) for water, equilibrated to the respective states.57 The 6-31G basis set was employed for all presented results; a comparison with the 6-31G** basis and the relativistically contracted atomic natural orbitals basis (ANO-RCC,58 both of double and triple ζ qualities, including scalar relativistic corrections through a Douglas−Kroll Hamiltonian) for the gas-phase calculations of 2-TP showed only small (≤0.08 eV) differences in the relative energies. Trajectory Analysis and Visualization. To analyze the structure of the solvatization, the radial distributions functions (RDFs) and spatial distribution functions (SDFs) were generated between key sites of the compound in its different forms and the aqueous solvent from the trajectories. RDF data, structural data, and molecular orbitals were obtained with visual molecular dynamics (VMD).59 An SDF is the threedimensional analogue of the well-known RDF, which describes the probability of finding a particle as a spatial (as opposed to
■
RESULTS AND DISCUSSION Electronic and Structural Stability. To better understand the solvatization structure presented in later sections, we first characterize all of the considered species and evaluate their
5557
DOI: 10.1021/acs.jpcb.9b03084 J. Phys. Chem. B 2019, 123, 5555−5567
Article
The Journal of Physical Chemistry B
Table 1. Excitation Character and Relative Energies (in Units of Electronvolt, w.r.t. 2-TP S0 for Each Level of Theory) of the S0 and T1 States of the 2-TP and 2-MP Molecular Formsa DFT state
2-TP(g)
2-TP(aq)
state
2-MP(g)
2-MP(aq)
S0 T1 (σnS → π*) T2 (π → π*)
0.00 2.14 2.75
0.00 ± 0.23 3.18 ± 0.16
S0 T1 (π →π*) T2 (π → π*)
0.10 3.36 7.53b
0.35 ± 0.10 3.82 ± 0.15
state
2-TP(g)
2-TP(PCM)
state
2-MP(g)
2-MP(PCM)
S0 T1 (π → π*) T2 (σnS → π*)
0.00 2.46 2.97
0.00 2.64 3.10
S0 T1 (π → π*) T2 (π → π*)
0.12 3.93 5.19
0.31 4.09 5.33
CASPT2
a
Results are presented for DFT gas-phase optimizations (g), sampling of the AIMD simulations in aqueous solution (aq), and single-state CASPT2 calculations in the gas phase (g) and with a polarizable continuum model (PCM) for water. Additionally, T2 states were calculated with the maximum overlap method for the DFT gas phase and with CASPT2. bThe calculation, depending on the initial conditions, converges to three different π → π* states with varying π-conjugations and relative energies of 7.53, 7.66, and 10.86 eV, which in all cases constitute a clear overestimation.
overlap method DFT, depending on the initial conditions, converges to three different solutions with varying πconjugation at clearly overestimated energies. Overall, it is clear that the excitation characters of low-lying triplet states in 2-TP are more difficult to assign than in 2-MP and depend not only on experimental conditions such as choice of solvent or phase but also on computational choices of methods and approximations. As several of the previous computational studies on this particular system have not reported the excitation character of their triplet states,8,10,11 they may therefore unknowingly even have considered “different” T1 states. In fact, performing as a test the same DFT gas-phase geometry optimization with the B3LYP functional, the 2-TP T1 state converged instead to (π → π*) excitation, whereas the 2MP T1 character was preserved. The T1 excitation character of 2-TP is even more ambiguous in the aqueous AIMD simulations than in the gas phase. To assign its character, we produced the molecular orbitals from a few different time steps toward the end of the trajectory; the results from two of the snapshots are shown in Figure 2d. At 57.11 ps, the excitation can be confidently assigned as σnS → π*, but at 57.10 ps it is not even possible to clearly distinguish between σ and π symmetric orbitals. As the highest doubly occupied molecular orbital and the first singly occupied molecular orbital (SOMO1) both have a strong S 3p character and nearly equal energies, they may easily rotate with each other when affected by the solvent-induced symmetry breaking. Consequently, the excitation characters of T1 and T2 cannot be rigorously separated over the course of the trajectory, and both states (as assigned from the gas phase) are likely to be present at some point in the simulation. While the molecular orbitals could not reasonably be saved and inspected over the whole trajectory, we can infer from the solvatization structure presented in Radial Distribution Functions section and Spatial Distribution Functions section that the simulated triplet state is at most times of σnS → π* character. This is not the transient species previously detected at ps time scales upon photoexcitation in aqueous solution, as excitation from the σnS orbital would not open the experimentally observed second transient N 1s X-ray absorption resonance.12 However, as the probe is spectroscopically insensitive to this excitation, and in light of the discussion above, it is not at all unlikely that triplet states of both characters are present at this or other times in the measurements,12 and it is in any case meaningful to
stability as expressed through energetic shifts of their electronic states and structural distortions of their molecular forms. Electronic States. To characterize the variations in the electronic structure between the different species, we turn to a molecular orbital picture. Figure 2a−c shows the frontier molecular orbitals obtained for the DFT geometry optimizations of the different forms in the singlet electronic state. While some structural and orbital relaxation occurs upon triplet excitation for all of the forms, the T1 electronic state can to a good degree be characterized as excitation from the highest occupied molecular orbital (HOMO) to the lowest unoccupied molecular orbital (LUMO) as presented. Interestingly, the triplet excitation character thereby differs between the forms: σnS → π* for 2-TP and 2-TP− and π → π* for 2-MP. In a simple orbital picture, this is directly rationalized by how the σnN orbital in 2-MP is 1.43 eV lower in energy (−6.11 eV) than the σnS orbital in 2-TP (−4.68 eV). Similar to the stabilization of valence excitations in thiobases,22 these orbital energies can in turn be explained by how the lone pair is strongly localized on an atom with a larger (N) or smaller (S) ionization potential. However, these results are in disagreement with our previously reported multiconfigurational calculations (corroborated below) that instead predict a π → π* character of the T1 state, both for vertical excitation from 2-TP S0 and for all molecular forms when structurally relaxed in their different electronic states.12 To ensure direct comparability to the current DFT calculations, state-specific CASPT2 calculations were performed at the DFT optimized geometries. In addition to the S0 and T1 states, the T2 state was also calculated with the maximum overlap method DFT65 and CASPT2 for the T1 optimized geometry. Excitation characters and relative energies are presented for both DFT and CASPT2 in Table 1 for 2-TP and 2-MP. For reasons that are described later in the section, 2-TP− is not included in this comparison. CASPT2, in contrast to DFT but in accordance with previous results,12 predicts a similar π → π* excitation character for the T1 state of both molecular forms. This can be understood as a multiconfigurational effect as the π* orbitals undoubtedly correlate better with the π orbitals than the σ symmetric lone pairs. Yet, both levels of theory predict a T2 state with the opposite character and only 0.5−0.6 eV higher energy than T1 for 2-TP. CASPT2 predicts a secondary π → π* character T2 state at roughly twice the energy separation for 2-MP, whereas the maximum 5558
DOI: 10.1021/acs.jpcb.9b03084 J. Phys. Chem. B 2019, 123, 5555−5567
Article
The Journal of Physical Chemistry B consider how the differences in electronic structure connect to differences in the solvatization for the compound. From the presented scrutinization, we conclude that the excitation character of 2-TP T1 is highly ambiguous and must be kept in mind for the analysis, but we will, in the discussion of our AIMD trajectories, for simplicity, mainly refer to the different species as the singlet and triplet states of our simulations for all molecular forms. To further evaluate the stability of the molecular forms and their electronic states, we next consider the solvent-induced shifts of the relative energies. The DFT gas-phase geometry optimizations are compared with sampled values from the AIMD simulations and, additionally, with the CASPT2 results obtained with and without the PCM model for water; all of the values are presented in Table 1. As a state-of-the-art electronic structure method, CASPT2 may be considered a benchmark for the gas-phase calculations to assess the overall performance of the DFT method. Oppositely, the explicit solvatization model provided by the AIMD simulations may be considered as a reference for the solvent-induced energy shifts to evaluate the applicability of the simplified PCM model. We note that PCM methods in principle measure differences in free energy, as opposed to the purely electronic energy of the AIMD simulations, but still find a comparison relevant as the two approaches are often comparably applied to describe solventinduced energy shifts in, e.g., simulations of spectroscopy. As also the case in a previous work,12 the relative energies of the protonated and deprotonated forms (i.e., 2-TP and 2-MP against 2-TP−) can unfortunately not be directly compared, as none of the calculations describes a solvent-mediated proton abstraction (i.e., for all of the calculations of 2-TP−, the proton is simply removed from the system, which results in a net negative charge for the gas phase and PCM calculations, whereas charge neutrality in the AIMD simulations is instead achieved with a negative counterion), which rules out a comparison of unshifted singlet energies. As seen in the Molecular Structure section, triplet 2-TP− does not remain stable in the AIMD simulations, which further rules out a meaningful comparison between the two electronic states of 2TP−. Yet, the energies obtained for both electronic states of 2-TP and 2-MP can be directly compared for all of the calculations. 2-TP is slightly stabler than 2-MP in the gas phase, as estimated from the singlet-state electronic energies at both levels of theory, but the energy difference is completely canceled (
