10.34 – Fall 2006 Homework #10 Due Date: Friday, December 1st, 2006 – 9 AM ** Note: Please read the entire problem set before starting, there is important information throughout, even at the very end. For this problem, you do NOT need to have the Matlab code generate all of the results for part A-E by running it once. However, it should be able to take a temperature and number of points and generate all of the desired plots for that set of inputs. The most popular way to experimentally test a proposed geometrical structure for a large molecule (such as a protein) is by X-ray crystallography. However, some proteins are hard to crystallize; for these proteins, proposed geometrical structures can be tested using nuclear magnetic resonance (NMR). NMR measures the through-space magnetic coupling between two atoms which are not directly bonded to each other; this magnetic coupling is proportional to <1/R6>, where R is the distance between the two atoms. The symbol < > means the Boltzmann average over all the possible molecular geometries; in the classical limit and neglecting some minor complications due to the integral over the kinetic energy we can write: 1 ⎡− ( 1, x2,..., xN ) 3 3 31 Vx ⎤1 = ∫∫∫ 6 ⋅ exp ⎢ ⎥ dxdx 1 2...dxNR6 ⎤ kT Q ⎡⎣Rx( 1, x2,..., xN )⎦ ⎣ B ⎦ where Q is the classical partition function: ∫∫∫ ⎢⎡−Vx( 1, kTx2,..., xN )⎥⎤ 3132 3 NQ = exp dxdx ...dx ⎣ B ⎦ This high-dimensional integral can be computed for a proposed structure using Monte Carlo techniques. Of course for a molecule with a large number of atoms this can be quite challenging. Here we instead ask you to compute this integral for a small molecule. Note that it is very easy to figure out the equilibrium geometry from this analytical expression for V (note V = 0 at the equilibrium geometry). We suggest you use Metropolis’s method, and start your Monte Carlo steps from the equilibrium geometry. Write a set of Matlab functions which use Monte Carlo integration to compute <1/RHH6> at a given Temperature, where RHH is the distance between the two H atoms in HOOH. A. Determine the equilibrium structure at 0 K of HOOH by minimizing the potential energy of HOOH. Plot the structure of the molecule in 3-D using the plot3 command in Matlab. State the 0 K equilibrium values for <1/RHH6> and <RHH> in Angstroms. Do you expect the value of <RHH> to be different for T = 300 K, why? Cite as: William Green, Jr., course materials for 10.34 Numerical Methods Applied to Chemical Engineering, Fall 2006. MIT OpenCourseWare (http://ocw.mit.edu), Massachusetts Institute of Technology. Downloaded on [DD Month YYYY].B. Use your code to solve for the value of <1/RHH6> at 300 K. Report the value obtained for <1/RHH6>, the value of <RHH>, and the number of Monte Carlo steps attempted and accepted. C. Plot the 3-D location of all of the MC points obtained in the above simulation using the plot3 command again. It may also be instructive to plot the equilibrium structure underneath the MC points. This can be done with something like (obviously, this syntax will need to be modified to your problem): plot3(equil,’-‘,’linewidth’,4.0); hold on; plot3(MC_points,’.’); hold off; D. Repeat this to generate plots for temperatures of 600 K, 1000 K, and 5000 K. Generate these plots with a minimum of 10000 MC steps. Generate a histogram of the <RHH> values for each temperature, using the same x-axis scale for all figures. Also create a histogram (50 bins) showing the distribution of potential energies that the molecule adopts for each temperature (you don’t need to use the same x-axis scale). Find the bin with the largest frequency and compare this energy value with the value of kBT. E. Generate a plot for each T showing the evolution of the <1/RHH6> as the number of MC points increases (ideally this curve will converge to the actual value of as N Æ ∞ ). F. For your answer in part B, give your best guess at the uncertainty in your predicted value of <1/RHH6>, and explain how you derived it. Assume this is the expression for the potential energy of HOOH: = R +V (R + 2 k ⋅(R −L )+ 2 k ⋅⎣θ −θ) ( θ 0VVOH (O H 11 ) OH O H 2 2 ) 1 OO OO 02 1 θ ⎡(HOO 02 + OOH −θ where: V ()= DOH ⋅(1exp ⎡⎣−α⋅(r − LH )⎤⎦)2 r −OH kJV =80 ⋅sin (1 θ )⋅sin (1 θ )⋅⎡cos φ−cos ( ) ⎤φ mole 2 HOO 2 OOH ⎣() φ0 ⎦ kJ −1 JDOH =360 mole LH =1.05 Å α=1.5 Å kOO =300 m2 −6 pJL0 = 1.6 Å kθ= 10 radian 2 θ0 = 1.8 radians φ0 = 1.7 radians The R’s are Cartesian distances between the atoms. θHOO is the angle defined by H1-O1-O2 and θOOH is the angle defined by O1-O2-H2 (you can compute these using law of cosines). An expression for the dihedral angle φis given below. )2 ⎤⎦+Vφ Cite as: William Green, Jr., course materials for 10.34 Numerical Methods Applied to Chemical Engineering, Fall 2006. MIT OpenCourseWare (http://ocw.mit.edu), Massachusetts Institute of Technology. Downloaded on [DD Month YYYY].Hint: Molecule fixed axes Molecular potentials V do not depend on the position of the molecule in space, nor on its orientation, but only on the relative position of the atoms. Hence one can usually cut 6 degrees of freedom (corresponding to the position of the molecule and its angular orientation (Euler angles)) out of molecular problems. In this particular problem, we suggest using molecule-fixed axes where the position of atom O1 sets the origin, atom O2 lies on the x axis, and atom H1 lies in the x-y plane. Then one can remove these 6 degrees of freedom from the problem: (xO1,yO1,zO1,yO2,zO2,zH1). (You can set them all equal to zero). When you remove the orientational degrees of freedom you pick up some Jacobian volume elements; including these the new expression for the integral (again approximating away some minor terms related to rotational kinetic energy) is: 1 1 ⎛ xO 22 ⋅ 6 = ∫∫∫⎜ yH 61 ⋅ exp ⎡⎢ −V ⎤⎥⎞⎟ dxO 2 dxH1 dyH 1 dxH 2 dyH 2 dzH 2 R Q ⎜ kT ⎟HH red ⎝[RHH ]⎣ B ⎦⎠ ⎛ ⎡−V ⎤⎞Qred = ∫∫∫⎜ xO 2 ⋅ ⋅ exp ⎢ ⎥⎟ dx O2 dx H1 dy H1 dx H 2 dy H 2 dz H 2yH kT⎝ 2 1 ⎣ B ⎦⎠ In this molecule-fixed axis system, the expression for the dihedral angle is: yycos ()φ = H1 H2 ⋅ yH 22 + zH 22yH1 One other point of interest how to move around this multi-dimensional space, and the step sizes to take. When you are moving atoms around in 3-D Cartesian space (not