Anisotropic Rotational Diffusion Studied by ... - ACS Publications

Laboratory Experiment pubs.acs.org/jchemeduc

Anisotropic Rotational Diffusion Studied by Nuclear Spin Relaxation and Molecular Dynamics Simulation: An Undergraduate Physical Chemistry Laboratory Michael M. Fuson* Department of Chemistry and Biochemistry, Denison University, Granville, Ohio 43023, United States S Supporting Information *

ABSTRACT: Laboratories studying the anisotropic rotational diffusion of bromobenzene using nuclear spin relaxation and molecular dynamics simulations are described. For many undergraduates, visualizing molecular motion is challenging. Undergraduates rarely encounter laboratories that directly assess molecular motion, and so the concept remains an abstraction. Spin relaxation and molecular dynamics simulations are complementary tools in studying molecular motion in solution. Simulations allow direct visualization of motion. Comparison of simulation results to experimental spin relaxation results enables discussion of the strengths and weaknesses of each technique. Performance of spin relaxation experiments in addition broadens undergraduate perceptions of the utility of NMR spectroscopy. Laboratories which enable the study of the motion of small molecules by both techniques are presented.

KEYWORDS: Upper-Division Undergraduate, Laboratory Instruction, Physical Chemistry, Kinetic-Molecular Theory, Molecular Mechanics/Dynamics, NMR Spectroscopy, Computational Chemistry, Liquids, Transport Properties, Hands-On Learning/Manipulatives

O

Nuclear spin relaxation rates depend directly on the rotational motion of the molecule, can be computed directly from a MD simulation, and can be determined experimentally. While there is an example in the chemical education literature of measuring 13C spin−lattice relaxation times (T1) and relating them to correlation times for CH vector reorientation (for segmental motion of 1-hexanol10), no examples have been published which exploit the connection to MD simulations. Here, laboratories are described which investigate the anisotropic molecular motion of bromobenzene dissolved in chloroform. Bromobenzene is chosen as the subject of investigation because it is a rigid molecule (aside from molecular vibrations), and this greatly simplifies discussion of its rotational motion when compared to that of 1-hexanol, for example. Further, the ability to measure T1 values associated with the motion of multiple CH spin pairs, which are rigidly oriented relative to each other, allows experimental characterization of the anisotropy of its rotational motion, which is also novel in this journal.

ne of the fundamental challenges facing undergraduate chemistry students as they seek to make the transition from novice to expert chemists is to understand and deeply internalize the fundamental basis of kinetic molecular theory: molecules in motion. For 100 years, chemists have used physical models to help visualize the molecular realm. Over the last 20 years, 3D graphics have become integrated into teaching materials. Both of these are enormously useful in training students to think of molecules as discrete objects with particular structures that are related to their properties and behavior. A full comprehension of the dynamic aspects of molecular behavior remains less common.1 Realistic molecular dynamics (MD) simulations have long been recognized as a potentially important tool for teaching but, until recently, remained too computationally demanding (expensive) to be readily available. With the development of sophisticated MD programs such as NAMD2 and Gromacs3 that are freely available and run efficiently on desktop computers, use of MD simulations in teaching is becoming more widespread. Published examples include simulations of water using Amber,4 of pentane using Tinker,5 and of a dipeptide using NAMD,6 a decapetide using Abalone,7 and a 12 amino acid peptide using Gromacs.8 While these examples are interested in comparing computed properties of the system to experimental quantities, they focus on structure and energy. Recently, a paper has appeared using MD to calculate translational diffusion rates of hexane using Amber.9 © XXXX American Chemical Society and Division of Chemical Education, Inc.

Received: September 11, 2016 Revised: February 7, 2017

A

DOI: 10.1021/acs.jchemed.6b00699 J. Chem. Educ. XXXX, XXX, XXX−XXX

Journal of Chemical Education

■

Laboratory Experiment

BACKGROUND

surrounded on all sides with replicas of itself). Under periodic boundary conditions, if a molecule diffuses out of one side of the box, it simultaneously appears on the opposite side of the box. After the system is constructed, typically, the next step is an energy minimization process to remove close contacts and to begin closing holes left by omitting partially overlapping molecules during the generation of the initial solution structure. The molecular dynamics simulation then generates a “trajectory” containing the positions of all the atoms as a function of time. This is begun by assigning each atom a random velocity whose root mean squared magnitude is consistent with the desired temperature. Based their positions, the forces and thus the accelerations on each atom are calculated and then new positions for the atoms after a brief time step (typically 1 or 2 fs) are calculated. This process is then repeated over and over to build up the trajectory. Typically, the initial part of the trajectory is discarded to allow for “equilibration” of the system, which is to allow any artifacts introduced by the choices of initial positions and velocities to be forgotten. The time required for the simulation (other than the efficiency of the software and hardware) is a strong function of the number of atoms in the system and also depends on the length of the trajectory. Even today, MD simulations longer than 1 μs are unusual. After the simulation is complete, the motion of the atoms can of course be directly visualized. In addition, a wide variety of properties can be computed from the trajectory, sometimes explicitly as a function of time and other times looking at probability distributions or averages.

NMR Spin−Lattice Relaxation

The fundamentals of NMR spin−lattice relaxation and its measurement have appeared in this journal previously.10 As shown in Figure 1, the inversion recovery method for

Figure 1. Schematic representation of the inversion recovery experiment.

measuring T1 (a) begins with magnetization at equilibrium, aligned with the external field; (b) after a 180° pulse, the spin magnetization is inverted; (c) waiting a variable amount of time (delay τ) allows for partial recovery toward equilibrium; (d) giving a 90° pulse rotates the magnetization into the transverse plane (perpendicular to the external field); (e) it precesses in the transverse plane where it may be sampled to acquire a FID; and finally, (f) the FID is Fourier transformed to give the spectrum. The intensity of the observed peak varies as a function of τ according to eq 1: M(τ ) = M(∞)(1 − (1 + α)e−τ / T1)

Connecting MD and Spin Relaxation

As has been stated above, the T1 values of a molecule contain information about the motion of the molecule. The source of this information can be shown by considering the mechanism causing spin−lattice relaxation.16,17 For a 13C with a directly bonded proton, the dominant relaxation mechanism is the dipole−dipole interaction of two spins placed in a static magnetic field. (Other mechanisms must be minimized which includes dipole−dipole interactions with more remote spins and chemical shift anisotropy, which can be important at high fields.) The strength of such an interaction turns out to be dependent on the orientation of the internuclear vector between spins i and j relative to the static field (the angle θ in Figure 2). As a molecule moves in solution, the angle θ will fluctuate, and the strength of the dipolar interaction will fluctuate. If this fluctuation has frequency components at or near the Larmor frequency (ω0) of the spin of interest, this will

(1)

where M(∞) is the magnetization after an infinite delay or the equilibrium magnetization and α is the fractional inversion (i.e., degree of incomplete inversion due to instrumental imperfections). Molecular Dynamics Simulations

Considering the importance of molecular modeling and molecular dynamics simulations to modern biochemistry, surprisingly little introduction to these methods is given in modern physical chemistry textbooks. A concise introduction is given in Grant and Richards’ Computational Chemistry from the Oxford Chemistry Primers series.11 To summarize briefly, the process begins by generating an initial structure either by direct computation or by importing a structure from an external resource (either a computational chemistry program such as Chem3D12 or Spartan13 or, for biomolecules, a database such as the Worldwide Protein Data Bank14). The energy associated with the structure is based on a force field from which an energy for each bond, bond angle, torsion angle, intramolecular nonbonded interaction, and intermolecular nonbonded interaction based on atom types, atom positions, and potential functions (for example, harmonic potentials for bond stretching or Lennard-Jones potentials for nonbonded interactions) can be computed. Here, we use the OPLS-AA force field.15 For a liquid-state simulation, the molecule of interest is surrounded by solvent molecules and the effect of a macroscopic system is achieved by imposing periodic boundary conditions (that is, a small “box” of molecules is

Figure 2. Orientation of a pair of nuclei, i and j, relative to the static magnetic field, B0. B




The diffusion rates used above can be computed from correlation times which in turn are computed directly from the MD simulations. The correlation time is the integral of the correlation function between the orientation of an appropriately chosen unit vector at time t, êi(t), and its orientation at a later time, t + τ:19

allow transfer of energy between the spin system and the motion of the molecules, leading to spin relaxation. Following Woessner (who generalized the analysis to account for anisotropic motion),18 we have 2⎛ ℏ⎞ 1 9 ⎛ μ ⎞ γγ i j = R1 = ⎜ 0 ⎟ ⎜⎜ 3 ⎟⎟ (J1(ω0) + J2 (2ω0)) T1 8 ⎝ 4π ⎠ ⎝ rij ⎠ 2

(2)

τc =

where μ0 is the permeability of free space, γi is the magnetogyric ratio of nuclei i, rij is the internuclear distance between i and j, and Jk(kω0) is the spectral density at frequency kω0 (k = 1 or 2). For the simplest case of isotropic motion and in the limit of extreme narrowing (ω02τc2 ≪ 1), J1(ω0) + J2(2ω0) = 4/3 τc, and thus the relaxation rate is directly proportional to the correlation time of the CH bond vector, τc. (Correlation times are discussed in more detail below. See eq 4.) For bromobenzene, we will model the rotational motion of the molecule as the rotational diffusion of a symmetric top with the symmetry axis collinear with the CBr bond (see Figure 3).

∫0

∞

1 ⟨3(eî (t ) ·eî (t + τ ))2 − 1⟩dτ 2

(4)

where the integrand is the correlation function and the angle brackets stand for an average over all possible starting times, t, and over all molecules. The dot product between the unit vectors, êi(t)·êi(t + τ), yields the cosine of the angle between the orientation of the vector at time t and its orientation at time t + τ. For a model of isotropic motion, the unit vectors are simply aligned with the relevant internuclear vectors. For anisotropic motion, unit vectors aligned with a molecule fixed frame of reference are chosen (see Figure 3). A short correlation time means rapid motion of that vector. Somewhat confusingly, the diffusion rate with respect to a given axis is the rate of motion about that axis as opposed to the rate of motion of that axis as characterized by the correlation time. Thus, the relationships between correlation times and diffusion rates are not simple reciprocal relations. Relations expressing these Cartesian correlation times to the diffusion rates are given by Fuson et al.19 Solving these equations for the diffusion rates yields Dx =

Dy = Figure 3. Axes of rotation defined with respect to bromobenzene. The origin of the coordinate system is at the center of mass, the z axis is aligned with the CBr bond, and the y axis is in the molecular plane.

Dz =

■

While the molecule’s C2v symmetry does not mean that the moments of inertia with respect to the two axes perpendicular to the symmetry axis are equal, we can qualitatively justify this approximation by noting that diffusion about axes which reorient the heavy Br atom (x and y) will be significantly slower than about the z axis. Examining this approximation is one of the advantages of pairing this experimental work with the MD simulation. Under the symmetric top approximation, D∥ = Dz is the diffusion rate around the symmetry axis and D⊥ = Dx = Dy is the diffusion rates perpendicular to the symmetry axis. We then can write17 4 J1(ω0) + J2 (2ω0) = (Aτ1 + Bτ2 + Cτ3) (3) 3

τy + τz − 5τx 6(τx2

+

τy2

+

τy2

+

τy2

+ τz2 − 2τxτy − 2τyτz − 2τzτx)

(5)

τx + τz − 5τy 6(τx2

+ τz2 − 2τxτy − 2τyτz − 2τzτx) τx + τy − 5τz

6(τx2

+ τz2 − 2τxτy − 2τyτz − 2τzτx)

EXPERIMENTAL SECTION

NMR T1 Experiments

Samples were prepared in advance and consisted of 50% by volume bromobenzene in CDCl3 sealed in 5 mm NMR tubes after degassing with three freeze−pump−thaw cycles. Proton decoupled 13C inversion recovery experiments were performed on a Bruker Avance 400 spectrometer using TopSpin 2.1 software. Twelve delay times were used with a maximum value of 75 s (approximately 5T1). Collecting eight transients per spectrum led to S/N of approximately 1000 for the fully relaxed spectra and required 2.5 h for the complete experiment. Rather than using TopSpin’s built-in T1 calculation, peak amplitudes were read manually and the data were fit in Mathematica20 using NonLinearModelFit. This is a pedagogical choice designed to make the students interact more closely with the data and to allow control over the type of error estimates obtained. D∥ and D⊥ were calculated from the T1 values of the meta and para carbons.

where 1/τ1 = 6D⊥, 1/τ2 = D∥ + 5D⊥, and 1/τ3 = 4D∥ + 2D⊥. Furthermore, the constants A, B, and C depend on the angle of the CH internuclear vector to the symmetry axis, φ,18 which differ for the meta and para CH vectors of bromobenzene (120 and 180°, respectively) and are defined as 1 A = 4 (3 cos2 φ − 1)2 , B = 3 cos2 φ(1 − cos2 φ)2, and

MD Simulations

3

C = 4 (cos2 φ − 1)2 . Measurement of relaxation times for these two CH pairs with differing values of φ allows determination of D∥ and D⊥ through solving two equations in two unknowns.

Simulations were run using Gromacs 5.1 running on Linux based personal computers. Structure files of bromobenzene and an equilibrated CDCl3 solvent box were prepared in advance. Students construct a 4.0 nm cube containing the appropriate C




mole ratio of solute and solvent (e.g., 183 bromobenzene and 240 CDCl3). After several equilibration steps, a 200 ps trajectory was calculated using the OPLS-AA force field11 with a 1 fs step size and NVT conditions (constant number of molecules, volume, and temperature). Lengths of equilibration steps and the final trajectory were chosen to minimize computation time while not affecting the final correlation times obtained by more than 5%. NVT conditions were chosen to most nearly replicate the environment of a sealed NMR tube under thermostatic control. (Strictly speaking, linear momentum is conserved only in simulations run under NVE conditions, and so simulations probing dynamics would be best run under these conditions. In simulations run under NVT or NPT conditions, linear momentum is not conserved but is instead coupled to “bath variables”, which reasonably leads to concerns about whether the dynamical properties of the system are effected. In practice, however, the selection of different ensembles often leads to identical results.) Atomic coordinates were stored every 100 fs. Correlation functions were calculated from the atomic coordinates of the trajectories using a Fortran program provided by the instructor. Correlation times can be obtained by direct numerical integration of the correlation function. The diffusion rates are then calculated from the correlation times using expressions given in eq 5. All input files required for the MD simulation are provided in the Supporting Information.

Typical results for the diffusion constants (or rates) obtained both from the spin relaxation results and from the MD simulation are given in Table 2. For the NMR data, Table 2. Rotational Diffusion Constants for Bromobenzene Dissolved in CDCl3a

a

HAZARDS CDCl3 and bromobenzene are both irritants for skin, eye, and respiratory system. CDCl3 is a suspected carcinogen. Strong magnetic fields are present in the NMR laboratory.

■

RESULTS AND DISCUSSION Typical results for the 13C T1 values for bromobenzene dissolved in CDCl3 (50% by volume) are given in Table 1. Table 1. 13C T1 Values for Bromobenzene Dissolved in CDCl3a location

a

T1/s 8.7 11.7 11.2 8.0

± ± ± ±

NMR

MD

Dx Dy Dz

47 ± 2 47 ± 2 104 ± 11

55 ± 4 32 ± 2 90 ± 6

Errors are 95% confidence intervals.

uncertainties are propagated from the T1 uncertainties. The uncertainties of the values derived from the MD runs are based on the 95% confidence intervals of 10 independent runs. (In the experiment as described, each student calculates only one trajectory and so cannot generate similar confidence intervals for the correlation times from the MD runs. If desired, confidence intervals could be estimated from pooled class data if students perform independent simulations. How to do this is discussed in the Instructor Notes.) Quite reasonable agreement can be seen in the magnitudes of the diffusion constants obtained by both techniques. This agreement can be difficult to obtain since the results from MD simulations are very sensitive to the quality of the force field. Larger diffusion constants mean more rapid diffusion with respect to the specified axis. A molecule “forgets” its initial orientation on a time scale on the order of 1/D. As predicted, we obtain larger values for Dz as compared to Dx and Dy. The values of Dx and Dy from the MD simulations are both significantly smaller than Dz, somewhat validating the approximation of the motion as a symmetric top, as discussed above. At the same time, they are significantly different from each other, illustrating the limits of the axially symmetric model. Translational diffusion constants are also calculated from the MD simulations. This allows students to compare rotational and translational diffusion and deepen their understanding of molecular motion. An interesting extension of this lab would be to add an experimental measurement of the translational diffusion constant using pulsed field gradient NMR methods.22 These laboratories have been used in physical chemistry and advanced topics classes with lab sections of 8−15 junior and senior undergraduate students. Four 3 h laboratory periods were used: one for spin relaxation data acquisition, one for spin relaxation analysis and calculations, one for running the MD simulations, and the final period for analysis of the MD simulations. However, it would be easy to reduce the number of laboratory periods required to two, as discussed in the Instructor Notes. These two laboratories represent two out of six laboratories that we do as the lab portion of a four credit hour lecture/lab course. Students who are familiar with 1H and 13C NMR data acquisition can run the inversion recovery experiment fairly independently, but often we take the approach of the instructor running the experiment with groups of students and offering extensive commentary about NMR methods as this is done. The spin relaxation calculations are sufficiently challenging that we give students a fairly extensive set of instructions. Likewise, for the MD simulations, we provide students an explicit set of instructions because the software is complex and oriented to the command line. Understanding student mistakes and being

■

C1 C2 (ortho) C3 (meta) C3 (para)

D/ns−1

1.6 0.3 0.4 0.3

Errors are 95% confidence intervals.

Errors are 95% confidence intervals. If our assumption that the rotational motion of bromobenzene can be modeled as a symmetric ellipsoid is correct, we would expect that T1 for the ortho and meta carbons would be nearly identical, whereas that for the para carbon would be different; we indeed see this pattern. (We might expect small differences in the relaxation of the ortho and meta carbons due to differing numbers of neighboring protons. Based on distances and the r−6 dependence of the interaction strengths, we would expect neighboring protons to have 1.6% of the impact on the relaxation time as the directly bonded proton. In our analysis, we have assumed that the relaxation is due solely to interactions with directly bonded protons.) Another point of interest is that while C1 has no directly bonded proton and thus little relaxation due to the dipole−dipole interaction, it is fairly efficiently relaxed by the uncommon mechanism of scalar coupling to the quadrupolar nuclei 79Br and 81Br.21 D




Molecular Dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781− 1802. (3) Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. GROMACS: A Message-Passing Parallel Molecular Dynamics Implementation. Comput. Phys. Commun. 1995, 91, 43−56. (4) Speer, O. F.; Wengerter, B. C.; Taylor, R. S. Molecular Dynamics Simulations of Simple Liquids. J. Chem. Educ. 2004, 81, 1330−1332. (5) Galembeck, S. E.; Caramori, G. F.; Romero, J. R. A New Exploration of the Torsional Energy Surface of n-Pentane Using Molecular Models and Molecular Modeling Software. J. Chem. Educ. 2005, 82, 1800−1804. (6) Carlotto, S.; Zerbetto, M. Computational Study of Environmental Effects on Torsional Free Energy Surface of N-Acetyl-N′-Methyl-LAlanylamide Dipeptide. J. Chem. Educ. 2014, 91, 96−102. (7) Spitznagel, B.; Pritchett, P. R.; Messina, T. C.; Goadrich, M.; Rodriguez, J. An Undergraduate Laboratory Activity on Molecular Dynamics Simulations. Biochem. Mol. Biol. Educ. 2016, 44, 130−139. (8) Rodrigues, J. P. G. L. M.; Melquiond, A. S. J.; Bonvin, M. J. J. Molecular dynamics characterization of the conformational landscape of small peptides: A series of hands-on collaborative practical sessions for undergraduate students. Biochem. Mol. Biol. Educ. 2016, 44, 160− 167. (9) Eckler, L. H.; Nee, M. J. A Simple Molecular Dynamics Lab To Calculate Viscosity as a Function of Temperature. J. Chem. Educ. 2016, 93, 927−931. (10) Gasyna, Z.; Jurkiewicz, A. Determination of Spin-Lattice Relaxation Time Using 13C NMR. J. Chem. Educ. 2004, 81, 1038− 1039. (11) Grant, G. H.; Richards, W. G. Computational Chemistry; Oxford University Press: New York, 1995; pp 32−58. (12) Chem3D; PerkinElmer Informatics, 940 Winter Street, Waltham, MA 02451. (13) Spartan; Wavefunction, Inc., 18401 Von Karmen, Suite 370, Irvine, CA 92715. (14) Protein Data Bank: Biological Magnetic Resonance Bank; http://www.wwpdb.org/ (accessed January 2017). (15) Kaminski, G.; Friesner, R. A.; Tirado-Rives, J.; Jorgensen, W. L. Evaluation and Reparametrization of the OPLS-AA Force Field for Proteins via Comparison with Accurate Quantum Chemical Calculations on Peptides. J. Phys. Chem. B 2001, 105, 6474. (16) Hore, P. J. Nuclear Magnetic Resonance; Oxford University Press: New York, 1995; pp 56−71. (17) Keeler, J. Understanding NMR Spectroscopy, 2nd ed.; Wiley: Chichester, UK, 2010; pp 241−274. (18) Woessner, D. E. Nuclear Spin Relaxation in Ellipsoids Undergoing Rotational Brownian Motion. J. Chem. Phys. 1962, 37, 647−654. (19) Fuson, M. M.; Brown, M. S.; Grant, D. M.; Evans, G. T. Cartesian Correlation Times for NMR AX2 Spin Systems. J. Am. Chem. Soc. 1985, 107, 6695−6698. (20) Mathematica; Wolfram Research, 100 Trade Center Drive, Champaign, IL 61820. (21) Levy, G. C. Carbon-13 Spin-Lattice Relaxation: CarbonBromine Scalar and Dipole-Dipole Interactions. J. Chem. Soc., Chem. Commun. 1972, 1972, 352−354. (22) Harmon, J.; Coffman, C.; Villarrial, S.; Chabolla, S.; Heisel, K. A.; Krishnan, V. V. Determination of Molecular Self-Diffusion Coefficients Using Pulsed-Field-Gradient NMR: An Experiment for Undergraduate Physical Chemistry Laboratory. J. Chem. Educ. 2012, 89, 780−783.

prepared to respond to them does require that instructors are quite familiar with the MD software. For both laboratories, student access to and familiarity with software such as Mathematica is required for nonlinear curve fitting, solving systems of equations, and some of the more complex error propagation. When offered as part of the physical chemistry sequence, at least as it is taught at our institution, these laboratories are largely independent of content covered in the lecture portion of the course. While they draw on fundamental ideas of kinetic molecular theory and NMR spectroscopy, spin relaxation and MD simulations are introduced only in the context of these laboratories and represent an opportunity to broaden the experience of students. In addition, such concepts as relaxation times, correlation functions, and anisotropic motion are new and engaging for students. Using these laboratories as part of an advanced topics class offers the chance for much closer alignment of the course content with the laboratories. In either context, students find these laboratories both challenging and rewarding. They are asked to use instruments in unfamiliar ways, to interact with unfamiliar software, and to think about new concepts. Many students respond with pride at having risen to the challenges and with fascination at the new insights on how to think about molecules in motion.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available on the ACS Publications website at DOI: 10.1021/acs.jchemed.6b00699. Instructor Notes (PDF, DOCX) Student laboratory handouts for molecular dynamics simulation of bromobenzene (PDF, DOCX) Student laboratory handouts for spin−lattice relaxation of bromobenzene (PDF, DOCX) Input files for MD simulations and Fortran program to calculate correlation functions from trajectory output (ZIP)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Michael M. Fuson: 0000-0002-0900-9250 Notes

The author declares no competing financial interest.

■

ACKNOWLEDGMENTS The work here is directly derived from an exercise used by my graduate advisor, James Prestegard, many years ago. Much of the initial development of these laboratories was tested by my students Randy Ruffner and Yuqing Liang. The work described here has been supported by Denison University and the National Science Foundation (MRI 0215955).

■

REFERENCES

(1) Jones, L. L. How Multimedia-Based Learning and Molecular Visualization Change the Landscape of Chemical Education Research. J. Chem. Educ. 2013, 90, 1571−1576. (2) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kalé, L.; Schulten, K. Scalable E


Anisotropic Rotational Diffusion Studied by ... - ACS Publications

Recommend Documents