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 V x1 x2 x N 3 1 1 1 3 3 exp 6 d x1 d x2 d x N 6 R Q k T R x1 x2 x N B where Q is the classical partition function V x1 x2 x N 3 Q exp d x1 d 3 x2 d 3 x N k BT 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 1 1 2 2 2 V VOH RO1H1 VOH RO2 H 2 kOO ROO L0 k HOO 0 OOH 0 V 2 2 where VOH r DOH 1 exp r LH 2 kJ sin 1 2 HOO sin 1 2 OOH cos cos 0 V 80 mole DOH 360 L0 1 6 kJ mole LH 1 05 1 5 1 k 10 6 0 1 8 radians 0 1 7 radians pJ radian 2 kOO 300 J m2 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 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 xO2 y H V 2 1 exp k T dxO 2 dxH 1 dyH 1 dxH 2 dy H 2 dz H 2 R 6 B HH V xO2 2 y H1 exp dxO 2 dxH 1 dy H 1 dxH 2 dy H 2 dzH 2 k BT 1 6 RHH Qred 1 Qred In this molecule fixed axis system the expression for the dihedral angle is cos y H1 y H 2 y H1 y H2 2 z H2 2 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 using internal coordinates you are not very well constrained along the normal modes of the atom i e you have a small chance of making a large jump from one low energy position to another For example consider a water molecule with the O centered at 0 0 one H1 at …