How to fix Error Link9999 in Gaussian

Question

I have a problem with geometry optimization in Gaussian software with M062X method. The structure was smoothly optimized using M062X/Gen. Nevertheless, I need very accurate results so afterwards the job was submitted once more, but this time with Opt=VeryTight and Int=UltraFine keywords. Nevertheless, I got the error:

Optimization stopped.
    -- Number of steps exceeded,  NStep= 126
    -- Flag reset to prevent archiving.

... and at the end:

Error termination request processed by link 9999.
Error termination via Lnk1e in /usr/local/gaussian/2016-C.01/g16/l9999.exe at Tue Aug 25 16:05:14 2020.

What can I do? The input is following and, as I told, geometry is taken directly from the previous step which was performed using the same method and basis but without Opt=VeryTight and Int=UltraFine keywords.

I'm stuck, don't know what to do. Any idea?

%chk=vacCFI_7_8.chk
#p M062X/Gen
opt=(Cartesian,maxcycles=200,verytight,restart) freq
int=ultrafine
geom=check guess=read
scf=xqc

vacCFI_7_8

0 1

H     0
S   3   1.00
     34.0613410              0.60251978D-02
      5.1235746              0.45021094D-01
      1.1646626              0.20189726
S   1   1.00
      0.32723041             1.0000000
S   1   1.00
      0.10307241             1.0000000
P   1   1.00
      0.8000000              1.0000000
P   1   1.00
      0.95774129632D-01      1.0000000
****
C     0
S   6   1.00
  13575.3496820              0.22245814352D-03
   2035.2333680              0.17232738252D-02
    463.22562359             0.89255715314D-02
    131.20019598             0.35727984502D-01
     42.853015891            0.11076259931
     15.584185766            0.24295627626
S   2   1.00
      6.2067138508           0.41440263448
      2.5764896527           0.23744968655
S   1   1.00
      0.57696339419          1.0000000
S   1   1.00
      0.22972831358          1.0000000
S   1   1.00
      0.95164440028D-01      1.0000000
S   1   1.00
      0.48475401370D-01      1.0000000
P   4   1.00
     34.697232244            0.53333657805D-02
      7.9582622826           0.35864109092D-01
      2.3780826883           0.14215873329
      0.81433208183          0.34270471845
P   1   1.00
      0.28887547253           .46445822433
P   1   1.00
      0.10056823671           .24955789874
D   1   1.00
      1.09700000             1.0000000
D   1   1.00
      0.31800000             1.0000000
D   1   1.00
      0.90985336424D-01      1.0000000
F   1   1.00
      0.76100000             1.0000000
****
O     0
S   6   1.00
  27032.3826310              0.21726302465D-03
   4052.3871392              0.16838662199D-02
    922.32722710             0.87395616265D-02
    261.24070989             0.35239968808D-01
     85.354641351            0.11153519115
     31.035035245            0.25588953961
S   2   1.00
     12.260860728            0.39768730901
      4.9987076005           0.24627849430
S   1   1.00
      1.1703108158           1.0000000
S   1   1.00
      0.46474740994          1.0000000
S   1   1.00
      0.18504536357          1.0000000
S   1   1.00
      0.70288026270D-01      1.0000000
P   4   1.00
     63.274954801            0.60685103418D-02
     14.627049379            0.41912575824D-01
      4.4501223456           0.16153841088
      1.5275799647           0.35706951311
P   1   1.00
      0.52935117943           .44794207502
P   1   1.00
      0.17478421270           .24446069663
P   1   1.00
      0.51112745706D-01      1.0000000
D   1   1.00
      2.31400000             1.0000000
D   1   1.00
      0.64500000             1.0000000
D   1   1.00
      0.14696477366          1.0000000
F   1   1.00
      1.42800000             1.0000000
****

Here is .xyz geometry:

O         -4.44721        1.17772       -0.00003
C         -3.68389        0.06486       -0.00000
O         -4.16201       -1.03886        0.00003
C         -2.24755        0.39366       -0.00001
C         -1.34975       -0.59256        0.00001
C          0.10468       -0.46469        0.00000
C          0.87630       -1.62398       -0.00001
C          2.26287       -1.55445       -0.00001
C          2.88570       -0.32236       -0.00001
O          4.23767       -0.13929       -0.00001
C          2.12513        0.85396        0.00001
O          2.73759        2.06167        0.00002
C          0.74994        0.77741        0.00001
H         -5.36979        0.88538       -0.00002
H         -1.98671        1.44280       -0.00003
H         -1.74093       -1.60620        0.00003
H          0.38777       -2.58968       -0.00001
H          2.86062       -2.45869       -0.00002
H          4.69415       -0.98608       -0.00002
H          3.69353        1.92694        0.00001
H          0.19248        1.70441        0.00002

Answer 1

TL;DR: There really isn't a straightforward way to fix this error, or a way to guarantee success. This usually involves a very detailed hands-on troubleshooting session, because everything is dependent on everything and even tiny screws could lead to a significant change and often the underlying issues are worse than toggling a switch (or adding a keyword).

What does the error Error termination request processed by link 9999. actually mean?

Gaussian often isn't very straight forward with its errors. There are multiple explanation sites of common errors on the internet because of this. The description of it is actually given many lines above. If you look up link 9999 on the gaussian online manual, then you'll find:

• L9999: Finalizes calculation and output

So that's the something-went-wrong-catch-all of Gaussian. That usually means Gaussian's internal procedures set a flag to prevent it from terminating normally.

In your specific case (and I experienced none other) it simply means the optimiser didn't find a stationary point within the allocated resources.

Anyone can reproduce this error quickly with the following minimal example:

%chk=error.chk
#P PM6 opt(maxcycle=3) 

L9999 error producing input

0,1
O 0. 0. 0.
H 1. 0. 0.
H 0. 1. 0.

! Blank line at the end

Obviously, the above example forces the error for demonstration purposes only. BTW, the ! is the comment character; these lines will be ignored. (The last line is there because SE decides to remove crucial blank lines automagically. So the example is copy-pastable with the comment.)

Now, Gaussian isn't the best at communicating. You'll find (as you did) the following:

[...]
 Optimization stopped.
    -- Number of steps exceeded,  NStep=   3
    -- Flag reset to prevent archiving.
                           ----------------------------
                           ! Non-Optimized Parameters !
                           ! (Angstroms and Degrees)  !
 --------------------------                            --------------------------
[...]

There it tells you that at NStep= 3 the number of steps is exceeded (more to that end later), and that'll reset some flag to prevent something. That's actually the internal procedure to cause the final message.

Steps and Cycles. Apparently Gaussian isn't very consistent in its terminology. Steps may be called cycles. So the keyword option maxcycle= sets the maximum NStep.

Unfortunately, the L9999 error isn't as much an error as it is a message to you: Basically the program is trying to tell you that it failed to reach what you have asked of it. The underlying issue is somewhere in the set-up of your calculation.

Fixing the 'error'

You need to perform a deep-dive into the nitty-gritty details of your calculation.

In the most common case, your starting guess simply was too far away from the stationary point. You can plot the convergence with each cycle to see what is going on. Most molecular viewers have that built in.
In the best case, you are a few cycles away from convergence. You can then simply restart the calculation. I recommend setting up a new calculation to prevent accidental overwriting. For the example above, this would be (in the same directory, requires > G09 D.01):

%oldchk=error.chk
%chk=error-cont.chk
#P PM6 opt(maxcycle=300, restart) 

L9999 error restart

0,1
O 0. 0. 0.
H 1. 0. 0.
H 0. 1. 0.

! Blank line at the end

So increasing the maxcycle value might already fix the 'error'. Note that the default is 128, so you need to set something higher than that if you restart. Another common way is to simply extract the last coordinates and start a new calculation.
In most cases on tight optimisation criteria and a reasonably well-behaved molecule, the default cycles will suffice. You can always save time by running pre-optimisations first. I recommend starting with semi-empirics, then increase to a pure functional with a split valence basis set. The new Gaussian versions in conjunction with a queueing system can handle chaining those calculations easily (use %oldchk), so you don't waste time.

Everything beyond that is magic.

Often, when you get to know your molecules, you will also know whether the default number of cycles will suffice and you can adjust your default route section.

---

Your example is a bit more than the default as you are trying to converge to verytight. The first question you need to ask yourself is why you are doing this and whether this is really necessary. In my personal and very subjective experience, it takes about thrice the time (and cycles) to converge from tight to verytight as it takes from loose to tight.

In order to get to convergence in fewer cycles, tighten the integration grid. The default is int(ultrafine) so you do not need to specify this. If you are running on a smaller grid, then your calculations are not reliable (Get Ready to Recalculate by Derek Lowe). If you are doing benchmark tests, and you want them to be benchmarks for the method not the grid, then use int(superfine). It might also be a good idea to force tighter SCF convergence with scf(conver=9).

The estimation or calculation of the step and step direction is directly influenced by the quality of the wave function/ density matrix. If that's too coarse, you likely circle around the minimum.

You can adjust the step size to force it to go slower, opt(maxstep=10), the default is 30. Whenever you do that, you should increase the cycles appropriately. An okay rule of thumb is by halving the step size to double the cycles. This is usually a good idea when trying to converge transition states.

Sometimes, the estimation of the Hessian is not good enough. You might want to start your calculation with something better than a guess, then use opt(calcFC). Alternatively, you can run a frequency calculation of your guess structure and analyse that before doing anything else. This might help identifying lower lying issues in the molecular structure. If you do, opt(RCFC) will be your friend.
If you really, really want it, and you know exactly what you are doing (or you don't care anymore), then calculating force constants at every step is the brute-force method, i.e. opt(calcall). I do not recommend this, it is very, very wasteful. There are intermediate options to do these things, so study the manual for more.

I always recommend using redundant coordinates. Often enough the default is good enough. In rare cases they won't. The evil ones are near 0 or 180° angles, sometimes rings, cages, strongly delocalised parts of the molecule, long bonds, dispersive interactions, etc., the list is long. You can build your own redundant coordinates and that might help in those cases.
If you are not yet doing it, you should switch to generalised internal coordinates (requires G16): geom=GIC. If you need examples, browse on Chemistry.SE.

If you really, really want to add more brute to your force, you can try optimising in cartesian coordinates. Note that you will not take advantage of many features of the Berny algorithm, you'll introduce more constraints, and you will need more cycles. (For a well-behaved molecule I'd add about 25% to be safe, but that is again as subjective as it gets.)

If you are still convinced that what you are looking for exists, you can try a different algorithm. I have never used it, but I heard about people using it.

The final failure

In most cases, if you get the L9999 error and it cannot easily be fixed (and I would include increasing the grid here), there is a deep-lying issue. You should check whether other methods produce sensible results. If your optimisations with very tight criteria succeed at BP86, PBE, TPSS, B3LYP, PBE0, TPSSh, etc., but not on M06-2X, then the problem might be that the Minnesota functional doesn't describe your system well enough. (I've only picked that because of the example. B3LYP is much more likely to fail producing sensible results. Statistically speaking.)
The usual way to go about that is to look for benchmarks studies in that field. Many things have been done already and you don't have to reinvent the wheel. (If only I could remember some of the quality review articles. I'm sure [or at least hopeful] someone will comment them below.)

The post-final Hoorah!

There are alternative programs to Gaussian. It might not be a terrible idea to cross-validate.