Imaginary frequency in CPCM solvent

Hello everyone,
I have a question regarding the ground state optimization in CPCM solvent model (acetone and toluene). Based on the following input, there is an imaginary frequency mostly related to the Boron and Fluorine atoms. I also used a different grid and smaller basis set (6-31G*), but there is still an imaginary frequency. I would appreciate it if you have any suggestions for this issue.

$molecule
0 1
C 2.7200343398 -0.5442187873 -0.9227170176
C 1.2520019676 -0.6846000582 -0.7228589964
C 0.5104592957 -1.3675387249 -1.6611317225
C -0.8823631026 -1.4698819173 -1.5012575241
C -1.5316549640 -0.8505654161 -0.4472884097
C -0.7798290535 -0.2527546328 0.6090802766
C 0.6337333168 -0.1453402398 0.4475033074
C -1.3621598827 0.1825700445 1.8294805008
C -0.5907263846 0.7471659753 2.8144044484
C 0.8029652695 0.9071220605 2.6305818174
C 1.4026604847 0.4637265594 1.4784612421
C -2.9814944582 -1.1014365850 -0.5776928817
C -3.0827381975 -2.3870900917 -1.3920411669
C -1.8634571281 -2.2442205408 -2.3411230133
C -4.0227409274 -0.2917480265 -0.2536748486
C -2.9626036621 -3.6464136140 -0.5266741492
C -3.9904047126 1.1581064930 0.0758328949
C -5.3039848987 1.5959757376 0.3426322606
C -6.2143031649 0.4567395383 0.1890417190
C -5.4614030276 -0.6692619423 -0.2098067855
C -2.9598148617 2.0985973469 0.0134214218
C -3.2371464247 3.4368817376 0.2915168133
C -4.5322173198 3.8509704289 0.6118570736
C -5.5793503764 2.9308702342 0.6240961804
C -7.5909065901 0.3809373226 0.3738761906
C -8.2341861255 -0.8374902486 0.1601911170
C -7.4985852743 -1.9614257984 -0.2201648559
C -6.1160383129 -1.8890030038 -0.3975908123
C 3.2142782908 0.6038542436 -1.5486738832
C 3.5755552030 -1.5611208460 -0.4898766149
C 2.5544407336 1.7702452517 -2.0436778878
C 3.5575858945 2.5828762555 -2.5564863778
C 4.7929999626 1.9277780336 -2.3774544620
N 4.5817109741 0.7500183184 -1.7744786626
N 4.9497921866 -1.4614601783 -0.7008549656
C 5.5426709350 -2.5467552153 -0.1860359403
C 4.5605439106 -3.3834628067 0.3825188304
C 3.3204570384 -2.7835209335 0.2044230441
C 7.0164613675 -2.7551645702 -0.2475702953
C 6.1533393013 2.3934772885 -2.7676653892
H 1.0074357822 -1.8239884374 -2.5195735528
H -2.4341607689 0.0565894156 1.9771015968
H -1.0545770568 1.0755047839 3.7458177944
H 1.4043629803 1.3693491788 3.4147942978
H 2.4814498885 0.5634271732 1.3482977417
H -4.0159190997 -2.4119180126 -1.9686605913
H -2.1452364591 -1.6733439453 -3.2412402485
H -1.4631854340 -3.2137131222 -2.6738447028
H -3.0127255562 -4.5521645272 -1.1503170659
H -3.7643152243 -3.6996528020 0.2244337250
H -2.0018169383 -3.6544932989 0.0123910090
H -1.9465368053 1.8101759898 -0.2585242827
H -2.4285449273 4.1686307341 0.2500301707
H -4.7288389081 4.9021573869 0.8299374807
H -6.6005783005 3.2558893301 0.8320608848
H -8.1592520868 1.2592398250 0.6857511296
H -9.3138864250 -0.9160788019 0.2988563505
H -8.0077760213 -2.9144589749 -0.3735906300
H -5.5739276013 -2.7921352425 -0.6697631230
C 1.0954533105 2.1057815544 -2.0217987085
H 3.4280423161 3.5590148239 -3.0185018169
H 4.7542813919 -4.3320680966 0.8780369271
C 2.0146546714 -3.3371435269 0.6834849890
H 7.3557777144 -2.7388845703 -1.2933142754
H 7.5361808431 -1.9339261160 0.2664558515
H 7.2922746651 -3.7112613992 0.2138638213
H 6.1014218470 3.3804065276 -3.2436950599
H 6.8043756468 2.4432179243 -1.8831837752
H 6.6146381014 1.6755953271 -3.4610486594
B 5.6756089483 -0.2796778967 -1.3935538334
F 6.5931587937 0.3040983443 -0.5042387079
F 6.3419193572 -0.7281844984 -2.5445507621
H 0.9227374304 3.0603813512 -2.5368747511
H 0.7208732627 2.1997602441 -0.9907523450
H 0.4882021515 1.3319281241 -2.5131261632
H 2.1864934910 -4.2575714501 1.2578348527
H 1.4819931033 -2.6228302870 1.3285278644
H 1.3406544941 -3.5750251250 -0.1527527085
$end

$rem
JOBTYPE = OPT
BASIS = cc-pVDZ
EXCHANGE = wB97X-D
GEOM_OPT_MAX_CYCLES = 5000
MAX_SCF_CYCLES = 5000
GUI = 2
SYM_IGNORE = TRUE
SOLVENT_METHOD = pcm
xc_grid = 000075000302
RPA = 0
$end

$pcm
Theory CPCM
Method SWIG
Solver Inversion
Radii Bondi
$end

$solvent
Dielectric 20.7
$end

@@@

$molecule
read
$end

$rem
JOBTYPE = FREQ
BASIS = cc-PVDZ
EXCHANGE = wB97X-D
GEOM_OPT_MAX_CYCLES = 2000
GUI = 2
MAX_SCF_CYCLES = 2000
SYM_IGNORE = TRUE
xc_grid = 000075000302
RPA = 0
SOLVENT_METHOD = pcm
$end

$pcm
Theory CPCM
Method SWIG
Solver Inversion
Radii Bondi
$end

$solvent
Dielectric 20.7
$end

Best wishes,

So there are a couple of options that can help resolve the imaginary frequency.

First being you can displace the geometry along that imaginary mode, then restart the geometry optimization. This will provide you with a better initial guess geometry which could help.

Second, since you have run a frequency calculation and if you have saved that jobs scratch folder you could restart the geometry optimization using the exact Hessian. The keyword that is needed is
GEOM_OPT_HESSIAN = READ. By default the initial Hessian is a guess (unit matrix) but the exact Hessian will always improve the geometry optimization algorithm. If you have not saved the Hessian that is fine, just use the last geometry from the optimization and use the following output.

$molecule
0 1
C xxx xxx xxx

(Final optimized geometry)
$end

$rem
JOBTYPE = FREQ
BASIS = cc-PVDZ
EXCHANGE = wB97X-D
GEOM_OPT_MAX_CYCLES = 2000
GUI = 2
MAX_SCF_CYCLES = 2000
SYM_IGNORE = TRUE
xc_grid = 000075000302
RPA = 0
SOLVENT_METHOD = pcm
$end

$pcm
Theory CPCM
Method SWIG
Solver Inversion
Radii Bondi
$end

$solvent
Dielectric 20.7
$end

@@@

$molecule
read
$end

$rem
JOBTYPE = OPT
BASIS = cc-pVDZ
EXCHANGE = wB97X-D
GEOM_OPT_MAX_CYCLES = 5000
MAX_SCF_CYCLES = 5000
GUI = 2
SYM_IGNORE = TRUE
SOLVENT_METHOD = pcm
xc_grid = 000075000302
RPA = 0
GEOM_OPT_HESSIAN = READ
$end

$pcm
Theory CPCM
Method SWIG
Solver Inversion
Radii Bondi
$end

$solvent
Dielectric 20.7
$end

@Kiana can you also paste the final geometry convergence check table, it should look something like this:

                 Maximum       Tolerance    Cnvgd?                    
  Gradient       0.000149      0.000300     YES                      
  Displacement   0.000386      0.001200     YES                      
  Energy change -0.000620      0.000001      NO

It should appear before this tag:

 ******************************                                                 
 **  OPTIMIZATION CONVERGED  **                                                 
 ****************************** 

Thanks for the quick reply.
Actually, I have tried the first option that you suggested, but it was not helpful.

Here is the final geometry table:

                         Maximum     Tolerance    Cnvgd?
     Gradient           0.000187      0.000300     YES
     Displacement       0.002470      0.001200      NO
     Energy change     -0.000001      0.000001     YES

Problem is unlikely to be the DFT grid (XC_GRID) but could (?) be the PCM surface discretization grid. You can adjust that with some additional keywords in your $pcm section, e.g.

$pcm
HPoints 310
HeavyPoints 590
$end

whereas the default values are 110 (for H atoms) and 194 (for heavy atoms), so this grid is quite a bit denser. Might be overkill. An intermediate suggestion would be

$pcm
HPoints 194
HeavyPoints 310
$end

Thanks! I will try these two options.

@ jherbert
I have used these options, however, 310 is not an acceptable value for the Lebedev grid. Can I use 302 instead of that?

Thanks,

yes, that was a typo on my part.

You can use any of the available Lebedev values that are listed in the manual (it’s not quite all of the Lebedev values that are available for XC_GRID, but there are quite a few choices anyway). In our early work on PCM, we used 590 points per atom, even for H atoms, to be absolutely converged to the grid limit. (Convergence tests here: https://doi.org/10.1063/1.3511297). This turns out to be overkill for most purposes and the defaults were later reduced, as users tried to do large molecules and complained. The current defaults are set to still be within ~1 kcal/mol of the converged limit for solvation energies, which I think is okay because the models aren’t that accurate anyway (so who cares so long as grid error is smaller than intrinsic accuracy of the model). However, it’s possible that for some applications the defaults could be too loose. I was trying to suggest some intermediate values but I misremembered or mistyped 302 as 310. There are some additional convergence tests here:
Redirecting.

@ jherbert
Thank you for the explanation. Unfortunately, those two suggested options did not work and there is still an imaginary frequency. I was thinking maybe using a different xc_grid in combination with different Lebedev values is helpful to get rid of this negative frequency. I would appreciate it if you have any other suggestions.

Thanks,

Have you tried Peter’s suggestion above to compute the exact Hessian at the first step?

Yes, I have tried, but no success.

So @epif is also running this job, he is trying to see if this job is related to a known issue with frequency jobs using SWIG-PCM. I do not know the status of his inquiry but your issue could be related to the frequency calculation and not the optimization. I am not familiar with PCM but could you use a different PCM method?

In addition, what were some the results of the optimization job that used the exact Hessian? Did it optimize to the same structure or a different one? Did it find a lower energy structure?

Analytic Hessian available only for C-PCM.

Correct, sorry probably should have been more clear. I was thinking of doing the optimization with a different method other than SWIG like ISWIG, but that may not change the overall issue.

Just so I can try on my end what version of Q-Chem are you using?

@Kiana can you also attach your output for the original job. That way I can look through the geometry optimization to see if I can spot anything.

There may be some issues with the force/hessian.

$molecule
-2 1
F
F 1 4.4
$end

Analytic force:
Gradient of SCF Energy
1 2
1 0.0000000 -0.0000000
2 0.0000000 -0.0000000
3 0.0144171 -0.0144171
Finite difference:
Atom X Y Z
1 2.932823579897104e-10 -3.008024184509850e-11 6.558327072939651e-04
2 -3.008024184509850e-11 3.760030230637312e-11 -6.558321733696724e-04

@pfmclaug
Thanks! The version of the Q-chem that I am using is 5.3.2. Actually, the structure of different optimizations is the same and also there is no too much difference in energy. Indeed, when I used the xc_grid the energy is lower compared to the jobs without xc_grid. Here is the output of one of my job.

Thanks for helping me!

C 2.7030573599 -0.5572687954 -0.9157540566
C 1.2409977597 -0.7274458221 -0.6978255480
C 0.5037981760 -1.4309579873 -1.6235874690
C -0.8887546292 -1.5390932006 -1.4571688000
C -1.5368319001 -0.9019023844 -0.4132334510
C -0.7860502306 -0.2927035465 0.6375788496
C 0.6264120775 -0.1815498192 0.4722400248
C -1.3673695336 0.1503893376 1.8560870264
C -0.5962047244 0.7243770475 2.8373653299
C 0.7975337638 0.8878157035 2.6491697652
C 1.3951453196 0.4382981407 1.4970530086
C -2.9874659910 -1.1382250064 -0.5526183006
C -3.0993331532 -2.4332953522 -1.3489967734
C -1.8668998896 -2.3263598914 -2.2862131413
C -4.0181371176 -0.3056916791 -0.2468590202
C -3.0108163111 -3.6808608456 -0.4622041119
C -3.9646818306 1.1483019979 0.0631058811
C -5.2755218968 1.6129740690 0.3057848693
C -6.2037631374 0.4875684763 0.1553636167
C -5.4635391047 -0.6575184624 -0.2159735970
C -2.9177922758 2.0735966729 0.0047236206
C -3.1785592006 3.4205878132 0.2622779428
C -4.4725765677 3.8603488241 0.5595709094
C -5.5354029705 2.9566370005 0.5677347255
C -7.5860072798 0.4408413126 0.3201261602
C -8.2475941371 -0.7703056813 0.1135154782
C -7.5248555826 -1.9135537939 -0.2399331686
C -6.1375555690 -1.8694704074 -0.3976281672
C 3.1643629380 0.6041198318 -1.5450487609
C 3.5887028630 -1.5546490003 -0.4980602051
C 2.4722211056 1.7533561849 -2.0332524025
C 3.4516636996 2.5933338010 -2.5550171375
C 4.7033917439 1.9716133271 -2.3900466887
N 4.5284580766 0.7840864266 -1.7837673133
N 4.9608383866 -1.4203251803 -0.7219943606
C 5.5845317521 -2.4954514302 -0.2099692076
C 4.6277314952 -3.3533853714 0.3649801487
C 3.3701062034 -2.7826204448 0.1968751863
C 7.0612160757 -2.6840246100 -0.2754593192
C 6.0415857666 2.4812867295 -2.8025440030
H 0.9973908282 -1.8950727066 -2.4796599766
H -2.4388221623 0.0207756588 2.0075683106
H -1.0583818562 1.0530779445 3.7697871370
H 1.3996006117 1.3556740391 3.4298017917
H 2.4738822890 0.5401547483 1.3670753136
H -4.0254100677 -2.4554912016 -1.9364196047
H -2.1299000753 -1.7781974317 -3.2055414096
H -1.4697127379 -3.3077037896 -2.5834788418
H -3.0761589207 -4.5908198503 -1.0773640529
H -3.8213309828 -3.7060044959 0.2810187783
H -2.0533523326 -3.7020263067 0.0825683259
H -1.9044076959 1.7684332367 -0.2484499388
H -2.3574539202 4.1390700282 0.2244131519
H -4.6546374454 4.9172925163 0.7630723328
H -6.5546549096 3.2991380826 0.7575737075
H -8.1420034621 1.3352070597 0.6089162850
H -9.3307804528 -0.8273918719 0.2356797975
H -8.0482069541 -2.8596392768 -0.3906138553
H -5.6105139267 -2.7871953834 -0.6505633368
C 1.0062779259 2.0555739349 -2.0063921090
H 3.2918606911 3.5659475177 -3.0155683757
H 4.8494507590 -4.2970934383 0.8586961197
C 2.0835377460 -3.3748322266 0.6810126574
H 7.4017647202 -2.6898407490 -1.3211169693
H 7.5773984493 -1.8543655844 0.2288637516
H 7.3444859687 -3.6301254139 0.2010701526
H 5.9450527036 3.4714918321 -3.2638886708
H 6.7126319849 2.5483772924 -1.9342462562
H 6.5090620679 1.7920576500 -3.5206396304
B 5.6357007047 -0.2270202145 -1.4270956561
F 6.5843253640 0.3726705791 -0.5661803270
F 6.2985908736 -0.6613852267 -2.5986575955
H 0.8121239864 3.0027851673 -2.5269753182
H 0.6340280808 2.1485420197 -0.9746531612
H 0.4159693976 1.2641086196 -2.4901670187
H 2.2886829016 -4.2829921648 1.2634562938
H 1.5265885681 -2.6714744475 1.3172199223
H 1.4227497499 -3.6433001315 -0.1564492644

The analytic gradient of pcm is correct. I was looking at the wrong section.


– total gradient after adding PCM contribution –

Atom X Y Z

1 -0.000000 -0.000000 0.000656
2 0.000000 0.000000 -0.000656

Thanks for the geometry Kiana. So I have been focusing on the optimization of this molecule and it appears you have the correct optimized structure. For your job there seems to be something going on with the frequency portion resulting in the imaginary frequency. So just to be sure, the frequency job was repeated for the optimized structure by calculating the Hessian using finite-difference with analytical gradient evaluations IDERIV = 1. This makes the frequency calculation quite expensive and should be used only to validate suspicious behavior. Which for your unique case with this finite-difference evaluation produced no imaginary frequency as seen below:

 Mode:                 1                      2                      3          
 Frequency:        16.09                  17.89                  23.47          
 Force Cnst:      0.0009                 0.0010                 0.0023          
 Red. Mass:       5.6350                 5.3766                 7.0236          
 IR Active:          YES                    YES                    YES          
 IR Intens:        0.399                  0.001                  2.060          
 Raman Active:       YES                    YES                    YES 

I am running some additional tests but for the time being if you add in IDERIV = 1 for the frequency job with the optimized structure you should be able to proceed. I will post an update with more information after my additional tests.

@pfmclaug
Your suggestion really helped me a lot. I used IDERIV = 1 for the frequency job and there is no imaginary frequency. Thanks a lot for all your efforts to resolve this issue.

Best wishes,