Effects of translational and rotational degrees of freedom on the properties of model water

Molecular dynamics simulations with separate thermostats for rotational and translational motions were used to study the effects of these degrees of freedom on the structure of water at a fixed density. To describe water molecules, we used the SPC/E model. The results indicate that an increase of the rotational temperature, $T_\textrm{R}$, causes a significant breaking of the hydrogen bonds. This is not the case, at least not to such an extent, when the translational temperature, $T_\textrm{T}$, is raised. The probability of finding an empty spherical cavity (no water molecule present) of a given size, strongly decreases with an increase of $T_\textrm{R}$, but this only marginally affects the free energy of the hydrophobe insertion. The excess internal energy increases proportionally with an increase of $T_\textrm{R}$, while an increase of $T_\textrm{T}$ yields a much smaller effect at high temperatures. The diffusion coefficient of water exhibits a non-monotonous behaviour with an increase of the rotational temperature.


Introduction
We dedicate this contribution to J.A. Barker and D. Henderson [1], authors of the important paper What is "liquid"? Understanding the states of matter, which was published in Reviews of Modern Physics, almost exactly 50 years ago. This excellent paper, which served as a textbook to many of us, is still a reference article for every person who wishes to study the theory of liquids. Our contribution here touches upon one particular, but extremely important liquid -water, from the standpoint of the individual degrees of freedom.
Complex molecules have, besides the translational, also rotational and vibrational degrees of freedom. They all make their own contribution to thermodynamic functions, which depend on temperature of the system. Temperature is accordingly considered to be the most important thermodynamic parameter. An interesting question is how the fundamental degrees of freedom contribute to the potential of mean force between water molecules, when their characteristic temperatures vary independently.
A biologically important situation arises when water molecules are irradiated with microwaves of appropriate frequency to excite their rotational motion, which affects the solvation and interaction between the solute molecules [2][3][4]. The interaction of microwaves with water has been previously studied using the non-equilibrium molecular dynamics simulation [5], where the electric field was modelled explicitly. It has been shown that microwaves induce heating by excitation of rotational motion in water molecules. The excess rotational kinetic energy is then transferred to translational degrees of freedom. Following previous studies [2][3][4], the direct interaction of microwaves with a solvent is replaced by a model system in which the rotational degree of freedom has a different temperature than the translational motion. Such a non-equilibrium steady-state picture, where different degrees of freedom can have different temperatures, has already been used to explain experimental phenomena, such as increased chemical reaction rates and altered polar solvent boiling points [2,6].
In the preceding contribution [4], we studied the effects in hydration of cations, anions and hydrophobes caused by independent variations of the rotational and translational temperatures. The main conclusions were: (i) an increase of the rotational temperature affects the hydration of cations in an opposite way to anions, (ii) an increase of the translational temperature always decreases the height of the first peak in the solute-water radial distribution function; (iii) an increase of the rotational temperature yields an increase in the first peak in the solute-water radial distribution function for hydrophobes and cations; while, (iv) in contrast to this, the solvation peak decreases around ions with a sufficiently large negative charge.
The focus of the preceding work was, therefore, on solvation of various solutes and not much attention was paid to the solvent itself. With this respect, the present work complements the previous one; here, we are primarily interested in the structure of the water under such non-equilibrium conditions. In order to explore the properties of the model water, we used molecular dynamics (MD) simulations, where the rotational and translational motions were coupled to separate thermostats. We used a rigid model of water (SPC/E) [7]; the effect of vibrational temperature variations was found to be small [3]. The focus of the present study is on the structure of water as reflected in various site-site distribution functions. The probability of water molecule to form a certain number of hydrogen bonds is calculated as a function of rotational and translational temperatures. The effects of the separate T R and T T variations on the excess internal energy of the system and on the dynamics of water molecules is examined.

Model and simulations
Molecular dynamics simulations were performed at a constant number of particles and volume (N ,V ) with 256 SPC/E water molecules in the simulation box. As before [4], the number density of water molecules was set to 1.0 g/cm 3 , and the time step used was 1.0 fs. Each MD simulation had at least 100 ps long MD equilibration run, while the statistics was collected during the 4900 ps long production run.
The translational motion of a single water molecule i was described by its center-of-mass position r i , velocity v i and the force f i acting on it. Similarly, the rotational motion of a single water molecule i was characterized by its orientation (using quaternion representation) q i , angular velocity ω i and torque τ i .
While the acceleration of the center-of-mass follows immediately from the force acting on it, the angular acceleration must be computed from Euler's equations of motion for rigid bodies. We employed the standard Verlet algorithm to integrate the equations of motion [8]. To hold T R and T T fixed at their respective values, we employed a simple velocity re-scaling scheme.

Radial distribution functions
We start the presentation of numerical results with the oxygen-oxygen, g OO , and oxygen-hydrogen distribution functions, g OH . First, in figure 1 we present the g OO distributions for various rotational (T R ) and translational (T T ) temperatures of the system.
In both cases, an increase of the temperature decreases the height of the first peak of the oxygenoxygen distribution function. However, there is a difference regarding how the position of the peak responds to the temperature variations. An increase of T R causes the first peak to move toward larger distances. In contrast to this, the position of the first peak in g OO does not change upon an increase of the translational temperature, T T . In addition, the peak becomes somewhat broader; the shape seems to reflect the fact that molecules with larger kinetic energy can easier "penetrate" into each other.
Liquid water is at ambient conditions characterized by the position of the second peak in this figure (figure 1). The peak is located at 4.5 Å rather than at 2σ W W ∼ 6 Å as it is at the absence of hydrogen bonds. The position of this peak is a consequence of the hydrogen-bonding of water molecules. The shift of the second peak toward ∼ 6 Å suggests that some hydrogen bonds break upon an increase of rotational temperature. The effects of temperature variations on the hydrogen bonds in the model water are most clearly shown in figure 2. Here, the oxygen-hydrogen distribution, g OH , is presented as a function of the respective temperature. The position of the first peak of g OH distribution is characteristic of a hydrogen bond. We see that the height of this peak is very sensitive to T R and much less to the T T variations. In both cases, however, the peak decreases in magnitude. From this distribution, the coordination number, n H , is calculated (for the definition see caption in figure 3). This quantity, as shown in figure 3, decreases almost linearly with T R . On the other hand, n H is much less sensitive to the T T variations at higher temperatures. In figure 4, we show the probability P (N HB ) that a water molecule forms N HB hydrogen bonds. Up to 500 K, the probability distribution, P (N HB ), is practically the same for the T R and T T dependence. Above T R =500 K, the peak of this distribution is gradually shifted toward smaller N HB values and there is an appreciable probability for only one hydrogen bond. The situation is different for an increase of T T where the P (N HB ) distribution does not change much upon an increase of the translational temperature above T T =500 K.
All the evidence gathered so far suggests that hydrogen-bonding is more sensitive to T R , than to T T variations. The reason for this has been discussed in reference [3] (page 024108-6). The formation of a hydrogen bond requires two water molecules to meet at an appropriate distance and a correct orientation. It is easy to understand that every increase in T R (faster rotation) will make favorable orientations less probable. At conditions of high T R , the model molecules will behave as independent rotators with the hydrogen bond interaction energy approaching zero.
The mechanism causing two molecules of water to lose the correlation upon the T T increase is different. In the latter case, the molecules should separate for a certain distance from each other, which is not that easy because each molecule is caged by the others. By increasing T T , having in mind that we perform simulations at a constant water density, we do not approach the ideal gas behaviour, even at high temperatures (T T =900 K) the system is under the effect of strong repulsive interaction. An initial increase of T T (up to T T = 500 K) indicates that some hydrogen bonds are broken and the radial order is disrupted to a certain extent.
A more detailed information on the water structure can be obtained if the radial distribution function, g OO , is divided into the separate contributions of the neighboring molecules coordinated on the chosen molecule. In our notation, g i is the contribution of the i -th neighbor to the total pair distribution function.
These results are shown in figure 5, where the left-hand panels belong to different rotational and the right-hand panels belong to different translational temperatures. In case of T R increase, we observe the shifts of the peaks toward larger distances. This applies to all four nearest neighbors. At this density water molecules are closely packed; while σ W W for this (SPC/E) model is 3.169 Å, the position of the first peak is actually at around 2.7 Å. Rising the T R causes for this distance to increase, and the molecules distribute at larger average distances. In contrast to this, an increase of the T T does not significantly change the position of the first peak. This is yet another consequence of the fact that faster rotation (increase in T R ) disrupts the hydrogen bond network more efficiently than an increase in T T . An interesting question, related to the changes of peak positions of various g i 's upon T R increase, has been posed by an anonymous reviewer. Figure 5 indicates that the position of the first peak, g 1 , is less affected than the positions of more distant peaks, such as g 3 and g 4 . Why? Considering that g 3 measures the contribution of the third neighbour to the oxygen-oxygen distribution function, we may conclude that correlation between the central molecule and the third shell molecules is significantly suppressed.

13004-4
Effects of translational and rotational freedom In other words, the effect of the central molecule does not propagate that far in the bulk as for T R = T T = 300 K. This is clearly shown in figure 1; the structure of the solution as measured by oxygen-oxygen pair distribution function, g OO , becomes more uniform upon an increase of either T R or T T . Furthermore, oscillations are not just smaller at elevated temperature, they are also "out of phase".

Probability of observing an empty cavity of a size of water
The simulations presented in this paper are performed for a constant number of particles in a given volume. An important information on a liquid system is associated with the density fluctuations. In this subsection we studied the probability distribution P (N ) for observing exactly N water molecules (actually their centers of gravity) in a spherical volume having diameter 3 Å. This probability distribution is approximately Gaussian [9] for small (molecular size) volumes. The width of the distribution reflects compressibility of the liquid: a wider distribution means a somewhat higher compressibility. Our simulations indicate that P (N ) curves change only marginally upon an increase of T R and T T . A notable exception is the P (0) value -i.e., the probability of observing an empty spherical cavity. Accordingly, only the dependencies of P (0) on T R and T T are shown in figure 6. An increase of T R from 300 to 900 K causes a decrease for about one order of magnitude for P (0) (note the natural logarithm plotted on y-axis in figure 6). The density fluctuations are suppressed due to a decreased number of hydrogen bonds. On the other hand, the variations of T T have little effect on P (0). Complementary to these calculations, we also performed the thermodynamic integration to simulate the free energy of a hydrophobe insertion, ∆F ins , into the model water (see figure 7). Hydrophobe was represented by the Lennard-Jones particle with the well depth equal to that of the model water molecule and σ SS = 3 Å. Lorentz-Berthelot mixing rule was used for cross interaction parameters. Under conditions where the hydrophobe-water coupling is weak, we expect for the ∆F ins to have qualitatively the same features as the same quantity for the hard-sphere solute.
While the P (0) values decrease with an increase of T R , the insertion free energy is only slightly affected under such conditions. On the other hand, T T , which has virtually no effect on the P (N ) distribution (see figure 6), has a relatively large effect on the ∆F ins . In both cases, the free energy for the particle insertion, ∆F ins , increases; it becomes increasingly more difficult to add a hydrophobe to the system.
Here, we offer a qualitative explanation for this result. The insertion free energy of the hard-sphere solute is related to the probability of observing an empty spherical cavity, P (0), in solvent [9] (as usual k B is the Boltzmann constant): As the solute-water interaction depends only on the positions of water molecules (and not on their orientations) and since equation (3.1) applies to any solvent (not necessarily possessing rotational degrees of freedom), we might expect that temperature in equation (3.1) refers to the translational temperature.
While the P (0) values vary substantially with T R , the ln P (0) only increases for about 20 % in magnitude, which is close to the relative change of ∆F ins with T R . Much stronger increase of ∆F ins is observed upon the rise of T T and it may be attributed to an increase of the pre-factor k B T .

13004-6
Effects of translational and rotational freedom

Excess internal energy
It is of interest to examine the variations of the excess internal energy (E ex ) per mol of water molecules as functions of T R or T T (figure 8). We can see that elevation of T R causes an almost linear increase of E ex , while the same quantity changes less and less if T T is increased. The reason for such a behaviour lies, as explained before, in different effects of the two degrees of freedom to hydrogen bonding. These graphs also provide some insights into the T R and T T effects on the excess heat capacities (at constant N and V ), C ex v,r and C ex v,t . As we see, the rotational contribution (derivative of the upper curve) is temperature independent, while for the translational contribution (derivative of the lower curve) the excess heat capacity goes to zero for high T T values.

Dynamics of model water molecule upon an increase of T R and T T
How the different degrees of freedom affect the dynamics of water molecules? We can see the answer in figure 9, showing the effect of T R and T T on the (self) diffusion coefficient of a model water molecule.
At first, an increase of T R causes a faster diffusion, which can be attributed to the breaking of hydrogen bonds. As the water molecules are more weakly bound, they tend to diffuse faster. However, an increase of T R also means a decrease of P (0). Therefore, a random walk should proceed in smaller steps, causing the diffusion to slow down. This effect seems to dominate at the highest T R studied here. On the other hand, increasing T T , as expected, causes a monotonous increase of the diffusion coefficient in this temperature interval.

Conclusions
It is well known that microwaves induce heating by excitation of rotational motion in water molecules and that the rotational kinetic energy is then transferred to the translational degrees of freedom. The microwave irradiation is important in biology: there are reports that microwaves enhance the folding and unfolding kinetics of globular proteins [10][11][12], as well as protein aggregation [13].
The effect of microwaves on the properties of liquids and solutions can be studied by non-equilibrium molecular dynamics simulations [5]. In such studies, the electric field is modelled explicitly. Following the ideas, proposed in references [2,3], we perform calculations in which the rotational degree of freedom has the temperature different from the translational one. In a previous paper [4], we examined the effect of such non-equilibrium conditions on the solvation of simple solutes. In the present contribution, we studied the properties of a model water under non-equilibrium conditions where the translational and rotational temperatures vary independently. The main conclusions of this work are: (i) an increase of the rotational temperature causes the first peak of g OO to move toward larger distances, while its position remains unchanged upon T T increase. (ii) The height of the first peak in the hydrogen-oxygen distribution function is much more sensitive on the T R than on the T T variations. (iii) At high rotational temperatures, the number of hydrogen bonds may drop as low as to one. (iv) Free energy of the hydrophobe insertion does only marginally depend on the rotational temperature. (v) Water diffusion coefficient measured via the mean square displacement is, in contrast to the T T dependence, a non-monotonous function of the rotational temperature. These and other effects observed upon varying the rotational temperature may significantly affect the properties of water as solvent, which has to a certain degree been examined in reference [4]. The work exploring such effects on polymer and polyelectrolyte conformations in aqueous solution is currently underway -preliminary results clearly point toward a high probability of polymer collapse in cases of rotational heating.
In the end, we shall conclude that, despite important developments following reference [1], there is still ample room for improvements in theory and simulation of liquids, of which, water is by far most important.