A Mixed Quantum-Classical Molecular Dynamics Study of anti-Tetrol

Article pubs.acs.org/JPCB

A Mixed Quantum-Classical Molecular Dynamics Study of anti-Tetrol and syn-Tetrol Dissolved in Liquid Chloroform II: Infrared Emission Spectra, Vibrational Excited-State Lifetimes, and Nonequilibrium Hydrogen-Bond Dynamics Kijeong Kwac and Eitan Geva* Department of Chemistry, University of Michigan, Ann Arbor, Michigan 48109-1055, United States ABSTRACT: The effect of vibrational excitation and relaxation of the hydroxyl stretch on the hydrogen-bond structure and dynamics of stereoselectively synthesized syntetrol and anti-tetrol dissolved in deuterated chloroform are investigated via a mixed quantum-classical molecular dynamics simulation. Emphasis is placed on the changes in hydrogenbond structure upon photoexcitation and the nonequilibrium hydrogen-bond dynamics that follows the subsequent relaxation from the excited to the ground vibrational state. The propensity to form hydrogen bonds is shown to increase upon photoexcitation of the hydroxyl stretch, thereby leading to a sizable red-shift of the infrared emission spectra relative to the corresponding absorption spectra. The vibrational excited state lifetimes are calculated within the framework of Fermi’s golden rule and the harmonic-Schofield quantum correction factor, and found to be sensitive reporters of the underlying hydrogen-bond structure. The energy released during the relaxation from the excited to the ground state is shown to break hydrogen bonds involving the relaxing hydroxyl. The spectral signature of this nonequilibrium relaxation process is analyzed in detail.

I. INTRODUCTION The structure and dynamics of hydrogen bonds (H-bonds) have received much experimental and theoretical attention due to their important role in determining the thermodynamical and dynamical properties of a large number of key chemical and biological systems.1−38 Recently, Vöringer and co-workers reported using time-resolved infrared (IR) spectroscopy to investigate the intramolecular H-bond structure and dynamics in stereoselectively synthesized tetrols.39,40 The latter provide a low-dimensional model of an H-bond network. More specifically, in the all-syn tetrol (syn-tetrol), the hydroxyls are oriented in a manner designed to promote the formation of Hbonds among all four hydroxyls, while in the all-anti tetrol (anti-tetrol), the hydroxyls are oriented in a manner designed to minimize H-bonding between the hydroxyls (see Figure 1). In a previous paper,38 we employed a mixed quantum-classical molecular dynamics simulation in order to elucidate the equilibrium H-bond structure of the above-mentioned syntetrol and anti-tetrol dissolved in liquid chloroform. We found that H-bonding in this system is very sensitive to the conformation of the tetrol molecules and presented a molecular model of these systems which reproduces the infrared (IR) absorption spectra of the hydroxyl stretch rather accurately. In doing so, we were able to elucidate the rather intricate relationship between the absorption IR spectrum and the underlying H-bond structure. Within our model, syn-tetrol is found to support a stable four-member H-bonded hydroxyl chain. In contrast, anti-tetrol was observed to exhibit a rather © 2013 American Chemical Society

Figure 1. Chemical structure of anti-tetrol (A) and syn-tetrol (B). The hydroxyl numbering convention is defined as shown. The three carbon and one oxygen atoms in the gray area defines the dihedral angle ψ.

wide distribution of partially H-bonded structures that differ with respect to which hydroxyls form H-bonds as well as the directionality of the H-bonds. In the present paper, we extend the molecular model and mixed quantum-classical methodology developed in ref 38 to understanding the H-bond structure in the excited state and its effect on the emission spectra, calculate the excited-state Received: August 29, 2013 Revised: October 24, 2013 Published: October 25, 2013 14457

dx.doi.org/10.1021/jp408580n | J. Phys. Chem. B 2013, 117, 14457−14467

The Journal of Physical Chemistry B

Article

lifetimes and explain them in terms of the underlying H-bond structure, and investigate the nonequilibrium H-bond dynamics that follows vibrational relaxation from the excited to the ground state, as well as its spectral signature. The structure of the remainder of this paper is as follows. The underlying theory and simulation techniques are outlined in section II. The results are presented and discussed in section III. Concluding remarks are made in section IV.

Here, i = 1, 2, 3, 4 denotes the four hydroxyls, {ψ1 = 60°, ψ2 = −60°, ψ3 = 180°} are the three values of the dihedral angle ψ that define stable conformations (see Figure 1),38 {pψ1, pψ2, pψ3} are the equilibrium probabilities of those conformations (assumed to be independent of vibrational state), and

II. METHODS A. Mixed Quantum-Classical Molecular Dynamics. The mixed quantum-classical methodology and molecular model used here are similar to those employed in ref 38 and will therefore only be outlined briefly. Simulations were performed on a system that consists of a single tetrol molecule and 247 deuterated chloroform (CDCl3) molecules. One of of the four hydroxyl stretches is treated as quantum-mechanical, while the remaining degrees of freedom (DOF) are treated as classical. For a given configuration of the classical DOF, Q, the vibrational energy levels, En(Q), and corresponding stationary wave functions, Ψn(q; Q), of the quantum hydroxyl stretch correspond to solutions of the adiabatic Schrödinger equation:

where

Had(q ̂, p ̂ ; Q)Ψn(q; Q) = En(Q)Ψn(q; Q)

Sψ(ij)(t ) = e−t /2T1⟨μ01(t )μ10 (0) exp[−i

⟨···⟩1, ψj =

∫

dqΨ*n (q , Q) ∇Q V (Q, q)̂ Ψn(q , Q)

(1)

(2)

3

4

j=1

j

∑∫ i=1

0

∞

dtSψ(ij)(t ) exp(iωt )

(5)

exp[−βH1(Q 0, P0) + Vψj(Q 0)] Z1, ψj

μij (t ) = μ′⟨Ψi(Q t)|q|̂ Ψj(Q t⟩

···

(7)

where

μ′ =

(3)

∂μ ∂q

q = q0

(8)

Non-Condon effects are accounted by taking advantage of the previously found linear mapping relations between μ′ and the electric field along the hydroxyl bond, which can be obtained via QM/MM ONIOM calculations for each one of the four hydroxyls.38 C. Excited State Lifetime and Vibrational Relaxation. The population relaxation rate constant from the first-excited to the ground state of the ith hydroxyl stretch, as calculated within the framework of Fermi’s golden rule, is given by14

As a result, the structure and dynamics of the classical DOF are affected by the state of the quantum DOF, and vice versa. This should be contrasted with the common practice of assuming equilibrium structure and dynamics of the classical DOF, which are thereby assumed to be unaffected by the state of the quantum DOF. It should be noted that the latter treatment can only be justified in the limit of weak coupling between the classical and quantum DOF. The fact that in our case the interaction between the quantum and classical DOF involves H-bonding suggests that this is not the case in the system under consideration here. Thus, a treatment that accounts for the feedback between the classical and quantum DOF is pursued. All the simulations in this paper were carried out using the AMBER10 program package,41 where the source code has been modified in order to implement the mixed quantum-classical algorithm. The reader is referred to refs 13−15 and 38 for a more detailed discussion of the algorithm and its implementation. B. IR Emission Spectra. Within the mixed quantumclassical treatment, the IR emission spectrum can be calculated on the basis of the following formula:42 I(ω) = Re ∑ pψ

dτω10(i)]⟩1, ψj

where μ01(t) is the transition dipole moment, Vψj (Q0) is the harmonic restraining potential for the dihedral angle ψ to the value ψj, H1(Q0, P0) = KB(P0) + E1(Q0), and Z1,ψj = ∫ dQ0 ∫ dP0 exp[−βH1(Q0,P0) + Vψj (Q0)] is the classical DOF Hamiltonian and partition function for a given conformation when the quantum hydroxyl is in its f irst excited vibrational-state. (Q0,P0) represents the initial coordinates and (i) momenta of the classical DOF, ω(i) 10 (t) ≡ ω10 (Qt) ≡ [E1(Qt) − E0(Qt)]/ℏ is the transition frequency, Qt represents the classical DOF configuration after evolving for a period of time t on the first-excited adiabatic potential energy surface. Finally, T1 is the lifetime of the first-excited state, which we set to T1 = 1.0 ps, which corresponds to the average lifetime reported in ref 40. The transition dipole moment, μ01(t) (see eq 5), is given by

is the adiabatic Hamiltonian. Here, q̂ and p̂ are the coordinate and momentum operators of the quantum DOF, n = 0, 1, ... denotes a vibrational state, Kq(p̂) is the vibrational kinetic energy of the quantum hydroxyl stretch, and V(q̂, Q) is the overall potential energy. Within our mixed quantum-classical approach, the classical DOF evolve according to the classical equation of motion with a force that depends on the state of the quantum system: F(Q) = −

t

(6)

where Had(q ̂, p ̂ ; Q) = Kq(p ̂) + V (q ̂, Q)

∫ dQ 0∫ dP0

∫0

k 01(⟨ω10(i)⟩1) = Q HS(⟨ω10(i)⟩1) ⟨exp(i

∫0

t

∞

∫−∞ dt exp(i⟨ω10(i)⟩1t )

(i) dτδω10(i)(Q τ))[d10 (Q 0) ·Q̇ 0]

(i) [d10 (Q t) ·Q̇ t]⟩1

(9)

The corresponding excited-state lifetime is given by T(i) 1 = (i) (i) k−1 01 (⟨ω10 ⟩1). Here, ⟨ω10 ⟩1 is the average emission frequency of (i) (i) the ith hydroxyl, δω(i) 10 (Qτ) = ω10 (Qτ) − ⟨ω10 ⟩1 is the deviation of the instantaneous emission frequency from the average value, (i) d10 (Q) = ⟨ψ1(i)(Q)|∇Q |ψ0(i)(Q)⟩

(10)

is the nonadiabatic coupling vector,43 and Q HS(ω) = [e βℏω /2β ℏω/(1 − e−βℏω)]1/2

(4) 14458

(11)

dx.doi.org/10.1021/jp408580n | J. Phys. Chem. B 2013, 117, 14457−14467

The Journal of Physical Chemistry B

Article

is the harmonic-Schofield quantum correction factor,44−46 which accounts for the quantum nature of the high frequency accepting modes.45 The choice of this quantum correction factor was motivated by its successful usage in the work of Skinner and co-workers on HOD/D2O47 and its agreement with experiment in the case of methanol/CCl4 mixtures.14 The probability for finding a hydroxyl stretch which is excited at time t = t0 in the excited state at a later time, t > t0, is given by15 P01(t ) = k10e−k10(t − t0)

50 ps constant temperature (300 K) and constant energy mixed quantum-classical production run was carried out with a constrained dihedral angle for each of the three conformations.

III. RESULT AND DISCUSSION A. Emission Spectra and Stokes Shift. Figures 2 and 3 show the calculated conformation-specific and conformationally

(12)

Thus, the probability for a transition from the excited to the ground state within the simulation time step δ is given by15 δP01(t ) =

t+δ

∫t

dt ′P01(t ′) = e−k10t (1 − e−k10δ)

(13)

Importantly, each vibrational relaxation event is accompanied by energy released to the classical DOF. Within the mixed quantum-classical treatment, this energy release manifests itself as a momentum jump, which is given by15 ̂ (Q)·P ̅ ) ΔP(Q) = M1/2{sgn(d10 ̂ (Q) ·P ̅ )2 + 2ℏω (Q) (d10 01 ̂ (Q)·P ̅ )}d ̂ (Q) − (d10 10

Here, M

1/2

is a diagonal matrix with (Mj)

(14) 1/2

Figure 2. The calculated conformation-specific and conformationally averaged absorption (green) and emission (red) spectra of anti-tetrol. Also indicated are the relative subpopulations of the three different conformations (in %).

along the diagonal, ̂ = d /|d̅ | is a = (d10,j/(Mj) ) and d10 10 10

P̅j = (Pj/(Mj) ), d̅10,j unit vector along d̅10. D. Simulation Details. The force fields parameters used in this work are the same as in ref 38. The AM1-BCC method was used for assigning partial charges to the atoms that constitute the tetrol molecule. It should be noted that the AM1-BCC partial charges were found to give a better overall agreement with the experimental IR absorption spectra of both syn-tetrol and anti-tetrol than RESP partial charges.38 Simulations were carried out on a cubic simulation cell with periodic boundary conditions containing a single anti-tetrol or syn-tetrol molecule and 247 deuterated chloroform (CDCl3) molecules. The simulation time step was 1.0 fs. The cutoff radius for the nonbonded interaction was 10 Å. The long-range electrostatic interaction is taken into account by using the particle mesh Ewald method.48,49 The criteria used in this work to define H-bonding were that the (1) O···H distance is less than 2.65 Å and (2) O···H−O angle is greater than 120° (”···” denotes the hydrogen bond). The volume of the simulation cell was determined by a 800 ps constant pressure (1 atm) and constant temperature (300 K) fully classical equilibration run using Langevin dynamics for the temperature coupling with the coupling constant γ = 1.0 ps and Berendsen algorithm50 for the pressure coupling with the pressure relaxation time of 0.2 ps. Initial configurations were then collected at 2 ps intervals during a 200 ps fully classical constant volume and constant temperature (300 K) equilibrium trajectory. Each initial configuration was used to start a 200 ps constant temperature (300 K) classical equilibration run, after reassigning initial velocities from the Maxwell distributions and applying a restraint in the form of a harmonic potential V = k(ψ−ψ0)2 to restrain the dihedral angle ψ (k = 30.0 kcal/mol/ rad2 and ψ0 = 60°, −60°, or 180°). This was followed by a 35 ps mixed quantum-classical constant temperature (300 K) equilibration run with a constrained dihedral angle. Finally, a 1/2

1/2

Figure 3. The calculated conformation-specific and conformationally averaged absorption (green) and emission (red) spectra of syn-tetrol. Also indicated are the relative subpopulations of the three different conformations (in %).

averaged absorption and emission spectra of anti-tetrol and syntetrol, respectively. For both syn-tetrol and anti-tetrol, and regardless of conformation, the emission spectra are seen to be red-shifted and broader in comparison to the absorption spectra. This can be traced back to a higher propensity for forming H-bonds in the excited state in comparison to the ground state. It should be noted that similar behavior was observed in methanol/CCl4 mixtures, where photoexcitation of 14459

dx.doi.org/10.1021/jp408580n | J. Phys. Chem. B 2013, 117, 14457−14467

The Journal of Physical Chemistry B

Article

excited sates are very similar, and the main difference between the two is that the frequency range is red-shifted and broader in the excited state compared to the ground state. In the case of Figure 5, the bond length is seen to be larger and more widely spread in the excited state compared to the ground state. Both trends can be traced back to the higher propensity for Hbonding in the excited state. It should also be noted that the mapping relations for the different hydroxyls are very similar, which can be traced back to the dominant role of H-bonding as opposed to location along the molecular backbone. The mapping relations were obtained by fitting polynomial functions (see Figures 4 and 5 and Appendix). Once established, these mapping relations can be implemented within a fully classical molecular dynamics (MD) simulation, where the electric field along the tagged hydroxyl is calculated at each time step, and then used to get the stretch frequency and bond length via the mapping relations. In so doing, we bypass the computationally expensive on-the-fly diagonalization of the adiabatic Hamiltonian.14 It should be noted that using the mapping relations involves approximating the forces in the following manner:

the methanol hydroxyl stretch was seen to give rise to increased propensity for forming H-bonds, accompanied by a similar spectral signature.14 The red shift between the absorption and emission spectra (the so-called Stokes shift) is 21.8 cm−1 and 39.7 cm−1, for antitetrol and syn-tetrol, respectively. The fwhm of the main spectral feature is increased from the absorption to the emission spectra by ∼18% for the anti-tetrol (see Figure 2) and ∼22− 25% for the syn-tetrol (see Figure 3). The larger red shift and broadening of the emission spectra in syn-tetrol in comparison to anti-tetrol is indicative of and can be traced back to the fact that the hydroxyls in syn-tetrol are oriented in a manner more favorable for forming H-bonds. B. Mapping Relations. Following a procedure similar to that employed in our previous work on a methanol/CCl4 mixture,13−15 we found that nonlinear relations can be established between the OH stretch frequency and bond length, and the electric field along the OH bond, for each of the four hydroxyls (see Figures 4 and 5, respectively, for anti-tetrol;

⟨ψi(Q)|F(q , Q)|ψi ⟩(Q) ≈ F(⟨ψi(Q)|q|ψi(Q)⟩, Q)

(15)

To demonstrate the quantitative accuracy of our mapping procedure, we compare the calculated absorption spectra obtained from a fully classical simulation, with the help of the mapping relations, with those obtained via mixed quantumclassical MD simulations based on on-the-fly diagonalization (see Figure 6). The spectra obtained via both methods essentially coincide, thereby testifying for the accuracy of the mapping relations.

Figure 4. Correlation plots between the hydroxyl frequency and the electric field along the hydroxyl bond axis. Also shown are the nonlinear mapping relations between these two quantities.

Figure 5. Correlation plots between the hydroxyl bond-length and the electric field along the hydroxyl bond axis. Also shown are the nonlinear mapping relations between these two quantities.

Figure 6. Comparison between the calculated absorption spectra obtained from a fully classical simulation, with the help of the mapping relations (blue), with those obtained via mixed quantum-classical MD simulations based on on-the-fly diagonalization (red).

similar mapping relations can be established for syn-tetrol). In the case of Figure 4, the mapping relations in the ground and 14460

dx.doi.org/10.1021/jp408580n | J. Phys. Chem. B 2013, 117, 14457−14467

The Journal of Physical Chemistry B

Article

In contrast, the fourth hydroxyl exhibits opposite trends. More specifically: (1) The excited-state lifetimes of the fourth hydroxyl are seen to be significantly longer in comparison to the other three hydroxyls. This can be traced back to the possibility of forming an H-bond between the fourth hydroxyl and the backbone oxygen, which is weaker than the H-bonds between hydroxyls. (2) The lifetime of the fourth hydroxyl in syn-tetrol is seen to be longer than in the anti-tetrol. This can be traced back to the fact that, in syn-tetrol, the fourth hydroxyl forms an H-bond with the backbone oxygen as a donor, whereas in anti-tetrol there is a significant subpopulation where the fourth hydroxyl forms an H-bond with the third hydroxyl as a donor. Since an H-bond with the backbone oxygen is weaker than an H-bond with a neighboring hydroxyl, the fourth hydroxyl is in fact H-bonded more strongly in anti-tetrol than in syn-tetrol, thereby giving rise to a shorter lifetime. In the next step we consider the dependence of the lifetime on the excitation frequency, in an attempt to compare with the corresponding experimental results reported in ref 40. To this end, it is important to note that the absorption spectra of the four hydroxyls overlap extensively,38 so that photoexciting the system at a given frequency would usually correspond to exciting all four hydroxyls. However, the fraction of each hydroxyl in the photoexcited subpopulation can change with the excitation frequency, and given the fact that the lifetime differs from one hydroxyl to another this can give rise to a dependence of the lifetime on the excitation frequency. To estimate this effect, we define the excitation-frequency-dependent population relaxation function as follows:

The results reported in Figures 7−11 were obtained from fully classical MD simulations with the help of the mapping

Figure 7. The dependence on the excitation frequency of the calculated excited state lifetimes for anti-tetrol (ψ = −60° conformation, red triangles) and syn-tetrol (ψ = 60° conformation, blue circles). Also shown are the corresponding experimental results (adopted from ref 40).

relations. The only exception is the calculation of the nonadiabatic coupling vector at the time of the momentum jump, which requires an explicit on-the-fly diagonalization of the adiabatic Hamiltonian in order to obtain the adiabatic wave functions (see eq 10). C. Excited State Lifetimes. The calculated excited state lifetimes, eq 9, evaluated at the average excited state transition frequencies of the four hydroxyls, are shown in Table I for anti-

4

P(ωex , t ) =

∑ p(i) (ωex ) exp[−k 01(i)(ωex )t ] i=1

Table I. Calculated Excited State Lifetimes (in Picoseconds) of the Four Hydroxyl Stretches in anti-Tetrol and synTetrola anti avg frequency OH1 OH2 OH3 OH4

3489.6 3443.6 3421.6 3473.8

Here, ωex is the excitation frequency, k(i) 01 is the vibrational relaxation rate constant of the ith hydroxyl calculated at the excitation frequency, and 0 < p(i)(ωex) < 1 is the fraction of photoexcited ith hydroxyl in the overall photoexcited subpopulation. The latter are calculated by assuming that p(i)(ωex) ∝ I(i)(ωex), where I(i)(ω) is the absorption spectrum of the ith hydroxyl (adopted from ref 38). The frequencydependent lifetime was then obtained by fitting P(ωex,t) to an exponential function. The resulting excitation frequency dependence of the lifetime is shown in Figure 7, side by side with the corresponding experimental results. While the agreement is not perfect, several of the experimental trends are reproduced, namely, (i) the lifetimes of syn-tetrol and anti-tetrol are relatively similar in the low-frequency domain (