Rogers Computational Chemistry Using the PC

CURVE FITTING

83

In the multivariate case, one slope vector is not enough; a square slope matrix must be generated with dimensions equal to the number of independent unknowns xi one wishes to determine. The slope matrix is

M ¼	m11	m12	ð3-66Þ
	m21	m22

for a problem in two unknowns and larger for n unknowns. This is done by repeating the procedure just given for obtaining an m vector n times with different vectors y1; y2; . . . ; yp to produce an n n slope matrix M. Once having M, one can determine unknown values of x

x ¼ M 1y

ð3-67Þ

The n-fold procedure ðn > 2Þ produces an n-dimensional hyperplane in n þ 1 space. Lest this seem unnecessarily abstract, we may regard the n n slope matrix as the matrix establishing a calibration surface from which we may determine n unknowns xi by making n independent measurements yi. As a final generalization, it should be noted that the calibration surface need not be planar. It might, for example, be a curved surface that can be represented by a family of quadratic equations.

An illustrative example generates a 2 2 calibration matrix from which we can determine the concentrations x1 and x2 of dichromate and permanganate ions simultaneously by making spectrophotometric measurements y1 and y2 at different wavelengths on an aqueous mixture of the unknowns. The advantage of this simple two-component analytical problem in 3-space is that one can envision the plane representing absorbance A as a linear function of two concentration variables

A ¼ f ðx1; x2Þ.

Application: Simultaneous Analysis by Visual Spectrophotometry

Simultaneous Cr2O27 and MnO4 determination (Ewing, 1985) is a multivariate spectrophotometric analysis that requires determination of a matrix of four calibration constants, one for each unknown at each of two wavelengths,

a11x1 þ a12x2 ¼ A1

ð3-68Þ

a21x1 þ a22x2 ¼ A2

The elements aij are absorptivities (or are proportional to absorptivities, depending

	on the concentration units	and	cell dimensions), x ¼		x1	is the unknown
					x2
	concentration vector, and y	A1	is the absorbance	vector, observed at wave-
		A2
	lengths l1 and l2.	¼

The absorbance A is proportional to x through Beer’s law (see Computer

Project 2-1). The analytical problem is to solve the matrix equation
Ax ¼ y	ð3-69Þ

84	COMPUTATIONAL CHEMISTRY USING THE PC

for x from measured values of y once we have determined the matrix A. Be careful to distinguish between A and A. (For accepted nomenclature, see Ewing, 1985.) The wavelengths li should be chosen so as to make the A matrix as ‘‘nonsingular as possible,’’ that is, wavelengths should be selected so that, insofar as possible, absorbance by one species dominates all the rest. The wavelengths selected for this problem are 440 nm for Cr2O27 and 525 nm for MnO4 .

Dichromate-permanganate determination is an artificial problem because the matrix of coefficients can be obtained as the slopes of A vs. x from four univariate least squares regression treatments, one on solutions containing only Cr2O27 at 440 nm, one on the same solution at 525 nm, and one on solutions containing only MnO4 at each of these two wavelengths. We did this for five concentrations of each absorbing species and obtained the matrix

	¼	3:85 10 4	4:24 10 2	ð		Þ
A		4:39 10 3	2:98 10 3		3-70

Elements in the slope matrix A are proportional to absorptivities and concentrations are in parts per million. We shall take this as the ‘‘true’’ slope matrix.

Part 1. Generate the slope matrix from Eq. (3-65)

To obtain this matrix by the multivariate method, we first generate two absorptivity vectors a1j and a2j from a known concentration matrix in parts per million

0	27:0	13:0	1
B	53:0	8:65	C

BC

X ¼ B 80:0		4:33 C		ð3-71Þ
B	0:0	17:3	C
@	106	0:0	A

and the measured absorbance vectors yj, one at the lower wavelength and one at the higher wavelength. From the measured absorbance vector y ¼ f0:251; 0:149; 0:361; 0:049; 0:456g at 440 nm, we obtained the vector a1j ¼ f4:32 10 3; 2:72 10 3g rounded to the appropriate number of significant figures. This is the top row of the matrix A.

Mathcad				1
	0	27:0	13:0
		53:0	8:65
X :	B	80:0	4:33 C		XT	¼		53 27 80 0	106
	¼ B	0:0		C				8:65 13 4:33 17:3	0
	B		17:3 C
	B			C
	B	106	0:0	C
	B			C
	@			A

CURVE FITTING

85

¼

1:156

103

561:861

XT

X

2:117

104

1:156

103

¼

2:005 10 3

1:094 10 4

0

XT

X

1

0

5:32

10 5

1:094 10 4

:149

1

:568

1

:251

:401

y :

B

:361 C

y1 :

¼

B

:209 C

XT X

1

XTy

4:316

10 3

¼ B

C

B

C

3

B

:049 C

B

:740 C

Þ

¼

2:723

10

B

C

B

C ð

B

:456 C

B

:042 C

B

C

B

C

@

A

@

A

Repeating the process with a new measured absorbance vector {0.401, 0.568, 0.209, 0.740, 0.042} at 525 nm leads to the a2j vector y2 ¼ f3:85 10 4; 4:29 10 2g, the bottom row of the slope matrix. Together, they yield A

	¼	3:85 10 4	4:29 10 2	ð		Þ
A		4:32 10 3	2:72 10 3		3-72

Matrix (3-72) is essentially the same as matrix (3-70), but it is not exactly the same because it was obtained by the multivariate method from a different data set.

Part 2. Analyze unknown mixtures using the A matrix and two measurements of the absorbance, one at 440 nm and the other at 525 nm, as the y vector.

Having combined the two absorbance vectors into the absorbance matrix A, we are in a position to use A to solve for unknown concentration vectors x. Because y ¼ Ax, it follows that

x ¼ A 1y

Suppose that we have measured absorbance vectors y for three different solutions each containing both Cr2O27 and MnO4 . Let us call them ya, yb, and yc,

ya :¼	:401		yb :¼	:354		yc :¼	:723
	:331			:156			:177

where the measurement at 440 nm is the y1 (top) element in each y vector and the 525 measurement is the y2 (bottom) element in y. The solution concentrations in parts per million follow easily.

Mathcad

A :¼

:000385

:0429

ya :¼

:401

yb :¼

:354

yc :¼

:723

:004316

:002723

:331

:156

:177

86	8:71	COMPUTATIONAL CHEMISTRY USING THE PC
A 1 ya ¼		A 1 yb ¼	7:97	A 1 yc ¼	16:58
	71:2		31:11		30:55

The fourth significant figures in the vectors xb and xc are artifacts of the calculation. The concentration vectors should be reported as xa ¼ f71:2; 8:71g; xb ¼ f31:1; 7:97g, and xc ¼ f30:5; 16:6g ppm.

Error Analysis

Subtracting the slope matrix obtained by the multivariate least squares treatment from that obtained by univariate least squares slope matrix yields the error matrix

	0:0	0:50 10 3	ð		Þ
	0:07 10 3	0:26 10 3		3-73

where the univariate results are taken as ‘‘true’’ values.

Normally, one does not have ‘‘true’’ values of the elements of the slope matrix M for comparison. It is always possible, however, to obtain ^y, the vector of predicted y values at each of the known xi from any of the slope vectors m obtained by the

multivariate procedure
^y	¼	Xm	ð	3-74	Þ

This permits error analysis of that vector. (Note that the order Xm is necessary for the matrix and vector to be conformable for multiplication.) Repeating the procedure for all m vectors leads to error analysis of the entire matrix M.

We wish to carry out a procedure that is the multivariate analog to the analysis in the section on reliability of fitted parameters. A vector multiplied into its transpose gives a scalar that is the sum of squares of the elements in that vector. The ^y vector leads to a vector of residuals

e	¼	^y		y	¼	Xm		y	ð	3-75	Þ

The product eTe is the sum of squares of residuals from the vector of residuals. The variance is

s2	eTe	ð3-76Þ
	¼ p n 1

where p is the number of measurements made to establish n components of the slope vector. Under the assumption that the residuals are normally distributed, the best estimator of the variance of the ith element in m is s2dii, where dii is the ith diagonal element of D ¼ XTX 1. These variances are analogous to those obtained by Eqs. (3-34).

CURVE FITTING

87

For our sample data set of Cr2O27 and MnO4 absorbances, we seek the first of the two ^y vectors of residuals at 440 nm.

0 1

53:0 8:65

BC

B27:0 13:0 C

BC

X :¼ B 80:0 4:33 C

BC

BC

0:0 17:3 A@

106 0:0

0 1

0:252

BC

B0:152 C

BC

X m1 ¼ B 0:357 C

BC

BC @ 0:047 A

0:458

m1 :¼ :00432 :00272

This gives an error (residual) vector of

01

y :

:251

0

3 10 3

1

:361

¼ B

:149

1:488

10 3

:049 C

y

X

m1

¼ B

3:622

10 3

C

B

C

1:944

10

3

B

:456 C

B

1:92

3

C

B

C

B

10

C

@

A

B

C

@

A

which leads to

SSE :

y

X

m1 T

ð

y

X

m1

Þ

SSE

¼ ð

3:18

10 5

Þ

¼ ð

Þ

and an s2 value of

s2

eTe

¼

SSE

¼

3:18 10 5

¼

1:6

10 5

¼ p n 1

5 2 1

2

matrix D

T

X

1

The required inverse

5

¼

X

4 is

¼

1:09 10 4

2:01 10 3

D

5:32

10

1:09

10

which gives the uncertainties in row 1 of the slope matrix as

¼

ð

3:85 10 4

4:29 10 2

A

4:32

0:03Þ 10 3

ð2:72

0:18Þ

10 3

ð3-77Þ

ð3-78Þ

The reader is asked to find the standard deviations of the slopes of row 2 of the A matrix in Problem 9 below.

88	COMPUTATIONAL CHEMISTRY USING THE PC

Student’s t statistics (Rogers, 1983) follow in the usual way as do the 95%

confidence limits on the computed slopes, fð4:32 0:12Þ 10 3; ð2:72 0:77Þ 10 3g at 440 nm and fð3:85 2:7Þ 10 4; ð4:29 0:17Þ 10 2g at 525 nm.

These are not the same as the standard deviations due to the t statistic. The relative uncertainty on element a21 is large because the parameter is an order of magnitude smaller than the other elements in the slope matrix.

COMPUTER PROJECT 3-5 j

Let the generalization of Eq. set (3-55)

X

yi ¼ mjxij

Calibration Surfaces Not Passing

Through the Origin

ð3-79Þ

each contain a term, call it m0, with the stipulation that xi0 is some constant value, take it to be 1.0 for simplicity. The normal equations and the solution for the m vectors follow just as they did in the previous section except that each equation in set (3-58) contains an additive constant m0. The constant m0, a minimization parameter along with the rest of the mj, is the best estimator of the y intercept for a function not passing through the origin; it is the unique point at which the calibration surface cuts the y-axis.

To set up the problem for a microcomputer or Mathcad, one need only enter the input matrix with a 1.0 as each element of the 0th or leftmost column. Suitable modifications must be made in matrix and vector dimensions to accommodate matrices larger in one dimension than the X matrix of input data (3-56), and output vectors must be modified to contain one more minimization parameter than before, the intercept m0.

Procedure. The method can be tested using the matrix of concentrations, in micromoles per liter (mmol L 1), of tryptophan and tyrosine at 280 nm suitably modified to take into account constant absorption at 280 nm of some absorber that is neither tryptophan nor tyrosine

01

1:0	47:6	116
B 1:0	125	147 C

BC

X ¼ B 1:0		23:7	109 C	ð3-80Þ
B			C
@	1:0	156	48:3
	1:0	272	15:1 A

Columns 2 and 3 of matrix (3-80) are analogous to the concentration matrix (3-71). Using the absorbance vector {0.846, 1.121, 0.776, 0.599, 0.559}, compute the solution vector, m280. The first element of the solution vector is the intercept due to background absorption, and the second two elements are the absorbancies. What are the absorbancies of tyrosine and tryptophan at 280 nm by this method? Compare your results with the accepted values a280ðtyrÞ ¼ 1:28 103 and a280ðtryÞ ¼ 5:69 103 L mol 1 cm 1 (Eisenberg and Crothers, 1979). What are the units of your results for tyrosine and tryptophan? What is the intercept of absorbance at 280 nm due to compounds that are neither tyrosine nor tryptophan?

CURVE FITTING	89
COMPUTER PROJECT 3-6	j Bond Energies of Hydrocarbons

Determination of bond energies in hydrocarbons is a nontrivial example of multivariate analysis because lone C C and C H bonds cannot be observed in stable hydrocarbons; they always appear in groups. One could take the C H bond energy to be one-fourth of the energy of atomization of methane, subtract six times that value from the energy of atomization of ethane to get the C C single bond energy, and proceed in a like way, using the atomization energies of an alkene and an alkyne to generate the carbon-carbon double and triple bonds. This strategy would be risky at best because the C H bonds in a single molecule, methane, would be taken to represent all C H bonds, ethane would be taken to represent all C C bonds, and so on. It would be better to draw bond energies from a basis set of data for several, preferably many, molecules on the reasonable assumption that the mean result is more reliable than any single result from the set. Because the bonds cannot be observed singly, the problem is multivariate, and because we wish to generate a few bond energies from many experimental results, the input matrix will be overdetermined.

We will generate the energies for the carbon-hydrogen bond BCH and the carboncarbon single bond BCC using the five linear alkanes from ethane through hexane as

the five-member data base. The equation to be used is
hBCH þ sBCC ¼ Ha	ð3-81Þ

where h is the number of C H bonds in each hydrocarbon, s is the corresponding number of C C single bonds, and Ha is the enthalpy of atomization. Enthalpies of atomization of carbon and hydrogen were taken as 716.7 and 218.0 kJ per mole of atoms produced (Lewis et al., 1961) and were combined with the appropriate enthalpy of formation (Cox and Pilcher, 1970; Pedley et al., 1986; www.webbook. nist.gov) to obtain the enthalpy of atomization of each hydrocarbon by the method shown in Fig. 3-3 (see also Computer Project 2-1). The enthalpy of formation of methane, f H298ðmethaneÞ ¼ 74:5 kJ mol 1, is the enthalpy necessary to go from the elements in the standard state (0 by definition) to the molecule in the

					C atoms, H atoms
	4(218.0) =
		872.0 kJ/mol			C atoms
								1663
	716.7 kJ/mol							kJ/mol
0							methane
				–74.3 kJ/mol

Figure 3-3 Enthalpy Diagram for the Atomization of Methane.

90	COMPUTATIONAL CHEMISTRY USING THE PC

standard state at 298 K. In this illustration and in the computer project calculations that follow, discrepancies of 2 kJ mol 1 are not uncommon because of experimental uncertainty and differences among the various sources. We shall go into the difference between energy and enthalpy in this context in a later chapter.

Computation. Decide on an appropriate input matrix of bond numbers

0 h2		s2	1
	h1	s1
B ..		..	C
B hn		sn	C
B	. .		C
@			A

for ethane through hexane. The enthalpy of formation vector is { 20.24, 24.83,30.36, 35.10, 39.92} in the same order, where the units are kilocalories per mole. To convert from kilocalories per mole to kilojoules per mole, multiply by 4.184. Calculate the 5-fold enthalpy of atomization vector and the 2-fold vector of bond enthalpies. Obtain an error vector by comparing your result with the accepted values (Atkins, 1994) of 412 and 348 for the C H and C C bonds. respectively.

More than a reasonable number of significant figures is carried through the calculation to be rounded off at the end. When properly rounded, the uncertainty of the computed result should be reflected in the significant figures such that the rightmost digit is uncertain but no more than one uncertain digit is included in the final result. This is, of course, an approximate indicator of uncertainty; if a rigorous indicator is desired, the standard deviation, variance, or confidence level should be reported with the computed result.

COMPUTER PROJECT 3-7 j Expanding the Basis Set

Add ethylene, 1-propene, 1-butene, acetylene, and 1-propyne to the basis set. To do this, you must calculate five new atomization enthalpies from cycles similar to the one in Fig. 3-3. Also extend the input matrix to a 4 10 matrix. Generate the C H, C C, C C, and C C bond energies. Comment on the magnitude of the bond enthalpies, particularly the enthalpies of the C C single, double, and triple bonds. Is there any relationship between bond strength and bond energy for the three carbon-carbon bonds? Look up a set of accepted values for these bond energies and calculate a 4-fold error vector. Does the error for C H and C C get larger or smaller for the extended basis set as compared with the smaller basis set used in the first part of this experiment? Discuss this result.

PROBLEMS

1. Obtain the normal equations [Eq. set (3-63)] from the minimization conditions

q P di2 ¼ q P di2 ¼ 0 qm1 qm2

CURVE FITTING

91

2.Multiply XTy from Eq. set (3-61) to show that it is equal to the right side of Eq. set (3-63).

3.Multiply XTXm from Eq. set (3-61) to show that it is equal to the left side of Eq. set (3-63).

4.Can a rectangular matrix be both premultiplied and postmultiplied into its own transpose, or must multiplication be either preor postfor conformability? If multiplication must be either one or the other, which is it?

5.Show that eTe is the sum of squares of elements in the vector e ¼ f1; 2; 3g.

6.Does eTe commute with eeT ?

7.What is the average enthalpy of atomization of the four C H bonds in methane? Compare this value with the accepted value of the C H bond enthalpy.

8.Calculate the bond enthalpy of the C C bond in ethane using only the enthalpies of atomization of methane and ethane. Compare this result with the accepted result.

9.Find the standard deviations of the slopes in matrix (3-78) for row 2, which refers to absorbances measured at 525 nm.

10.When 6 moles of the substrate analog PALA combine with the enzyme ACTase, two things happen at the same time. The enzyme T unfolds to a more active form R

T ! R

and 6 moles of PALA bind to the enzyme. The measured enthalpy of both reactions together is

H ¼ 6 PALH þ T!RH ¼ 209:2 kJ mol 1

(Klotz and Rosenberg, 2000).

We would like to know the binding energy per mole of PALA and the enthalpy of the transformation T ! R, but we do not have enough information. Independent studies show that partial unfolding of ACTase occurs on binding

of less PALA, in particular, 1.8 mol of PALA cause 43% unfolding and 4.8 mol cause 86%. The enthalpy changes are 63.2 and 184.5 kJ mol 1 respec-

tively, leading to

H ¼ 1:8 PALH þ 0:43 T!RH ¼ 63:2 kJ mol 1 and

H ¼ 4:8 PALH þ 0:86 T!RH ¼ 184:5 kJ mol 1

Use all three equations to find the enthalpies of binding and of unfolding for this enzyme.

C H A P T E R

4

Molecular Mechanics:

Basic Theory

The first molecular modeling technique we shall look at is molecular mechanics

(MM). MM is a very fast method of determining the geometry, molecular energies, vibrational spectra, and enthalpies of formation of stable ground-state molecules. Because of its speed, it is widely used on large molecules such as those of biological or pharmaceutical importance that are currently beyond the reach of more computer-intensive molecular orbital methods. MM is an empirical method, relying on a large number of parameters, drawn from experimental data, called, collectively, the force field parameters. The major drawback of MM is encountered when one or more of the parameters necessary to solve a problem is not known. Because they are parameterized using data from molecules in the ground state, neither MM nor semiempirical methods are not as useful for modeling transitionstate chemistry as ab initio methods.

The Harmonic Oscillator

The harmonic oscillator (Fig. 4-1) is an idealized model of the simple mechanical system of a moving mass connected to a wall by a spring. Our interest is in very small masses (atoms). The harmonic oscillator might be used to model a hydrogen atom connected to a large molecule by a single bond. The large molecule is so

Computational Chemistry Using the PC, Third Edition, by Donald W. Rogers ISBN 0-471-42800-0 Copyright # 2003 John Wiley & Sons, Inc.

93