6
$\begingroup$

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
------------------------------------------------------------------------------

Reading the OPT-File                    .... done
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
Current trust radius                    ....  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.

$\endgroup$
6
  • 1
    $\begingroup$ 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. $\endgroup$
    – Tyberius
    Commented May 15 at 1:52
  • 1
    $\begingroup$ @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. $\endgroup$
    – S R Maiti
    Commented May 15 at 7:45
  • $\begingroup$ 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). $\endgroup$
    – Tyberius
    Commented May 15 at 13:45
  • $\begingroup$ I'm a little unclear as to what the problem is. Are you having trouble converging on a TS structure? $\endgroup$ Commented May 16 at 5:24
  • $\begingroup$ @isolatedmatrix Yes, the TS optimisation never converges. $\endgroup$
    – S R Maiti
    Commented May 16 at 9:49

1 Answer 1

2
$\begingroup$

I want to preface this by saying I'm not certain if the suggestions below will solve your problem, but some of these are generally useful to know (and it got too long for a comment).

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.

TS geometry optimizations are somewhat notorious for being difficult to converge. There are a few things you can try in a TS calculation to potentially improve convergence:

  • 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.

Another option to try if you have access and aren't tied to exclusively to Orca is AutoTS. It tries to provide a more black box experience for finding transition states given just the reactants and products without you needing to tweek various TS search parameters. It can also account for different conformations of the reactants/products/TS, which can potentially have a large impact on the computed barrier. (Disclaimer: I work on the team that develops AutoTS).

$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .