How does microwave irradiation affect the mechanism of water reorientation?

How does microwave irradiation affect the mechanism of water reorientation?

Journal of Molecular Liquids 302 (2020) 112522 Contents lists available at ScienceDirect Journal of Molecular Liquids journal homepage: www.elsevier...

2MB Sizes 0 Downloads 5 Views

Journal of Molecular Liquids 302 (2020) 112522

Contents lists available at ScienceDirect

Journal of Molecular Liquids journal homepage: www.elsevier.com/locate/molliq

How does microwave irradiation affect the mechanism of water reorientation? Tomaž Mohoriˇc a , Urban Bren a, b, * a b

Faculty of Chemistry and Chemical Technology, University of Maribor, Smetanova 17, Maribor SI-2000, Slovenia National Institute of Chemistry, Hajdrihova 19, Ljubljana SI-1001, Slovenia

A R T I C L E

I N F O

Article history: Received 12 November 2019 Received in revised form 6 January 2020 Accepted 16 January 2020 Available online 20 January 2020 Keywords: Non-equilibrium microwave effects Water reorientation H-bond exchange Molecular dynamics

A B S T R A C T We employed non-equilibrium molecular dynamics simulations to investigate the hydrogen-bond dynamics in bulk water upon microwave heating. Microwave irradiation has been taken into account implicitly by performing simulations at a rotational temperature exceeding the translational one (Tr > Tt ). We have shown that the effect of an increased Tr results in a modified geometry of the average path for hydrogenbond switch and in decreased characteristic decay times for the reorientation of a water molecule. We also confirmed the validity of the extended jump model to describe hydrogen-bond dynamics under the studied non-equilibrium conditions. © 2020 Elsevier B.V. All rights reserved.

1. Introduction Many of the important properties of liquid water stem from its ability to form dynamic and labile H-bond networks [1-5]. Connectivity of these networks changes permanently and in a large extent due to the reorientation of water molecules. This constantly rearranging H-bond network influences a wide range of phenomena, such as proton transfer [6, 7], SN 2 reactions in aqueous solutions [8], proton transport in bulk water [9-11] as well as in water wires through biological channels, where water reorientation probably represents the rate limiting step [10, 12]. Water reorientation plays an important role in lability of hydration shells of proteins and other biological macromolecules [13-18], hence influencing the molecular recognition of ligands [2, 19, 20]. It has been well established that the mechanism for water molecule reorientation consists of largeamplitude angular jumps during which a central water molecule breaks a former and forms a new H-bond [21-23]. This notion provided a starting point to examine the nonequilibrium microwave effects on bulk water, in particular its influence on geometrical and dynamical properties of H-bond switches. Upon microwave irradiation, water molecules tend to align their dipole moments with the alternating external field, thereby

* Corresponding author at: Faculty of Chemistry and Chemical Technology, University of Maribor, Smetanova 17, Maribor SI-2000, Slovenia. E-mail address: [email protected] (U. Bren).

https://doi.org/10.1016/j.molliq.2020.112522 0167-7322/© 2020 Elsevier B.V. All rights reserved.

enhancing their rotational motion [24-26]. Under continuous microwave irradiation the excess rotational kinetic energy cannot fully dissipate into the remaining degrees of freedom (translations and vibrations) and builds up in the form of faster rotations. The net effect results in the formation of “rotationally excited” water molecules, which may be described by a non-equilibrium state, where the rotational temperature exceeds the translational and vibrational one. It should be noted that the term “temperature” usually refers to the equilibrium temperature, thus a state where the rotational temperature exceeds the translational one represents the non-equilibrium effect of microwaves. The effects of microwaves have been extensively examined using non-equilibrium molecular dynamics with explicitly included electric (and possibly magnetic) fields. Such simulations have been successfully applied to study liquid water [27, 28], aqueous solutions of salt [28-30], proteins [31-33], ionic liquids [32] and some other systems [34, 35]. We have proposed an alternative approach in which the rotational and translational degrees of freedom are controlled by separate thermostats [36-38]. Using such an approach we are able to keep the rotational temperature (Tr ) above the translational one (Tt ). The resulting non-equilibrium state serves as a replacement model to study the effects of microwave irradiation. It should be pointed out that Auerbach and coworkers [39-41] have indeed provided convincing experimental evidence that under continuous microwave irradiation the rotational temperature in a system may be substantially higher than the translational one. Using quasielastic neutron scattering measurements on microwave

2

T. Mohoriˇc and U. Bren / Journal of Molecular Liquids 302 (2020) 112522

irradiated zeolite-guest systems they: (i) demonstrated the selective heating of polar substances versus the non-polar ones; (ii) determined the effective rotational and translational temperatures, where the former substantially exceeded the latter one [41]. These experiments support our approach to study the non-equilibrium microwave effects using molecular dynamics with separate thermostats for translational and rotational degrees of freedom. Moreover, the non-equilibrium microwave effects have been also reported to affect DNA hybridization [42], enzymatic catalysis [43] as well as electric conductivity of aqueous solutions [44]. It would be also beneficial to study the non-equilibrium state where the rotational temperature exceeds the translational one using analytical theories, such as [45], which was shown to successfully describe anomalous properties of liquid water. Such a theoretical approach would offer several advantages: beside the instantaneous computation of system’s properties it also provides an interpretation in terms of molecular physics. However, the issue with analytical approaches is that they are often limited to equilibrium cases, since integration over degrees of freedom usually employs only a single temperature. Although the main goal of the current study is to examine the effect of the increased rotational temperature on the geometric and dynamic properties of H-bond switches in bulk water, we additionally explored the effects of a corresponding increase in translational temperature. In Section 2 we provide the model and simulation details, which is followed by the presentation of the influence of elevated Tr and Tt on geometric, Subsection 3.1, and dynamic, Subsection 3.2, features of H-bond exchange.

2. Model and simulation details We performed molecular dynamics simulations using a constant number of particles and volume. The system consisted of 256 rigid SPC/E water molecules [46]. The volume was adjusted so that the number density of water molecules corresponded to a density of 1.0 g/cm3 (yielding a cube box with edge dimension of 19.7 Å). Intramolecular geometrical constraints (bond lengths and angles) were satisfied by employing the Euler’s equations for rigid body dynamics [47]. A radial cutoff of 9 Å was applied for the non-bonded Lennard Jones interactions, while for the long-range effect of electrostatic interactions the damped shifted force method was used [48]. Simulations were initiated with 500 ps long equilibration runs followed by 5 ns long production runs during which the data were collected. The timestep used was 0.5 fs and trajectory was sampled every 5 fs. Like in our previous studies we used an in-house developed code to separately solve the equations of translational and rotational motion [49-55]. This approach allowed us to assign translational and rotational temperatures, Tt and Tr , to separately control the kinetic energies of translational and rotational degrees of freedom. To keep these temperatures constant, we employed the Woodcock thermostat (which is identical to the Berendsen thermostat with its time constant tB equal to the simulation timestep Dt) [56]. Although Woodcock thermostat constrains the system’s kinetic energy, the energy of every single particle is still allowed to fluctuate. We have also verified (for details see the following paragraph) that this restraint does not influence dynamic properties which are the subject of this study (i.e. average jump path, orientational time correlation functions). Moreover, it is well known that Woodcock thermostat may lead to artifacts since in general system’s kinetic energy does not decrease. However, in our nonequilibrium microwaved case the rotational kinetic energy is actually being reduced through energy dissipation from rotational into translational degrees of freedom. Woodcock thermostat is indeed used to prevent the rotational temperature from decreasing to the equilibrium value (the heat flow according to the laws of thermodynamics

occurs from the higher Tr to the lower Tt ). In previous studies we explored alternative options: the use of the Andersen thermostat was found to be inefficient, while the Nose–Hoover thermostat, designed to correctly sample the canonical ensemble, was not well suited to tackle our problem. Woodcock thermostat turned out to be particularly suitable for our investigation, since it allowed us to strictly maintain the preset temperatures. Although in real systems the energy would dissipate from the degrees of freedom with higher temperature to those with lower temperature, this is the exact reason which makes it difficult to study non-equilibrium microwave effects. Hence, the purpose of Woodcock thermostat is not to replicate the real system in all its details, but mainly in enabling us to answer questions that are otherwise difficult to tackle. Notice that we studied the non-equilibrium but steady-state situation where Tr exceeds Tt , which serves as a replacement model for the direct incorporation of external electromagnetic fields into nonequilibrium molecular dynamics simulations. In order to maintain this steady-state situation we have to counter the natural phenomenon of continuous energy dissipation from excited rotational into colder translational degrees of freedom by enforcing rather strict thermostating conditions. Namely, we rescale particles’ velocities at every timestep such that they fully comply with the prescribed temperatures. This is ensured by setting the thermostat’s constant equal to the timestep (tB = Dt). In Berendsen thermostat the system’s temperature relaxes towards the preset value exponentially with a time constant tB . If we increase the thermostat’s time constant such that tB > Dt, the system is no longer forced to keep its temperature at a preset value, i.e. we allow the temperatures to drift and fluctuate. We found that average jump path and orientational time correlation functions from different thermostating conditions (at least up to tB = 20Dt) to completely overlap. Hence, dynamical properties are not significantly affected by the damped temperature fluctuations under velocity rescaling regime. To explore the effects of microwave–excited rotational motion we kept Tt fixed at 300 K, while Tr assumed values of 300, 400, 500 and 600 K. In addition, we also investigated the opposite case where only Tt values increased (to 400, 500 and 600 K), while keeping Tr at 300 K. While the first case mimicked the effect of faster rotational motion of water on H-bond exchange, the latter served as a model for the effect of accelerated translational motion of water. System at the equilibrium, Tr = Tt = 300 K, represented a reference state for exploring the effects of faster rotational or translational motion. To justify the choice of the temperature range studied, we refer to one of our previous studies [38], where the microwave power required to maintain one mole of water molecules at a given elevated rotational temperature above the translational temperature of 300 K was calculated. It was shown that a microwave power of 361, 795, and 1230 W is needed to maintain the rotational temperature at 400, 500, and 600 K, respectively [38]. Therefore, in a typical microwave reactor with a power of 850 W an increase in rotational temperature of up to 200 K above the translational temperature could in principle be achieved [38]. Hence, the studied temperature range corresponds well to situations that would typically occur in laboratory practice. Within our model of microwave heating we ignore the vibrational motion of water molecules for two main reasons. (i) In one of our previous publications we have demonstrated that the effect of an increased vibrational temperature, Tv , on the structural and dynamic properties of bulk and shell water is negligible up to at least Tv = 700 K [38]. (ii) Based on the fact that actual kinetic energy needed to excite bending motion of a water molecule is about 5 kcal/mol [57], which corresponds to a vibrational temperature of about 2500 K, we can safely assume that vibrational motion at the temperatures studied in the current work does not significantly differ from the equilibrium state.

T. Mohoriˇc and U. Bren / Journal of Molecular Liquids 302 (2020) 112522

3

Fig. 1. A) Schematic representation of the adopted definition of a H-bond. B) Schematic representation of the main geometric parameters describing the H-bond exchange. The O∗ H∗ covalent bond of the central water molecule is rotating in a clock-wise direction (designated by the gray curved arrow) from its initial H-bond acceptor, Oa , to the final one, Ob . Angle h is depicted as the transparent yellow area, while angle hOa O∗ Ob is represented by the blue area. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 2. Time evolution of H-bond coordination numbers. Top row: average number of received H-bonds by the initial acceptor Oa (red line), the final acceptor Ob (blue line) and the central O∗ atom (black line). Bottom row: average number of donated H-bonds by the initial acceptor Oa (red line), the final acceptor Ob (blue line), the exchanging O∗ H∗ bond (black line) and the non-exchanging O∗ H bond (green line). Different columns correspond to different thermostatting conditions: equilibrium at Tr = Tt = 300 K (left), non-equilbrium Tr = 600 K (middle) and non-equilibrium Tt = 600 K (right). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

4

T. Mohoriˇc and U. Bren / Journal of Molecular Liquids 302 (2020) 112522

3. Results and discussion The main focus of the current study is on the mechanism of water molecule reorientation in liquid phase. It has been well established that the main mechanism for water reorientation consists of large-amplitude angular jumps during a H-bond exchange. To detect H-bond exchange events we employed the approach outlined in previous studies [21, 22]. From the calculated trajectories we extracted all events for which the O∗ H∗ covalent bond formed a H-bond to the oxygen Oa for t < 0 and to the oxygen Ob for t > 0. We adopted tight geometric conditions for the H-bond existence [22], namely: rO∗ O < 3.1 Å, rH∗ O < 2.0 Å and hH∗ O∗ O < 20◦ - also see Fig. 1A) for schematic representation of these geometric parameters. For each of these switching events we examined the short time interval preceding and following it, as long as no other H-bond exchange occurred, and in Subsection 3.1 we investigate the time evolution of several H– bond parameters along the exchange trajectory under equilibrium and non-equilibrium conditions. In a similar manner we devote the Subsection 3.2 to closely inspect time correlation functions for the reorientation of water molecules. 3.1. Geometric considerations of the H-bond exchange We start the discussion by looking at the trends in the coordination number of water molecules. Changes in the coordination of water molecules involved in the H-bond switch are depicted in Fig. 2. The top row shows the time evolution of the number of received H-bonds by each of the oxygen atoms, while the bottom row analogously presents the evolution of donated H-bonds by the initial (Oa ) and final (Ob ) H-bond acceptor as well as both OH bonds on the central water molecule (O∗ H∗ and O∗ H). In equilibrium at Tr = Tt = 300 K the H-bond exchange starts with the initial H-bond acceptor becoming overcoordinated and the final

H-bond acceptor becoming undercoordinated. This difference in Hbond coordination facilitates the exchange of partners for the central covalent O∗ H∗ bond. Immediately after the switch, the H–bond coordination numbers of two acceptor oxygens are reversed and with the diffusion away of the initial H-bond acceptor they all tend towards their bulk values (≈1.8). Meanwhile the numbers of donated Hbonds by the Oa and Ob , as well as the coordination number of the O∗ H do not change significantly. Of course, the lack of H-bond coordination of the central O∗ H∗ covalent bond during the jump is clearly reflected in a decrease of donated H-bonds at t ≈ 0. In non-equilibrium states at Tr = 600 K or Tt = 600 K the coordination numbers are generally shifted towards lower values due to the smaller average number of H-bonds per water molecule, while the main characteristics of the curves are retained. Next we look at the time evolution of the angles and distances characteristic for a H-bond exchange. The main geometric parameters of the H-bond switch, which are delineated in the text below, are illustrated in Fig. 1B). Fig. 3 shows the time evolution of the oxygenoxygen distances of the three water molecules involved in the Hbond exchange in the form of a probability distribution. Note that red lines in Fig. 3 represent the time evolution of the average value of each geometric parameter. In equilibrium at Tr = Tt = 300 K the initial H-bond between the O∗ H∗ and Oa is characterized by a short distance rO*Oa (top row, left column), which steeply increases around t = 0, indicating a H-bond breakage, and is followed by a further increase in rO*Oa but with a smaller slope, as non-H-bonded water molecules diffuse apart. The probability distribution of the rO*Oa parameter is strongly peaked at ≈2.9 Å before the jump occurs and widens substantially after the breakage of the initial H-bond. Time evolution of a distance between the O∗ and Ob , rO*Ob (middle row), is a mirror image of the rO*Oa (top row): it starts with slowly decreasing values characteristic of non-H-bonded water molecules and around time t = 0 steeply decreases due to the H-bond formation, followed

Fig. 3. Probability distributions of the distance to the initial H-bond acceptor, rO*Oa (top row), distance to the final H-bond acceptor, rO*Ob (middle row), and of the distance between the initial and final H-bond acceptors, rOa Ob (bottom row) along the average jump path, where the system crosses the transition state at time t = 0. Different columns correspond to different thermostatting conditions: equilibrium at Tr = Tt = 300 K (left column), non-equilbrium Tr = 600 K (middle column) and non-equilibrium Tt = 600 K (right column). The red line plotted on top of each probability distribution represents the time evolution of the average value. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

T. Mohoriˇc and U. Bren / Journal of Molecular Liquids 302 (2020) 112522

by a more or less constant value of about rO*Ob ≈ 2.9 Å, indicative of H-bonded water molecules. Contrary to the described behaviour of rO*Oa and rO*Ob is the distance between the O atoms of the initial and final H-bond acceptor, rOa Ob (bottom row), which is almost constant throughout the examined time interval, except for a slight decrease around the time origin t = 0, where the H-bond exchange takes place. Furthermore the probability distribution of rOa Ob is rather wide (from ≈3 Å to ≈6 Å and is mainly peaked at ≈3 Å with a less pronounced second peak at ≈4.5 Å. Note that around t = 0 the second peak in the probability distribution of rOa Ob also converges towards 3 Å. Fig. 3 also depicts the time evolution of the same geometric parameters for rotationally (middle column) and translationally (right column) excited water molecules. First we note that an increase in Tr or Tt does not affect the rO*Oa and rO*Ob distances when the two molecules are H-bonded to O∗ . However, an increase in Tr inhibits the diffusion of non-H-bonded water molecules that drift apart more slowly - see the average rO*Oa after and the average rO*Ob before the switch (red lines). This can be attributed to the diminished empty space and reduced fluctuations in water number density upon a Tr increase [50]. All in all, an increase in Tr does not appear to affect the probability distributions of rO*Oa and rO*Ob to a significant extent. On the other hand an increased Tt drastically accelerates the diffusion of non-H-bonded water molecules as evidenced by a strong spread of the probability distributions and a steep increase of the average rO*Oa after and the average rO*Ob before the switch (red lines). The effects of increased Tr and Tt are markedly different for the rOa Ob parameter. A faster rotational motion effectively diminishes the first peak in the rOa Ob probability distribution – designating much smaller likelihood of the initial and final H-bond acceptors to be directly H-bonded and increasing the likelihood of them to be further apart. This is also reflected in a higher Tr shifting the average distance between the initial and final H-bond acceptor, rOa Ob , to

5

larger values (compare the average values as indicated by red lines). As revealed later, this is already indicative of an increase in the jump angle. On the other hand, a faster translational motion through faster diffusion substantially widens the probability distribution of rOa Ob away from the time origin, while retaining a greater likelihood for the two water molecules to be directly H-bonded near the switching event. These changes are reflected in a higher Tt increasing the average rOa Ob , but mainly on the account of faster translational diffusion of water molecules, while the initial and final acceptor still approach each other more than in the case of a corresponding Tr increase. Next we examine the time evolution of the various angles. Fig. 4 shows time dependence of probability distributions of angles between the O∗ H∗ covalent bond and the OO vector to the initial, hOa O∗ H∗ (top row), or final, hOb O∗ H∗ (middle row), H-bond acceptor. In equilibrium at Tr = Tt = 300 K the probability distribution of angle hOa O∗ H∗ between two H-bonded water molecules (top middle row, left column) is highly localized around 15◦ for t < 0 due to the librational motion of H-bond and becomes rather wide after the initial H-bond breakage (t > 0). The average value of hOa O∗ H∗ (red line, top row, left column) steeply increases to about 70◦ , while hOb O∗ H∗ (red line, middle row, left column) steeply decreases from about 70◦ to 15◦ due to the new H-bond formation. The effect of a higher Tr causes the distribution to widen, indicating a faster librational motion resulting in a greater equilibrium angle hOO∗ H∗ of H-bonded water molecules, while the shape of the distribution does not change qualitatively. On the other hand, an increase in Tt has no effect on the librational motion of H-bonded water molecules, thus leaving the equilibrium angles hOa O∗ H∗ and hOb O∗ H∗ almost intact. Bottom row in Fig. 4 presents the time evolution for the probability distribution of angle between the three water molecules involved in a H-bond switch, hOa O∗ Ob . The equilibrium distribution is very wide (from ≈30◦ to ≈120◦ ) and is peaked at ≈50◦ throughout the trajectory. Near the H-bond switch (t = 0) this probability distribution

Fig. 4. Probability distributions of the angle between the O∗ H∗ covalent bond and the OO vector to the initial H-bond acceptor, hOa O∗ H∗ (top row), angle between the O∗ H∗ covalent bond and the OO vector to the final H-bond acceptor, hOb O∗ H∗ (middle row), and of the angle between OO vectors to the initial and final H-bond acceptor, hOa O∗ Ob (bottom row) along the average jump path, where the system crosses the transition state at time t = 0. Different columns correspond to different thermostatting conditions: equilibrium at Tr = Tt = 300 K (left column), non–equilbrium Tr = 600 K (middle column) and non-equilibrium Tt = 600 K (right column). The red line plotted on top of each probability distribution represents the time evolution of the average value. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

6

T. Mohoriˇc and U. Bren / Journal of Molecular Liquids 302 (2020) 112522

exhibits a shoulder at ≈60◦ . The average value (red line) is about 70◦ and remains almost unchanged throughout the H-bond exchange. This indicates that the final H-bond acceptor (Ob ) approaches the central O∗ H∗ covalent bond radially in such a way, that the angle hOa O∗ Ob remains constant. Note also that the average hOa O∗ H∗ angle after the switch also approaches the hOa O∗ Ob angle. Expectedly, the orientation of the newly formed H-bond does not differ significantly from the orientation of the O∗ Ob vector. With an increase in Tr the angle hOa O∗ Ob is shifted towards larger values (about 85◦ ) – in agreement with the shift towards larger rOa Ob distances. This indicates that water molecules with a greater rotational kinetic energy are more likely to shift for a larger angle when switching the H-bonds. In nonequilibrium state with Tr = 600 K (middle column) the shape of distribution changes substantially. The peak at ≈50◦ diminishes and the distribution becomes rather uniform, except for t ≈ 0, where it exhibits a weaker and wider peak in the range of 40◦ to 80◦ . In line with previous observations a corresponding increase in Tt (right column) retains the shape of the distribution at least near the time origin. An increase in Tr has a greater influence on the hOa O∗ Ob than a corresponding increase in Tt – again, in agreement with the observed effects on rOa Ob . To continue our discussion about the angles related to the H-bond exchange, in Fig. 5 we present the time evolution of the probability distribution of the angle between the O∗ H∗ covalent bond and  the bisector of the Oa O∗ Ob angle, h (top row), as well as of the out of-plane angle between the O∗ H∗ covalent bond and the Oa O∗ Ob plane, 0 (bottom row). The latter distribution is localized around 0 = 0 with the average being equal to zero (red line) throughout the switch, confirming that the O∗ H∗ covalent bond fluctuates around the plane of the three oxygen atoms. Expectedly, thermostatting conditions exhibit no effect on the time dependence of 0. The former average angle, h, exhibits a large average amplitude jump (red line) during the switch: from ≈−33◦ before to ≈+33◦ after the

H-bond exchange – corresponding roughly to the hOa O∗ Ob angle at t = 0. Note also the small peaks in the time evolution of h immediately before and after the jump, indicating an oscillation of the covalent O∗ H∗ bond. Both, an increase in Tr (middle column) and in Tt (right column) magnify the angular jump in accordance with the previously seen increase in hOa O∗ Ob – the former effect being more pronounced than the latter. However, there is yet another subtle difference, between the jumps in rotationally and translationally excited water. Water molecules with higher rotational kinetic energy exhibit substantially more pronounced oscillations of the covalent O∗ H∗ bond immediately before and after the jump, while these oscillations tend to diminish in translationally excited water. The same effect can be observed in Fig. 4, where the amplitude of librational motion increases with Tr . 3.2. Reorientation times In this section we focus on the determination of the reorientation times of water molecules. To describe water reorientation kinetics we rely on the time-correlation function (TCF) of a vector describing the orientation of a water molecule: Cn (t) = Pn [u(0) • u(t)],

(1)

where Pn represents Legendre polynomial of rank n, and u(t) is the vector describing orientation of a water molecule at time t. Fig. 6 presents such TCFs for water molecule’s OH covalent bond vector, C2OH (black lines). It has been well established that for bulk water at room temperature C2OH exhibits a fast but limited decay on a subpicosecond timescale due to the fast librational motion – see the inset of Fig. 6. On a picosecond timescale this decay is exponential, as evidenced by a constant slope in Fig. 6. The reorientation of a water molecule on a picosecond timescale has been attributed

 Fig. 5. Probability distributions of the angle between the O∗ H∗ covalent bond and the bisector of the plane Oa O∗ Ob , h (top row), as well as of the out-of-plane angle between  the O∗ H∗ covalent bond and the plane Oa O∗ Ob , 0 (bottom row), along the average jump path, where the system crosses the transition state at time t = 0. Different columns correspond to different thermostatting conditions: equilibrium at Tr = Tt = 300 K (left column), non-equilbrium Tr = 600 K (middle column) and non-equilibrium Tt = 600 K (right column). The red line plotted on top of each probability distribution represents the time evolution of the average value. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

T. Mohoriˇc and U. Bren / Journal of Molecular Liquids 302 (2020) 112522

7

where uOH is the unit vector of the OH covalent bond in a laboratory frame, uOO is the unit vector of the oxygen–oxygen axis between the two H-bonded water molecules in a laboratory frame and u˜ OH is the unit vector of the OH covalent bond in the OO frame. As each of these TCFs is assumed to be exponential, the decay time tn of the Cn (t) can be given as: 1 1 1 = J + F, tn t tn n

(3)

J

where tn is the decay time associated with the jump mechanism and tnF with the frame reorientation. The latter decay time can be extracted from the COO n,nj , which are for different thermostatting conJ

ditions plotted in Fig. 7. The former decay time, tn , can be estimated using the Ivanov jump model, which takes into account the average geometry of the H-bond jump and the frequency of the jumps, 1/t0 : Fig. 6. Time-correlation functions of water’s OH covalent bond, C2OH (black lines), and cross correlation functions, Cab (red lines), at different thermostatting conditions: equilibrium at Tr = Tt = 300 K (solid lines), non-equilbrium Tr = 600 K (dashdotted lines) and non-equilibrium Tt = 600 K (dotted lines). Inset: Time-correlation functions of water’s OH covalent bond, C2OH , on a shorter time scale. The effects of variations in Tr (left) and Tt (right) are depicted with different colors. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

to two mechanisms: switches of H-bonds and frame reorientation (tumbling). The extended jump model(EJM), which was extensively described in [21, 22], defines contributions of each of these two mechanisms to the total TCF. EJM assumes independence of the jump and frame motions, and decomposes CnOH into the TCFs within a laboratory frame and within the local OO frame: CnOH (t) = Pn [uOH (0) • uOH (t)] = = Pn [u˜ OH (0) • u˜ OH (t)] × Pn [uOO (0) • uOO (t)],

(2)

1 J tn

=

1 t0

 1−

 sin[(n + 1/2)Dh] 1 , 2n + 1 sin(Dh/2)

(4)

where t0 represents the jump time between two consecutive jumps and Dh is the jump amplitude. While the latter parameter can be estimated from Fig. 5, the former can be obtained as a characteristic decay time of the cross correlation function for H-bond switching: Cab = na (0)nb (t).

(5)

Here na equals to 1 when H* is H-bonded to the oxygen atom Oa and to 0 otherwise, and nb equals to 1 when H* is H-bonded to the oxygen atom Ob and to 0 otherwise. In addition, absorbing conditions are assumed for nb so that when a new H-bond has formed, it remains intact. We thus discard any contribution from the back reaction or from the subsequent rupture of the O∗ Ob H-bond. Hence, the extracted time, t0 , is associated with the forward rate constant for

Fig. 7. Time-correlation functions of oxygen–oxygen vector of two H-bonded water molecules between the two consecutive H-bond switches, COO n,nj , for n = 1 (left column), n = 2 (middle column) and n = 3 (right column). The top row depicts the effects of variations in Tt , while the bottom row presents the results for different Tr .

8

T. Mohoriˇc and U. Bren / Journal of Molecular Liquids 302 (2020) 112522

Table 1 Average angular jump amplitudes Dh and characteristic decay times (in picoseconds): t0 is characteristic decay time of Cab cross correlation function (obtained from simulations), J tnF represents frame reorientation times obtained from COO n,nj time-correlation function (Fig. 7), tn represent jump reorientation times obtained from Ivanov jump model (Eq. EJM

(4) using Dh and t0 as determined in simulation), tn represent characteristic decay times of the CnOH time-correlation functions calculated from Eq. (3), and tnMD represent characteristic decay times of the CnOH time-correlation functions as determined from the simulations (Fig. 6, inset). J

J

J

EJM

EJM

EJM

Tt [K]

Tr [K]

Dh [◦ ]

t0

t1F

t2F

t3F

t1

t2

t3

t1

t2

t3

t1MD

t2MD

t3MD

300 300 300 300 300 400 500 600

300 400 500 600 300 300 300 300

70 75 81 86 70 74 76 77

2.5 1.7 1.1 0.7 2.5 1.0 0.6 0.5

12.9 10.6 6.7 5.2 12.9 4.6 2.9 1.8

4.7 3.3 2.2 1.8 4.7 1.5 0.9 0.6

2.4 1.7 1.2 0.9 2.4 0.8 0.5 0.3

5.7 3.4 2.0 1.1 5.7 2.1 1.2 1.0

2.6 1.6 1.0 0.6 2.6 1.0 0.6 0.5

2.0 1.4 0.9 0.6 2.0 0.8 0.5 0.4

4.0 2.6 1.5 0.9 4.0 1.4 0.8 0.6

1.7 1.1 0.7 0.4 1.7 0.6 0.3 0.2

1.1 0.8 0.5 0.4 1.1 0.4 0.2 0.2

4.2 2.5 1.5 0.8 4.2 1.6 1.0 0.7

2.3 1.5 0.9 0.5 2.3 0.8 0.5 0.3

1.5 1.1 0.6 0.5 1.5 0.5 0.3 0.2

the H-bond switching. As has been already observed in the past and confirmed in Fig. 6 (red lines), this correlation function exhibits a monoexponential decay with a decay time t0 varying depending on the thermostatting conditions. Generally, an increase in Tt reduces the characteristic decay time of Cab or C2OH more than a corresponding increase in Tr . We have calculated the characteristic decay times at different thermostatting conditions, in order to distinguish the effects of faster translational or rotational motion on the reorientation of bulk water. These results are collected in Table 1 and plotted in Fig. 8. Note that t0 represents the characteristic decay time from the cross correlation function Cab (Eq. (5) and red lines in Fig. 6), tnF are frame reorientation J times extracted from the time-correlation functions COO n,nj (Fig. 7), tn represent jump reorientation times computed using the Ivanov jump EJM model (Eq. (4), using Dh and t0 from simulations), tn are characterOH istic decay times of the Cn time-correlation function as calculated

from Eq. (3), and tnMD represent characteristic decay times of the CnOH time-correlation function determined from the simulation (Fig. 6, inset). As expected, the characteristic decay times tend to diminish with an increase in Tt or Tr , due to the faster dynamics at elevated temperatures. As deduced from the Fig. 8, the characteristic decay times computed from the MD simulations (i.e. from the CnOH directly) generally fall in between the predictions of Ivanov jump model and EJM. In most cases the latter more closely agrees with the results of the MD simulations than the former. Nevertheless, this difference is not big and in most cases it diminishes with an increase in Tt or Tr , which confirms that the EJM model may be used also under the studied non-equilbrium conditions. Note the difference in the scales of the top and bottom rows in Fig. 8. On one hand, we obtain a linear relationship between the inverse translational temperature T−1 t and log tn , while on the other

Fig. 8. Characteristic decay time tn of the time-correlation function for OH covalent bond reorientation, CnOH , as a function of inverse translational (top row) and rotational (bottom row) temperatures. Colors depict results obtained with the jump model of Ivanov (J; black triangles), extended jump model (EJM; red squares) and molecular dynamics (MD; blue circles). Different tn values are plotted in different columns: n = 1 (left column), n = 2 (middle column) and n = 3 (right column). Lines are only guides to the eye. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

T. Mohoriˇc and U. Bren / Journal of Molecular Liquids 302 (2020) 112522

hand the inverse rotational temperature T−1 is linearly correlated r with the tn . This also suggests a way to observe non-equilibrium effects of microwave heating in the experiment. For example: one could experimentally determine characteristic decay times tn (e.g. using the NMR technique) of a microwaved sample in a steady state regime with strictly controlled translational temperature (note that the conventionally measured temperature corresponds to the translational temperature). Changes in characteristic decay times could then be attributed to the rotational excitation of water molecules. It would be interesting to examine the relationship between the microwave power and tn . 4. Conclusions To summarize, we performed non-equilibrium molecular dynamics simulations of bulk water. We investigated the mechanism of Hbond exchange at increased rotational or translational temperature, the former situation mimicking the effects of microwave heating. We reproduced the main characteristics of H-bond switches that were reported in the past. We confirmed that the mechanism remains generally the same even at non-equilibrium conditions under study. We examined in detail not only the average path of a H-bond switch, but also the corresponding distributions. An increase in Tr resulted in more pronounced librational motion and in a wider distribution of jump angles for H-bond switches than a corresponding increase in Tt . However, while the effect of increasing Tr affects the geometry of H-bond switches more than a corresponding increase in Tt , the effect on reorientational dynamics is reversed. There an increase in Tt lowers the characteristic decay times of timecorrelation functions more than the corresponding increase in Tr . In this respect we also confirmed that the EJM describes the H-bond dynamics in bulk water well even at non-equilbrium conditions. Interestingly, we also found that the characteristic decay time is proportional to the inverse rotational temperature T−1 and to the r exp(1/Tt ). CRediT authorship contribution statement Tomaž Mohoriˇc: Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Writing - original draft, Visualization. Urban Bren: Funding acquisition, Project administration, Supervision, Writing - review & editing. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgments This study was supported by the Slovenian Research Agency (ARRS) funds through Program 0103–0201 as well as through Projects J1-6736 and J1-5448.

9

References [1] D. Eisenberg, W. Kauzmann, The Structure and Properties of Water, Oxford Clarendon Press, London, 1969. [2] B. Bagchi, Chem. Rev. 105 (2005) 3197. [3] I. Ohmine, H. Tanaka, Chem. Rev. 93 (1993) 2545. [4] C.J. Fecko, J.D. Eaves, J.J. Loparo, A. Tokmakoff, P.L. Geissler, Science 301 (2003) 1698. [5] R. Rey, K.B. Møller, J.T. Hynes, J. Phys. Chem. A 106 (2002) 11993. [6] K. Ando, J.T. Hynes, J. Mol. Liq. 64 (1995) 25. [7] K. Ando, J.T. Hynes, J. Phys. Chem. B 101 (1997) 10464. [8] B.J. Gertner, R.M. Whitnell, K.R. Wilson, J.T. Hynes, J. Am. Chem. Soc. 113 (1991) 74. [9] D. Marx, M.E. Tuckerman, J. Hutter, M. Parrinello, Nature 397 (1999) 601. [10] N. Agmon, Chem. Phys. Lett. 244 (1995) 456. [11] S. Cukieman, Biochim. Biophys. Acta 1757 (2006) 876. [12] R. Pomes, B. Roux, Biophys. J. 71 (1996) 19. [13] P. Ball, Sci, Proc. Natl. Acad. USA 114 (2017) 13327. [14] M.D. Fayer, Acc. Chem. Res. 45 (2012) 3. [15] D. Laage, T. Elsaesser, J.T. Hynes, Chem. Rev. 117 (2017) 10694. [16] M.K. Kong, S. Erramilli, J. Biol. Phys. 38 (2012) 1. [17] Y. Itoh, T. Aida, Nat. Chem. 9 (2017) 934. [18] M.-Cl. Bellissent-Funel, A. Hassanali, M. Havenith, R. Henchman, P. Pohl, F. Sterpone, D. van der Spoel, Y. Xu, A.E. Garcia, Chem. Rev. 116 (2016) 7673. [19] Y. Levy, J.N. Onuchic, Annu. Rev. Biophys. Biomol. Struct. 35 (2006) 389. [20] M. Chaplin, Nat. Rev. Mol. Cell. Bio. 7 (2006) 861. [21] D. Laage, J.T. Hynes, Science 311 (2006) 832. [22] D. Laage, J.T. Hynes, J. Phys. Chem. B 112 (2008) 14230. [23] D. Laage, G. Stirnemann, F. Sterpone, R. Rey, J.T. Hynes, Annu. Rev. Phys. Chem. 62 (2011) 395. [24] H.M. Kingston, J.S. Haswell, Microwave-Enhanced Chemistry, American Chemical Society, Washington DC, 1997. [25] G. Keglevich, Milestones in Microwave Chemistry, Springer International Publishing, Switzerland, 2016. [26] V. Polshettiwar, R.S. Varma, Aqueous Microwave Assisted Chemistry Royal Society of Chemistry, 2010. [27] N.J. English, J.M.D. MacElroy, J. Chem. Phys. 118 (2003) 1589. [28] M. Tanaka, M. Sato, J. Chem. Phys. 034509 (2007) 126. [29] W.Y. Tian, K.M. Huang, L.J. Yang, Z.N. Guo, F.H. Liu, Chem. Phys. Letters 607 (2014) 15. [30] J. Liu, G. Jia, J. Mol. Liq. 227 (2017) 31. [31] N.J. English, D.A. Mooney, J. Chem. Phys. 091105 (2007) 126. [32] N.J. English, G.Y. Solomontsev, P. O’Brien, J. Chem. Phys. 131 (2009) 035106. [33] G.Y. Solomentsev, N.J. English, D.A. Mooney, J. Chem. Phys. 133 (2010) 235102. [34] N.J. English, J.M.D. MacElroy, J. Chem. Phys. 120 (2004) 10247. [35] J.A. Garate, N.J. English, J.M.D. MacElroy, J. Chem. Phys. 131 (2009) 114508. [36] U. Bren, A. Kržan, J. Mavri, J. Phys. Chem. A 112 (2008) 166. [37] M. Bren, D. Janežiˇc, U. Bren, J. Phys. Chem. A 114 (2010) 4197. [38] U. Bren, D. Janežiˇc, J. Chem Phys. 137 (2012) 024108. [39] C. Blanco, S.M. Auerbach, J. Am. Chem. Soc. 124 (2002) 6250. [40] C. Blanco, S.M. Auerbach, J. Phys. Chem. B 107 (2003) 2490. [41] H. Jobic, J.E. Santander, J.C. Conner, G. Wittaker, G. Giriat, A. Harrison, J. Ollivier, S.M. Auerbach, Phys. Rev. Lett. 106 (2011) 157401. [42] W.F. Edwards, D.D. Young, A. Deiters, Org. Biomol. Chem. 7 (2009) 2506. [43] D.D. Young, J. Nichols, R.M. Kelly, A. Deiters, J. Am. Chem. Soc. 130 (2008) 10048. [44] K. Huang, X. Yang, W. Hua, G. Jia, L. Yang, New J. Chem. 33 (2009) 1486. [45] T. Urbic, K.A. Dill, J. Am. Chem. Soc 140 (2018) 17106. [46] H.J.C. Berendsen, J.R. Grigera, T.P. Straatsma, J. Phys. Chem. 91 (1987) 6269. [47] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Clarendon Press. 1989. [48] C.J. Fennell, J.D. Gezelter, J. Chem. Phys 124 (2006) 234104. [49] T. Mohoriˇc, B. Hribar-Lee, V. Vlachy, J. Chem. Phys. 140 (2014) 184510. [50] T. Mohoriˇc, B. Hribar-Lee, V. Vlachy, Cond. Matt. Phys. 18 (13004) (2015) 19. [51] T. Mohoriˇc, U. Bren, V. Vlachy, Acta Chim. Slov. 62 (2015) 489. [52] T. Mohoriˇc, U. Bren, V. Vlachy, J. Chem. Phys. 244510 (2016) 143. [53] T. Mohoriˇc, U. Bren, J. Chem. Phys. 146 (2017) 044504. [54] T. Urbic, T. Mohoric, J. Chem. Phys. 146 (2017) 094505. [55] T. Mohoriˇc, U. Bren, J. Mol. Liq. 266 (2018) 218. [56] L.V. Woodcock, Chem. Phys. Lett. 10 (1971) 257. [57] J. Petersen, K.B. Moller, R. Rey, J.T. Hynes, J. Phys. Chem. B 117 (2013) 4541.