Imaginary modes, CPCM

Hi,

At our computational facility we have a user whose geometry optimization exhibits persistent imaginary modes. I have a script to iteratively displace the geometry along the imaginary mode and reoptimize (which I’m happy to share if it’s of interest to others), and the opt+freq cycle continues indefinitely. Gas-phase optimization works fine. I’ve tried changing xc_grid, and separately changing every relevant CPCM keyword:

SwitchThresh to 12
SASradius to 0.05
SurfaceType to SES
HPoints 302 and HeavyPoints 590,

all with the same outcome. I’ve tried reducing GEOM_OPT_TOL_GRADIENT to 3, which has the same result; here are the final energies of the optimizations:

ropt5_iter_10.out: Final energy is -759.095968138437
ropt5_iter_11.out: Final energy is -759.095966551625
ropt5_iter_12.out: Final energy is -759.095967984947
ropt5_iter_13.out: Final energy is -759.095966394464
ropt5_iter_14.out: Final energy is -759.095967838025
ropt5_iter_15.out: Final energy is -759.095966250109
ropt5_iter_1.out: Final energy is -759.095967228929
ropt5_iter_2.out: Final energy is -759.095968663064
ropt5_iter_3.out: Final energy is -759.095967107375
ropt5_iter_4.out: Final energy is -759.095968544787
ropt5_iter_5.out: Final energy is -759.095966979221
ropt5_iter_6.out: Final energy is -759.095968414041
ropt5_iter_7.out: Final energy is -759.095966843379
ropt5_iter_8.out: Final energy is -759.095968280314
ropt5_iter_9.out: Final energy is -759.095966700833

At one point I was thinking the problem is limited to the freq calulation, but these energy variations, although small, but considering the tight convergence tolerance, suggest that the problem possibly extends to the CPCM gradients as well.

I see a similar problem has been discussed previously: Imaginary frequency in CPCM solvent. So it appears the problem is still present in new versions (I’m using 6.02).

Here’s an example input file:
$molecule
-1 1
C -3.15855 2.38634 -0.145622
C -2.47795 1.06058 -0.371107
O -1.14541 1.17615 -0.189907
O -3.04181 0.0267281 -0.676997
C -0.341712 -0.0272315 -0.238338
C 0.912436 0.289266 -1.06424
C 1.81806 1.3372 -0.501815
N 2.96214 1.03343 0.0945355
O 3.76283 1.94341 0.523864
O 3.31118 -0.224597 0.269313
C -0.162378 -0.488376 1.22216
N 1.03208 -1.00047 1.56233
O -1.12867 -0.433838 1.9969
C 1.25944 -1.47831 2.91831
H -2.63465 3.18215 -0.681251
H -4.19719 2.33053 -0.471105
H -3.12401 2.62704 0.922086
H -0.926312 -0.804251 -0.741514
H 1.47584 -0.637634 -1.20484
H 0.551382 0.599781 -2.05315
H 1.60363 2.39369 -0.584651
H 1.8759 -0.831314 0.980104
H 2.29522 -1.81404 2.99421
H 1.08537 -0.689666 3.66052
H 0.595334 -2.31606 3.15661
$end

$rem
JOBTYPE = opt
BASIS = 6-31G**
METHOD = B3LYP
DFT_D = D3
SCF_ALGORITHM = DIIS_GDM
SOLVENT_METHOD = PCM
SYMMETRY = FALSE
SYM_IGNORE = TRUE
MEM_TOTAL = 8000
$end

$pcm
THEORY CPCM
METHOD SWIG
Radii BONDI
$end

Any suggestions will be appreciated, thanks.

It’s a known problem, there’s an open ticket about this issue but I don’t have any updates.

Ok, thanks, we’re going to proceed for the time being with gas-phase optimizations followed by PCM single points.

Also at UCB here, and had similar issue using Orca, same theory level, similar solvent effects, but without D3 corrections. I was able to solve the issue with a finer grid (DEFGRID3) and extreme convergence criteria (VeryTightOpt and ExtremeSCF). I would suggest trying Orca since this is a common issue with Q-chem. It is very easy to set up from a q-chem input file. Feel free to look me up on the email directory if you want to chat more.

Note that those adjustments could be made in Q-Chem also, we just don’t give cute names to thresholds (which I personally think has the negative effect of putting distance between users and the numerics of the procedure). You can tighten XC_GRID (to 2, 3, or 000099000590, for example, and then there are three different stopping criteria for geometry optimizations, GEOM_OPT_TOL_* (* = GRADIENT, DISPLACEMENT, or ENERGY); see the manual for the numerical values. It is also possible to tighten the PCM discretization grid using HPoints and HeavyPoints keywords in the $pcm section.

1 Like

That is a great point, look up Orca’s manual for the specific numerical criteria, and you can reproduce them in Q-Chem. In the system I was looking at, using a finer grid has been the most effective in removing negative frequencies for numerical frequency calculations.

Thanks for your reply, Pedro. We indeed use ORCA frequently at the facility and I had earlier confirmed that it works for this molecule. But this user’s workflow is qchem-based, and so she’s hoping to continue with using qchem. I had tried tightening all of those keywords (and others, see above), even to 1e-6 for the gradient threshold, but the problem persists. So I think this is likely bug related.

1 Like