# Geometry optimisation with low frequency modes

I am currently running a transition state (TS) optimisation in ORCA. One of the issues that I am seeing is that whenever the Hessian (updated or calculated) has an eigenvalue with small magnitude (i.e. low frequency mode) ORCA will replace the eigenvalue with a larger number of the same sign. (There is of course a large imaginary mode corresponding to the TS)

I can understand why this is done, since the partitioned-RFO step that is used in the optimiser will take huge steps along the low eigenvalue modes. Since the step is calculated as $$f_i/(b_i-\lambda)$$ along the $$i$$th mode with $$b_i$$ being the eigenvalue and $$f_i$$ being gradient component along that eigenvector. So if $$b_i$$ is almost zero for some mode, then the shift parameter $$\lambda$$ will have to be larger, otherwise the step will reach infinity. So the steps along other modes will be suppressed.

Note that this is after projecting out the translation and rotation modes from the Hessian, so the low frequency mode is actually a vibrational mode. I suspect there is a inflection point near this TS which is why that eigenvalue is near zero. My question is - what is the best way to handle near-zero modes (non translational or rotational) in geometry optimisation (specifically quasi-Newton or RFO like methods)?

I could not unfortunately find much information on this on any program manual.

My input file looks like this (NEB-TS) calculation:

! PBE0 D3BJ def2-SVP def2/J TightSCF defgrid3 RIJCOSX NEB-TS TightOpt

%pal nprocs 8 end

%neb
product "reaction-2_PROD_pbe0.xyz"
maxiter 1000
end

%geom
trust -0.1
maxstep 0.1
maxiter 300
end

%maxcore 4000

*xyzfile 1 1 reaction-2_RCT_pbe0.xyz


The problem arises in the TS optimisation part of NEB-TS. (NEB-TS first does a CI-NEB calculation and starts a TS opt from the highest point.)

I cannot share the geometry or the full output as this is still not published. But here is a snippet of the geometry optimisation part from ORCA output showing the issue.

------------------------------------------------------------------------------
ORCA GEOMETRY RELAXATION STEP
------------------------------------------------------------------------------

Getting information on internals        .... done
Copying old internal coords+grads       .... done
Making the new internal coordinates     .... (new redundants).... done
Validating the new internal coordinates .... (new redundants).... done
Calculating the B-matrix                .... done
Calculating the G,G- and P matrices     .... done
Transforming gradient to internals      .... done
Projecting the internal gradient        .... done
Number of atoms                         ....  51
Number of internal coordinates          .... 443
Current Energy                          .... -1572.221647268 Eh
Current gradient norm                   ....     0.029414124 Eh/bohr
Maximum allowed component of the step   ....  0.100
Updating the Hessian (Bofill)          .... Diagonalizing the Hessian               .... done
Dimension of the hessian                .... 443
Lowest eigenvalues of the Hessian:
-0.002171712 -0.000031447  0.007609177  0.011020000  0.012391733
WARNING! Eigenvalue   2 too small, replaced by -0.000100000
Hessian has   2 negative eigenvalues
Taking P-RFO step
Searching for lambda that maximizes along the lowest mode


What happens after that is that the first eigenvalue becomes smaller in magnitude but remains negative, while the second eigenvalue is shifted to positive or negative depending on which side of 0 it is on. The P-RFO calculation never converges and the RMS gradient fluctuates quite a lot even if given a large number of steps.

• I believe what you are seeing is the so called quasi-harmonic approximation where (in one common form) frequencies below some cutoff like 100 cm^-1 are bumped up to that cutoff. For a TS optimization, you would either not apply this or would have it recognize that some mode of interest should not be changed. Usually this mode could be identified by it's large negative/imaginary frequency, but you are suggesting that your mode of interest has a small magnitude frequency, so it's difficult to recognize automatically.
– Tyberius
Commented May 15 at 1:52
• @Tyberius Sorry, I think I didn't write the question clearly. I have an imaginary mode corresponding to the TS mode. But I also have a very small negative eigenvalue that is not the TS mode. So ORCA shifts that mode up to a larger negative eigenvalue - and now the TS mode and that mode have similar negative eigenvalues and P-RFO constantly keeps switching between which mode to follow up. Commented May 15 at 7:45
• That is odd, I would have thought it would have shifted the small negative value to a small (but not as small) positive value. It may help to include your input/output file (with the actual structure redacted if you can't share this).
– Tyberius
Commented May 15 at 13:45
• I'm a little unclear as to what the problem is. Are you having trouble converging on a TS structure? Commented May 16 at 5:24
• @isolatedmatrix Yes, the TS optimisation never converges. Commented May 16 at 9:49

The Hessian eigenvalue lower limit is controlled by Hess_MinEV, which seems to default to 0.0001. This value (with the appropriate sign) replaces eigenvalues with a smaller absolute value. Some QM codes provide a different lower cutoff for eigenvalues that should be set to 0, as their value is likely just numerical noise and could cause numerical instabilities. I looked in the ORCA 4.2.1 docs and didn't find such a keyword, but maybe there is a more recent version that has this sort of keyword.
• Recalculate the true Hessian more frequently: the keyword Recalc_Hess n will recompute the Hessian every nth. It's possible this additional small Hessian eigenvalue is being caused by the approximate scheme (in your case the Bofill algorithm) used to update the Hessian at each step.
• Pick a better coordinate to follow: by default, the TS optimization will try to optimize along the most negative frequency mode. But if you have a better since of what the transition coordinate should be (e.g. a particular bond breaks or angle bends), you can specify this with TS_Mode <coord_type> <indices>, where coord_type is B, A, or D for bond, angle, dihedral and indices are the zero indexed atom indices.
• Try a different eigenvalue modification: as you noticed, the Hessian eigenvalues are shifted when using P-RFO. You could instead try Hess_Modification EV_Reverse. I admittedly don't know exactly what this does, but its possible the shift is what's making it have a hard time distinguish between the real TS mode and the small eigenvalue.