Rogers Computational Chemistry Using the PC

SEMIEMPIRICAL CALCULATIONS ON LARGER MOLECULES

275

In the general case of an electronic Hamiltonian for atoms or molecules under the Born–Oppenheimer approximation,

^	1 n	2	n ZI		n n 1
	X		X		X X
HiðriÞ ¼ 2 i¼1 ri i¼1				r1	þ i¼1 J<I	rij		ð9-37Þ
use of the variational method with antisymmetrized orbitals
ð	ð		; nÞ X H^ iðriÞcð1; 2; ; nÞdt
E ¼ c ð1; 2; . .								ð9-38Þ

produces very many integrals, but most of them drop out as they did in Exercises 9-1 and 9-3. Two classes of integrals arise from two groupings of terms in the Hamiltonian for a many-particle system, one from a sum of one-electron terms

X H^ iðriÞ ¼ X	21ri2		Z
			ri

and the other from a sum of two-electron terms

X ^ X 1

HiðriÞ ¼

rij

In the first case, permutations P and Q must be identical for nonzero terms in the antisymmetrized sum just as they were in Exercise 9-3, leaving only

ð fiðriÞ 21ri2 ri fiðriÞdt		ð9-39Þ
	Z

as the integral that contributes to the energy of the system. This integral is given the symbol Ii.

In the second case

E ¼	ð	ð c ð1; 2; . . . ; nÞ X rij	cð1; 2; . . . ; nÞdt	ð9-40Þ
		1

the operator is only a premultiplicative factor 1=rij, so evaluation of the multiple integral is again similar to the normalization in Exercise 9-3. Terms remain under the condition P ¼ Q, leaving double integrals in place of multiple integrals

ð ð fi ðriÞfj ðrjÞ rij fiðriÞfjðrjÞdt		ð9-41Þ
1

where integration is over all space ri and rj. We give this integral the name Jij.

276	COMPUTATIONAL CHEMISTRY USING THE PC

The arbitrary labels i and j that we use in our equations have no influence on the physics of the real system, so the labels fj and fi are just as valid as the labels fi and fj for a two-electron interaction. Thus, in addition to the set of Jij integrals for which P ¼ Q, we get a nonzero set of double integrals contributing to the basis set in which P and Q are related by a one electron exchange

ð ð

1
fiðriÞfjðrjÞ rij	fiðrjÞfjðriÞdt	ð9-42Þ

This integral is called Kij. If we sum over doubly occupied orbitals, the I and J integrals contribute twice to the energy of each orbital, once for each electron but the Kij integral contributes only once because there can be only one exchange of two electrons; hence the ratio of 2:1 for integrals I and J versus the integral K. Also the Kij integral is negative because it arises from a single exchange.

In summary, the Hartree–Fock equation for antisymmetrized orbitals is written

N	N	N
X	X X					ð9-43Þ
E ¼ 2 Ii þ	i¼1		ð2Jij KijÞ
		j¼1
where
Ii ¼ fi ðriÞ 21ri2			Z
			ri		fiðriÞdt	ð9-44Þ
Jij ¼ ð ð fi ðriÞfj ðrjÞ rij fiðriÞfjðrjÞdt						ð9-45Þ
			1
and	Þfj ðrjÞ rij fiðrjÞfjðriÞ dt
Kij ¼ ð ð fi ðri						ð9-46Þ
		1
Note that integrals	can	be referred to as energies in				this context because of
Eq (9-38).

The Fock Equation

By this time, we have introduced so many approximations and restrictions on our wave function and energy spectrum that is no longer quite legitimate to call it a ‘‘Schroedinger equation’’ (Schroedinger’s initial paper treated the hydrogen atom only.) We now write

Fc ¼ ec

ð9-47Þ

SEMIEMPIRICAL CALCULATIONS ON LARGER MOLECULES

277

as the Hartree–Fock equation. One-electron orbitals obey the equation

^	¼ eifi									ð9-48Þ
Fifi
^
where Fi is called the Fock operator.
The Fock operator
F^ ¼ ^fi þ Xj				ð2J^j K^jÞ						ð9-49Þ
is made up of a one-electron part
^	1	2			Z
fi ¼ 2ri					ri					ð9-50Þ
and two two-electron parts
J^jðr1	Þ ¼	ð fj ðr2					Þ r12 fjðr2Þdt			ð9-51Þ
								1
K^jðr1	Þ ¼	ð fj ðr2					Þ r12 fiðr2Þdt			ð9-52Þ
							1
	^						^
Note that Jjðr1Þ and Kjðr1Þ are operators that go to make up the Fock operator.
They operate on functions. One often sees the notation
ð fj ðr2Þ r12 fjðr2Þdtfiðr1Þ										ð9-53Þ
		1
and										ð9-54Þ
ð fj ðr2Þ r12 fiðr2Þdtfjðr1Þ
		1

^ ð Þ

or something similar, used to show their operator nature. The notation Jj r1 and

^ ð Þ

Kj r1 emphasizes that the Fock operators are functions of the (probable) locations of the electrons, which are known only through their orbitals. The orbitals, in turn, are obtained through the secular matrix of Fock operators. Once again, we arrive at a circular calculation starting with a set of assumed orbitals, calculating the elements of the F matrix leading to a new set of orbitals, and so on to selfconsistency.

Having the Slater atomic orbitals, the linear combination approximation to molecular orbitals, and the SCF method as applied to the Fock matrix, we are in a position to calculate properties of atoms and molecules ab initio, at the Hartree– Fock level of accuracy. Before doing that, however, we shall continue in the spirit of semiempirical calculations by postponing the ab initio method to Chapter 10 and invoking a rather sophisticated set of approximations and empirical substitutions

278	COMPUTATIONAL CHEMISTRY USING THE PC

that permit calculation of molecular properties with great computational efficiency. The semiempirical methods so arrived at are probably the most widespread research-level molecular orbital calculations carried out in both academic and industrial laboratories.

The Roothaan--Hall Equations

Application of the variational self-consistent field method to the Hartree–Fock equations with a linear combination of atomic orbitals leads to the Roothaan–Hall equation set published contemporaneously and independently by Roothaan and Hall in 1951. For a minimal basis set, there are as many matrix elements as there are atoms, but there may be many more elements if the basis set is not minimal.

The LCAO approximation for the wave functions in the Hartree–Fock equations

X

fi ¼ aijwj ð9-55Þ

which is essentially Eq. (7-22), gives the Roothaan equations

XX

Fi

aijwj ¼ ei

aijwj

ð9-56Þ

where the wj are LCAO basis functions. In the more general case the basis functions need not be atomic orbitals. The Roothaan equations are simultaneous equations in the minimization parameters aij. The normal equations are

ðF11 S11e1Þa11 þ ðF12 S12e1Þa12 ðF21 S21e2Þa21 þ ðF22 S22e2Þa22

.

.

.

ðFn 1 Sn 1enÞan 1 þ ðFn 2 Sn 2enÞan 2

þðF13

þðF23

þðFn 3

S13e1Þa13 þ þ ðF1 n

S23e2Þa23 þ þ ðF2 n

.

.

.

Sn 3enÞan 3 þ þ ðFnn

S1 ne1Þa1 n ¼ 0

S2 ne2Þa2 n ¼ 0

SnnenÞann ¼ 0

ð9-57Þ

These are just the secular equations shown in equation set (7-2) with F in place of H and the ‘‘stacked matrix’’ Eq. (7-6) of eigenvectors in place of a single eigenvector. In matrix notation

0

ðF21

S21ej

Þ

Þ

10a21

a22

a2n

1

ðF11

S11ej

Þ

ðF12 S12ejÞ ðF1n

S1nej

a11

a12

a1n

B

Sn1 j

. . .

Fnn

Snn j

CB ..

. .

an 1n C

¼

ð Þ

B Fn1

CBan1

ann 1

ann

C

0

9-58

B

e

e

CB

.

.

C

@ð

Þ

ð

Þ A@

A

that is,

ðF SeÞA ¼ 0

ð9-59Þ

SEMIEMPIRICAL CALCULATIONS ON LARGER MOLECULES	279
or
FA ¼ SAe	ð9-60Þ

which is the same as Eq. (7-17) except that the Fock matrix replaces the Huckel matrix. Given the Fock operator [(Eq. (9-49)] and a basis set, we can calculate the

Fock matrix elements Fij ¼		^			fifjdt. An
		fiFfjdt and the overlap elements Sij ¼
initial approximate	eigenvalue spectrum		e	, which might come from a PPP-SCF
	Ð				Ð

calculation, gives us everything we need to calculate ðF SeÞ; hence, we can solve Eq. (9-59) for the matrix of eigenvectors A. Matrix A gives us the coefficients for new linear combinations fi; fj, new Fock operators, a Fock matrix for the next iteration, which leads to an improved orbital energy spectrum, and so on.

The Roothaan–Hall equation set (9-57) is often written in the notation

N

X

ðFmn eiSmnÞcni ¼ 0

ð9-61Þ

n¼1

leading to the one-orbital energies ei where Fmn is a Fock matrix element and Smn is an element in the overlap matrix. These equations become zero when the

determinant Fmn eiSmn becomes zero because, in general, cn i 6¼0. (For a brief, readable account of the development of this equation, see Zerner, 2000.)

The Semiempirical Model and Its Approximations:

MNDO, AM1, and PM3

If we assume that S ¼ I (which is not true in general), the matrix form of the Fock equation can be written

F A ¼ AE

ð9-62Þ

where A, is the matrix of molecular orbital coefficients, having elements aij in the basis set expansion. The assumption that S ¼ I is not necessary, but it saves on computer resources (time and memory). This is the neglect of diatomic differential overlap (NDDO) simplification. If the NDDO approximation is made for twocenter terms in the Fock matrix elements, it must hold for 3- and 4-center terms. This approximation is made in all three MOPAC methods, MNDO, AM1, and PM3. (The names are trivial: Modified Neglect of Differential Overlap, Austin (TX) Method 1, and Parameterized Method 3.)

Once the format of the Fock matrix is known, the semiempirical molecular problem (and it is a considerable one) is finding a way to make valid approximations to the elements in the Fock matrix so as to avoid the many integrations

necessary in ab initio evaluation of equations like Fij ¼			^
			fiFfjdt. After this has
been done, the matrix equation (9-62) is solved by	self-consistent methods not
		Ð

unlike the PPP-SCF methods we have already used. Results from a semiempirical

280	COMPUTATIONAL CHEMISTRY USING THE PC

calculation include or may include the optimized molecular geometry, the energy values of all the quantum levels in the system, charge densities and bond strengths, electronic and vibrational spectral transitions, and derived information like ionization potentials and dipole moments.

For the purpose of approximating and parameterizing the numerous integrals necessary to obtain Fij, the energies representing electronic and nuclear interactions are broken up into categories. Valence electrons are separated from nonvalence electrons. Nonvalence electrons are taken to be part of a core of nuclei and electrons influencing valence electrons through their classical Coulombic force field. Further semiempirical approximations to the elements Fij may be made in many ways.

A widely used protocol (Thiel, 1998) is as follows:

1.A specific carbon atom attracts ‘‘its own’’ electron by what is called a onecenter, one-electron interaction.

2.Valence electrons are repelled by other electrons in valence orbitals of the same carbon atom, a one-center, two-electron repulsion. These interactions are often parameterized with spectroscopic transition energies.

3.Two-center, one-electron resonance (bonding) is treated essentially as in PPP theory; two-center electronic repulsion is treated classically (as a Coulombic repulsion). Two-center repulsion integrals represent the energy of interaction between the charge distribution on different atoms. Based on a classical model

of charge interaction, Dewar has obtained a semiempirical function of distance f ðRijÞ between point charges i and j where the distance Rij is determined from the internuclear distance between atoms and the function is fitted to give correct values in the limits of R ¼ 0 and R ¼ 1.

4.One-electron-core integrals are parameterized classically.

5.Core-core interactions are also parameterized classically.

The result is an essentially classical model with the exception of the two-center, one-electron resonance integral, which is of quantum mechanical origin. Although the parameterization is dominated by classical interactions, it is used within the quantum mechanical framework of the Fock matrix. Thus MNDO, AM1, and PM3 are legitimate quantum mechanical molecular orbital methods with a strong influx of classical empirical parameterization. Typically, a molecular orbital package (MOPAC) contains the parameters necessary for each class of calculation, MNDO, AM1, or PM3. Such a package is said to contain the MNDO, AM1, and PM3 Hamiltonians (strictly, Fock operators), which are called up with the appropriate keyword as the first line of the input program. The point here is that the three methods are really the same except for different parameterizations. They produce the same information, but the computed numerical values they arrive at are different.

Beyond these approximations as to what the actual integrals in the F matrix are, in some more recent semiempirical methods further adjustment of any or all of the elements of the F matrix is used to bring the calculated results into the best possible agreement with standard thermochemical results, largely f H298. The number of parameters per element is 5–7 in MNDO and 18 in PM3 (Thiel, 1998).

SEMIEMPIRICAL CALCULATIONS ON LARGER MOLECULES

281

Methods parameterized with thermochemical f H298 produce, of course, f H results at 298 K, as distinct from ab initio results, which are total energies at 0 K, E0. Results at 298 K are the data most useful in practical applications, for example, in the determination of standard free energy changes of reaction and equilibrium constants. Having the data in immediately useful form can be merely a convenience or, in some cases, a real advantage, especially if one does not believe that calculated heat capacity data necessary for correction from f E0 to f H298 are reliable.

Exercise 9-4

Use MNDO, AM1, and PM3 (MOPAC, ccl.net) to determine the ionization potential of the hydrogen atom

H ! Hþ þ e

Note that this exercise refers to the standard MOPAC implementation, which does not have a graphical user interface (gui).

Solution 9-4

The input file for this calculation consists of only four lines, two of which may be blank

MNDO doublet

h

The first line specifies the method and gives the spin multiplicity for one electron ðn þ 1Þ ¼ 2. The second line may be blank or it can be used for an identifying message like

Semiempirical treatment of the hydrogen atom

This identifier will be echoed in the output file. The third line is a spacer or a second comment line, and the fourth line identifies the atom (in either upper or lower case).

The MOPAC executable can be run from DOS by using the command mopac. Respond to the prompt asking for an input file with the full input filename, including the file extension if any. For example, this might be h.txt if you used a text editor to create the input file. Some systems show only the filename without the extension but you still need the extension for MOPAC. The dir command in DOS will give you the full filename. Alternatively, you can use rename h.txt h to obtain the input file h with no .txt extension. After editing an input file, for example, changing the keyword from MNDO to AM1, you will need a new filename, say, h1, to avoid redundancy with the output files created during the MOPAC run. Redundant files are not normally erased by a new run.

The MNDO output from this four-line input file contains the ionization energy along with other information (Fig. 9-3). The results for the three methods are MNDO: 11.91, AM1: 11.40, and PM3: 13.07 eV. The experimental value is 13.61 eV.

282	COMPUTATIONAL CHEMISTRY USING THE PC
	HEAT OF FORMATION	=	52.102000	KCAL
	ELECTRONIC ENERGY	=	–11.906076	EV
	CORE–CORE REPULSION	=	.000000	EV
	DIPOLE	=	.000000	DEBYE
	NO. OF FILLED LEVELS	=	0
	AND NO. OF OPEN LEVELS	=	1
	IONIZATION POTENTIAL	=	11.906276	EV
	MOLECULAR WEIGHT	=	1.008
	SCF CALCULATIONS	=	2
	COMPUTATION TIME	=	.880 SECONDS

Figure 9-3 Partial Output from the MNDO Calculation of the Ionization Potential of Hydrogen.

Exercise 9-5

Find the ground-state energy and the equilibrium bond distance (length) for the hydrogen molecule H2 with the Arguslab implementation of MOPAC (arguslab.com) and the AM1 Hamiltonian. The Arguslab implementation of MOPAC has a gui.

Solution 9-5

Build C2 by clicking on File!new and following the steps in the help tutorial (right click, drag, control right click). Getting used to a new molecular structure package takes some time. Don’t be discouraged if you have to return to the tutorial frequently at first.

When your C2 molecule appears correctly in the window, go to selection mode (arrow). Change the default sp3 C atoms to H by right clicking on one atom followed by a left click on change atom!H[s]!H_hydrogen. Repeat for the other atom. Left click calculation!energy!AM1!OK, followed by calculation!run. If a Save a Molecule screen opens up, save under a unique filename like your initials. You will probably get something like E ¼ 0:7 au (hartrees), but this depends on the bond length, which hasn’t been optimized yet.

Go through the same run routine, clicking calculation!optimize geometry etc. Your initial geometry may not converge. If it doesn’t, use the geometry cleaning tool at the upper right of the arguslab window. Optimize the geometry again. (The geometry cleaning tool isn’t exact.) After any run, you can get bond information by right clicking on the bond in the molecular diagram, which gives a series of options including Bond Info. Left clicking on Bond Info shows you the bond length ¼

˚

0.7081 A. The complete information file on the run is obtained by left clicking

Edit!latest output file.

Repeat the energy calculation (calculation!energy!AM1!OK, followed by calculation!run) at the optimum geometry. You should get 1.011 hartrees (Fig. 9-4).

Exercise 9-6

Repeat Exercise 9-5 using the Arguslab package and the MNDO and PM3 Hamiltonians. Determine the total energy of atomization to the separated stationary atoms and electrons in hartrees.

SEMIEMPIRICAL CALCULATIONS ON LARGER MOLECULES				283
************** Final Geometry *************
H	1.18325673	0.00000000	0.00000000	1
H	1.89132627	0.00000000	0.00000000	1

Figure 9-4 Cartesian (x; y; z) Geometric Output from Arguslab AM1 Calculation. The difference between the coordinates on the x-axis is 0.708 angstroms.

Solution 9-6

	MNDO	AM1	PM3	experiment
˚	0.663	0.708	0.699	0.741
r, A
E, hartrees	1.040	1.011	1.148	1.174

The Programs

MOPAC is available from many commercial and freeware sources. Among the commercial sources are Serena Software and Arguslab (See Appendix Sources). A freeware source is ccl.net (http://ccl.net/cca/software/MS-DOS/mopac_for_dos/ index.shtml). Quantum Chemistry Program Exchange (QCPE) is between commercial and freeware sources because it charges a small fee for handling and software storage. MOPAC underwent a long evolution up to MOPAC 6.0 under government financial support; hence, it is public domain software. We shall distinguish between standard MOPAC, which does not have a gui, and commercial MOPAC, which usually does. (What you pay for is the gui. Writing your own gui is beyond the capabilities of most nonspecialists.)

The standard MOPAC 6.0 implementation has an identifier similar to the following:

**************************************************************************

** FRANK J. SEILER RES. LAB., U.S. AIR FORCE ACADEMY, COLO. SPGS., CO. 80840 **

**************************************************************************

AM1 CALCULATION RESULTS

**************************************************************************

*

MOPAC: VERSION 6.00

CALC’D.

-

-

etc.

********************************************************************

Commercial PCMODEL (Serena) does not include MOPAC as part of the package, but Serena makes it available as a collateral program. PCMODEL has a gui that permits MM optimization of a molecular geometry followed by a save option in mopac format that saves the file as filename.mop in the correct format for input to standard MOPAC 6.0. This option is virtually essential for molecules larger than four or five heavy (nonhydrogen) atoms, which would be daunting to input by hand.

284	COMPUTATIONAL CHEMISTRY USING THE PC
COMPUTER PROJECT 9-1	j	Semiempirical Calculations on
		Small Molecules: HF to HI

A. Repeat Exercises 9-5 and 9-6 for the hydrogen halides HF to HI using the Arguslab package. To identify an atom, right click on it and left click on Atom Info. To change an atom, do the same thing except that you left click on change atom. You will find that several calculations give you the error message Unsupported Element, indicating that the parameters have not been included for that particular calculation. Merely put a * in your table when this happens. You will get at least one calculation (PM3) for each hydride. You should have 4 tables with 3 columns each.

B. Make tables similar to the one in Part A for the dihydrides of the group 6 elements, starting with oxygen, H2O and going to selenium. Include the simple bond angle of the dihydride, for example, H O H ¼ 107.6 (PM3), with the geometric result. To find a simple angle, hold down shift, left click on atoms H, O, and H, left click on Monitor!Angle.

C. Find the MNDO, AM1, and PM3 estimates of the dipole moment and f H298 of formaldehyde H2O O. The calculation is run just like a geometric optimization, but be sure the Dipole moment box is checked before the program run.

COMPUTER PROJECT 9-2 j Vibration of the Nitrogen Molecule

Using MOPAC and the MNDO Hamiltonian, calculate the energy and equilibrium bond length of N2 to 4 significant figures. The input file is

mndo nitrogen

N
N	1.1	1	1	0	0

where the value 1.1 in the fifth line of the file is a starting, approximate, value for the equilibrium bond length. The third and fourth entries in that line are a designator that the second N atom is connected to atom 1 (also N) and a switch saying do optimize.

˚

Increase the bond length of the molecule by 0.001 A and recalculate the energy at this fixed value. Fixing a bond length at the value given in the second entry of the fifth line of the input file (preventing optimization to the equilibrium bond length) is brought about by changing the switch from 1 to 0 in the fourth entry of that line. Leaving a blank in place of 1 would accomplish the same purpose, but it is good practice to enter 0 as a placeholder. Increase the bond length from 1.100 to 1.110 in

˚

steps of 0.001 A. Plot the energy E, calculated by MNDO, as a function of bond length r. What is the mathematical function you observe? Is this consistent with what you know about the harmonic oscillator? According to MNDO, is it reasonable to regard the N N bond as a harmonic (Hooke’s law) spring? Print out the input data to your E vs. r plot and use the file to answer the questions below.