Rogers Computational Chemistry Using the PC

154	COMPUTATIONAL CHEMISTRY USING THE PC

but both have been accounted for in the empirical bond energy scheme by the methyl group increments; therefore, no TORS correction is added to the BE sum to obtain f H298.

In n-butane,

H3C CH2CH2CH3

the central C C bond contributes a torsional energy that is not carried into f H298 by a contiguous methyl group. Hence, in n-butane, 1 unit of torsional energy for the 2-3 bond is added to the f H298 sum. In n-pentane, 2 units are added, one for each internal C C bond. In isobutane, there is no internal C C bond contribution, isopentane has 1, and so on. Ethane is an exceptional case. There is one torsional motion at the C C bond, but there are two methyl groups so the molecule has been overcorrected. One torsional energy parameter must be subtracted to obtain the correct f H298.

The value of the torsional energy increment has been variously estimated, but TORS ¼ 0:42 kcal mol 1 was settled on for the bond contribution method in MM3. In the full statistical method (see below), low-frequency torsional motion should be calculated along with all the others so the empirical TORS increment should be zero. In fact, TORS is not zero (Allinger, 1996). It appears that the TORS increment is a repository for an energy error or errors in the method that are as yet unknown.

COMPUTER PROJECT 5-2 j The Heat of Hydrogenation of Ethylene

In this experiment, we shall run two MM3 files, one for ethylene and one for ethane. Each run gives a value for the ground-state enthalpy of formation without torsional corrections of the target compound. After correcting ethane for the torsional increment of two methyl groups and only one torsional motion, calculate the ‘‘heat’’ (strictly, enthalpy) of hydrogenation hH298 for the reaction

CH2 ¼ CH2 þ H2 ! CH3 CH3

Procedure. 1) MM3. Atom types for carbon are 1 and 2 for sp3 and sp2, respectively, and hydrogen is type 5 as previously seen. Be sure the switch is on in column 65 or you will not get an enthalpy of formation. If you develop your ethane input file from the output file of ethylene (which is probably the best way to do it if you are running without a GUI) remember that you can use a nonzero z-coordinate to take into account the nonplanar nature of the ethane molecule. One should be aware of the three-dimensional nature of molecules in drawing an input geometry because, starting from a planar input file, it is possible for atoms that should move out of the input plane to become frozen in it, leading to the wrong geometry and energy (see false minima below). An example is tetrahedral methane, which can become frozen in a false planar configuration.

H

H C H

H

MOLECULAR MECHANICS II: APPLICATIONS

155

Activate MM3 with the command mm3. Answer questions: file? ethene.mm3, parameter file? Enter (default) line number 1, option 2. The default parameter set is the MM3 parameter set; don’t change it. The line number starts the system reading on the first line of your input file, and option 2 is the block diagonal followed by full matrix minimization mentioned at the end of the section on the Hessian matrix. You will see intermediate atomic coordinates as the system minimizes the geometry, followed by a final steric energy. End with 0, output Enter, coordinates Enter.

Your output is stored as TAPE4.MM3, and your final geometry is TAPE9.MM3. Read TAPE4.MM3 (cat TAPE4.MM3 in UNIX). If you are using a UNIX or LINUX system, this is case sensitive, that is, tape4.mm3 won’t work because it is lower case. Once you get into Program MM3, procedures are identical from one system to another but file handling and editing are another matter. They are usually system specific, and you may need the help of your local system guru to establish a routine for your computer.

2)PCMODEL. Draw the carbon skeleton (a single horizontal line for ethylene) and add a bond using Add_B. Add hydrogens H/AD select Force field MM3 and Compute!minimize. Record Hf, which is f H298 uncorrected for torsional energy for the molecule. Use the Edit!erase option to clear the screen and repeat the entire procedure for ethane omitting only the Add_B step because we want a single bond between the carbon atoms.

3)Output. Upon successful execution, you will obtain an output file from which you can follow the geometry change during iterative minimization. As the atoms approach their respective potential energy minima, they are moved less and less until the criterion of minimum geometry change is met (see also PART 2, File 4-3).

It is wise to keep an archive of all input and output files using Save. At the end of the output, various energies are calculated along with the ‘‘heat’’ of formation that we seek in order to obtain hydH298. In calculating hydH298, remember that the enthalpy of formation of an element in the standard state is zero.

Pi Electron Calculations

It has long been known that the enthalpy		of hydrogenation of benzene
ð		1
	49:8 kcal mol 1; Conant and Kistiakowsky, 1937) is not the same as three times

the enthalpy of hydrogenation of cyclohexene ð3 28:6 kcal mol Þ. Evidently, the double bonds that we write in the Kekule structure of benzene

are very different from the double bond in cyclohexene. We ascribe a resonance stabilization to benzene and similar molecules (toluene, naphthalene, etc.) to account for the difference. Molecules possessing a benzenoid resonance stabilization

156	COMPUTATIONAL CHEMISTRY USING THE PC

energy are said to be aromatic. Stabilization is observed to a much lesser degree in a molecule like 1,3-butadiene when its enthalpy of hydrogenation is compared with two times the enthalpy of hydrogenation of 1-butene. Double bonds occupying alternant positions in a molecular structure are called conjugated double bonds. In the case of conjugated double bonds, we call the energy discrepancy a conjugation energy. Both resonance energy and conjugation energy are quantum mechanical in origin. Neither is treated by pure molecular mechanical calculations because MM is a classical mechanical theory as distinct from a quantum mechanical theory.

Each C C bond in a conjugated or aromatic system has two electrons, designated p electrons, which are largely responsible for its quantum mechanical features. One approach to the mechanics of a system of alternant double bonds is to carry out a relatively simple quantum mechanical treatment of the p system called a valence electron self consistent field (VESCF) calculation. Self-consistent field calculations bring about a change in carbon-carbon bond orders that are no longer simply single or double. Force constants depend on bond order; hence, they are changed by the VESCF calculation. The VESCF calculation is followed by a MM minimization, which causes the atom coordinates to change, necessitating a new VESCF optimization, which necessitates a new MM calculation, and so on. This alternation between classical and quantum mechanical calculations is cut off at some point when repeated calculations do not produce a significant change in the energy. The procedure is not as difficult as it sounds, and with a fast workstation or microcomputer conjugated and aromatic systems can be treated rather easily.

Exercise 5-4

Using graph paper, construct a simple hexagon of carbon atoms and place hydrogen atoms

˚

outside the hexagon about 1 A from the carbon atoms. The x, y, z-coordinates of all 12 atoms plus atom type designators constitute a minimal MM3 starting geometry for benzene.

Solution 5-4

One of many possible solutions is

0.8 1. 0. 2

0.81. 0. 2

1.7	0.	0.	2
.8	1.	0.	2
0.8	1.	0.	2
1.7	0.	0.	2
1.5	2.	0.	5
1.5	2.	0.	5

3.0. 0. 5

1.52. 0. 5

1.5 2. 0. 5

3. 0. 0. 5

MOLECULAR MECHANICS II: APPLICATIONS

157

Exercise 5-5

Add the necessary control lines to obtain the full MM3 minimal input file and run the file under the MM3 force field to obtain the enthalpy of formation of the aromatic molecule benzene.

Solution 5-5

minbenz											1	12
T T T T T T
1					6							1
1	2	3	4	5	6	1
1	7	2	8	3	9	4	10	5	11	6	12

Placing the 5-line control block above the geometry specification block of Exercise 5-4 gives the complete minimal input file for benzene, which we can call minbenz.mm3 (or anything else you like with the extension .mm3). Aside from the geometry block, there are two important differences between minbenz.mm3 and the file minimal.mm3 for ethylene in File 4-1a. One is the switch in column 61 of the first line, the other is the set of switches TTTTTT that constitutes the entire second line. The first switch tells the system that there are alternant sp2 carbon atoms in the system, and the series of T designators responds to the logical statement ‘‘atoms 1 through 6 are alternant sp2 carbon atoms.’’ The input line says that this statement is ‘‘True for each of the first 6 carbon atoms.’’ Remember that in digital logic, T for True is 1 and False is 0, so we are justified in calling these letters logical switches. Other differences between minbenz.mm3 and minimal.mm3 are the absence of a 4 in column 67 of line 1, which changes the amount of output, and the enthalpy of formation switch in column 65 of line 3 in the minbenz.mm3. The computed enthalpy of formation of benzene found by execution of this file in Program MM3 is 20:36 kcal mol 1.

COMPUTER PROJECT 5-3 j The Resonance Energy of Benzene

Compute the resonance energy of benzene.

Procedure. We already have the computed f H298 of benzene from the previous paragraph. After adding four hydrogen atoms in the appropriate places of the sketch you have made for benzene and modifying the minimal input file for the new molecule, compute the f H298 of cyclohexene. Don’t forget to change the control block of your input file by changing the added atoms list, the total atoms entry, and so on. Repeat the procedure with cyclohexane.

Alternatively, any or all three files for benzene, cyclohexene, and cyclohexane can be generated using the draw option of PCMODEL. Either way, the cyclohexene file is

Cyclohexene											0	16	3	0	0
0	1	0.00000			10	0	0	0	0	0	0	1
	1	2	3	4	5	6	1	0	0	0	0	0	0	0
	1	7	2	8	3	9	3	10	4	11	4	12	5	13	6 14
	6	15	6	16

158		COMPUTATIONAL CHEMISTRY USING THE PC
0.05627	1.36350	0.99040	2
1.21221	0.74441	0.70886	2
1.29552	0.51251	etc.	1
0.00069	etc.		1
etc.			1
			etc.

File 5-2. An Input File For Cyclohexene.

If one remembers that zeros and blanks are the same, this is essentially a minimal file for cyclohexene. Once having obtained f H298ðbenzeneÞ ¼ 20:36 kcal mol 1,f H298 (cyclohexene), and f H298 (cyclohexane), you have enough information to obtain hydH298 (benzene), hydH298 (cyclohexene), and the resonance energy of benzene. Construct an enthalpy-level diagram showing how the resonance energy of benzene is obtained from the enthalpies of hydrogenation you have calculated.

Strain Energy

A useful computed property of some molecules is the strain energy, included in the MM3 output, which is the difference between the energy (enthalpy) of a target molecule and the energy of the same molecule calculated with a parameter set derived from reference molecules that are supposed to be strain free. (The rigorous thermodynamic distinction between energy and enthalpy drops out in the comparison.) Note that the steric energy is not the strain energy because most normal molecules have a nonzero steric energy but not all molecules have strain. Cyclohexane has a steric energy of 7:72 kcal mol 1 but a strain energy (SE) of only 1:4 kcal mol 1 in the MM3 force field. Relative values of strain energies are useful in qualitative chemical thinking, but absolute values are questionable because of the arbitrary selection of a set of strain-free reference molecules.

False Minima

In trying alternate ways of building up input files, one is bound to discover paths that lead to different enthalpies of formation for the same molecule. In general, MM minimizations seek the extremum closest to the starting geometry. Not all extrema are minima, nor is any one minimum necessarily the global minimum. Molecular geometries may come to rest at a local minimum, which is a conformer of the minimum geometry or (rarely) a saddle point on the potential energy surface. A conformer is distinguished from an isomer by the fact that one can go from one conformer to another without breaking chemical bonds, for example, by rotation of the 2-3 bond in butane whereas a change from one isomer to another requires breaking chemical bonds, for example, the isomerization from cisto trans-2-butene.

False minima may entrap the unwary; a structure may be mistaken for the ground state that does not represent the most stable conformer. If so, the calculated

MOLECULAR MECHANICS II: APPLICATIONS

159

enthalpy of formation will be higher than the thermodynamic enthalpy of formation. The error may be small, as in the case of a gauche conformer of butane mistaken for the ground-state anti conformer, or it may be tens or hundreds of kilocalories per mole resulting from a geometry that is a far cry from the desired equilibrium geometry.

There is no known computational method for finding the absolute molecular energy minimum. In the case of small molecules, one can determine the enthalpy of formation for all of the limited number of minima that can exist simply by starting from many different geometries and convincing oneself that all reasonable alternatives have been tried. Depending on molecular rigidity and symmetry, the number of conformational choices may increase rapidly for larger molecules.

The opposite approach to judicious choice uses a program that perturbs the starting geometry in a random way and then allows it to relax (Saunders, 1987). After many random ‘‘kicks’’ of this kind, the output files are examined for the lowest energy. Entry into Program MM3 is just as it was for the single minimizations leading to f H298. Activating the random or stochastic search routine in MM3 is by choosing option 8 or 9 from the routine menu presented at the start of the MM3 run. The system needs to know how many random kicks are requested as well as the kick size. This information is provided by adding a line at the bottom of the input file with an integer right justified in a field of 5 followed by a floating point number anywhere in the next field of 10 columns. Using 0 as a place marker for the integer input, the line 00005 2.0 as the very last line of a file like File 5-2 would

˚

specify five kicks of size 2.0 A.

One sees much more output as the computer ‘‘kicks’’ or ‘‘pushes’’ the molecular geometry, carries out a new minimization, and stores the energy on each kick. You are free to specify a different number of kicks or a different kick size. The kick size

˚

of 2.0 A has been found by experience to be a pretty good compromise between kicking the molecule so hard that it ends up in a bizarre geometry and not kicking it hard enough to go over conformational energy barriers.

The stored output of a stochastic search, STOCHASTIC.MEM, for the various local minima of n-pentane using Program MM3 is accessed by the following four lines

select

MM3 COORDINATE FILE (INPUT FILE)? NPENTA.MM3

LINE NUMBER ¼ 1

SEARCH OUT PUT FILE [STOCHASTIC.MEM]? STOCHASTIC.MEM

in which lines 2, 3, and 4 are responses to program prompts (be careful of case in line 2; it may be npenta.mm3). The computer responds with File 5-3.

npenta

4 CONFORMATIONS FOUND

	ITERATION		# OF TIMES	NEGATIVE
	FIRST FOUND	ENERGY	FOUND	FREQUENCY
1	1	4.1418	8	0

160		COMPUTATIONAL CHEMISTRY USING THE PC
2	2	4.9948	21	0
3	3	5.7613	13	0
4	4	7.4594	8	0

File 5-3 Output From a Stochastic Search of the n-Pentane Molecule in MM3. The input structure was subjected to 50 ‘‘kicks.’’

The search reported in File 5-3 there found four energetically distinguishable conformers. Steric energies 2, 3, and 4 are of racemic pairs. They are degenerate. Strictly speaking, energy 4 corresponds to a pair of structural conformers that are somewhat twisted relative to one another; hence, though they are degenerate, they are not exactly a racemic pair (Mencarelli, 1995).

The energy differences among conformers relative to the ground state are 0.0, 0.85, 1.62, and 3:32 kcal mol 1. The relative populations of the states, judged by the number of times they were found in a random search or 50 trials, are 0.16, 0.21, 0.15, and 0.08 when degeneracy is taken into account. In the limit of very many runs, a Boltzmann distribution would lead us to expect a ground state that is much more populous than the output indicates, but this sample is much too small for a statistical law to be valid.

The last column in File 5-3 shows that no imaginary frequencies are found in this example. In general, imaginary frequencies are found when optimization settles on a saddle point for a transition from one conformer to another. Because the force constant on a saddle point is negative, it has an imaginary root, leading to an imaginary frequency. Searching the potential energy surface of propene sometimes reveals a saddle point with the double bond eclipsed by one hydrogen of the methyl group as it rotates from one staggered conformation to another.

Dihedral Driver

A dihedral or ‘‘twist’’ angle A-B-C-D that A makes with D can be driven around the B-C axis in arbitrarily chosen angular increments. At each step during the drive, an energy is recorded so that the locus of energies as a function of drive angle gives a profile of the potential energy hill a molecule goes over on being driven away from a minimum. For example, a switch 1 in column 80 of line 2 of an input deck with a line like

10001000020000300004000000180:000:0010:0

as the last line of an input deck, rotates 1 the dihedral angle connecting atoms 1 through 4, 0001000020000300004, (skip 5 spaces, 00000), from the anti angle 0180. to the eclipsed angle 000.0 in increments of 010.0 degrees. (If you know some FORTRAN this incantation will make sense to you.) The zeroes are not necessary parts of the input file; they serve as place keepers. If a zero is omitted, however, a blank must take its place because this file is in strict MM3 format.

MOLECULAR MECHANICS II: APPLICATIONS

161

	9
	8
Energy	7
	6
Steric
	5
	4
	3
	2											Figure 5-9	Steric Energy as
												Determined using the Dihedral
	0	20	40	60	80	100	120	140	160	180	200
												Driver Option of MM3. Com-
						ω
												pare with Figure 4-12.

Exercise 5-6

Using Program MM3, the MM3 force field, SigmaPlot, and the comments in this section on dihedral driver, plot the steric energy of n-butane as a function of the dihedral angle from the anti conformation to the eclipsed conformation. Do not use a full matrix option as your optimization method; use the block diagonal option 1. Out of curiosity, you might want to try the other options to see what happens, just don’t regard the results as your solution. The answer to any ‘‘What would happen if I...’’ question is, ‘‘Try it.’’

Solution 5-6

Your plot should resemble Fig. 5-9. Note that the angles of the dihedral driver rotate through negative values but it doesn’t make any difference because the plot comes out the same, except for the change of o ¼ ð 180 þ dihedral angleÞ. The steric energies are most conveniently read from the geo file, which is TAPE9.MM3. PCMODEL contains a convenient dihedral driver and graphing combination (see PCMODEL user’s manual).

One precaution is that, especially with congested molecules, these potential energy loci should not be taken too literally because rotated atoms or groups (within the model) can ‘‘stick’’ during rotation, then suddenly ‘‘snap into place’’, giving a potential energy discontinuity that has no counterpart in the real molecule.

Full Statistical Method

Very early force fields were used in an attempt to calculate structures, enthalpies of formation, and vibrational spectra, but it was soon found that accuracy suffered severely in either the structure-energy calculations or the vibrational spectra. Force constants were, on the whole, not transferable from one field to another. The result was that early force fields evolved so as to calculate either structure and energy or spectra, but not both.

162	COMPUTATIONAL CHEMISTRY USING THE PC

To determine geometry, energy, and vibrational spectral frequencies, one must be able to determine the location, depth, and shape of the potential well corresponding to each bond in the model. Off-diagonal elements correspond to interactions like stretch-bend interactions that are reflected in the shape of the well. Early structure-energy force fields contained few off-diagonal elements, largely because of machine limitations. Machine limitations have become less severe during the later evolution of MM methods so that contemporary force fields can be used to calculate geometry and f H and also give good results for vibrational spectral lines. MM3 is one example of a newer, more versatile MM force field. The partial MM3 output file in Fig. 5-10 shows good agreement between calculated and experimental values.

no	Frequency	Symmetry		A(i)
1.	3721.5	(B2	)	vs
2.	3666.3	(A1	)	vs
3.	1593.7	(A1	)	vs

Figure 5-10 Partial MM3 Output as Related to the Vibrational Spectrum of H2O. The experimental values of the two stretching and one bending frequencies of water are 3756, 3657, and 1595 cm 1. The IR intensities are all very strong (vs).

Entropy and Heat Capacity

A force field that can produce vibrational spectra has a second advantage in that thef H calculations can be put on a much more satisfactory theoretical base by calculating an enthalpy of formation at 0 K as in ab initio procedures and then adding various thermal energies by more rigorous means than simply lumping them in with empirical bond enthalpy contributions to f H298. The stronger the theoretical base, the less likely is an unwelcome surprise in the output.

Calculations by the more rigorous procedure yield, in MM3, a sum of (a) bond energies, (b) steric energy, (c) vibrational zero point and thermal energies, and

(d) structural features POP and TORS. Energies (a), (b), and (d) are calculated as before. Bond energy parameters appear to be quite different from those of the default MM3 calculations carried out so far because zero point and thermal energies are not included in the parameters but are added later.

Energies (c) are calculated from the harmonic oscillator model (Polyatomic Molecules, Chapter 4). The quantized harmonic oscillator has a ladder of equidistant energy levels within a parabolic potential energy well. With a force field that provides good vibrational spectra, we can calculate a molecular energy at 0 K by summing bond energies from the constituent atoms and add energies (b) and (d) and then add (c), a half-quantum of energy called the zero point energy, and a vibrational contribution to the thermal energy that is calculated from rigorous statistical mechanical principles. The zero point and vibrational energy

MOLECULAR MECHANICS II: APPLICATIONS					163
	ENERGY	ENTHALPY	ENTROPY	FREE ENERGY	HEAT CAPACITY
	(Kcal/mol)	(Kcal/mol)	(eu)	(Kcal/mol)	(cal/mol/deg)
Translational	.889	1.481	34.593	–8.833	4.967
Rotational	.889	.889	10.375	–2.205	2.980
Vibrational	12.842%	12.842%	.008*	12.839%	.054*
Potential	.000#	.000#	.000	.000#	.000
Mixing	.000	.000	.000	.000	.000
Total	14.619	15.211	44.976	1.801	8.002

Figure 5-11 Partial MM3 Output Showing Entropy and Heat Capacities for Water. Experimental values are 45.13 cal K 1 mol 1 for the absolute entropy and 8.03 cal K 1 mol 1 for the heat capacity.

contributions are known as accurately as the vibrational spacings are known, that is, rather well if enough accurate off-diagonal elements are involved in the Hessian matrix (Fig. 5-11).

Fig. 5-11 shows that, for water, entropy and heat capacity are summations in which two terms dominate, the translational energy of motion of molecules treated as ideal gas particles, and rotational, energy of spin about axes having nonzero moments of inertia terms (see Problems).

Free Energy and Equilibrium

Having calculated the standard values f H and S for the participants in a chemical reaction, the obvious next step is to calculate the standard Gibbs free

energy change of reaction rG and the equilibrium constant from

X

X

5-51

H

¼

f H

H

ð

5-50

Þ

r

prod

f react

r

S

¼ Xr

Xr

S

ð

Þ

r

¼

prod

react

ð

Þ

G H T S

5-52

and

rG ¼ RT ln Keq

ð5-53Þ

or, writing Eq. (5-53) in an equivalent form

Keq ¼ e rG =RT

ð5-54Þ

Equation (5-54) makes clear a difficulty that will bedevil us throughout computational chemistry: Although the accuracy of computational chemistry is extremely high, the demands placed on our results may be even higher. In the present case, the equilibrium constant is dependent on the exponential of the standard free energy