Currently, I have several transition state structures for a reaction of the kind A + B -> AB. These TS structures were calculated using Gaussian's own Berny-Optimization procedure and validated as first-order saddle-points by performing frequency calculations. Each of the structures has one, and only one, negative frequency. (It really doesn't matter for this question, but these are DFT calculations and they don't take too long as the system is moderatly small).

To validate and analyze these transition-states in the realm of the bigger picture of the whole PES some more, I would like to perform optimizations starting from them to yield the two minima connected by the transition-state.

Using the Gaussian software the most primitive idea to do this is to perform a regular optimization. I of course did this and, oddly, for every transition state the structure of the product of the reaction is yielded. Conceptually, a true transition-state on the PES should have zero gradient. This means that starting a primitive, gradient-based optimization from a hypothetical, 100% converged transition state should not yield another minimum anyway. I assume that the Gaussian software is in this case too "smart" and "user-friendly": The optimization might not be performed along the imaginary mode or using only the gradients, but most likely all convergence-accelearting black magic available is used. For most cases, this is a reasonable default scenario as it speeds up the optimization towards a truely relevant structure. However, I am also interested in the "path" of the optimization downward the imaginary mode. So this does not work for me.

It seems like the "EigenvectorFollowing" (EF) option for Opt in Gaussian is exactly what is needed to solve my problem. However, I cannot seem to get it to work for a regular optimization starting from the TS structure. Keywords that have failed include:

#p opt(z-matrix,EF) freq b3lyp/aug-cc-pVDZ
#p opt(cartesian,EF) freq b3lyp/aug-cc-pVDZ
#p opt(EF) freq b3lyp/aug-cc-pVDZ
#p opt(calcall,EF) freq b3lyp/aug-cc-pVDZ
#p opt(calcfc,EF) freq b3lyp/aug-cc-pVDZ
#p opt(readfc,EF) freq b3lyp/aug-cc-pVDZ

In every case, the error is the same:

 (Enter /apps/gaussian/g09_D/g09/l113.exe)



Searching for this error on the internet yields many results, however they are all realted to MP4 / CCSD(T) optimizations and not to my scenario. I have also tried to specify the molecule as Z-Matrix or Cartesian coordinates.

Note that - speaking plainly - the system in concideration is two organic molecules and previous investigations showed that it is quite boring from a quantum-chemical viewpoint. By this I mean that the structures are "well-behaved" and not biradicals or other quantum-chemical-corner-cases where DFT fails.

I hope that this issue can be resolved by an experienced gaussian user. I am using Gaussian 09 Rev. E. The first part of my question therefore is how can I use the EF option of Gaussian to perform an optimization along an imaginary mode starting from a transition structure?

Secondly, I still do not have any idea how to reach the "second" minimum connected by the TS (assuming the reaction is concerted, this would be the seperated educts). The only possible approach to tackle this which I can think of currently is to perturb the system manually along the imaginary mode and perform EigenVector-Following optimizations from these perturbed structres. Is there a more elegant way, circumventing manual structural modifications, to yield the two minima connected by the TS?


how can I use the EF option of Gaussian to perform an optimization along an imaginary mode starting from a transition structure?

I have never used the EF mode of Gaussian 09, but...

Is there a more elegant way, circumventing manual structural modifications, to yield the two minima connected by the TS?

I believe the standard way of proving that your TS is the TS you want is by an intrinsic reaction coordinate calculation:


Make sure to check out the "options" tab. A decent write-up on a small system is here:


Basically, you need to restart the calculation using the checkpoint (chk) file of your TS calculation, and as you calculated the vibrational spectrum, you can use "irc=rcfc" to read the force constants from the restart file. By default the IRC calculation does 10 points in both directions of the "reaction coordinate" starting from the TS, you can change this using MaxPoints=n, and you can also specify which direction you want to study using other keywords listed in the "options" tab on the Gaussian website.

Opening the logfile in Gaussview, you can see the IRC profile - the same way you can look at the optimization profile of an "opt" calculation. I'd follow up on the IRC by optimizing each of the end-points of the profile using a standard "opt" calculation, just to make sure.

Hope that helps.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.