A Simulation protocol

Dynamical and orientational structural crossovers in low-temperature glycerol


Mean square displacements of hydrogen atoms in glass-forming materials and proteins, as reported by incoherent elastic neutron scattering, show kinks in their temperature dependence. This crossover, known as the dynamical transition, connects two approximately linear regimes. It is often assigned to the dynamical freezing of subsets of molecular modes at the point of equality between their corresponding relaxation times and the instrumental observation window. The origin of the dynamical transition in glass-forming glycerol is studied here by extensive molecular dynamics simulations. We find the dynamical transition to occur for both the center of mass translations and the molecular rotations at the same temperature, insensitive to changes of the observation window. Both the translational and rotational dynamics of glycerol show a dynamic crossover from the structural to a secondary relaxation at the temperature of the dynamical transition. A significant and discontinuous increase in the orientational Kirkwood factor and in the dielectric constant is observed in the same range of temperatures. We, however, do not find any indications of a true thermodynamic transition to an ordered low-temperature phase. We therefore suggest that all observed crossovers are dynamic in character. The increase in the dielectric constant is related to the dynamic freezing of dipolar domains on the time-scale of simulations.

87.14.E-, 87.15.N-, 87.15.Pc

I Introduction

Displacements of atoms and molecules induced by thermal agitation generally increase with temperature. A linear growth of the mean-squared displacement (MSD) with increasing temperature is predicted by the Nyquist (fluctuation-dissipation) theorem Kubo (1966); Hansen and McDonald (2003). The MSD is experimentally extracted from either the intermediate scattering function of the neutron scattering experiment Gabel et al. (2002) or from the fraction of recoilless -ray emission of the nucleus in the Mössbauer experiment Parak (2003); Young et al. (2011). The Nyquist theorem was found to be violated for a number of glass-forming materials, where a kink in the MSD vs. temperature is often observed at the laboratory glass transition Angell (1995). More complex behavior, with several kinks Roh et al. (2005); Khodadadi et al. (2010); Hong et al. (2011), was observed for proteins in partially hydrated powders or in the polycrystalline form Doster et al. (1989); Zaccai (2000).

A typical temperature dependence of the protein MSD starts with the linear increase in accord with the Nyquist theorem and the corresponding vibrational density of states Zaccai (2000); Achterhold et al. (2002). It is followed by one or two low-temperature crossovers and, finally, with a much stronger increase above the temperature of the dynamical transition K Doster (2008). This latter temperature depends on a number of factors, including the resolution of the spectrometer, i.e., effectively the time period over which the atomic displacements are recorded Magazù et al. (2011); Fenimore et al. (2013). This phenomenology has attracted significant attention since enhanced flexibility and, therefore, the ability to perform biological function can develop at Schirò et al. (2015).

A somewhat unexpected observation came recently from Capaccioli et al Capaccioli et al. (2012), who presented two key observations based on the analysis of a large database of neutron scattering data accumulated so far: (i) the MSD measured in 50:50 lysozyme-glycerol mixture can be nearly seamlessly overlaid with corresponding measurements for the pure glycerol and (ii) there are two crossover temperatures common to lysozyme-glycerol and glycerol systems, at and 276 K.

The first observation is significant for assigning the modes of the protein-solvent system responsible for the protein’s extended flexibility at high temperatures. High protein flexibility is required for its biological action Smith (1991); Zaccai (2003); Parak (2003), and this perspective connects protein function with specific physical modes and fluctuations of the protein-solvent system Matyushov (2015). Frauenfelder and co-workers suggested that the solvent mode coupled to the protein atomic displacements has to be attributed to the hydration shell Fenimore et al. (2004); Lubchenko et al. (2005). They also noted that this mode is decoupled from the -relaxation of the bulk solvent (structural or collective relaxation with the longest relaxation time and usually connected to the liquid viscosity). The relaxation time of the hydration shell is both faster than -relaxation and is Arrhenius, with the activation energy usually smaller than that of -relaxation. Taken together, these features point to its -character in the established classification of glass science Angell et al. (2000); Donth (2001). Since secondary -relaxation processes exist also in the bulk solvent, the fluctuations localized in the hydration shell of the protein are classified as -relaxation and are expected to carry the dynamics distinct from the bulk Chen et al. (2008). The dynamical transition then occurs when the -relaxation of the hydration shell slows sufficiently down, with lowering temperature, to become longer than the instrumental time-scale (dynamical freezing) Khodadadi et al. (2008); Frauenfelder et al. (2009).

The observation of a near-equivalence of MSDs recorded by neutron scattering in lysozyme-glycerol and pure glycerol systems puts under question the hydration-shell hypothesis, or at least the part of it attributing -relaxation specifically to the shell, in contrast to a faster relaxation mode of the bulk (of presumably -character). The question posed by this observation is whether the modes of the solvent coupled with protein flexibility are hydration-shell specific or generic to the bulk material. Furthermore, since the dynamical transition is a general phenomenon common to glass-forming materials, including molecular liquids and biopolymers Angell (1995), the question here is what are the modes that experience dynamical freezing at and whether the instrumental resolution must necessarily be a part of the explanation. Addressing some of these mechanistic questions is a goal of this study.

In order to avoid the complexities of protein solutions, we address these basic questions by focusing solely on bulk glycerol, for which we report here extensive molecular dynamics (MD) simulations. The temperature dependence of hydrogen MSDs is analyzed in terms of separate contributions of the center of mass translations and rotations relative to the center of mass of the molecule. Both translational and rotational MSDs show a crossover at the same temperature K consistent with experimental data. The temperature of translational and rotational dynamical transitions does not change when the observation time is significantly altered. We also find that the same temperature characterizes the dynamic crossover from to relaxation as measured by glycerol’s diffusivity and rotational dynamics.

The consistent picture arising from our observations is that a structural crossover occurs in glycerol at K, which affects both the MSDs and relaxation times. However, there is no indication from our data that this crossover should be identified with a true thermodynamic transition. We therefore suggest that all observed crossovers are dynamical in character. In particular, the structural crossover to a low-temperature state of glycerol, characterized by long-ranged dipolar correlations, becomes possible because these collective correlations cannot relax on the limited observation time. The dynamical transition in the MSD recorded by neutron scattering is not the result of crossing of the time-scale of single-particle translational/rotational diffusion with the observation time-scale, but rather the crossing of the latter with the time-scale of multi-body relaxation of polarized domains. A corresponding significant increase in the orientational Kirkwood factor and the jump in the dielectric constant at low temperatures are caused, in our simulations, by the crossing of the relaxation time of dipolar domains and the observation (simulation) time. This phenomenology is similar to that of relaxor ferroelectrics where dynamic freezing of ferroelectric domains is responsible for the high dielectric constant of the low-temperature phase Samara (2003).

Ii Incoherent neutron scattering

The experimental MSDs are extracted from incoherent elastic neutron scattering. The reported signals are affected by the instrumental resolution function convoluting with the dynamic structure factor , for which we assume the scattering momentum directed along the -axis of the laboratory frame. The function is the time Fourier transform of the intermediate scattering function


where is the displacement of a hydrogen atom and the sum runs over hydrogen atoms in the system; denotes an ensemble average.

In what follows we will consider all hydrogens in the system identical, although we will separate two groups of hydrogens of glycerol: 3 hydroxyl hydrogens and 5 hydrogens bonded to carbon atoms. Correspondingly, experimental results for partially deuterated glycerol Wuttke et al. (1995) (g-d3) and (g-d5) will be analyzed by considering the corresponding groups of hydrogen atoms not substituted by deuteration.

The intensity of the elastic scattering function at gives access to the MSD Gabel et al. (2002); Wuttke et al. (1995). The corresponding function , depending on the resolution window of the spectrometer , can be approximated by , where the resolution time is related to the resolution window of the spectrometer. According to Doster et al Doster et al. (2003), the connection is , where is the width at half maximum of the resolution function.

The intermediate scattering function in Eq. (1) can be estimated in the Gaussian approximation Yi et al. (2012), which leads to


If the time autocorrelation function , decays sufficiently to zero on the resolution time , the second term in Eq. (2) disappears and one gets an estimate of the mean square fluctuation (MSF) from the linear slope of vs Wuttke et al. (1995); Schirò et al. (2012). Otherwise one obtains half of the MSD from the slope of vs .

Figure 1: for g-d5 (upper panel) and g-d3 (lower panel). The experimental data obtained from IN13 spectrometer for correspondingly deuterated glycerol Capaccioli et al. (2012) are compared to MD simulations. The simulated MSDs are separated into displacement of the glycerol center of mass (“Trans”) and the displacements of hydrogens relative to the center of mass (“Rot”). The dashed lines are the linear regressions drawn through the corresponding points from MD simulations.

The data presented here were obtained from extensive MD simulations of glycerol described by the OPLS-AA force field Caleman et al. (2012a) as explained in Appendix A. Our main purpose in the analysis of the intermediate scattering function is to extract the relative contributions to the observed MSD arising from center of mass translations and molecular rotations relative to the center of mass. The question that we address here is whether the dynamical transition, if observed, occurs at the same temperature for these two modes. In addition to general mechanistic insights that such an analysis can produce, this question is relevant to testing the idea of dynamical freezing of a subset of molecular motions as the reason for the experimentally observed kink in the dependence of the MSD on temperature Doster et al. (1989); Doster (2008); Zaccai (2000); Magazù et al. (2011), identified with . If the kink is caused by reaching the equality between the relaxation time and the instrumental observation window Magazù et al. (2011), the dynamical transition temperature should be different for translations and rotations having their distinct relaxation times, unless they happen to be close. This is not what we observe from our simulations: the dynamical transition temperatures are the same for rotations and translations when calculated from fitting the intermediate scattering function to Eq. (2) (Fig. 1).

The separation of the center of mass translations and rotations relative to the center of mass assumes the factorization of the intermediate scattering function into the translational, , and rotational, , components


We therefore calculated and separately and produced the linear fits of the corresponding functions vs with ps for both g-d3 and g-d5 liquids. No deuteration was actually performed in simulations and only the corresponding groups of hydrogen atoms were selected to produce the intermediate scattering functions.

The accuracy of translation/rotation factorization in Eq. (3) was tested previouly and is usually found to hold Teixeira et al. (1985); Chen et al. (1997); Liu et al. (2002). Indeed, one expects this separation to be accurate in the Gaussian limit since translations and rotations carry different symmetry. If one separates into the center of mass displacement and the rotation relative to the center of mass , the MSD becomes the sum of two self terms and the translational-rotational cross term


Figure 2 shows an example of the analysis of the three correlation components in Eq. (4) from MD simulations. The translational and rotational components of the MSD are close in magnitude, while the cross-correlation is negative and is much smaller.

The translational and rotational MSDs are shown separately in Fig. 1 to indicate the common point of the kink at K. The same temperature of the dynamical transition is reported experimentally Wuttke et al. (1995); Capaccioli et al. (2012). However, the absolute values of MSDs from experiment (closed diamonds in Fig. 1) are below the simulation results, which is easy to see from the plot since the overall MSD follows from adding up the translational and rotational components (Eq. (4)). The most probable explanation of this discrepancy is that fitting the experimental neutron scattering data in a limited range -values used in the measurements Wuttke et al. (1995) allows one to probe only a limited subset of motions Chen et al. (1999); Liu et al. (2005), presumably the translational diffusion. Indeed, the agreement between simulations and experiment for the center of mass MSD is quite good. We also note that the agreement between the calculated coefficient of self-diffusion of glycerol and the results of measurements by NMR Chen et al. (2006a) is also reasonable (Fig. 6 below).

Figure 2: vs time at K. The overall MSD (long-dashed line) is separated into the center of mass (solid line), rotational (dash-dotted line), and cross (dashed line) components (Eq. (4)).
Figure 3: Center of mass MSD, , of glycerol at different temperatures indicated in the plot.

The time dependence of MSDs shown in Fig. 2 also helps to understand the physical origin of MSDs recorded by neutron scattering. Both the translational and rotational components of the MSDs are characterized by two distinct regimes: a fast ( ps) growth due to ballistic motions in the liquid’s cage (localized diffusion Mamontov (2012)), followed by a much slower, long-range diffusion with . The main observation here is that most of the MSD on the resolution time-scale ps is caused by the ballistic displacement associated with a secondary relaxation and not by the diffusional motion associated with the primary relaxation process. This conclusion holds both below and above (Figs. 2 and 3). The increase of the observation window from 25 ps to 135 ps makes the time spent by the particle on the linear, diffusion portion of the MSD longer (Fig. 2) and thus increases the slope of the high temperature part of the MSD curve (Fig. 4). It is important to realize that fast cage dynamics, resulting to the main portion of the observed MSD, are much faster than the resolution time and in fact become even faster with lowering temperature because of a greater rigidity of the low-temperature glycerol. It is the amplitude of the ballistic displacement which gets larger with increasing temperature, resulting in the observed temperature dependence of the MSD. The crossing of the resolution time of the spectrometer (25 ps) and the relaxation time of these ballistic motions never occurs and, therefore, the kink in the MSD vs temperature cannot be attributed to the finite resolution time.

Figure 4: for g-d5 measured on ps (open points) and ps (closed points). The center of mass (“T”, squares) and rotational (”R”, circles) contributions are shown separately. The dashed and dash-dotted lines are linear regressions drawn through the low-temperature and high-temperature points.

The change of the form of the MSD vs with the changing observation window is shown in Fig. 4. It adds additional evidence to the suggestion that the kink in the MSD’s temperature dependence is not caused by the equality between the relaxation time and the observation window. While the high-temperature portion of the MSD has a steeper slope for a higher , in agreement with experiment Capaccioli et al. (2012), the temperature of the dynamical transition has little sensitivity to . In addition, the equality between the dynamical transition temperatures for the translational and rotational MSDs is preserved between ps and ps. If one assumes that the consistency in for ps shown in Fig. 1 is a mere coincidence, it is hard to see how it can be preserved at ps. One has to accept the conclusion that the kink in the MSD is not related to the observation window Schirò et al. (2012); Fomina et al. (2014) and, instead, should be attributed to the softening of the liquid cage, with increasing temperature, in which a glycerol molecule finds itself for a relatively short time of ps. The rattling inside the cage is followed by an escape and the onset of long-range diffusion, but this component simply adds to the main displacement achieved by the ballistic cage rattling. The next question is whether structural distinctions of the entire liquid producing the difference between the low-temperature rigid cage and the high-temperature soft cage can be identified.

Iii Dynamic crossover

An explanation alternative to the instrumental resolution effect for the appearance of the kink in the proton MSD involves the dynamic crossover, i.e., a corresponding kink in the dependence of the system relaxation time on the inverse temperature Chen et al. (2006b). This phenomenology, known as the fragile-to-strong transition in glass science Angell et al. (2000), represents the crossover from the structural -relaxation at high temperatures above the crossover to a secondary -relaxation at low temperature below the crossover. Correspondingly, the activation barrier of the high-temperature -relaxation is higher than the activation barrier of the low-temperature -relaxation. We show below that this phenomenon is not connected to the kink in the MSD reported by neutron scattering and, at least in our simulations, has a trivial explanation of slower dynamics exceeding in its relaxation time the observation window (simulation time in the case of MD).

The problem of dynamic crossover in confined water has been extensively studied Zanotti et al. (2005); Liu et al. (2005); Cupane et al. (2014) and it has been established that the temperature of the dynamic crossover of confined water is generally consistent with of proteins Chen et al. (2006b); Schirò et al. (2013). The temperature was also found to be independent of the protein hydration level Schirò et al. (2013); Fomina et al. (2014); Mallamace et al. (2015) even though the relaxation times themselves are strongly affected by hydration. This latter observation points to the connection between and some sort of structural change in confined water.

The dynamic crossover results for water are necessarily limited to confined systems since bulk water is unstable to nucleation below K Mallamace et al. (2006); Cupane et al. (2014). Since our present simulations apply to bulk glycerol, it would be of significant interest to establish a phenomenology similar to that found for confined water for a material available in bulk phase both in simulations and in the laboratory experiment.

Figure 5: Average relaxation time obtained from rotational correlation function (“rot”) and from the electric field correlation function (“field”). The solid line refers to the average relaxation time Richert (2014) obtained by fitting the dielectric loss spectrum to the Cole-Davidson function Menon et al. (1992). The dashed line is a regression drawn through the MD points obtained from the electric field correlation function. (dotted line) indicates the crossover temperature.

It is useful to start off with an estimate of how the dynamic crossover in the relaxation time can potentially affect the MSD measured on the resolution time . This can be illustrated for the rotational MSD, which can be rewritten in terms of the rotational MSF and the normalized autocorrelation function of rotations




The generic form of is the initial ballistic (Gaussian) decay, followed by exponential collective relaxation: (or, alternatively, multi-exponential or stretched exponential term) Hansen and McDonald (2003). In the entire temperature range studied for glycerol we find that falls between the time of ballistic relaxation and the time of collective exponential relaxation : . One therefore gets


In the limit of , the relaxation time is not expected to affect the MSD. Its magnitude is mostly determined by the amplitude of the Gaussian component of the relaxation dynamics, in agreement with the arguments presented in relation to Figs. 2 and 3. Therefore, if the dynamic crossover and the kink of the MSD occur at the same temperature Schirò et al. (2013) one has to relate this coincidence to a structural change and not to a direct effect of the relaxation time on the MSD. The hypothesis that the crossover in the relaxation time affects the MSD is, therefore, not supported by our simulation results.

Figure 6: Diffusion coefficients recorded experimentally by NMR (“Exp”, Chen et al. (2006a)) and obtained from the simulations (“MD”). The dashed line is a regression drawn through the MD points.

The results for the average rotational relaxation time for all protons in glycerol are shown in Fig. 5. It is calculated by integrating the time correlation function


where corresponds to the normalized time correlation function in Eq. (6). These results are shown by the open points in Fig. 5.

We have additionally calculated the time correlation function based on the dynamic variable of the electric field produced by the rest of the glycerol liquid at the center of mass of a given target molecule (). The microscopic electric field is therefore a fluctuating local field producing a torque on the glycerol’s dipole moment. The results for the average relaxation times obtained from the corresponding time correlation functions through Eq. (8) are shown by the closed points in Fig. 5. There is a good agreement between and suggesting that the electric field fluctuations are caused by molecular rotations, as one would anticipate from the standard Debye model of dielectric relaxation Fröhlich (1958); Richert (2014).

The average relaxation times from MD simulations are compared in Fig. 5 with the average relaxation time calculated from the Cole-Davidson fit of glycerol’s loss spectrum reported by broad-band dielectric spectroscopy Menon et al. (1992) (solid line). There is a very good agreement between the simulations and experimental dielectric data at high temperatures, suggesting that the adopted force field Caleman et al. (2012a) (see Appendix A) is well parametrized for glycerol rotations. There is a less satisfactory agreement between the diffusion coefficient calculated from MD and measured by NMR (Fig. 6). Differences between quasi-elastic neutron scattering (QENS) and NMR/viscosity data for glycerol self-diffusion have been documented in the past Cuello et al. (1998); Mamontov (2012) and might contribute to the discrepancy.

Figure 7: The Kirkwood factor (a) and dielectric constant (b) of glycerol calculated from MD (circles) and measured in bulk samples experimentally Matyushov and Richert (2016) (squares). The dashed lines are linear fits to the corresponding subsets of data to guide the eye. The Kirkwood factors in (a) were obtained both in NVE and NVT separate simulation runs.

The dynamic crossover occurs in the range of temperatures when the -relaxation time becomes comparable to the length of the simulation trajectory ns. In fact, the time window on which the time correlation function is calculated from the simulation trajectories is always shorter, . We therefore stop observing the slow relaxation in simulations when the -relaxation time becomes longer than . The relative weight of the fast relaxation in increases and we observe this as a dynamical crossover.

What our data do not seem to address is why the kinks in the rotational and translational MSDs and the corresponding dynamical crossovers in the rotational relaxation times and translational diffusion (Figs. 5 and 6) all occur in the same range of temperatures. A possible scenario to explain this coincidence might include a structural transition resulting in a drop of the configurational entropy Matyushov and Angell (2007). According to the general arguments based on the Adam-Gibbs relation Angell et al. (2000), this would result in a much slower main relaxation process, which would sharply disappear from the observation window of our numerical experiment. While our results presented below do support alteration of glycerol’s orientational structure, we do not have a direct evidence for a discontinuous change in the configurational entropy.

In order to identify possible structural changes, we have looked at the temperature dependence of the Kirkwood factor reflecting orientational correlations in the liquid


Here, are the unit vectors of molecular dipoles (4.6 D in the force field used in our simulations). The Kirkwood factor was in turn used in the Kirkwood-Onsager relation Fröhlich (1958) to calculate the dielectric constant (the glycerol force field is non-polarizable and the refractive index is equal to unity). The results of these calculations are shown in Fig. 7.

Figure 8: Projection of the pair distribution function of glycerol on the orientational invariant calculated from MD simulations at different temperatures.

The Kirkwood factor shows a discontinuous increase at K, which results in the corresponding increase of the dielectric constant calculated from MD simulations. The increase in is caused by the emergence of long-range orientational correlations of glycerol dipoles at low temperatures, as is illustrated in Fig. 8. We show there the projection of the pair correlation function of glycerol , depending on the distance between two molecules and their orientations and , on the rotational invariant of the scalar product between the unit vectors of the dipole moments . The corresponding pair distribution function Hansen and McDonald (2003)


at different temperatures in shown in Fig. 8.

Figure 9: Orientational order parameters and in Eq. (11) calculated from MD trajectories at different temperatures.

It is clear that a long-range oscillatory pattern, reflecting preferential parallel alignments of the dipoles, appears at low temperatures. The dipolar alignments are responsible for an increase in the low-temperature Kirkwood factor, , is the number density. Despite these long-range orientational correlations, the low-temperature phase does not show any specific orientational order, as confirmed by calculations of the first and second orientational order parameters Ayton and Patey (1996); Vertogen and de Jeu (1988) (Fig. 9) as explained below. No translational order is observed either: the radial pair distribution functions are nearly identical at low and high temperatures (Fig. 10). We therefore can conclude that the low-temperature phase is a disordered liquid.

Figure 10: Pair distribution functions of the center of mass of glycerol calculated at the temperatures indicated in the plot. The calculated functions nearly coincide on the scale of the plot.

The orientational order can be detected by orientational order parameters typically defined for liquid crystals Vertogen and de Jeu (1988). The order parameter is the average th order Legendre polynomial


relative to the liquid director ; is the number of molecules in the liquid. The director is identified as the eigenvector corresponding to the largest eigenvalue of the tensor


where and are the Cartesian projections and is the Kronecker delta function. The results of calculations for the first and second order parameters () are shown in Fig. 9. No orientational order can be identified at low temperatures from these calculations.

The jump in the simulated dielectric constant is in stark disagreement with the linear dielectric experiment Matyushov and Richert (2016) where no discontinuities were observed (squares in Fig. 7b). The results of simulations are in fair agreement with experiment at high temperatures, but the increase in the Kirkwood factor at lower temperatures (Fig. 7a) makes the dielectric constant much higher than observations. Since the crossover temperature for the dielectric constant is roughly consistent with the kinks in the rotational and translational MSDs, we conclude that restricting the observation window not only makes changes to the observable relaxation dynamics, but also does not allow certain orientational correlations to relax. As a result, we observe a long-range orientational order frozen on the observation time-scale. This implies that both the low-temperature Kirkwood factor and the corresponding dielectric constant shown in Fig. 7 are non-equilibrium quantities. A similar, about five times compared to the bulk (Fig. 9 in Ref. Kasina et al., 2015), increase in the dielectric constant was observed for ultrathin films of glycerol obtained by vapor deposition Capponi et al. (2012). Subsequent combined dielectric and calorimetry measurements have suggested the existence of rigid polar clusters, which relax as a whole, with an enhanced cluster dipole moment Kasina et al. (2015). There is also recent evidence of an unrelaxed orientational order in organic glasses obtained by surface deposition Dalal et al. (2015).

Figure 11: at ps (black lines) and ns (blue lines) calculated from NVT simulations of glycerol at different temperatures indicated in the plot. The red line indicates NVE simulation at 270 K.

The existence of highly correlated clusters should be seen in the heterogeneity of binary correlations expressed in terms of fourth-order correlation functions Berthier and Biroli (2011). In order to test this hypothesis, we made the next step of calculating the distance- and time-dependent correlations between binary dipolar orientational correlations expressed through the instantaneous Kirkwood factors. Specifically, the quantity


was constructed at each point of the simulation trajectory to reflect the instantaneous binary correlations of the chosen dipole moment with all remaining dipoles in the liquid. Obviously, one has . We then constructed the distance- and time-dependent correlation between the local binary correlations as follows


where the average is taken along the simulation trajectory and is the liquid volume. The normalization of relates it to the Kirkwood factor


Similarly to in Fig. 8, but significantly more pronounced, we observe the rise of long-range heterogeneous correlations at low temperatures (Fig. 11).

Iv Discussion and implications for the protein dynamical transition

We obtained here, by computer simulations, both a kink in the temperature dependence of the MSD (dynamical transition) and the dynamical crossovers in the relaxation times. Both effects have been observed experimentally and a link between them has been suggested through some sort of structural transition in the liquid Chen et al. (2006b); Cupane et al. (2014); Mallamace et al. (2015). The answer to the ongoing discussion of whether a purely dynamical crossover or a structural transition explains the data might be that both are present. However, in contrast to the scenarios involving thermodynamic liquid-liquid transitions, both the structural and relaxation time crossovers have a dynamic origin. The structural crossover is caused by the inability of certain structural correlations to relax on the observation window. There is nothing in our data that connects the appearance of such structural correlations to a thermodynamic transition between two phases of a bulk material. This distinction becomes, however, less loaded with physical meaning in the low-temperature state. When the relaxation time of the “orientationally correlated liquid” becomes much longer than any conceivable experimental time, one has to distinguish this state of the material as an “orientationally correlated glass”, with all relevant properties distinct from the “ordinary” glass. One arrives at polyamorphism of the glass state Angell et al. (2000) caused by long-ranged orientational correlations.

The observation of an increase in the dielectric constant of glycerol below the dynamical transition, here by simulations and for vapor deposited glasses experimentally Capponi et al. (2012); Kasina et al. (2015), adds a structural component to the standard picture of ergodicity breaking of glass science. The standard paradigm is that the glass does not have the ability to relax, but maintains the structure of the liquid. This is indeed true for the positional structure of the glycerol molecules. However, the inability of dipolar orientations to relax causes orientational heterogeneity represented by correlated dipolar clusters, which do not relax on the observation time-scale. The long-sought growth of the structural order of glass-formers on approach to the laboratory glass transition might be, therefore, best discovered by experiments probing the heterogeneity of orientational multipolar correlations.

The conclusion that no thermodynamic transformation is at work in creating dipolar domains does not make our observations less “interesting”. In particular, this scenario is relevant to the role of dynamics and structure of protein’s hydration shells in the protein function. About anything related to the protein structure and function has to be described as metastable. Protein itself is unstable to either hydrolysis or association, both bringing it to a thermodynamically more stable state Scheraga (2015). The function of proteins as enzymes catalyzing specific biochemical reactions is even more affected by the notion of a finite “observation window” Matyushov (2015). This idea implies that any dynamical or structural information related to the protein itself or to its hydration shell has to be considered from the perspective of a finite observation window provided by the reaction rate, i.e., the characteristic time on which the reactants climb the activation barrier separating them from the products. A dynamic process slower than the rate becomes dynamically frozen and does not contribute to the fluctuation spectrum of the bath driving the reaction.

The ability of the solvent to preserve a specific structure distinct from its thermodynamic state on a given observation window immediately implies that an enzymetic reaction will “see” different solvents, with potentially dramatically different properties (such as polarity), depending on the reaction rate. Figure 7 provides a dramatic confirmation of this possibility showing the ability of glycerol to possess a very high dielectric constant due to its inability to relax its long-range orientational correlations on a given observation window. A related example, with a similar phenomenology, is the appearance of polarized (ferroelectric) domains in the hydration shells of proteins observed on the time-scale of simulations Martin and Matyushov (2015). Similarly to our present results for glycerol, these domains might well equilibrate to zero overall dipole on longer time-scales, but a non-zero net dipole of the shell will be recorded by any kinetic process occurring faster than the domain relaxation dynamics.

Bulk glycerol studied by linear dielectric spectroscopy does not display the features indicative of domain formation. There is a general agreement that linear dielectric spectroscopy does not directly probe heterogeneity of a bulk material Richert (2014). However, it might still be illuminating to ask why the relaxation of oriented domains in the bulk is not observed by dielectric spectroscopy. One possible answer to this is that the lifetime of a domain is smaller than its rotational relaxation time. The domains dissolve before there is a chance to probe their rotational relaxation. Increasing the lifetime of domains, as potentially achieved by surface vapor deposition Capponi et al. (2012); Kasina et al. (2015), might create conditions for observing the large dipole of the correlated domain.

The identification of the MSD crossover with the cage dynamics, in the combination with nearly identical behavior of MSD of glycerol and lysozyme-glycerol Capaccioli et al. (2012), puts under question the need for a special relaxation process of the hydration shell Fenimore et al. (2004); Lubchenko et al. (2005); Frauenfelder et al. (2009) to explain these data. It appears that fast secondary relaxation of bulk glycerol ( in the standard classification of glass science Ediger et al. (1996); Lunkenheimer and Loidl (2002); Lunkenheimer et al. (2010)) is sufficient to describe the glycerol-protein system. It does not necessarily mean that the same situation repeats itself for a hydrated protein, or applies equally well to the Mössbauer experiment with a much longer resolution time of ns Parak (2003). Some experimental data indeed claim the existence of independent relaxation processes of the protein hydration shells with significantly slower relaxation times Pal and Zewail (2004); Bhattacharyya (2008). The resolution of this claim, however, depends on the water mode probed by the observations. There is a relatively insignificant slowing down of water’s single-molecule rotational dynamics in hydration shells Laage et al. (2011). An attempt to find a separate dynamic process in density fluctuations (translations) probed by depolarized light scattering resulted in the realization that cross protein-water correlations, instead of a separate dynamic process, can explain the data Martin and Matyushov (2014). However, the collective variable of the shell dipole moment can be characterized as a separate dynamic process, which is both significantly slower and is spatially extended into the bulk Matyushov (2012). From a general perspective, a strong perturbation of the forces existing in the bulk is required for a new dynamic process to appear. If a significant alteration of the hydrogen-bond network is achieved in the solvation layer, one can expect a separate dynamic process to show up. The extent of such network perturbation is where the distinction between glycerol and water might be found.

This research was supported by the National Science Foundation (CHE-1464810) and through XSEDE (TG-MCB080116N). We are grateful to Ranko Richert for his help with dielectric data for glycerol and for helpful suggestions on the manuscript. Discussing glass science with Austen Angell has inspired many of the questions raised in this paper.

Appendix A Simulation protocol

Molecular dynamics (MD) simulations were performed for twelve different temperatures (147, 168, 179, 195, 214, 239, 255, 275, 287, 302, 312, 334 K) in a cubic box consisting of 1000 glycerol molecules using the OPLS-AA (Optimized Potentials for Liquid Simulation - All Atoms) force field Caleman et al. (2012b) within the Gromacs Hess et al. (2008) simulation package. The initial box with 1000 glycerol molecules was downloaded from virtual chemistry website Frenkel et al. (2000); Caleman et al. (2012b). After the initial NPT and NVT equilibration runs, the raw box was used to produced 50 ns equilibration trajectories in the NVE ensemble with no constraints. The output was then used to generate the NVE trajectories.

Each system was initialized with a 300 ps NVT run using a Nose-Hoover thermostat with H-bonds constrained followed by a 300 ps run with no constraints. Then a 1-3 ns NVE run was generated to check for stability before doing the 50 ns production run for each temperature. The time step for all production runs was 0.5 fs, with all atoms (including hydrogens) free to move according to the OPLS-AA force field parameters. The group cutoff-scheme was used with an update time of 5 ns and a cutoff distance of 1.1 nm for the shifted Lennard-Jones and electrostatic interactions with a group list distance of 13 Å renewing every 10 simulation steps. Long-ranged electrostatic interactions were calculated with the particle mesh Ewald method. Additional trajectories (tens of ns) were generated for a separate raw box under the NVT ensemble as a means to compare the results between NVE and NVT ensembles. These NVT simulations with the H-bond constraints were equilibrated with an initial 3 ns run and an additional 3 ns with the constraints removed, assuring minimization and equilibration before production runs for each temperature were executed. The 50 ns NVT simulations were carried out for temperatures 230, 240, 250, 260, 270, 280 K. In this case, the Verlet cutoff-scheme was implemented and a Nose-Hoover thermostat was used. All simulations were carried out using periodic boundary conditions.


  1. R. Kubo, Rep. Prog. Phys. 29, 255 (1966).
  2. J. P. Hansen and I. R. McDonald, Theory of Simple Liquids (Academic Press, Amsterdam, 2003).
  3. F. Gabel, D. Bicout, U. Lehnert, M. Tehei, M. Weik, and G. Zaccai, Quat. Rev. Biophys. 35, 327 (2002).
  4. F. G. Parak, Rep. Prog. Phys. 66, 103 (2003).
  5. R. D. Young, H. Frauenfelder, and P. W. Fenimore, Phys. Rev. Lett. 107, 158102 (2011).
  6. C. A. Angell, Science 267, 1924 (1995).
  7. J. H. Roh, V. N. Novikov, R. B. Gregory, J. E. Curtis, Z. Chowdhuri, and A. P. Sokolov, Phys. Rev. Lett. 95, 038101 (2005).
  8. S. Khodadadi, A. Malkovskiy, A. Kisliuk, and A. P. Sokolov, Biochim. Biophys. Acta 1804, 15 (2010).
  9. L. Hong, N. Smolin, B. Lindner, A. P. Sokolov, and J. C. Smith, Phys. Rev. Lett. 107, 148102 (2011).
  10. W. Doster, S. Cusack, and W. Petry, Nature 337, 754 (1989).
  11. G. Zaccai, Science 288, 1604 (2000).
  12. K. Achterhold, C. Keppler, A. Ostermann, U. van Bürck, W. Sturhahn, E. E. Alp, and F. G. Parak, Phys. Rev. E 65, 051916 (2002).
  13. W. Doster, Eur. Biophys. J. 37, 591 (2008).
  14. S. Magazù, F. Migliardo, and A. Benedetto, J. Phys. Chem. B 115, 7736 (2011).
  15. P. W. Fenimore, H. Frauenfelder, S. Magazù, B. H. McMahon, F. Mezei, F. Migliardo, R. D. Young, and I. Stroe, Chem. Phys. 424, 2 (2013).
  16. G. Schirò, Y. Fichou, F.-X. Gallat, K. Wood, F. Gabel, M. Moulin, M. Härtlein, M. Heyden, J.-P. Colletier, A. Orecchini, et al., Nat. Comm. 6, 6490 (2015).
  17. S. Capaccioli, K. L. Ngai, S. Ancherbak, and A. Paciaroni, J. Phys. Chem. B 116, 1745 (2012).
  18. J. C. Smith, Quat. Rev. Biophys. 24, 227 (1991).
  19. G. Zaccai, J. Phys.: Condens. Matter 15, S1673 (2003).
  20. D. V. Matyushov, J. Phys.: Condens. Matter 27, 473001 (2015).
  21. P. W. Fenimore, H. Frauenfelder, B. H. McMahon, and R. D. Young, Proc. Natl. Acad. Sci. USA 101, 14408 (2004).
  22. V. Lubchenko, P. G. Wolynes, and H. Frauenfelder, J. Phys. Chem. B 109, 7488 (2005).
  23. C. A. Angell, K. L. Ngai, G. B. McKenna, and S. W. Martin, J. Appl. Phys. 88, 3113 (2000).
  24. E. Donth, The Glass Transition: Relaxation Dynamics in Liquids and Disordered Materials (Springer, Berlin, 2001).
  25. G. Chen, P. W. Fenimore, H. Frauenfelder, F. Mezei, J. Swenson, and R. D. Young, Phil. Mag. 88, 33 (2008).
  26. S. Khodadadi, S. Pawlus, J. H. Roh, V. G. Sakai, E. Mamontov, and A. P. Sokolov, J. Chem. Phys. 128, 195106 (2008).
  27. H. Frauenfelder, G. Chen, J. Berendzen, P. W. Fenimore, H. Jansson, B. H. McMahon, I. R. Stroe, J. Swenson, and R. D. Young, Proc. Natl. Acad. Sci. USA 106, 5129 (2009).
  28. G. A. Samara, J. Phys.: Condens. Matter 15, R367 (2003).
  29. J. Wuttke, W. Petry, G. Coddens, and F. Fujara, Phys. Rev. E 52, 4026 (1995).
  30. W. Doster, M. Diehl, R. Gebhardt, R. E. Lechner, and J. Pieper, Chem. Phys. 292, 487 (2003).
  31. Z. Yi, Y. Miao, J. Baudry, N. Jain, and J. C. Smith, J. Phys. Chem. B 116, 5028 (2012).
  32. G. Schirò, F. Natali, and A. Cupane, Phys. Rev. Lett. 109, 128102 (2012).
  33. C. Caleman, P. J. van Maaren, M. Hong, J. S. Hub, L. T. Costa, and D. van der Spoel, J. Chem. Theory Comp. 8, 61 (2012a).
  34. J. Teixeira, M.-C. Bellissent-Funel, S. H. Chen, and A. J. Dianoux, Phys. Rev. A 31, 1913 (1985).
  35. S. H. Chen, P. Gallo, F. Sciortino, and P. Tartaglia, Phys. Rev. E 56, 4231 (1997).
  36. L. Liu, A. Faraone, and S.-H. Chen, Phys. Rev. E 65, 041506 (2002).
  37. S. H. Chen, C. Liao, F. Sciortino, P. Gallo, and P. Tartaglia, Phys. Rev. E 59, 6708 (1999).
  38. L. Liu, S.-H. Chen, A. Faraone, C.-W. Yen, and C.-Y. Mou, Phys. Rev. Lett. 95, 117802 (2005).
  39. B. Chen, E. E. Sigmund, and W. P. Halperin, Phys. Rev. Lett. 96, 145502 (2006a).
  40. E. Mamontov, Chem. Phys. Lett. 530, 55 (2012).
  41. M. Fomina, G. Schirò, and A. Cupane, Biophysical Chemistry 185, 25 (2014).
  42. S.-H. Chen, L. Liu, E. Fratini, P. Baglioni, and E. Mamontov, Proc. Natl. Acad. Sci. 103, 9012 (2006b).
  43. J.-M. Zanotti, M.-C. Bellissent-Funel, and S.-H. Chen, Europhys. Lett. 71, 91 (2005).
  44. A. Cupane, M. Fomina, I. Piazza, J. Peters, and G. Schirò, Phys. Rev. Lett. 113, 215701 (2014).
  45. G. Schirò, M. Fomina, and A. Cupane, J. Chem. Phys. 139, 121102 (2013).
  46. F. Mallamace, C. Corsaro, D. Mallamace, S. Vasi, C. Vasi, H. E. Stanley, and S.-H. Chen, J. Chem. Phys. 142, 215103 (2015).
  47. F. Mallamace, M. Broccio, C. Corsaro, A. Faraone, U. Wanderlingh, L. Liu, C.-Y. Mou, and S. H. Chen, J. Chem. Phys. 124, 161102 (2006).
  48. R. Richert, Adv. Chem. Phys. 156, 101 (2014).
  49. N. Menon, K. P. O’Brien, P. K. Dixon, L. Wu, S. R. Nagel, B. D. Williams, and J. P. Carini, J. Non-Crystal. Solids 141, 61 (1992).
  50. H. Fröhlich, Theory of dielectrics (Oxford University Press, Oxford, 1958).
  51. G. J. Cuello, F. J. Bermejo, R. Fayos, R. Fernandez-Perea, A. Criado, F. Trouw, C. Tam, H. Schober, E. Enciso, and N. G. Almarza, Phys. Rev. B 57, 8254 (1998).
  52. D. V. Matyushov and R. Richert, J. Chem. Phys. 144, 041102 (2016).
  53. D. V. Matyushov and C. A. Angell, J. Chem. Phys. 126, 094501 (2007).
  54. G. Ayton and G. N. Patey, Phys. Rev. Lett. 76, 239 (1996).
  55. G. Vertogen and W. H. de Jeu, Thermotropic Liquid Crystals, Fundamentals (Springer-Verlag, Berlin, 1988).
  56. A. Kasina, T. Putzeys, and M. Wübbenhorst, J. Chem. Phys. 143, 244504 (2015).
  57. S. Capponi, S. Napolitano, and M. Wübbenhorst, Nat. Comm. 3, 1233 (2012).
  58. S. S. Dalal, D. M. Walters, I. Lyubimov, J. J. de Pablo, and M. D. Ediger, Proc. Natl. Acad. Sci. USA 112, 4227 (2015).
  59. L. Berthier and G. Biroli, Rev. Mod. Phys. 83, 587 (2011).
  60. H. A. Scheraga, Quart. Rev. Biophys. 48, 117 (2015).
  61. D. R. Martin and D. V. Matyushov, J. Phys. Chem. Lett. 6, 407 (2015).
  62. M. D. Ediger, C. A. Angell, and S. R. Nagel, J. Phys. Chem. 100, 13200 (1996).
  63. P. Lunkenheimer and A. Loidl, Chem. Phys. 284, 205 (2002).
  64. P. Lunkenheimer, U. Schneider, R. Brand, and A. Loid, Contemporary Physics 41, 15 (2010).
  65. S. K. Pal and A. H. Zewail, Chem. Rev. 104, 2099 (2004).
  66. K. Bhattacharyya, Chem. Commun. pp. 2848–2857 (2008).
  67. D. Laage, G. Stirnemann, F. Sterpone, R. Rey, and J. T. Hynes, Annu. Rev. Phys. Chem. 62, 395 (2011).
  68. D. R. Martin and D. V. Matyushov, J. Chem. Phys. 141, 22D501 (2014).
  69. D. V. Matyushov, J. Chem. Phys. 136, 085102 (2012).
  70. C. Caleman, P. J. van Maaren, M. Hong, J. S. Hub, L. T. Costa, and D. van der Spoel, J. Chem. Theory Comput. 8, 61 (2012b).
  71. B. Hess, C. Kutzner, D. van der Spoel, and E. Lindahl, Journal of Chemical Theory and Computation 4, 435 (2008).
  72. M. Frenkel, X. Hong, Q. Dong, X. Yan, and R. Chirico, Densities of Halohydrocarbons (Springer, Berlin / Heidelberg, 2000).
This is a comment super asjknd jkasnjk adsnkj
The feedback cannot be empty
Comments 0
The feedback cannot be empty
Add comment

You’re adding your first comment!
How to quickly get a good reply:
  • Offer a constructive comment on the author work.
  • Add helpful links to code implementation or project page.