Structural dynamics of supercooled water from quasielastic neutron scattering and molecular simulations

One of the outstanding challenges presented by liquid water is to understand how molecules can move on a picosecond time scale despite being incorporated in a three-dimensional network of relatively strong H-bonds. This challenge is exacerbated in the supercooled state, where the dramatic slowing down of structural dynamics is reminiscent of the, equally poorly understood, generic behavior of liquids near the glass transition temperature. By probing single-molecule dynamics on a wide range of time and length scales, quasielastic neutron scattering (QENS) can potentially reveal the mechanistic details of water's structural dynamics, but because of interpretational ambiguities this potential has not been fully realized. To resolve these issues, we present here an extensive set of high-quality QENS data from water in the range 253-293 K and a corresponding set of molecular dynamics (MD) simulations to facilitate and validate the interpretation. Using a model-free approach, we analyze the QENS data in terms of two motional components. Based on the dynamical clustering observed in MD trajectories, we identify these components with two distinct types of structural dynamics: picosecond local (L) structural fluctuations within dynamical basins and slower interbasin jumps (J). The Q-dependence of the dominant QENS component, associated with J dynamics, can be quantitatively rationalized with a continuous-time random walk (CTRW) model with an apparent jump length that depends on low-order moments of the jump length and waiting time distributions. Using a simple coarse-graining algorithm to quantitatively identify dynamical basins, we map the newtonian MD trajectory on a CTRW trajectory, from which the jump length and waiting time distributions are computed. The jump length distribution is gaussian and the rms jump length increases from 1.5 to 1.9 A? as the temperature increases from 253 to 293 K. The rms basin radius increases from 0.71 to 0.75 A? over the same range. The waiting time distribution is exponential at all investigated temperatures, ruling out significant dynamical heterogeneity. However, a simulation at 238 K reveals a small but significant dynamical heterogeneity. The macroscopic diffusion coefficient deduced from the QENS data agrees quantitatively with NMR and tracer results. We compare our QENS analysis with existing approaches, arguing that the apparent dynamical heterogeneity implied by stretched exponential fitting functions results from the failure to distinguish intrabasin (L) from interbasin (J) structural dynamics. We propose that the apparent dynamical singularity at ～220 K corresponds to freezing out of J dynamics, while the calorimetric glass transition corresponds to freezing out of L dynamics.     

Energetics and dynamics in MbCN: CN--vibrational relaxation from molecular dynamics simulations

The dynamics of the cyanide anion bound to sperm-whale myoglobin is investigated using atomistic simulations. With density-functional theory, a 2D potential energy surface for the cyanide-heme complex is calculated. Two deep minima with a stabilization energy of approximately 50 kcal/mol corresponding to two different binding orientations (Fe-CN and Fe-NC) of the ligand are found. The Fe-CN conformation is favored over Fe-NC by several kcal/mol. Mixed quantum mechanics/molecular mechanics calculations show that the binding orientation affects the bond strength of the ligand, with a significantly different bond length and a 25 cm-1 shift in the fundamental CN-frequency. For the molecular dynamics (MD) simulations, a 3-center fluctuating charge model for the Fe-CN unit is developed that captures polarization and ligand-metal charge transfer. Stability arguments based on the energetics around the active site and the CN- frequency shifts suggest that the Fe-CN conformation with epsilon-protonation of His epsilon 64 are most likely, which is in agreement with experiment. Both equilibrium and nonequilibrium MD simulations are carried out to investigate the relaxation time scale and possible relaxation pathways in bound MbCN. The nonequilibrium MD simulations with a vibrationally excited ligand reveal that vibrational relaxation takes place on a time scale of hundreds of picoseconds within the active site. This finding supports the hypothesis that the experimentally observed relaxation rate (3.6 ps) reflects the repopulation of the electronic ground state.     

Energy relaxation of intermolecular motions in supercooled water and ice: a molecular dynamics study

We investigate the energy relaxation of intermolecular motions in liquid water at temperatures ranging from 220 K to 300 K and in ice at 220 K using molecular dynamics simulations. We employ the recently developed frequency resolved transient kinetic energy analysis, which provides detailed information on energy relaxation in condensed phases like two-color pump-probe spectroscopy. It is shown that the energy cascading in liquid water is characterized by four processes. The temperature dependences of the earlier three processes, the rotational-rotational, rotational-translational, and translational-translational energy transfers, are explained in terms of the density of states of the intermolecular motions. The last process is the slow energy transfer arising from the transitions between potential energy basins caused by the excitation of the low frequency translational motion. This process is absent in ice because the hydrogen bond network rearrangement, which accompanies the interbasin transitions in liquid water, cannot take place in the solid phase. We find that the last process in supercooled water is well approximated by a stretched exponential function. The stretching parameter, β, decreases from 1 to 0.72 with decreasing temperature. This result indicates that the dynamics of liquid water becomes heterogeneous at lower temperatures.     

Local order and dynamics in supercooled water: a study by IR spectroscopy and molecular dynamic simulations

Micron-sized water droplets in a cryogenic flow tube were probed by IR spectroscopy. The analysis of the IR spectra suggests that there is a relative increase of about 30% in the fraction, f(L), of low density domains in water on cooling over the temperature range between 300 and 240 K. The results derived from the experiments agree qualitatively with those of molecular dynamics (MD) simulations in terms of the increase in the f(L) values. The MD simulations show that the intensities of the mode at about 100 cm(-1) for the molecules in the low density domains are reduced in comparison to the average, while the intensities and frequencies of the librational mode at 600 cm(-1) are increased. Furthermore, the reorientations (dielectric relaxation times) in these domains are found to be somewhat slower, pointing to the fact that these low density "cages" live longer than the average local molecular environments in supercooled water.     

Ab initio molecular dynamics simulations of local structure of supercooled Ni

We report results of first-principles molecular dynamics simulations for stable and undercooled nickel liquids. The calculated structure factors as a function of temperature are discussed with respect to recent experimental measurements. In addition, structural analysis using bonding orientational order and three-dimensional pair analysis techniques have been performed in detail and the effect of undercooling on the microstructure has been analyzed. More particularly, we show the importance of fivefold symmetry local structures.     

Rotational dynamics of benzene and water in an ionic liquid explored via molecular dynamics simulations and NMR T1 measurements

The rotational dynamics of benzene and water in the ionic liquid (IL) 1-butyl-3-methylimidazolium chloride are studied using molecular dynamics (MD) simulation and NMR T(1) measurements. MD trajectories based on an effective potential are used to calculate the (2)H NMR relaxation time, T(1) via Fourier transform of the relevant rotational time correlation function, C(2R)(t). To compensate for the lack of polarization in the standard fixed-charge modeling of the IL, an effective ionic charge, which is smaller than the elementary charge is employed. The simulation results are in closest agreement with NMR experiments with respect to the temperature and Larmor frequency dependencies of T(1) when an effective charge of ±0.5e is used for the anion and the cation, respectively. The computed C(2R)(t) of both solutes shows a bi-modal nature, comprised of an initial non-diffusive ps relaxation plus a long-time ns tail extending to the diffusive regime. Due to the latter component, the solute dynamics is not under the motional narrowing condition with respect to the prevalent Larmor frequency. It is shown that the diffusive tail of the C(2R)(t) is most important to understand frequency and temperature dependencies of T(1) in ILs. On the other hand, the effect of the initial ps relaxation is an increase of T(1) by a constant factor. This is equivalent to an "effective" reduction of the quadrupolar coupling constant (QCC). Thus, in the NMR T(1) analysis, the rotational time correlation function can be modeled analytically in the form of aexp (-t/τ) (Lipari-Szabo model), where the constant a, the Lipari-Szabo factor, contains the integrated contribution of the short-time relaxation and τ represents the relaxation time of the exponential (diffusive) tail. The Debye model is a special case of the Lipari-Szabo model with a = 1, and turns out to be inappropriate to represent benzene and water dynamics in ILs since a is as small as 0.1. The use of the Debye model would result in an underestimation of the QCC by a factor of 2-3 as a compensation for the neglect of the Lipari-Szabo factor.     

Rotational dynamics of a discotic liquid crystal studied by deuteron spin relaxation

We report on measurements of deuteron quadrupolar splitting and spin-lattice relaxation times T_1Q and T_1Z in the columnar phase of a ring-deuteriated hexaoctyloxyrufigallol at 46 MHz as a function of temperature. To describe small-step rotations of these molecules within each column in the columnar phase, a space-fixed frame is used to diagonalize the molecular diffusion tensor. The principal diffusion constants in this so called 'anisotropic viscosity' model D_alpha and D_beta are for rotations of a molecule around and perpendicular to the columnar axis, respectively. A global target analysis of the spectral densities at seven temperatures in a minimization procedure was carried out. We found that D_alpha > D_beta, which is consistent with the picture that the motion towards or away from the local director tends to disrupt the packing of molecules within the column.     

Dissociative water potential for molecular dynamics simulations

A new interatomic potential for dissociative water was developed for use in molecular dynamics simulations. The simulations use a multibody potential, with both pair and three-body terms, and the Wolf summation method for the long-range Coulomb interactions. A major feature in the potential is the change in the short-range O-H repulsive interaction as a function of temperature and/or pressure in order to reproduce the density-temperature curve between 273 K and 373 at 1 atm, as well as high-pressure data at various temperatures. Using only the change in this one parameter, the simulations also reproduce room-temperature properties of water, such as the structure, cohesive energy, diffusion constant, and vibrational spectrum, as well as the liquid-vapor coexistence curve. Although the water molecules could dissociate, no dissociation is observed at room temperature. However, behavior of the hydronium ion was studied by introduction of an extra H+ into a cluster of water molecules. Both Eigen and Zundel configurations, as well as more complex configurations, are observed in the migration of the hydronium.     

Rotational motion in liquid water is anisotropic: a nuclear magnetic resonance and molecular dynamics simulation study

Experimental NMR measurements of the deuterium and (17)O T(1) relaxation times in deuterium-enriched liquid water have been performed from 275 to 350 K. These relaxation times can yield rotational correlation times of appropriate molecule-fixed unit vectors if the quadrupole coupling constants and asymmetry parameters are known. We determine the latter from ab initio studies of water clusters and experimental chemical shift measurements. We find that the rotational correlation time for the OD bond vector in D(2)(16)O varies from 5.8 ps at 275 K to 0.86 ps at 350 K, and that the rotational correlation time for the out-of-plane vector of dilute D(2)(17)O in D(2)(16)O varies from 4.4 ps at 275 K to 0.64 ps at 350 K. These results indicate that the rotational motion of water is anisotropic. Molecular dynamics simulations of liquid water are in good agreement with these experiments at the higher temperatures, but the simulation results are considerably faster than experiment at the lower temperatures.     

Void-induced dissolution in molecular dynamics simulations of NaCl and water

To gain a better understanding of the interaction of water and NaCl at the surface during dissolution, we have used molecular dynamics to simulate the interface with two equal-sized slabs of solid NaCl and liquid water in contact. The introduction of voids in the bulk of the salt, as well as steps or pits on the surface of the NaCl slab results in a qualitative change of system structure, as defined by radial distribution functions (RDFs). As an example, the characteristic Na-Na RDF for the system changes from regularly spaced narrow peaks (corresponding to an ordered crystalline structure), to a broad primary and smaller secondary peak (corresponding to a disordered structure). The change is observed at computationally short time scales of 100 ps, in contrast with a much longer time scale of 1 mus expected for complete mixing in the absence of defects. The void fraction (which combines both bulk and surface defects) required to trigger dissolution varies between 15%-20% at 300 K and 1 atm, and has distinct characteristics for the physical breakdown of the crystal lattice. The void fraction required decreases with temperature. Sensitivity studies show a strong dependence of the critical void fraction on the quantity and distribution of voids on the surface, with systems containing a balanced number of surface defects and a rough surface showing a maximum tendency to dissolve. There is a moderate dependence on temperature, with a 5% decrease in required void fraction with a 100 K increase in temperature, and a weak dependence on water potential model used, with the SPC, SPC/E, TIP4P, and RPOL models giving qualitatively identical results. The results were insensitive to the total quantity of water available for dissolution and the duration of the simulation.     

Vesicle shapes from molecular dynamics simulations

Lipid bilayer membranes are known to form various structures such as large sheets or vesicles. When the two leaflets of the bilayer have an equal composition, the membrane preferentially forms a flat sheet or a spherical vesicle. However, a difference in the composition of the two leaflets may result in a curved bilayer or in a wide variety of vesicle shapes. Vesicles with different shapes have already been shown in experiments and diverse vesicle shapes have been predicted theoretically from energy minimization of continuous curves. Here we present a molecular dynamics study of the effect of small changes in the phospholipid headgroups on the spontaneous curvature of the bilayer and on the resulting vesicle shape transformations. Small asymmetries in the bilayers already result in high spontaneous curvature and large vesicle deformations. Vesicle shapes that are formed include ellipsoids, discoids, pear-shaped vesicles, cup-shaped vesicles, as well as budded vesicles. Comparison of these vesicles with theoretically derived vesicle shapes shows both resemblances and differences.     

Estimating entropies from molecular dynamics simulations

While the determination of free-energy differences by MD simulation has become a standard procedure for which many techniques have been developed, total entropies and entropy differences are still hardly ever computed. An overview of techniques to determine entropy differences is given, and the accuracy and convergence behavior of five methods based on thermodynamic integration and perturbation techniques was evaluated using liquid water as a test system. Reasonably accurate entropy differences are obtained through thermodynamic integration in which many copies of a solute are desolvated. When only one solute molecule is involved, only two methods seem to yield useful results, the calculation of solute-solvent entropy through thermodynamic integration, and the calculation of solvation entropy through the temperature derivative of the corresponding free-energy difference. One-step perturbation methods seem unsuitable to obtain entropy estimates.     

Dynamic properties of water around a protein-DNA complex from molecular dynamics simulations

Formation of protein-DNA complex is an important step in regulation of genes in living organisms. One important issue in this problem is the role played by water in mediating the protein-DNA interactions. In this work, we have carried out atomistic molecular dynamics simulations to explore the heterogeneous dynamics of water molecules present in different regions around a complex formed between the DNA binding domain of human TRF1 protein and a telomeric DNA. It is demonstrated that such heterogeneous water motions around the complex are correlated with the relaxation time scales of hydrogen bonds formed by those water molecules with the protein and DNA. The calculations reveal the existence of a fraction of extraordinarily restricted water molecules forming a highly rigid thin layer in between the binding motifs of the protein and DNA. It is further proved that higher rigidity of water layers around the complex originates from more frequent reformations of broken water-water hydrogen bonds. Importantly, it is found that the formation of the complex affects the transverse and longitudinal degrees of freedom of surrounding water molecules in a nonuniform manner.     

Microstructure and dynamics in lyotropic liquid crystals. Principles and applications of nuclear spin relaxation

Nuclear spin relaxation of quadrupolar nuclei provides access to a wide range of properties of lyotropic liquid crystals, ranging from the molecular ordering and dynamics at the interface to the macroscopic viscoelastic behaviour. We emphasize here the unique capability of the spin relaxation method to provide detailed geometric and dynamic information relating to the microstructure of lyotropic liquid crystals, i.e. the metric, curvature, and fluctuations of the dividing interface that separates polar and non-polar regions. This information is conveyed to the spin system via the translational diffusion of surfactants or counterions over the interface. The general principles of the spin relaxation method, as applied to lyotropic liquid crystals, are described, with emphasis on the model-independent information content of the relaxation observables and on the relation to microstructure. Specific results for lamellar, hexagonal, cubic, and nematic phases are also described.     

Microstructure and dynamics in lyotropic liquid crystals. Principles and applications of nuclear spin relaxation

Abstract

Nuclear spin relaxation of quadrupolar nuclei provides access to a wide range of properties of lyotropic liquid crystals, ranging from the molecular ordering and dynamics at the interface to the macroscopic viscoelastic behaviour. We emphasize here the unique capability of the spin relaxation method to provide detailed geometric and dynamic information relating to the microstructure of lyotropic liquid crystals, i.e. the metric, curvature, and fluctuations of the dividing interface that separates polar and non-polar regions. This information is conveyed to the spin system via the translational diffusion of surfactants or counterions over the interface. The general principles of the spin relaxation method, as applied to lyotropic liquid crystals, are described, with emphasis on the model-independent information content of the relaxation observables and on the relation to microstructure. Specific results for lamellar, hexagonal, cubic, and nematic phases are also described.     

Micellar crystals in solution from molecular dynamics simulations

Polymers with both soluble and insoluble blocks typically self-assemble into micelles, which are aggregates of a finite number of polymers where the soluble blocks shield the insoluble ones from contact with the solvent. Upon increasing concentration, these micelles often form gels that exhibit crystalline order in many systems. In this paper, we present a study of both the dynamics and the equilibrium properties of micellar crystals of triblock polymers using molecular dynamics simulations. Our results show that equilibration of single micelle degrees of freedom and crystal formation occur by polymer transfer between micelles, a process that is described by transition state theory. Near the disordered (or melting) transition, bcc lattices are favored for all triblocks studied. Lattices with fcc ordering are also found but only at lower kinetic temperatures and for triblocks with short hydrophilic blocks. Our results lead to a number of theoretical considerations and suggest a range of implications to experimental systems with a particular emphasis on Pluronic polymers.     

Hydrogen nuclear spin relaxation in hydrogen-ice clathrate

H2 in D2O ice clathrate has been studied by hydrogen NMR. In a previous report, the H2 line shape was shown to be due to incompletely averaged intramolecular dipolar interactions. Here the relaxation times T1, T1rho, and T2 are reported. T1 passes through a minimum at 10 K, indicating that the rotational transition rate Gamma between the three sublevels of J = 1 passes through the resonance frequency at this temperature. On the cold side, T1 varies as T(-2.6); on the hot side, the rate T1(-1) varies as T(-2) plus a constant (due to paramagnetic impurities). These indicate a two-phonon process drives the rotational transitions Gamma. The spin-echo T2 is nearly independent of temperature and in reasonable agreement with the Van Vleck intermolecular H2-H2 second moment. T1rho deviates from the expected T1rho = T1 behavior above 85 K, revealing an additional slow-motion source of relaxation. The deviation is driven by the hopping of H2 between large cages. Ortho-para conversion is measured to be much slower in the clathrate than in the bulk solid, reflecting the greater distances between the H2 molecules.     

Kinetic aspects of the thermostatted growth of ice from supercooled water in simulations

In experiments, the growth rate of ice from supercooled water is seen to increase with the degree of supercooling, that is, the lower the temperature, the faster the crystallization takes place. In molecular dynamics simulations of the freezing process, however, the temperature is usually kept constant by means of a thermostat that artificially removes the heat released during the crystallization by scaling the velocities of the particles. This direct removal of energy from the system replaces a more realistic heat-conduction mechanism and is believed to be responsible for the curious observation that the thermostatted ice growth proceeds fastest near the melting point and more slowly at lower temperatures, which is exactly opposite to the experimental findings [M. A. Carignano, P. B. Shepson, and I. Szleifer, Mol. Phys. 103, 2957 (2005)]. This trend is explained by the diffusion and the reorientation of molecules in the liquid becoming the rate-determining steps for the crystal growth, both of which are slower at low temperatures. Yet, for a different set of simulations, a kinetic behavior analogous to the experimental finding has been reported [H. Nada and Y. Furukawa, J. Crystal Growth 283, 242 (2005)]. To clarify this apparent contradiction, we perform relatively long simulations of the TIP4P/Ice model in an extended range of temperatures. The temperature dependence of the thermostatted ice growth is seen to be more complex than was previously reported: The crystallization process is very slow close to the melting point at 270 K, where the thermodynamic driving force for the phase transition is weak. On lowering the temperature, the growth rate initially increases, but displays a maximum near 260 K. At even lower temperatures, the freezing process slows down again due to the reduced diffusivity in the liquid. The velocity of the thermostatted melting process, in contrast, shows a monotonic increase upon raising the temperature beyond the normal melting point. In this case, the effects of the increasing thermodynamic driving force and the faster diffusion at higher temperatures reinforce each other. In the context of this study, we also report data for the diffusion coefficient as a function of temperature for the water models TIP4P/Ice and TIP4P/2005.