1 Rotational Diffusion Propagator of the ... - ACS Publications

Subscriber access provided by READING UNIV

Article

Rotational Diffusion Propagator of the Intramolecular ProtonProton Vector in Liquid Water: A Molecular Dynamics Study Wisman Acharige Monika Madhavi, M. S. Samantha Weerasinghe, and Konstantin I Momot J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.7b07551 • Publication Date (Web): 14 Nov 2017 Downloaded from http://pubs.acs.org on November 16, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 53

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Rotational Diffusion Propagator of the Intramolecular Proton-Proton Vector in Liquid Water: A Molecular Dynamics Study

W. A. Monika Madhavi1,2, M. S. S. Weerasinghe3, Konstantin I. Momot1*

1

School of Chemistry, Physics and Mechanical Engineering, Queensland University of Technology (QUT), GPO Box 2434, Brisbane, Qld 4001, Australia 2

Department of Physics, University of Colombo, Colombo 03, Sri Lanka

3

Department of Chemistry, University of Colombo, Colombo 03, Sri Lanka

*Corresponding Author:

Konstantin I. Momot School of Chemistry, Physics and Mechanical Engineering Queensland University of Technology (QUT) GPO Box 2434 Brisbane, Qld 4001 Australia

E-mail address:

[email protected]

Tel.:

+61 7 3138 1173

Fax:

+61 7 3138 9079

1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 53

ABSTRACT

Rotational motion of water molecules plays the dominant role in determining NMR spin relaxation properties of liquid water and many biological tissues. The traditional theory of NMR spin relaxation predominantly uses the assumption that the reorientational dynamics of water molecules is described by a continuous-time rotational-diffusion random walk with a single rotational diffusion coefficient. However, recent experimental and theoretical studies have demonstrated that water reorientation occurs via large, discrete angle jumps superimposed on a continuous-time rotational diffusion process. We have investigated the rotational diffusion propagator of the proton-proton (H−H) vector of water molecules in liquid water at 298 K using molecular dynamics simulations. Analysis of the MD-simulated reorientational trajectories reveals that reorientation of the intramolecular H–H vector occurs through a combination of the two mechanisms, rotational diffusion proper and discrete largeangle jumps. We demonstrate that, empirically, the rotational diffusion propagator of the intramolecular H−H vector in liquid water can be described in terms of multiple rotational diffusion coefficients. A model with two rotational diffusion coefficients was found to provide a reasonable (albeit imperfect) fit of the MD-simulated propagator on the time scales relevant to NMR spin relaxation near room or physiological temperature (picoseconds to nanoseconds). We report the apparent values of the two rotational diffusion coefficients determined from the propagator analysis at 298 K and discuss their physical meaning.


Page 3 of 53

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


INTRODUCTION Molecular dynamics of liquid water, including rotational dynamics, plays a crucial role in fundamental biological processes such as ion transport and protein folding.1–4 Rotational dynamics of water molecules also provides insights into water’s hydrogen bonding network, which in turn is considered to be the source of most of the vital physical and chemical properties of water.5–8 Experimental methods used to study rotational motion of water include nuclear magnetic resonance spin relaxation measurements,9–11 time-resolved vibrational spectroscopy12 and quasi-elastic neutron scattering.13–15 While the temporal or the spatial resolution of these techniques is not sufficient to probe the detailed dynamics of individual molecules, they can yield various ensemble-averaged motional correlation functions characterising the rotational motion. It is generally acknowledged that molecular dynamics of water is complex, and its comprehensive understanding requires consideration of a wide range of temporal and spatial scales.4,16,17 The original theoretical model of reorientational dynamics of liquid water is the continuous rotational diffusion model proposed by Debye and later developed by Perrin.18,19 In this model, rotational diffusion is treated as a continuous-time Brownian process representing the stochastic rotation of a sphere in a uniform viscous fluid. More recent studies, however, have suggested that the Debye model is not capable of fully explaining the water rotational dynamics, and there is both experimental and theoretical evidence that water reorientation proceeds via a combination of Debye continuous rotational diffusion and large, discrete angular jumps of the water molecules.4,10,20,21 Jump model developed by Ivanov22 proposes that water reorientation takes place in finite-angle discrete jumps rather than as a continuous rotational diffusion process. Relatively recent studies have shown that, picosecond-scale reorientation dynamics of the OH bond,4,21 dynamical structure factor of QENS data,23 water dynamics in hydration shells of hydrophilic and hydrophobic aqueous solutions,24 NMR 3 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 53

relaxation data in supercooled water10 and molecular dynamics simulations results25 are better explained by jump-like models of water reorientation. Reorientation of water molecules through large angular jumps is observed in MD simulations by analysing the single molecular trajectories at high time resolution (femtosecond time scale).8 Nuclear magnetic resonance (NMR) spin relaxation has been extensively used to investigate water dynamics under a wide range of physical, chemical and biological conditions.1,4,9–11,26– 28

The phenomenon of NMR spin relaxation originates in the modulation of the spin

Hamiltonian by random thermal motions of molecules in a liquid or gas. These motions typically encompass rotational tumbling, translational motion or chemical exchange.25,29 The longitudinal NMR spin relaxation rate constant (1/T1) of 1H nuclei in water (S = ½) is almost entirely determined by intramolecular and intermolecular dipolar interactions between the 1H spins, with the intramolecular contribution typically being the dominant one. As the intramolecular H−H distance is on average constant on the time scale of the Larmor precession frequency, the intramolecular component of (and therefore the main contribution to) the relaxation rate arises entirely from the rotational motion of the molecules. Despite the advances in the theoretical and experimental understanding of water rotational dynamics over the last 1-2 decades, in the NMR literature the treatment of molecular reorientation in liquids remains based largely on Debye’s original continuous rotational diffusion model.30 Conversely, most of the latest studies of the rotational dynamics of molecules in liquids are not NMR-specific and are focused on general motional autocorrelation functions.8,9,21,31,32 The rotational diffusion propagator, which serves as the key to the description of NMR spin relaxation, has not been investigated in detail outside of the Debye rotational-diffusion regime. Therefore, in this study we have aimed to investigate the implications of “nonDebye” dynamics for the spin autocorrelation functions and the full rotational diffusion


Page 5 of 53

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


propagator dependent on both spatial and temporal dimensions that determine 1H NMR spin relaxation rates. Molecular dynamics (MD) simulations provide a platform to study the detailed dynamics of individual water molecules, while also enabling the calculation of measurable macroscopic quantities related to spin relaxation. Therefore, MD is a computational tool widely used to understand the dynamics of liquid water.9–11,25,28,33–36 In turn, the dynamics of water and the resulting spin-relaxation properties enable the interrogation of the structure and organisation of materials and biological tissues.37–44 As the long-term aim of this research is to understand the NMR spin relaxation properties of water in detail, the study has focused on the dynamics of the H−H vector of the water molecule on the time scales relevant to 1H spin relaxation at high magnetic field (1−20 T). In this paper, we present our study of the rotational dynamics of the intramolecular proton-proton vectors of liquid water using molecular dynamics simulations at non-extreme temperatures and pressures. We have tested several popular allatom rigid water models developed for MD simulations and compared their ability to reproduce the experimental NMR 1H longitudinal spin relaxation data. We then used the bestperforming water model (TIP4P/2005) for the detailed study of the rotational diffusion propagator of the intramolecular H−H vector.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 53

THEORY Theory of spin relaxation The theoretical model accounting for the NMR relaxation of nuclear spins was proposed by Bloembergen, Purcell, and Pound45 and further developed by Redfield.30 This model is semiclassical, whereby the spins are treated quantum-mechanically and the molecular motion is treated classically. In water with the natural-abundance mix of the hydrogen isotopes, almost every water molecule contains two 1H nuclei (S = ½). As discussed above, intermolecular and intramolecular dipole-dipole interactions between neighbouring 1H nuclei represent the dominant source of 1H spin relaxation in liquid water. For any given 1H−1H pair, the dipolar spin Hamiltonian can be represented as a combination of stochastically modulated irreducible spherical-tensor operators of rank 2 and the orders ranging from -2 to +2. The dipolar spin Hamiltonian scales with the internuclear distance, r, as 1/r3; for this reason, the relatively short-range intramolecular interactions dominate the longitudinal spin relaxation in liquid water, while the longer-range intermolecular interactions typically account for approximately 1/3 of the measured 1/T1.30 The rate constants of spin relaxation due to dipolar interaction can be represented as a combination of the contributions from spectral densities of motion at the multiples of the Larmor frequency, = :

{

}

1 3 4 2 = γ h I ( I + 1) J (1) (ω0 ) + J ( 2) (2ω0 ) T1 2

(1)

1 15 3 3  = γ 4h2 I ( I + 1) J (0) (0) + J (1) (ω0 ) + J ( 2) (2ω0 ) T2 4 8 8 

(2)

where J (ω) is the spectral density associated with the irreducible spherical-tensor operator (q)

of the order q sampled at the frequency ω.


Page 7 of 53

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Theory of rotational diffusion The standard model used to describe the rotational motion of molecules in the context of NMR spin relaxation in liquids is the Debye model, which considers molecular reorientation as analogous to the continuous-time rotational diffusion of a sphere in a viscous fluid. The process is described by a rotational diffusion propagator, , which is found as a solution of the rotational diffusion equation:30 ∂Ψ = D∆ s Ψ ∂t

(3)

where Δ is the Laplacian operator on the surface of a sphere and D is the rotational diffusion coefficient in rad2 s−1. The latter can be estimated, under the assumption of non-slip boundary conditions and a large hydrodynamic radius a, using the Stokes’ formula:

D=

kT 8π a3η

(4)

Rotational diffusion propagator is the probability density of finding the molecule in the orientation Ω at time t when it is known that the molecule was oriented at Ω0 at the zero time. The propagator can be represented as a spherical-harmonics expansion: ∞

Ψ (Ω, Ω 0 , t ) = ∑

l

∑ Υl*,m (Ω 0 )Υl ,m (Ω )e

−

t

τl

(5)

l =0 m=−l

where Ω0 is the initial orientation of the molecule (at t = 0); Ω is its orientation at time t; and the rotational correlation times τl are related to the rotational diffusion coefficient:

1

τl

= Dl(l + 1)

(6)

In the coordinate frame where the initial H−H vector was aligned with the z axis (i.e. θ0 = 0), integration of the propagator over the azimuthal angle (φ) gives the probability density that at



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 53

time t the H−H vector has undergone a given angular displacement α from its initial orientation. In the general coordinate system, if the H−H vector of a water molecule was at the polar angle θ0 at t = 0 and found at θ at time t, and the angle between the two orientations is α, the orientational correlation functions 〈 , 〉 (where is the l-th rank Legendre polynomial) can be calculated by considering the orthogonality property of spherical harmonics. In the case of a single rotational diffusion coefficient and the propagator given by Eq. (5), this correlation function is a single exponential function: Pl ( cos α , t ) =

1 4π

2π π 2π π

∫ ∫ ∫ ∫ P ( cos α ) Ψ ( Ω , Ω , t ) sin θ sin θ l

0

0

dθ d φ d θ 0 d φ0

0 0 0 0

(7) =e

− l ( l +1) Dt


Page 9 of 53

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


METHODS

Selection of a water model for MD simulations: Comparison between simulations and experiment As the first step, the performance of several all-atom rigid water models in reproducing experimental spin-relaxation data was investigated. Five commonly used water models were selected: TIP3P, SPC, SPCE, TIP4P and TIP4P/2005. For each one, CHARMM-compatible46 topology and parameter files were prepared based on the parameter values found in the literature.47–49 Parameter files for the TIP4P/2005 water model49 developed by J. Faraudo were taken from http://catalan.quim.ucm.es/frames_en.html.

MD Simulations For each water model investigated, a cubic simulation box containing 10,000 water molecules was constructed. Prior to the MD production run, each system was energy minimised for 40 ps in 2 fs time steps in a NVT ensemble at the target temperature and equilibrated in a NPT ensemble at 1 atm. The equilibration was achieved in one or more equilibration runs, each of length 100 ps in 2 fs time steps. At the end of each run, fluctuation of the total energy of the system over the course of the run was measured. The equilibration run was extended in increments of 100 ps until the total energy of the system reached an asymptotic minimum with a relative standard deviation less than 1% in order to ensure that each system reached both the thermal equilibrium and the minimum-energy configuration prior to the production MD run. Therefore, the total duration of equilibration run differed depending on the model. TIP3P and SPC models equilibrated in 200 ps while SPCE, TIP4P and TIP4P/2005 models needed 240, 300 and 500 ps, respectively, for equilibration.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 53

Production runs of the MD simulations were performed in a NPT ensemble at 1 atm pressure for the duration of 10 ns. The Newton equations of motion were solved with a time step of 1 fs, and electrostatic interactions were updated with a 10 fs time step. The atomic coordinates were saved every 0.1 ps. All bonds between heavy atoms and hydrogen atoms were maintained rigid during all the simulations using SHAKE.50 The cut-off distance for nonbonding Lennard-Jones interactions was 1.2 nm, with a switching function starting at 1.0 nm. The Ewald summation method was employed for electrostatic interactions with a grid size of 1 Å. Constant temperatures were maintained using a Langevin thermostat51 with a relaxation constant of 1 ps−1 , and constant pressure was attained using a Langevin piston52 with the decay constant 50 ps−1. For each water model, the simulations were performed at five equidistantly spaced temperatures ranging from 278 K to 318 K. As a test to ensure the absence of simulation artifacts due to the thermostat and/or the piston, a molecular dynamics simulation of 10,000 TIP4P/2005 water molecules at 298 K was performed without the temperature and pressure controls (NVE ensemble). Additionally, in order to probe the rotational dynamics at a higher time resolution, a 100 ps MD simulation of TIP4P/2005 water at 298 K temperature and 1 atm pressure was performed with the atomic coordinates saved every time step (1 fs). All MD simulations were performed using NAMD53 version 2.10. The simulations were executed on a supercomputing cluster containing 3780 processors; the simulation jobs utilised 64 parallel processors. The typical CPU time for each production run was 400 processor-hours and the typical clock time was 24 hours.


Page 11 of 53

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Measurement of the intramolecular contribution to the spin relaxation rates Out of the simulated 10,000 water molecules, a subset of 500 molecules was selected as tracer molecules. These molecules were distributed approximately uniformly in the simulation volume. This was done by an in-house tcl script written using VMD, the visualisation component of the NAMD software.54 The initial configuration of the simulation box was divided into 512 (8×8×8) cells and out of that 500 cells were randomly selected. One water molecule from each selected cell was taken as a tracer molecule. For each tracer molecule, the Redfield stochastic functions F(q), defined as the coefficients of the rank-2 irreducible spin tensor operators comprising the dipole-dipole interaction Hamiltonian,30 were calculated at each 0.1 ps time step. The spectral densities J(q) used in Eqs. (1) and (2) were calculated using the Wiener-Khinchin theorem55 as the squared-absolute-values of the spectrum of the stochastic functions F(q)(t). The spectral densities were averaged over the 500 tracer molecules to obtain the ensemble-averaged spectral densities of reorientational motion. The corresponding longitudinal (1/T1) and transverse (1/T2) spin relaxation rates were then calculated using Eqs. (1) and (2).

Intermolecular contribution to the spin relaxation rate constants For TIP4P/2005 water at 298 K, the contribution from intermolecular interactions to the spin relaxation rates was investigated. 200 hydrogen atoms (belonging to 200 different water molecules), which were approximately uniformly distributed within the simulation volume, were selected as the tracer atoms using the same tcl code as for the treatment of the intramolecular relaxation contribution. For each tracer atom, all the hydrogen atoms within a sphere of 5 Å radius from that atom were considered; at each 0.1 ps time step, the H–H vectors were calculated as were the values of the Redfield stochastic functions F(q) for the H–



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 53

H vectors. The net value of each function F(q) corresponding to a selected tracer atom was calculated by summation of all the individual stochastic functions corresponding to the individual H−H vectors considered. The spectral densities of motion and the spin relaxation rates were calculated according to the procedure described in the previous paragraph. The resulting value of the intermolecular contribution to the spin relaxation rate constant was compared with the intramolecular contribution.

Fitting the simulated diffusion propagator with theoretical diffusion propagator models MD simulation trajectories of TIP4P/2005 water at 298 K sampled at the time resolution of 0.1 ps were used to investigate the numerical representation of the rotational diffusion propagator of water. Due to the isotropic nature of liquid water, the full 3D rotational diffusion propagator given by Eq. (5) was reduced to the uniangular propagator Ψ′ , , given by Eq. (8), which depends only on the absolute angular displacement α between the initial time and time t. The value of α at t = 0 is by definition zero. Therefore, only the zerothorder spherical harmonics (m = 0) contribute to the uniangular propagator. The uniangular propagator was obtained by integration of Eq. (5) and has the form ∞

Ψ ' (α , t ) = 2π ∑ Υl*,0 (0 )Υl , 0 (α )e

−

t

τl

(8)

l =0

For each water molecule, the angular displacement α was defined as the angle that the H−H vector at time t makes with that vector’s direction at t = 0, as illustrated in Figure 1.


Page 13 of 53

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1: Snapshots of the trajectory of a water molecule at t = 0, t = 5 ps and t = 10 ps. The path of the tip of the H−H vector is shown in dots, and the orientation of the H−H vector at 5 ps and 10 ps is shown as the blue arrows. The angle α was defined as the angle that the H−H vector makes at t = tk with the direction of the same vector at t = 0.

The orientations of the H−H vectors of all the 10,000 water molecules in the simulation box were sampled from the MD simulation trajectories at equidistant times (∆t = 0.1 ps) for the first 100 ps of the simulation. The angular displacements of the intramolecular H−H vectors were calculated for each time step (k), each molecule (m) as

αm ( tk ) = arccos[nm (t0 ) ⋅ nm (tk )]

(9)

where nm(t0) and nm(tk) are the unit vectors describing the orientation of the intramolecular H−H vector of the m-th molecule at the initial time and time tk, respectively.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 53

At each time step, the resulting 10,000 α values were grouped into 50 equally-sized bins ranging from 0 to π, and the histograms were constructed. The histograms were normalised such that, at any given time point tk, the properly normalised continuous probability density function of molecular orientations formed an envelope of the respective histogram. The 10 ns long MD simulation trajectory was split into blocks of 100 ps length. For the purposes of propagator analysis, each block was considered as a separate trajectory, and the time evolution of the H–H vector and the diffusion propagator histograms were obtained for each individual block. The full uniangular propagator was then constructed by averaging all the histograms obtained from the individual 100 ps blocks. The uniangular rotational diffusion propagator Ψ′ , obtained from the MD simulations was least-squares fitted using Levenberg-Marquardt method56,57 with the theoretical propagator model given by Eq. (8) using the in-built Mathematica (Wolfram Inc) function “NonLinearModelFit”. The model was truncated at l = 30 in order to reduce the computational cost while also minimising the short-time oscillations.

Fitting of the angular and temporal cross-sections of the uniangular propagator The full uniangular propagator Ψ′(α, t) sampled from the MD trajectories formed a threedimensional plot, with α and t being the independent variables and Ψ being the dependent variable. Its slices in the angular and the temporal dimensions were selected and least-squares fitted with the uniangular propagator model containing a single diffusion coefficient. The corresponding plots of the fit residuals were visually examined to evaluate the quality of the fit. The number of rotational diffusion coefficients was gradually increased, and the fitting repeated, until an apparent near-random distribution of the residuals was obtained.


Page 15 of 53

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Fitting of the full uniangular propagator Prior to least-squares fitting the MD-simulated data with the uniangular diffusion propagator model given by Eq. (8), the accuracy of the fitting method was tested with a synthetic numerically simulated data set. A set of trajectories of virtual particles undergoing a Brownian random walk on the surface of a sphere with a single rotational diffusion coefficient was numerically generated. The numerical rotational diffusion propagator was extracted from the trajectories of these virtual particles. The time dependence of the ensemble-averaged Legendre polynomial P1(t) given by Eq. (7) (l = 1) and the uniangular diffusion propagator Ψ′(α, t) were calculated from the synthetic data. These data were leastsquares fitted with the models given by Eqs. (7) and (8), respectively. When tested on the synthetic data, the fitting of the uniangular propagator Ψ′(α, t) has exhibited poor convergence and failed to reproduce the known input value of the rotational diffusion coefficient when the entire block of data (α = 0…π and t = 0…tmax) was used for the fitting. However, the fitting returned the correct value of the diffusion coefficient when only the data points corresponding to large angular displacements which were not too close to π were considered: a block of data corresponding α ranging from π/2 to 2π/3 was selected, and the segment was least-squares fitted with the model of Eq. (8). Fitting of the data within this angular displacement range was capable of reproducing the correct value of the rotational diffusion coefficient; this type of block-fitting was therefore used as the starting point for the fitting procedure for the actual MD data. The data points corresponding to larger and smaller angles were progressively introduced into the data set being fitted as described below, while checking the uniformity of the fit residuals. The sensitivity of the fit parameters to the initial guess values was then investigated. Initially, an initial guess remote from the optimum was used, and the fit parameters were obtained; these were not necessarily physically meaningful. Remote guess values were adjusted until a 15 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 53

physically meaningful fit was obtained. The fitting was then repeated with several different initial guesses in order to ascertain the initial-guess range where the optimum at which the fit converged was independent of the initial guess. To examine the quality of the fit, a “non-uniformity of the fit” was defined. For each data set being fitted, the space of the independent variables (t and α) was divided into four equal segments {(t < α <

! !

), (t >

& α <

! !

& α >

), (t

! !

), (t >

&

)}, and the mean and the standard deviation of

the residuals was calculated in each segment. The difference between the largest and the smallest mean, divided by the standard deviation of the entire set of residuals, was taken as the “non-uniformity of the fit” (ζ). The fitting of the MD-simulated uniangular diffusion propagator was performed in iterative steps as follows. Initially, only the block of data corresponding to large values of the angle α (from π/2 to 2π/3) was selected, and the segment was least-squares fitted with the model containing two D values. The fit was considered to have converged when the relative change in the sum of the squares of the residuals between iterative steps fell below the precision goal specified (10−8). Data points corresponding to the smaller and larger values of the angle α were progressively added to the fitted data set, and the fitting was repeated until the convergence criterion was met. This procedure was extended until the entire data range (0 < t < 100 ps and 0 < α < π) was included in the fitted data set. In each instance, values of the fit parameters obtained in the previous fitting step were used as the initial guess for the fitting of the current dataset. Visual examination of the distribution of the fit residuals, as well as the value of ζ, were used to evaluate the quality of the fit. When the plot of the residuals started to exhibit noticeable systematic patterns, the number of the diffusion coefficients was gradually increased and the 16 ACS Paragon Plus Environment

Page 17 of 53

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


fitting was repeated until no systematic patterns in the residuals plot were visually observed. The number of the diffusion coefficients was incremented in one of two ways: by fixing the values of the diffusion coefficients already found and allowing their corresponding amplitudes to vary, and by allowing all parameters to vary.

Fitting of the ensemble-averaged Legendre polynomials In order to further characterise the rotational dynamics of the water molecules, we have sampled the Legendre polynomials Pl (α m (tk )) for each molecule m, at each time point tk (in 0.1 ps timesteps) for l = 1 − 4. The ensemble-averaged Legendre polynomials were then calculated by averaging the values of Pl (α m (tk )) over the set of 10,000 molecules (index m) for given values of l, k. The procedure was repeated for the one-hundred 100 ps blocks of the MD trajectory, as described above for the analysis of the uniangular diffusion propagator; the Legendre polynomial functions were then additionally averaged over the one-hundred 100 ps blocks. The four Legendre polynomial functions (l = 1 − 4) were individually least-squares fitted with a mono-exponential decay function. Plots of the fit residuals were used to evaluate the quality of each fit. The ratios of the time constants of the four exponential decay functions were compared with the theoretical ratios given by Eq. (7). The procedure was repeated for the data obtained from the NVE simulation.

Investigation of individual molecular trajectories The time dependence of angular displacement α (defined according to Eq. (9) and Fig. 1) of individual water molecules was sampled for all the molecules in the simulation box over a 17 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 53

100 ps time course at the time resolution of 0.1 ps. This data was used for the analysis of angular jumps exhibited by the molecules. A “single jump” was defined as a continuous sequence of time steps during which the value of α for a given molecule changed in the same direction. The magnitudes of all such jumps were recorded for each individual molecule. A combined histogram of the values of jump angles was then constructed for all the molecules in the simulation box. Events where large angular jumps were followed within 0.5 ps by a reversal of the change of α were excluded from the histogram.

Decay time constants of the Legendre polynomials In order to analyse the effect of discrete angular jumps on the decay time constants of the Legendre polynomials (see Eq. (7)), sets of synthetic data were simulated as follows. 5,000 virtual particles underwent a combined reorientational stochastic process whereby Debye rotational-diffusion random walk with a known diffusion coefficient was superimposed with discrete jumps of a given amplitude and probability per time step. Reorientational trajectories of these virtual particles were sampled, and the ensemble-averaged Legendre polynomials (l = 1 – 4) were calculated as a function of time. This simulation was repeated for a series of jump amplitudes and probabilities, as shown in Supporting Information. Following each simulation, each Legendre polynomial decay curve was least-squares fitted individually with a three-parameter monoexponential decay function. The decay time constants and their ratios were tabulated for each simulation (see Supporting Information).


Page 19 of 53

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


RESULTS Spin relaxation rate constants of water As the first step, the total longitudinal spin relaxation rate constant (comprising both the intermolecular and the intramolecular contributions) was calculated from the TIP4P/2005 MD data at T = 298 K in the limit of low Larmor frequency (ω0

1 Rotational Diffusion Propagator of the ... - ACS Publications

Recommend Documents