CHEM10600 PROGRAMMING IN MATLAB

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: [

] (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)

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:

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.

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

) 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)
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.
WX：codehelp