I optimized a structure, then calculated the frequencies for it. The frequency calculation showed the structure was not converged even though the optimization completed. Is my structure reliable?
The structure you have is likely not reliable and can lead to incorrect results, including incorrect frequencies. If the frequency calculation does not say “Stationary point found.”, then the structure is not a true stationary point: minimum, transition structure, or saddle point. This is true even if the number of imaginary frequencies is correct (e.g., 0 for a minimum). In summary, the short answer is No. This question addresses what this means in detail, and tells you what you can do about it.
Convergence Disagreements Between Optimizations and Frequency Calculations
Occasionally, the convergence checks performed during the frequency step will disagree with the ones from the optimization step. After running an optimization, best practice is to also run a frequency job to verify the nature of the stationary point as a minimum or transition structure.
You can run the frequency calculation automatically using Opt Freq in the route section. This allows you to verify that the final structure is a true stationary point. If it is not, then it is not a true minimum (or transition structure); though, it is likely to be very close to one.
The following example illustrates the discrepancy. This first output section is from the optimization, performed using a DFT method:
Item Value Threshold Converged? Maximum Force 0.000038 0.000450 YES RMS Force 0.000014 0.000300 YES Maximum Displacement 0.000635 0.001800 YES RMS Displacement 0.000367 0.001200 YES
This second part is from the frequency calculation at the “optimized” structure (performed using the same model chemistry):
Item Value Threshold Converged? Maximum Force 0.000038 0.000450 YES RMS Force 0.000014 0.000300 YES Maximum Displacement 0.005385 0.001800 NO RMS Displacement 0.002819 0.001200 NO
As you can see, the frequency calculation indicates that the structure is not converged even though the optimization completed successfully.
In general, a stationary point is found when one of two criteria is met for a point in an optimization calculation:
-
All four values listed in the output are smaller than the indicated thresholds.
-
The Maximum Force and RMS Force are two orders of magnitude smaller than the thresholds shown, regardless of the values of the displacements.
In this case, the displacements are different between the optimization and the frequency calculation because they rely on the Hessian: the matrix of second derivatives of the energy with respect to displacement of the atoms. In a frequency calculation, this is calculated analytically for methods for which second derivatives are implemented (as is true for the DFT example here). In contrast, an estimated Hessian is used in a geometry optimization unless you explicitly request a computed one using the CalcFC or CalcAll option to the Opt keyword.
In this example, the results of the convergence tests differ, leading to convergence in one case, but not the other. Clearly, the exact Hessian from the frequency calculation is more accurate than the estimated Hessian from the optimization, so the frequency results are the ones which indicate the true status of this molecular structure: very near a minimum but not a true minimum. As is well known, frequency and thermochemistry results are based on a harmonic analysis that is only valid at true stationary points. Accordingly, some results will be incorrect at non-stationary points.
You will need to continue the optimization to get a reliable structure: one that is at a true stationary point. Fortunately, the partially optimized structure you already have is an excellent starting point, and you can take advantage of having calculated the Hessian from the frequency calculation to further speed the final optimization. Generally, the optimization will complete in very few steps.
Simply start a new calculation using the checkpoint file from the previous one, using a route section like the following:
%Chk=myfile # method/basis OPT=ReadFC Freq Geom=AllCheck Guess=Read
These changes tell Gaussian to read the geometry and initial guess orbitals from the checkpoint file, to read in the force constants (Hessian) calculated in the frequency job, and then to do an optimization followed by a frequency calculation. You should check again at the end of the second calculation to make sure that a stationary point is found.
How Much Difference Does It Make?
You will not always see large changes in the results of completing the optimization in these cases. Indeed, in many cases, the changes in predicted quantities will be extremely modest. However, you cannot know in advance which cases will lead to significantly different results and which ones will not. Thus, prudence requires that every optimization be performed to a true and valid final stationary point as its final structure.
What If the The Second Frequency Job Still Fails to Converge?
If there still NOs for one or more convergence criteria after reoptimization and a second frequency calculation, then there is a problem. The optimization may be proceeding in a region of the potential energy surface (PES) that is quite flat. Geometry optimizations to a stationary point that is located in a relatively flat region are typically more sensitive to numerical accuracy than regions where the potential is steep. In the particular case of DFT methods, the calculation of the exchange-correlation terms involve numerical integrations over a grid of points, which can be source of numerical “noise” and instability. A denser (finer) integration grid will reduce the numerical instability in such calculations.
In Gaussian 16, the default grid is the UltraFine grid (a 99,590 grid), which can help obtain smoother convergence to the stationary point in cases where the potential energy surface is very flat. It is also useful for reducing numerical instability in general (e.g., DFT calculations in solution), making such calculations generally more well behaved at a modest increase in computational cost. For these reasons, it was chosen as the default grid for Gaussian 16.
Gaussian 09’s default grid is the smaller Fine grid; some Gaussian 16 users elect to use this grid as well. If this is the grid you are using, try rerunning the calculation with Int=UltraFine.
Since modification of the numerical integration grid will bring changes in the predicted total energy, it must be considered an integral part of the model chemistry. Thus, it is necessary to the same integration grid in all the calculations in the same study in order for the computed energies and molecular properties to be comparable.
See the discussion of the Integral keyword for more information about numerical integration grids.
Last updated on: 1 November 2021.
