代写CHEM10600、MATLAB程序语言代做

日期：2023-10-27 10:48

EXPERIMENT 20

PROGRAMMING IN MATLAB

Aims and Intended Learning Objectives

Use the MATLAB software to manipulate matrices, plot graphs and write programs.

By the end of this experiment you will be able to:

Perform matrix multiplication in MATLAB, using as an example the rotation of a

molecule in space, and determining the normal modes

Plot graphs in MATLAB, using as an example the analysis of the kinetics of a reaction

that is first order in each reactant

Skills & Techniques Developed

During this experiment you will utilise the following skills and/or techniques:

MATLAB software

Matrix multiplication, eigenvalues and eigenvectors

Graph plotting

Programming and problem solving

Aligned Lecture Course

CHEM10212 (kinetics)

CHEM20611 (molecular vibration)

CHEM20212 (molecular mechanics)

This experiment builds on previous teaching lab experiments:

CHEM10600 (MATLAB)

Safety

No CRA form is required as this experiment is entirely computational.

Make sure that you take regular breaks from the screen.

Preparation

Read the Introduction and Theory pdf. This includes going through chapters 4-6 of the

online ebook.

Pre-lab Quiz: Complete the pre-lab quiz (ensure you have looked over the experiment

before attempting).

Copy and save all relevant output, m files and graphs to a single word file named

E20_Firstname_Lastname. Write down relevant working, notes etc in your lab book, as

you would for a practical experiment.

2

EXPERIMENTAL

Part 1. Molecular Rotation

Molecules are dynamic and can undergo translation, rotation and vibration. Here we will use

matrix multiplication to calculate the coordinates of atoms of carbon dioxide upon rotation.

In two dimensions the coordinates of atoms of the CO2 molecule (shown below) may be

represented as three vectors: [

−1

0

], [

0

0

] and [

1

0

] (assuming a bond length of 1 for simplicity).

1.1 In Matlab, open a new m-file called rot.m and do the following:

Create a vector v for the oxygen atom at [

1

0

] (Chapter 2, page 13)

Set the variable theta equal to π/4 radians (i.e. 45°) (Hint: pi is already a variable in

Matlab).

Create a rotation matrix R in terms of the variable theta (see Introduction and Chapter

2, page 14)

Multiply the matrix and vector to generate the new coordinates vnew = R v (note the

order of multiplication matters). Your output vector should be [

0.7071

0.7071].

1.2 Repeat the above process, but rotate the molecule 2π (i.e. 360°) in steps of π/6 (30°). To

do this make the following m-file called rotco2.m (Chapter 3, page 31) with the text below, fill

in the blanks (...), and add axis labels to the graph (see Introduction for more information

about plotting). Copy the resulting graph (it will be saved as figure1.png) into your word

document and add an appropriate figure caption to describe it.

max = 12; %max is the number of points to calculate, 2 in steps of /6

v = ...; % v is the initial vector of the atom coordinates

figure(1); % create a figure to plot to

for i = 0:max % Loop (until ‘end’), i is the iterator & steps from 0 to max

theta = ...; % theta is the angle,

R = ...; % R is the rotation matrix,

vnew = ...; % vnew is the set of coordinates at angle theta

fprintf ('Angle %3.1f, Coordinates %5.2f %5.2f\n', theta, vnew);

hold on; % this allows data to be plotted each cycle of the loop

plot(vnew(1), vnew(2), ‘ko’); % plot the point as a black(k) circle(o)

end

print(‘figure1.png’, ‘-dpng’) % Save as an image in the current folder

1.3 Now rotate all three atoms in CO2 by making the following modifications to your code:

Set v = [

−1 0 1

0 0 0

] which now contains the coordinates of all three atoms.

Add four more fields %5.2f to the fprintf statement after “Coordinates”.

Execute the completed m-file and copy your code and the output to your word

document. [You do not need to generate any graphics output for this part.]

3

Part 2: Normal modes

In three dimensions, linear molecules containing N atoms such as HCl and CO2 have 3N

degrees of freedom in total, 3 of which are translations (x, y and z) and 2 are rotations. This

leaves 3N−5 internal degrees of freedom, otherwise known as vibrational modes.

We start by considering HCl in one dimension.

2.1 How many translational, rotational and vibrational modes will there be for HCl in two and

one dimension?

As outlined in the introduction the mass-weighted Hessian, Hm, matrix for HCl is

𝐇𝐦 = [

𝑘⁄𝑚H − 𝑘⁄√𝑚H𝑚Cl

− 𝑘⁄√𝑚H𝑚Cl 𝑘⁄𝑚Cl

]

(1)

2.2. In Matlab generate a new script called HCl.m, define variables for 𝑘 = 480 N m−1 and the

masses of H and Cl in atomic mass units (g mol−1

), then create the array for 𝐇𝐦 shown above.

Find the eigenvalues 𝜆 and eigenvectors of 𝐇𝐦 (Chapter 4, page 53).

2.3 Either separately, or by adding to your script, calculate the vibrational frequency of each

normal mode in wavenumbers (𝜈̅, cm−1)

𝜈̅=

1

2𝜋𝑐

√𝜆 × 1000𝑁A

(2)

where c is the speed of light (cm s−1) and 1000 NA converts from g mol−1 to kg.

2.4 Describe the normal mode specified by each eigenvector.

2.5 Compare your values of 𝜈̅with the vibrational frequency of HCl calculated using

𝜈̅=

1

2𝜋𝑐 √

𝑘

𝜇

(3)

where 𝜇 is the reduced mass, and

1

𝜇

=

1

𝑚H

+

1

𝑚Cl

.

2.6 By changing the mass of the elements in your script HCl.m calculate the vibrational

frequencies of D35Cl and D37Cl. (D = deuterium = 2H). Report these values in your word file.

Three atoms vibrating in one dimension:

Using CO2 (O1=C=O2) as an example, the Hessian and Mass matrices are:

𝐇 = [

𝑘 −𝑘 0

−𝑘 2𝑘 −𝑘

0 −𝑘 𝑘

] 𝐌 = [

𝑚O 0 0

0 𝑚C 0

0 0 𝑚O

]

(4)

2.7 In Matlab, create a CO2.m script. Define variables for 𝑘 = 1600 N m−1 and the masses of

O and C in atomic mass units, then create arrays for 𝐇 and M. Generate the mass-weighted

Hessian matrix, Hm = M-½HM-½ and obtain the eigenvalues 𝜆 and eigenvectors of 𝐇𝐦.

Calculate the vibrational frequency of each normal mode in wavenumbers (cm−1) with

Equation (2), describe each normal mode and explain which will be infra-red active. Report

these values in your word file and copy your m script.

4

Part 3. Graph Plotting: Chemical Kinetics

The esterification of ethanol (A) by acetic anhydride (B) to produce ethyl acetate (P) in carbon

tetrachloride proceeds as

A + B → P + H2O

This reaction is first order with respect to each reactant i.e.

𝑣 = 𝑘[A][B] (5)

where 𝑣 = −𝑑[A]/𝑑𝑡 is the rate of reaction, 𝑘 is the rate constant in units of M-1

s

-1

, and 𝑡 is

time. The integrated rate equation is:

[A]

[B]

=

[A]0

[B]0

exp(−𝑘𝑡([B]0 − [A]0

))

(6)

The reaction was performed with starting concentrations [A]0 = 0.096 M and [B]0 = 0.137 M

at 55C. Observations of [A] were made hourly for 16 hours and are recorded in Table 1.

Table 1. Concentration of ethanol [A] over time during esterification with acetic anhydride [B].

Time / h Ethanol [A] / M

0 0.096

1 0.079

2 0.067

3 0.057

4 0.050

5 0.044

6 0.039

7 0.034

8 0.031

9 0.028

10 0.025

11 0.023

12 0.021

13 0.019

14 0.017

15 0.016

16 0.015

3.1 Create a new m file and load the data for [A] from Table 1 into an array a (Double check

that you have the correct number of values!)

i.e. a = [0.096, 0.079, 0.067, 0.057, 0.050, 0.044, 0.039, ...];

3.2 Create an array h of time values in hours. (Hint: to create an array h of integers from 0 to

10 (inclusive) you could use: h = 0:10; ).

3.3 Use [A]0, [A] and [B]0 to calculate the array of [B] values in an array (Hint: [P] = [A]0 −

[A] = [B]0 − [B] ).

3.4 Create a figure and plot of [A] and [B] against time in hours as circles (not as a line) on the

same axis. We plot data like these as discrete points because each point is a measurement.

[Plotting your data as a line would imply that there were an infinite number of measurements

taken over your time range – this isn’t true!] Your plot should look similar to Figure 1 below.

Save the plot and insert it into your word document with an appropriate caption.

5

Figure 1. Reactant concentrations against time for the reaction of ethanol (A, red circles) and acetic

anhydride (B, blue circles).

3.5 The integrated rate Equation (6)states that [A]/[B] should have an exponential decay over

time if it is first order with respect to both reactants. Add lines to your script to calculate

[A]/[B] as a function of time and assign these values to the array ab (Hint: you need to divide

the two arrays in an element-by-element fashion, i.e. ab = rdivide(a, b); ).

3.6 Add to your script to create a new figure and plot [A]/[B] against time with your data

points as circles, not as a line. Make sure to add axis labels, then save the plot and insert it

into your word document with an appropriate caption. Comment on whether the values

decay exponentially by inspection.

3.7 The rate constant 𝑘 for this reaction can be calculated by taking the logarithm of Equation

(6) to give

ln (

[A]

[B]

) = ln (

[A]0

[B]0

) − 𝑘𝑡([B]0 − [A]0

)

(7)

The resulting equation is of the form 𝑦 = 𝑚𝑥 + 𝑐, where x = t, 𝑚 = −𝑘([B]0 − [A]0

) is

the gradient, and 𝑐 = ln([A]0/[B]0

) is the intercept with the y axis. Add code to calculate

ln([A]/[B]) and assign it to the vector logab, then convert your time values in h to seconds

and assign the result to the vector t. Plot ln([A]/[B]) against time in seconds. Confirm that

your data is (roughly) linear, save the labelled plot in your word document with an appropriate

caption.

3.8 You should see that the data is not perfectly linear as a function of time - experimental

data is rarely perfect! To obtain a value for 𝑘, we need to fit a straight line to the data – you

might know this as the line of best fit. Matlab can obtain the gradient and y-axis intercept of

this line using the polyfit function.

p = polyfit(t,logab,1);

This returns a vector p, the first element of which is the gradient of the line of best fit, and

the second element is the y-axis intercept. Use this gradient in Equation (7) to calculate a

6

value for 𝑘 and print it to the screen with two significant figures. Confirm the intercept with

the y-axis is ln([A]0/[B]0

).

3.9 To check the quality of this fit, we should plot it on top of the experimental values. First

calculate the values of the line of best fit from the straight line equation 𝑦 = 𝑚𝑡 + 𝑐, using

your values of 𝑡 from experiment, and the values of 𝑚 and 𝑐 from the polyfit function. Assign

the resulting values to the array sline. Now replot your experimental data as before along with

the values of sline. For clarity, plot logab as circles and sline as a dashed line (see Introduction

for details). Confirm there is a good fit which provides further evidence that the reaction

obeys rate equation (6). Make sure to add axis labels, and a legend with appropriate names

for the two lines. Save the plot and insert it into your word document with an appropriate

caption. Copy the final complete script into your document too.

Part 4: Hydrogen Atomic Orbitals

The atomic orbitals for a hydrogen atom can be determined by solving the Schrodinger

equation (8). In this case the eigenvectors and eigenvalues describe the form of the orbital

and its energy.

The Schrodinger wave equation for the hydrogen atom is given in its simplest form by:

𝑯𝝍 = 𝐸𝝍 (8)

The eigenfunctions can be written as two components:

𝜓𝑛,𝑙,𝑚𝑙

(𝑟, 𝜃, 𝜑) = 𝑅𝑛,𝑙(𝑟)𝑌𝑙

𝑚𝑙

(𝜃,𝜑) (9)

For a hydrogen atom (Z = 1), and with r, the distance from the nucleus, expressed in atomic

units (1 a0 = 52.9 pm) the calculated radial functions, R(r) are as follows.

Orbital R(r)

1s 2𝑒

−𝑟

2s 1

2√2

(2 − 𝑟)𝑒

−𝑟/2

2p 1

2√6

𝑟𝑒

−𝑟/2

3s 2

81√3

(27 − 18𝑟 + 2𝑟

2

)𝑒

−𝑟/3

3p 4

81√6

(6𝑟 − 𝑟

2

)𝑒

−𝑟/3

These are converted into radial distribution density plots using the formula

𝑟

2𝑅(𝑟)

2

(9)

7

4.1 Write a Matlab script to plot the 1s radial distribution function for distances of r = 0 to 20

a0 in steps of 0.1 a0. The plot generated should be similar to that shown below

4.2 Add to your script so that the rdd for the 1s atomic orbital and at least one other of the

orbitals listed above are displayed on the same graph. Include in your write-up a copy of your

script and a suitably annotated figure of the resulting plot.

4.3 Write a script that asks for values of the n and l quantum numbers of the above orbitals,

and in response plots the corresponding RDD for that atomic orbital.

Assessment

After this experiment you will be assessed on:

Correctness of code.

Correct and properly formatted graphs and numerical answers.

Understanding, accuracy and interpretation of the work.

Use of your lab book to plan, prepare and maintain a record of the experiment.