Scanning Potential Energy Surfaces
Transition states can be localized in many cases through scanning one structural degree of freedom. This will, of course, only be successful if the reaction pathway can be described by essentially ONE structural parameter. In order to cover the relevant part of a potential energy surface, a series of calculations must be performed in which one of the structural parameters is fixed to a certain value, while all other parameters are optimized to their most favorable values. This relaxed potential energy surface scan can be performed automatically in either a Z-Matrix or redundant internal coordinate system. Both options will be illustrated using the rotational transition state in hydrogen peroxide (H-O-O-H) as an example:
1) Coordinate Driving in Internal Coordinates
The following is a Z-Matrix describing the hydrogen peroxide molecule together with its structural variables:
#P HF/6-31G(d) opt=Z-Matrix
H2O2 rotational potential 0. to 180., HF/6-31G(d) level
first step at d4=0.0
0 1
H1
O2 1 r2
O3 2 r3 1 a3
H4 3 r2 2 a3 1 d4
r2=1.0
r3=1.3
a3=110.
d4=0.0 F
|
 |
The last line of the Z-Matrix describes a value of 0.0 degrees for the H/O/O/H dihedral angle d4, the tailing character F indicating that this variable is frozen and not to be varied during the geometry optimization. Once the (partial) geometry optimization has completed, a series of additional partial geometry optimizations can be performed, fixing the H/O/O/H dihedral angle to larger and larger values. A complete rotational potential can thus be calculated through a series of separate, constrained geometry optimizations, varying the dihedral angle d4 from 0.0 to 180.0.
The same final result can be achieved in a single job using the following input:
#P HF/6-31G(d) opt=Z-Matrix nosymm
H2O2 rotational potential 0. to 180., HF/6-31G(d) level
internal coordinates
0 1
H1
O2 1 r2
O3 2 r3 1 a3
H4 3 r2 2 a3 1 d4
r2=1.0
r3=1.3
a3=110.
d4=0.0 S 18 +10.0
|
|
The last line of the Z-Matrix again describes an initial value of 0.0 degrees for the H/O/O/H dihedral angle d4 but also specifies a Scan of 18 steps, in each of which the dihedral angle d4 is varied by +10.0 degrees. In order to avoid problems caused through changes in the point group along the pathway (C2v at d4=0.0, C2 for d4=+10.0 - +170.0, C2h at d4=180.0) the nosymm keyword has been added. The choice of internal coordinates ensures, however, that both O-H bond distances as well as both H-O-O bond angles are identical all along the pathway.
Gaussian performs a series of constrained optimizations, writing the results of all of these optimizations to the standard output together with a summary of the overall results:
Summary of Optimized Potential Surface Scan
1 2 3 4 5
EIGENVALUES -- -150.75020-150.75054-150.75153-150.75305-150.75495
r2 0.94896 0.94907 0.94923 0.94948 0.94975
r3 1.40327 1.40328 1.40223 1.40043 1.39830
a3 106.69296 106.59654 106.35646 105.98433 105.53225
d4 0.00000 10.00000 20.00000 30.00000 40.00000
6 7 8 9 10
EIGENVALUES -- -150.75704-150.75910-150.76097-150.76251-150.76366
r2 0.94999 0.95014 0.95017 0.95008 0.94989
r3 1.39619 1.39440 1.39316 1.39262 1.39283
a3 105.03625 104.53256 104.04527 103.58323 103.14392
d4 50.00000 60.00000 70.00000 80.00000 90.00000
11 12 13 14 15
EIGENVALUES -- -150.76439-150.76474-150.76477-150.76457-150.76425
r2 0.94962 0.94934 0.94908 0.94887 0.94872
r3 1.39376 1.39528 1.39721 1.39932 1.40138
a3 102.72086 102.31073 101.91674 101.54825 101.21876
d4 100.00000 110.00000 120.00000 130.00000 140.00000
16 17 18 19
EIGENVALUES -- -150.76389-150.76357-150.76336-150.76328
r2 0.94863 0.94859 0.94857 0.94857
r3 1.40320 1.40461 1.40549 1.40549
a3 100.94334 100.73592 100.60734 100.60734
d4 150.00000 160.00000 170.00000 180.00000
Largest change from initial coordinates is atom 4 1.554 Angstoms.
In this particular case the energetically most favorable structure is the one at 120.0 degrees with an energy of -150.76477 Hartree. The structure at 180.0 degrees is only slightly less favorable at -150.76328 (+3.9 kJ/mol), while the structure at 0.0 degrees is substantially less favorable at -150.75020 Hartree (+38.3 kJ/mol). Given the symmetry of the molecule, the structures at 0.0 and 180.0 degrees are therefore transition states with respect to rotation around the central O-O bond. The following figure gives an overview over this part of the potential energy surface:
2) Coordinate Driving in Redundant Internal Coordinates
The Z-Matrix used in the previous input file for internal coordinates can also be used for the redundant internal coordinate definition. All that is required to perform constrained geometry optimizations in redundant internals is the modification of the opt keyword:
#P HF/6-31G(d) opt=ModRed
H2O2 rotational potential 0. to 180., HF/6-31G(d) level
redundant internals, structure at d4=0.0
0 1
H1
O2 1 r2
O3 2 r3 1 a3
H4 3 r2 2 a3 1 d4
r2=1.0
r3=1.3
a3=110.
d4=0.0
1 2 3 4 0.0 F
The modified keyword opt=ModRed leads the program to read additional input after specification of the structure of the system. This information is given on separate lines, one constraint per line. In the current example the dihedral angle specified through the centers 1, 2, 3, and 4 is set to 0.0 degrees and frozen to this value during the geometry optimization. That the dihedral angle defined through atoms 1/2/3/4 is indeed constrained to one value is visible in the list of redundant internal coordinates at the beginning of the output file:
----------------------------
! Initial Parameters !
! (Angstroms and Degrees) !
-------------------------- --------------------------
! Name Definition Value Derivative Info. !
--------------------------------------------------------------------------------
! R1 R(1,2) 1.0 estimate D2E/DX2 !
! R2 R(2,3) 1.3 estimate D2E/DX2 !
! R3 R(3,4) 1.0 estimate D2E/DX2 !
! A1 A(1,2,3) 110.0 estimate D2E/DX2 !
! A2 A(2,3,4) 110.0 estimate D2E/DX2 !
! D1 D(1,2,3,4) 0.0 Frozen !
--------------------------------------------------------------------------------
As in the Z-Matrix example before, a complete rotational potential can be constructed by performing a series of constrained optimizations with different values for dihedral angle 1/2/3/4. A complete relaxed rotational potential can be calculated in redundant internals in one calculation using the following input:
#P HF/6-31G(d) opt=AddRed nosymm
H2O2 rotational potential 0. to 180., HF/6-31G(d) level
redundant internals
0 1
H1
O2 1 r2
O3 2 r3 1 a3
H4 3 r2 2 a3 1 d4
r2=1.0
r3=1.3
a3=110.
d4=0.0
1 2 3 4 0.0 S 18 +10.0
The keyword options opt=AddRed and opt=ModRed produce identical results and are synonymous. The calculations performed in this case are very similar to those performed before in the internal coordinate system. That the dihedral angle defined through atoms 1/2/3/4 will be scanned is visible in the list of redundant internal coordinates at the beginning of the output file:
----------------------------
! Initial Parameters !
! (Angstroms and Degrees) !
-------------------------- --------------------------
! Name Definition Value Derivative Info. !
--------------------------------------------------------------------------------
! R1 R(1,2) 1.0 estimate D2E/DX2 !
! R2 R(2,3) 1.3 estimate D2E/DX2 !
! R3 R(3,4) 1.0 estimate D2E/DX2 !
! A1 A(1,2,3) 110.0 estimate D2E/DX2 !
! A2 A(2,3,4) 110.0 estimate D2E/DX2 !
! D1 D(1,2,3,4) 0.0 Scan !
--------------------------------------------------------------------------------
Again a summary of the potential energy surface scan is given at the end of the output file containing all geometrical parameters as well as the energy for each of the optimized points.
3) Driving a Bond Angle to Study Hydrogen Migration
Driving a bond angle can also be used to study the potential energy surface of reactions, the isomerization of HCN to CNH being a nice example. An input file for studying this system in internal coordinates at the HF/6-31G(d) level of theory is:
#P HF/6-31G(d) opt=Z-Matrix
HCN to CNH isomerization pathway, HF/6-31G(d)
0 1
N1
C2 1 r2
H3 2 r3 1 a3
r2=1.2
r3=1.1
a3=170.0 S 16 -10.0
|
 |
The scan starts in this case at 170.0 degrees in order to avoid problems related to the linear arrangement of atoms and the higher symmetry of the reactant and product structures at a3=180.0 and 0.0 degrees. For the same reason the scan only covers 17 steps (including the starting structure) and ends at a final value of a3=10.0 degrees. With the input file listed above the final summary provided by Gaussian is:
Summary of Optimized Potential Surface Scan
1 2 3 4 5
EIGENVALUES -- -92.87387 -92.86994 -92.86356 -92.85492 -92.84430
r2 1.13299 1.13427 1.13656 1.13974 1.14389
r3 1.05931 1.06020 1.06201 1.06478 1.06896
a3 170.00000 160.00000 150.00000 140.00000 130.00000
6 7 8 9 10
EIGENVALUES -- -92.83218 -92.81932 -92.80700 -92.79716 -92.79218
r2 1.14912 1.15491 1.16076 1.16562 1.16868
r3 1.07517 1.08410 1.09711 1.11650 1.14549
a3 120.00000 110.00000 100.00000 90.00000 80.00000
11 12 13 14 15
EIGENVALUES -- -92.79398 -92.80213 -92.81317 -92.82513 -92.83728
r2 1.17021 1.17303 1.17840 1.18177 1.17707
r3 1.19050 1.26788 1.41607 1.63724 1.85070
a3 70.00000 60.00000 50.00000 40.00000 30.00000
16 17
EIGENVALUES -- -92.84706 -92.85322
r2 1.16672 1.15763
r3 2.01158 2.10773
a3 20.00000 10.00000
According to the energies listed in this summary the transition state for this rearrangement is located at a bond angle of around 80.0 degrees. In order to pin down the location of the transition state more precisely, a second scan should be performed in which the bond angle a3 is varied between 90.0 and 70.0 degrees in smaller step sizes (such as 2.0 degrees). In this way one can locate the transition state to +/- 2.0 degrees. A more precise localization of the transition state structure must, however, rely on a gradient optimization algorithm.
4) Driving a Bond Distance to Study Hydrogen Exchange Reactions
Driving a bond distance is one possibility to study reactions in which bonds are broken. The bonds chosen for stepwise variation are obviously those that participate in the bond breaking and making processes. As there is usually more than one bond that breaks during a reaction, selecting the one that is most descriptive of the reaction coordinate is not always straight forward. The example chosen here to illustrate the situation is the identity reaction of hydrogen radical H with molecular hydrogen H2. Using an internal coordinate system and restricting the overall system to be linear (for the relevant part of the PES this is not much of a simplification) there are only two geometrical variables r3 and r4 left decribing the distances between the central hydrogen atom and the two outer hydrogen atoms:
Optimization of both distance variables at the HF/6-31G(d) level of theory leads to localization of a ground state reactant complex GS(R) with r3=72.9971 pm and r4=349.1454 pm. From this stationary point an attempt can be made to find the transition state for the exchange reaction through either stretching the short bond described by r3 or by shrinking the long bond described by r4 to shorter distances. While one could imagine these two options to be more or less equivalent in an identity reaction, we will see in the following that this is not so.
Choosing r3 as our scanning variable and selecting larger and larger values for this variable starting from the optimized value of r3=72.991 pm corresponds to pushing the covalent bond between H1 and H3 apart. An input file for scanning this variable to r3=272.991 pm in 20 steps can be found here. The result of this effort is presented in the upper part of following figure. It can clearly be seen that elongation of the bond connecting H1 and H3 leads to a continous increase in energy. The magnitude of the energy increase is far larger than expected for a radical substitution reaction and reaches values at r3=270 pm that are far greater than those for typical bond cleavage processes (here > 700 kJ/mol). Even though this simple system poses quite some problems for a number of theoretical methods (for an overview see E. Proyonv, H. Chermette, D. R. Salahub, J. Chem. Phys. 2000, 113, 10013), energies exceeding the experimentally measured reaction barrier of 40.6 kJ/mol (300K) by more than a factor of two may indicate that we are not heading for the transition state region in this type of scan.
Choosing r4 as the scanning variable and selecting smaller and smaller values for this variable starting from the optimized value of r3=349.1454 pm corresponds to pulling the covalent bond between H1 and H4 together. An input file for scanning this variable to r4=79.1454 pm in 27 steps can be found here. The result of this scan is featured in the lower part of the figure above. Aproaching from larger distances the energy now rises much more slowly until reaching a maximum at about 100 pm and just over 70 kJ/mol. Further contraction of r4 leads to a sudden decrease in energy suggesting a discontinous or at least very irregularly shaped potential energy surface. Despite the fact that the energies are much more reasonable in this second scan, it is not clear whether the energy maximum now corresponds to a true transition state.
How can the results obtained in both scans be rationalized? Plotting the bond distances r3 and r4 against each other gives a better view of the relevant part of the potential energy surface. The following figure represents a birds eye view of the PES and also contains the true transition state for the hydrogen exchange reaction at r3 = r4 = 93.4598 pm and an energy of +75 kJ/mol relative to the reactant complex.
We can now see that pulling in r4 is a well chosen scanning coordinate directly aiming at the transition state in terms of its structure as well as its energy, and that it is only after passing beyond the transition state region that r4 is not a good description of the reaction pathway anymore. Pushing apart r3, on the other hand, represents the wrong direction right from the start and effectively ruptures the existing H-H bond without forming a new one.