Rogers Computational Chemistry Using the PC

104	COMPUTATIONAL CHEMISTRY USING THE PC

Y, Angstroms

5

4

3

2

1

1

2

3

4

X, Angstroms

Figure 4-5 The Geometry of Ethylene with all Atoms at their Equilibrium Positions in the MM3 Force Field.

the program run when the operator declines or fails to specify a value in the input file. The second block of information in PART 1 is the INITIAL STERIC ENERGY. which we expect to be high because the molecule is not very close to its equilibrium geometry. The principal contribution to the initial steric energy is compression energy in the x-dimension brought about by the abnormally short C C bond in the input. There is a significant contribution from the H C H bond angle, which is abnormally small in the input structure. The dipole moment of this symmetrical structure is zero.

PART 1

CHEMICAL FORMULA : C( 2) H( 4)

FORMULA WEIGHT : 28.032

DATE : 09/14/2001

TIME : 11:53:17

THE COORDINATES OF 6 ATOMS ARE READ IN.

CONFORMATIONAL ENERGY, PART 1: GEOMETRY AND STERIC ENERGY OF INITIAL CONFORMATION.

CONNECTED ATOMS 1- 2-

ATTACHED ATOMS

1- 3, 1- 4, 2- 5, 2- 6,

INITIAL ATOMIC COORDINATES
ATOM	X	Y	Z	TYPE
C( 1)	2.00000	3.00000	0.00000	( 2)
C( 2)	3.00000	3.00000	0.00000	( 2)
H( 3)	1.00000	4.00000	0.00000	( 5)
H( 4)	1.00000	2.00000	0.00000	( 5)
H( 5)	4.00000	4.00000	0.00000	( 5)
H( 6)	4.00000	2.00000	0.00000	( 5)

MOLECULAR MECHANICS: BASIC THEORY		105
DIELECTRIC CONSTANT ¼ 1.500
INITIAL STERIC ENERGY IS 251.2407 KCAL.
COMPRESSION	218.0511
BENDING	31.9368
BEND-BEND	0.8293
STRETCH-BEND	0.2189
VANDERWAALS
1,4 ENERGY	0.0439
OTHER	0.0000
TORSIONAL	0.0000
TORSION-STRETCH	0.0000
DIPOLE-DIPOLE	2.3450
CHARGE-DIPOLE	0.0000
CHARGE-CHARGE	0.0000
DIPOLE MOMENT ¼	0.000 D

CPU time for initial calculation is 3.13 seconds.

PART 2 tracks the progress of iterative geometry minimization until the average

5 ˚

change in atomic position is less than the default limit of 1 10 A. In this case, the geometry is optimized by the end of the second cycle. Note that both the average and maximum atomic movements diminish as the iterations proceed. The root mean square (RMS) gradient is given in units of energy per unit distance.

CONFORMATIONAL ENERGY, PART 2: GEOMETRY MINIMIZATION

GEOMETRY OPTIMIZATION IS CONTINUED UNTIL ATOM

MOVEMENT CONVERGES WITHIN FOLLOWING VALUES:

AVERAGE MOVEMENT 0.00007 A

MAXIMUM MOVEMENT 0.00073 A

<<<<<<<<<<<<<<<<<C Y C L E 1>>>>>>>>>>>>>>>>>>

(CH)-MOVEMENT ¼ 1

ITER 1 AVG. MOVE ¼ 0.07057 A ITER 2 AVG. MOVE ¼ 0.05886 A ITER 3 AVG. MOVE ¼ 0.01364 A ITER 4 AVG. MOVE ¼ 0.00272 A ITER 5 AVG. MOVE ¼ 0.00031 A ITER 6 AVG. MOVE ¼ 0.00004 A

(MAX MOVE : ATOM 3 0.40396 A) (MAX MOVE : ATOM 3 0.43342 A) (MAX MOVE : ATOM 3 0.06326 A) (MAX MOVE : ATOM 3 0.01899 A) (MAX MOVE : ATOM 5 0.00196 A) (MAX MOVE : ATOM 3 0.00025 A)

<<<<<<<<<<<<<<<<<C Y C L E 2>>>>>>>>>>>>>>>>>

(CH)-MOVEMENT ¼ 0

ITER 7 AVG. MOVE ¼ 0.00001 A (MAX MOVE : ATOM 5 0.00003 A)

* * * * * * * * * GEOMETRY IS OPTIMIZED * * * * * * * * *

106 COMPUTATIONAL CHEMISTRY USING THE PC

WITHIN : AVERAGE MOVEMENT 0.00001 A (LIMIT 0.00007 A)

MAXIMUM MOVEMENT 0.00003 A (LIMIT 0.00073 A)

GRADIENT
RMS GRADIENT	0.000940 KCAL/MOL/A
MAX GRADIENT
X-DIRECTION : ATOM # 1	0.002595 KCAL/MOL/A
Y-DIRECTION : ATOM # 2	0.001302 KCAL/MOL/A
Z-DIRECTION : ATOM #***	0.000000 KCAL/MOL/A

DELTA TIME ¼ 0.00 SEC. TOTAL CPU TIME ¼ 3.18 SEC.

PART 3 is a repetition of PART 1 but with the final geometry (identical to the geo file). Both compression and bending contributions to the final energy of the ethylene molecule are small. Very small negative energies are sometimes encountered (e.g., the BEND-BEND energy) because steric energy is calculated relative to an arbitrary zero point. Dipole-dipole interactions are major contributors to the final steric energy. van der Waals repulsion and dipole-dipole interactions are discussed below. Although in this case relaxation of compression energy of the C C bond and relaxation of the bending energy of the H C H group are the major contributors to the decrease in initial steric energy from 251.2 kcal mol 1 to the final steric energy of 2.6 kcal mol 1, one should beware of making fine geometric distinctions drawn from steric energies ascribed to different modes of motion. Force field parameters are to some degree composite, and a single steric energy may represent more than one mode of motion.

Once we have the atomic coordinates relative to any origin, they can be translated so that the origin is at the center of mass, permitting calculation of the moments of inertia about the x-, y-, and z-axes. A single mass rotating in a plane at a fixed distance r from a center of rotation has a moment of inertia I ¼ mr2. For a collection of masses mi, each rotating at fixed ri, the moments of inertia are additive

X

Isystem ¼ miri2 ð4-12Þ

In the case of a polyatomic molecule, rotation can occur in three dimensions about the molecular center of mass. Any possible mode of rotation can be expressed as projections on the three mutually perpendicular axes, x, y, and z; hence, three moments of inertia are necessary to give the resistance to angular acceleration by any torque (twisting force) in x, y, and z space. In the MM3 output file, they are denoted IX, IY, and IZ and are given in the nonstandard units of grams square centimeters.

CHEMICAL FORMULA : C( 2) H( 4)

FORMULA WEIGHT : 28.032

DATE : 09/14/2001

TIME : 11:53:18

MOLECULAR MECHANICS: BASIC THEORY

107

* M M 3 * * * * * * * * * * * 1 9 9 2 *
*		* PARAMETER RELIABILITY
*		*
*		* BLANK ¼GOOD
*		*	B ¼ FAIR
*		*	C ¼ TRIAL
*		*	R ¼ READ IN
*		*	S ¼ SUBSTITUTED
*		*	U ¼ PARAMETER NOT EXIST
*		*
*		* THIS IS SHOWN AT THE END
*		* OF THE LINE.
* * * * * * * * * * * * * * * * * * * * * *
CONFORMATIONAL ENERGY, PART 3: GEOMETRY AND STERIC
ENERGY OF FINAL CONFORMATION.
CONNECTED ATOMS
1- 2-
ATTACHED ATOMS
1- 3, 1- 4, 2- 5, 2- 6,
FINAL ATOMIC COORDINATE
ATOM	X	Y	Z	TYPE
C( 1)	1.77425	3.00009	0.00000	( 2)
C( 2)	3.11180	3.00003	0.00000	( 2)
H( 3)	1.20884	3.94637	0.00000	( 5)
H( 4)	1.20876	2.05386	0.00000	( 5)
H( 5)	3.67729	3.94627	0.00000	( 5)
H( 6)	3.67722	2.05375	0.00000	( 5)
DIELECTRIC CONSTANT ¼ 1.500
FINAL STERIC ENERGY IS			2.6017 KCAL.
COMPRESSION		0.0190
BENDING		0.0160
BEND-BEND		0.0005
STRETCH-BEND		0.0020
VANDERWAALS
1,4 ENERGY		0.4905
OTHER		0.0000

TORSIONAL 0.0000

TORSION-STRETCH 0.0000

DIPOLE-DIPOLE 2.0748

CHARGE-DIPOLE 0.0000

CHARGE-CHARGE 0.0000

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

108				COMPUTATIONAL CHEMISTRY USING THE PC
COORDINATES TRANSLATED TO NEW ORIGIN WHICH IS
CENTER OF MASS
C( 1)	0.66877	0.00003		0.00000 ( 2)
C( 2)	0.66877	0.00003		0.00000 ( 2)
H( 3)	1.23419	0.94631		0.00000 ( 5)
H( 4)	1.23427	0.94620		0.00000 ( 5)
H( 5)	1.23426	0.94620		0.00000 ( 5)
H( 6)	1.23419	0.94631		0.00000 ( 5)
MOMENT OF INERTIA WITH THE PRINCIPAL AXES
(1) UNIT 10**(–39) GM*CM**2
IX ¼ 0.5994 IY ¼ 2.8021			IZ ¼ 3.4015
(2) UNIT ¼AU A**2
IX ¼ 3.6103 IY ¼ 16.8762			IZ ¼ 20.4865
DIPOLE MOMENT ¼ 0.000 D
COMPONENTS WITH PRINCIPAL AXES
	X ¼ 0.0000	Y ¼ 0.0000		Z ¼ 0.0000

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

End of

Total cpu time is 3.18 seconds.

This job completed at 11:53:18 (09/14/2001)

File 4-3 Conclusion of the Output File for Minimal Ethylene.

TINKER

For many individuals and chemistry departments the financial demands of commercial research level software packages are a burden. An extensive package of powerful MM and related programs called TINKER is available and can be downloaded from dasher.wustl.edu/tinker/ as freeware. These programs, from J. W. Ponder’s group at Washington University School of Medicine, St. Louis are not as easy to set up and operate as commercial programs, which are written to appeal to a wide audience of specialists and nonspecialists, but with an occasional hint from the local computer software guru, they can be up and running in a few hours. As with most molecular modelers, you will probably pick a few programs that do the job you want done and feel no guilt at all about ignoring the others. My favorite MM program is MM3. It is well to have a general idea what is out there, however, in case you run into a problem that your favorite program does not do or does poorly. The input file for a TINKER geometry minimization of ethylene is quite similar to minimal.mm3.

The name of the TINKER input file in File 4-4 is ethylene.xyz, where the .xyz indicates that the geometry is given in Cartesian coordinates. (There are other

MOLECULAR MECHANICS: BASIC THEORY

109

coordinate systems for input files.) Comparison with the MM3 input file should enable you to describe the function of every integer and floating-point number in File 4-4.

6Ethylene

1	C	2.000000	3.000000	0.000000	2	2	3	4
2	C	3.000000	3.000000	0.000000	2	1	5	6
3	H	1.000000	4.000000	0.000000	5	1
4	H	1.000000	2.000000	0.000000	5	1
5	H	4.000000	4.000000	0.000000	5	2
6	H	4.000000	2.090000	0.000000	5	2

6Ethylene

1	C	1.831255	3.008674	0.000000	2	2	3	4
2	C	3.168740	3.021325	0.000000	2	1	5	6
3	H	1.256891	3.949547	0.000000	5	1
4	H	1.274778	2.057116	0.000000	5	1
5	H	3.725220	3.972888	0.000000	5	2
6	H	3.743115	2.080450	0.000000	5	2

Files 4-4a and b. The Initial (top) and Final (bottom) Geometries of Ethylene Calculated by TINKER Using the MM3 Force Field.

Version 3.9 of TINKER offers no fewer than 18 different force fields, some still under construction or in testing, and many oriented toward biochemical and medicinal applications, as appropriate for a system developed at a medical school. In running the programs, both MM3 and TINKER prompt the operator to provide information, including the force field desired, before the program run. Comparing the geometric output for MM3 and TINKER using the MM3 force field, we see that the specific values of the coordinates differ slightly (the entire molecule has moved during minimization) but the molecular geometry itself is the same for the two

: : ˚

output files. The distance between carbon atoms is 1 338 0 001 A by both calculations.

Note the distinction between programs and force fields. It is possible to carry out a calculation using the MM3 force field and the MM3 program, or one can run the TINKER program with any one of its 18 resident force fields. It is also possible to modify a force field or to create one’s own force field. This is a difficult advanced task and is not recommended. Good force fields are the result of decades of testing by competent scientific teams. Addition of a bad force field to the literature is a disservice to the science.

Within reasonable limits, different starting coordinates can be given for the initial geometry; optimization leads to the same final geometry. Changing the input

				˚
	bond distance of the C	C bond to 1.5 A and the y-coordinates of two of the
	˚			˚
	hydrogen atoms by 0.2 A gives a final C			C bond length of 1.337 A, in agreement

110 COMPUTATIONAL CHEMISTRY USING THE PC

			˚
	with 1.338 A obtained from the previous starting geometry. The slight difference of
	˚		˚
	0.001 A		1 mA is found because the geometry criterion for exit from the loop can

never be zero but must be some small finite value. Thus two different calculations may arrive at slightly different locations very near the bottom of a potential energy well, both of which satisfy the geometric criterion for exit from the iterative loop.

COMPUTER PROJECT 4-1 The Geometry of Small Molecules

We have just seen how to construct a TINKER input file for ethylene. We shall now construct several new models and study their geometries.

Procedure. a) Using the procedure shown in constructing the ethylene input file, construct an approximate input file for H2O. The atom type for oxygen is 6. The approximate input geometry can be taken as in File 4-5, where all of the z- coordinates are set at 0. Go to the tinker directory in the MS-DOS operating system and create an input file for your H2O calculation. Be sure the extension of the input file is .xyz. Rename or edit the file as necessary.

y

O (0,1)

H (–1,0)

(1,0) H

x

3Water

1	H	1.	0.	0.	21	2
2	O	0.	1.	0.	6	1	3
3	H	1.	0.	0.	21	2

File 4-5. An Approximate Input File for the Water Molecule. A hydrogen attached to an oxygen has the special atom designator 21 as distinct from the designator 5 in hydrocarbons.

Run the file using program TINKER and force field MM3 to determine the H O bond lengths and the H O H bond angle. The program runs on the command minimize. In responding to the prompt requesting the name of the input file, include the .xyz extension. Respond to the parameter request with mm3. Take the default gradient (hit Enter). The output file is stored under the name of your input file with a tilde 1, that is, the input file h2o.xyz produces the output file h2o 1.xyz.

MOLECULAR MECHANICS: BASIC THEORY

111

Search the literature for the experimental results for the H O bond lengths and the H O H bond angle, and include a discussion of the comparison in your report.

Unlike program MM3, input format is not strict. The output file is formatted by TINKER, but the input file does not have to resemble it. After a successful run on H2O, try cutting down on the number of spaces between elements in the input file until you have arrived at File 4-4. Do the more compact files run? Does File 4-4 run?

3	Water
1	H	1.	0.	0.	21	2
2	O	0.	1.	0.	6	1 3
3	H	1.	0.	0.	21	2

File 4-6. TINKER Input File for the Water Molecule in Free Format

b)Construct a TINKER input file for ethane and determine its geometry. The numerical designator for an sp3 carbon atom is 1, and the designator for hydrocarbon H is 5. You will have some nonzero z-coordinates.

c)Using a piece of graph paper, plot the approximate coordinates of ethylene from File 4-4b, the minimized or ‘‘optimized’’ structure of the model for ethylene in the MM3 force field. Replace one of the hydrogens in File 4-4b with a methyl group. Change the coordinates of File 4-4b by replacing atom 4 (hydrogen) with the four atoms of a methyl group at the approximate geometry that you have found from your graph paper sketch. Renumber the atoms 1 to 9 as necessary (File 4-7). You now have an approximate input geometry for propene. Rename your file propene.xyz and minimize to obtain the final geometry for the propene model. Many variations are possible but the following file runs.

9Propene

1	C	1.831225	3.000000	0.000000	2	2	3	6
2	C	3.168775	3.000000	0.000000	2	1	4	5
3	H	1.265783	3.946259	0.000000	5	1
4	H	3.734217	3.946259	0.000000	5	2
5	H	3.734217	2.053741	0.000000	5	2

6	C		1.3		2.0	0.000000	1	1 7 8 9
7	H	0.2	2.0	0.0	5 6
8	H	0.5	1.5	1.0	5 6
9	H	1.2	1.2	0.0	5 6

File 4-7. TINKER Input File for Propene in Mixed Format. Mixed format can be used when one is modifying or editing an output file from a previous calculation.

Note that the strict format of the ethylene output file was not followed in adding new atoms. Be careful of your connected atom list to the right of the input file; it is a rich source of potential errors. Use your graph to keep the numbering straight.

112		COMPUTATIONAL CHEMISTRY USING THE PC
Table 4-1 Bond Lengths for Carbon-Carbon Bonds in Alkanes and Alkenes
Bond	C C	C C	C H
Bond length	—	—	—

d) Write out a table of bond lengths for the three bonds, C C, C C, and C H that you have studied using the MM3 force field (Table 4-1).

The GUI Interface

Up to this point, we have used a numerical input file to stress the fact that computers work on numbers, not diagrams. MM3 and TINKER work from numerical input files that are similar but not identical. Both can be adapted to work under the command of a graphical user interface, GUI (pronounced ‘‘gooey’’). Before going into more detail concerning MM, we shall solve a geometry optimization using the GUI of PCMODEL (Serena Software). The input is constructed by using a mouse to point and click on each atom of the connected atom list or skeleton of the molecule. This yields Fig. 4-6 (top).

After a rough estimate of the connected atom geometry has been entered, the hydrogen atoms can be added by exercising the H/AD option and the branch atom, by default a carbon atom at the top of Fig. 4-6, can be changed into an oxygen atom with the periodic table (PT) option. In this way, the skeleton at the top of Fig. 4-6 is converted to an approximate input structure at the bottom of the figure. The GUI translates the relative positions in the diagram into atomic coordinates to be input to the computer. The GUI also places the hydrogen atoms at default positions that, by the nature of the methyl group, for example, have nonzero z-coordinates (Fig. 4-7). As pictorial as all of this is, don’t forget that the computer processes only numbers, and binary numbers at that.

The parameter sets for PCMODEL V 8.0 are MMX, a derivative of the MM2 parameter set of Allinger’s group, MM3, MMFF94, AMBER, and Oplsaa. We shall continue using MM3. After minimization, the model, in this case that of propan- 2-ol, has bond distances and angles at the equilibrium values as determined by the force field. Once one knows the Cartesian coordinates of all the atoms in the model, one knows everything that can possibly be known about its geometry including all

Figure 4-6 Connected Atom Skeleton for Propan-2-ol (Top). Input Diagram with Hydrogen Atoms Added and the Oxygen Atom

Indicated (Bottom).	O

MOLECULAR MECHANICS: BASIC THEORY

113

Figure 4-7 Output of the Optimized Geometry of Propan-2-ol Depicted as a Pluto Model. Pluto is one of several pictorial options provided in PCMODEL.

bond lengths, simple angles, dihedral angles, and nonbonded interatomic distances. Each of these features is automatically calculated by PCMODEL using its QUERY option. During minimization, the image of the molecule on the CRT screen may change a lot or a little depending on how much the input geometry is changed to obtain the final geometry. After minimization, two energies (enthalpies) are displayed prominently at the right of the screen, the enthalpy of formation and the strain energy. They are related to the steric energy and are discussed below.

Parameterization

There is at present no methodical way of obtaining precise values of parameters necessary to construct a useful MM force field. Bond lengths and bond angles can be determined from rotational spectroscopy, X-ray diffraction, neutron diffraction, or electron diffraction experiments, but, unfortunately, experimental values from these sources are not strictly comparable (Burkert and Allinger, 1982). This is not because of any fundamental flaw in either the methods or the experiments but because the methods measure slightly different things. For example, X rays are diffracted by the electrons surrounding a nucleus and electrons are diffracted by the nucleus itself. If the electron probability density is centered at the nucleus these results are essentially the same, but if the electron probability density is distorted away from the nucleus toward the bond the results will be different. An example of distortion away from a nucleus toward a bond is the distortion of electrons away from the proton in a C H bond.

˚

In general, we know bond lengths to within an uncertainty of 0.005 A 0.5 pm. Bond angles are reliably known only to one or two degrees, and there are many instances of more serious angle errors. In addition to experimental uncertainties and inaccuracies due to the model (lack of coincidence between model and molecule), some models present special problems unique to their geometry. For example, some force fields calculate the ammonia molecule, NH3, to be planar when there is abundant experimental evidence that NH3 is a trigonal pyramid.

N

H H

H