Eigenvector Following with the Berny Algorithm
In the region of correct negative curvature, that is, in the vicinity of the actual transition state, one can use a local searching method to optimize the transition state structure. In order to work at all, it is therefore mandatory to provide a reasonable starting structure, which might either be obtained by guessing well or by a series of constrained optimizations (relaxed potential energy scans). Taking the isomerization of HCN to CNH again as an example, we could take the structure of highest energy along the hydrogen migration pathway and search from there using a local gradient method:
#P HF/6-31G(d) opt=(Z-Matrix,ts,calcfc,noeigen)
HCN to CNH isomerization ts opt, HF/6-31G(d)
0 1
N1
C2 1 r2
H3 2 r3 1 a3
r2=1.16868
r3=1.14549
a3=80.0
|
 |
The transition state optimization (ts) is specified here in terms of an internal coordinate system (Z-Matrix) and is based on a Hessian matrix calculated at the first point of the optimization process (calcfc). As more than one negative eigenvalue might appear in the Hessian matrix during the optimization procedure, we turn off checking the eigenvalues with noeigen. Explicit calculation of the Hessian matrix at the beginning of the optimization procedure is expensive, but necessary in transition state optimizations. That the second derivatives are actually calculated for all structural parameters can be seen at the beginning of the output file:
----------------------------
! Initial Parameters !
! (Angstroms and Degrees) !
---------------------- ----------------------
! Name Value Derivative information (Atomic Units) !
------------------------------------------------------------------------
! r2 1.1687 calculate D2E/DX2 analytically !
! r3 1.1455 calculate D2E/DX2 analytically !
! a3 80.0 calculate D2E/DX2 analytically !
------------------------------------------------------------------------
After each geometry optimization step the output file contains the current Hessian matrix, its eigenvalues, as well as the corresponding eigenvector:
Eigenvalues --- -0.21059 0.32009 1.24219
Eigenvectors required to have negative eigenvalues:
r2 r3 a3
1 -0.01659 -0.21851 0.97569
In this very simple case containing only three structural variables, there is indeed only one large and negative eigenvalue, the corresponding eigenvector being dominated by the angle a3. Transition state searches will usually only work if there is at least ONE negative eigenvalue. Multiple negative eigenvalues will not be a problem as long as one of these is significantly larger than all the others. The optimization algorithm follows the most negative (largest) eigenvalue in the optimization process.
With the derivative information in hand, the Berny algorithm steps into the supposedly correct direction uphill, at the same time lowering the energy gradient. For the new structure, a new Hessian is obtained using the previous Hessian and gradient information of the last points. This updating scheme is usually sufficient to lead to a successful transition state optimization within a small number of optimization cycles. If, however, the optimization is not successful after 10-15 optimization cycles while the structure of the system has changed considerably, a completely new Hessian should be generated by restarting the optimization with calcfc.
In the current case, the optimization takes only two steps until the default convergence criteria are fulfilled and the transition state is found with a bond angle of A(H-C-N)=77.4424 degrees, bond distances R(C-N)=116.92 pm and R(H-C)=115.46 pm, and an HF/6-31G(d) total energy of -92.7919533 Hartree. Relative to the respective HCN ground state with R(C-N)=113.25 pm, R(H-C)=105.90 pm, and a total energy of -92.8751975 Hartree, this represents an activation energy of +218.6 kJ/mol.
This system can also be used to demonstrate the importance of a good starting point for the transition state search. If the search is initiated from a bond angle of a3=90.0 instead of a3=80.0, the starting point is located further away from the actual transition state at a3=77.4424 degrees. Using this starting point, the transition state is located after four optimization cycles. Choosing larger initial bond angles of 100.0 or 110.0 degrees leads to even longer transition state searches of five or six steps.
In many cases an approximate transition state structure (and thus a good starting point for the transition state optimization) cannot be generated using the scanning strategies. This may either be too costly or due to the fact that the reaction path is described by more than one structural variable. In this situation it is desirable to generate a starting structure for the transition state optimization using a constrained optimization. The following example shows how the hydrogen migration transition state in HCN can be located following this strategy.
%chk=/scratch/test1.chk
#P HF/6-31G(d) opt=(Z-Matrix)
HCN to CNH isomerization ts opt, HF/6-31G(d)
first do partial optimization
0 1
N1
C2 1 r2
H3 2 r3 1 a3
r2=1.16868
r3=1.14549
a3=80.0 F
--Link1--
%chk=/scratch/test1.chk
#P HF/6-31G(d) opt=(Z-Matrix,ts,calcfc,noeigen,nofreeze)
geom=check guess=read
HCN to CNH isomerization ts opt, HF/6-31G(d)
now do full ts optimization
0 1
|
|
The first part of this compound job will perform a partial geometry optimization in the Z-Matrix coordinates given in the input file, freezing the bond angle a3 to a value of 80.0 degrees. In the second job step the preoptimized geometry is retrieved from the checkpoint file and the constrained bond angle is freed up using the nofreeze option of the opt keyword. This is reflected in the output file at the point of geometry retrieval:
Z-Matrix taken from the checkpoint file:
ztest12.chk
Charge = 0 Multiplicity = 1
N
C,1,r2
H,2,r3,1,a3
Variables:
r2=1.16868
r3=1.14549
Constants:
a3=80.
Recover connectivity data from disk.
Any frozen variables have been unfrozen.
Should there be more than one frozen variable then all of these constraints will be lifted at this point. The other keywords specify a transition state optimization after initial calculation of the Hessian. The transition state located in this fashion is identical to the one found before. The same sequence of steps can be performed in redundant internals using the following compound job:
%chk=/scratch/test1.chk #P HF/6-31G(d) opt=AddRed HCN to CNH isomerization ts opt, HF/6-31G(d) first do partial optimization 0 1 N1 C2 1 r2 H3 2 r3 1 a3 r2=1.16868 r3=1.14549 a3=80.0 1 2 3 80.0 F --Link1-- %chk=/scratch/test1.chk #P HF/6-31G(d) opt=(ts,calcfc,noeigen,AddRed) geom=check guess=read HCN to CNH isomerization ts opt, HF/6-31G(d) now do full ts optimization 0 1 1 2 3 80.0 A |
|