How to troubleshoot artifical negative frequencies

How To
#1

I am doing rather routine calculation: TD-DFT, S1 state, optimization and frequencies. Everything converges smoothly, but then I end up with one imaginary frequency, which I am sure is an artefact. I tried several things at random with no effect. Question to the experts: What are standard steps of troubleshooting this sort of a problem? What are the first thing to tighten in terms of thresholds/convergence/grids?

Here is my basic input, using defaults for everything:

Blockquote
$comment
OPT& Frequency calculation for
S1 structure optimized with
method = wB97X-D
basis = 6-311+G**
Also, use old optimzier
Start with slightly distorted structure to see if there is lower-energy
minimum compared to D2 struture
$end

$molecule
0 1
C 1.4029162007 1.3106017388 0.0000590696
C 0.7000528271 1.4568212458 1.2118385596
H 1.2444817885 1.4462938585 2.1496228374
C -0.7002534075 1.4565647257 1.2117792977
H -1.2447241991 1.4459115552 2.1495365571
C -1.4029162007 1.3106017388 -0.0000590696
C -0.7000528271 1.4568212458 -1.2118385596
H -1.2444817885 1.4462938585 -2.1496228374
C 0.7002534075 1.4565647257 -1.2117792977
H 1.2447241991 1.4459115552 -2.1495365571
C 2.8207494939 0.7871832610 0.0003491932
H 3.3544685496 1.1572220091 0.8803024505
H 3.3551773248 1.1579436494 -0.8789103039
C -2.8207494939 0.7871832610 -0.0003491932
H -3.35 1.1572220091 -0.8803024505
H -3.36 1.1579436494 0.8789103039
C 2.8207494939 -0.7871832610 -0.0003491932
H 3.3544685496 -1.1572220091 -0.8803024505
H 3.3551773248 -1.1579436494 0.8789103039
C 1.4029162007 -1.3106017388 -0.0000590696
C 0.7000528271 -1.4568212458 -1.2118385596
H 1.2444817885 -1.4462938585 -2.1496228374
C -0.7002534075 -1.4565647257 -1.2117792977
H -1.2447241991 -1.4459115552 -2.1495365571
C -1.4029162007 -1.3106017388 0.0000590696
C -0.7000528271 -1.4568212458 1.2118385596
H -1.2444817885 -1.4462938585 2.1496228374
C 0.7002534075 -1.4565647257 1.2117792977
H 1.2447241991 -1.4459115552 2.1495365571
C -2.8207494939 -0.7871832610 0.0003491932
H -3.3544685496 -1.1572220091 0.8803024505
H -3.3551773248 -1.1579436494 -0.8789103039
$end

$rem
jobtype = opt
geom_opt_driver optimize !use old algorithm
method = wB97X-D
basis = 6-311+G**
cis_n_roots = 1
cis_state_deriv = 1
cis_triplets=false
cis_singlets=true
mem_total = 10000
!XC_GRID = 3 !use SG3 grid.
$end

@@@
$molecule
read
$end

$rem
jobtype = freq
method = wB97X-D
basis = 6-311+G**
!XC_GRID = 3 !use SG3 grid
scf_guess = read
cis_n_roots = 1
cis_state_deriv = 1
cis_triplets=false
cis_singlets=true
mem_total = 10000
$end

#2 
GEOM_OPT_TOL_GRADIENT 300
GEOM_OPT_TOL_DISPLACEMENT 1200
GEOM_OPT_TOL_ENERGY 100

Those are the default values, I would reduce by factor of 2 in the situation you describe. Maybe also check DFT grid. SG-2 is the default for that functional, SG-3 is a better grid and if you need to increase from there, use XG_GRID = 000099000590 for EML (99,590).

#3

You can also re-run the optimization from your lowest energy structure and start with an exact initial Hessian.

Also using the new optimizer has a more rigorous topology generation which can better describe the system leading to smoother optimization.

In addition, as @jherbert mention convergence criteria can be changed, but with the old optimizer you are limited to integers.

#4

Thanks for the advice! Here is the report:

I tried both optimizers – no difference. I tried to reduce thresholds by a factor of 2 - it did not help. Then I tried these tight thresholds with SG-3 grid-- again, no difference, the negative frequency actually increased (62 instead of original ~40). I think I will submit this job as a ticket.
And I will try another functional to see if this is a problem with wB97X-D.

#5

As a last resort you can always perturb along the imaginary mode. No automated way to do that (used to be a script floating around but the script was wrong), have to do it by hand using the normal mode eigenvectors printed by the frequency job.