Energy is Converged, but RMS Convergence Seems 'Stuck'

sbembenek · May 5, 2024, 1:27am

Although the energy has converged, the RMS isn’t getting any smaller, and seems to be ‘stuck’ at a constant value.

How can I improve the RMS?

What parameter do I change to get it to just accept the current RMS and move onto the next geometry iteration?

While I cannot post the structure in it’s entirety, let me describe the system as completely as I can:
It’s a dimer system consisting of two non-bonded monomer units. Each monomer unit is a (garden variety) bi-aryl that has a protonated amine at one end, and a de-protonated acid at the other end (classic push/pull system); the overall system is neutral. Also note that a single monomer unit does converge (even with just DIIS under version 1.5)

The starting geometry for psi4 came from a semi-empirical calculation via XTB (GFN2-xTB). I have tried two different XTB starting geometries and saw the same RMS issue with both.

I am using Psi4 1.5 (and ma attempting to upgrade, but am having issue)

Here’s the input file parameters (I’ve tried several basis sets at this point and have also used b3lyp-d3bj) →

set scf_type df
set basis 6-31G(d_p)
set reference rks

mol.update_geometry()
mol.symmetrize(1e-3)

e, wfn = optimize('b3lyp', return_wfn=True)
oeprop(wfn, "MULLIKEN_CHARGES")
mol.print_out()
print_variables()

Here’s the last steps of the given iteration →

 ==> Iterations <==
                           Total Energy        Delta E     RMS |[F,P]|
   @DF-RKS iter  35: -1565.58429333171784   -3.05383e-05   1.18306e-03 DIIS
   @DF-RKS iter  36: -1565.58428661289395    6.71882e-06   1.18307e-03 DIIS
   @DF-RKS iter  37: -1565.58427872381435    7.88908e-06   1.18308e-03 DIIS
   @DF-RKS iter  38: -1565.58428963786287   -1.09140e-05   1.18307e-03 DIIS
   @DF-RKS iter  39: -1565.58428269595538    6.94191e-06   1.18308e-03 DIIS
   @DF-RKS iter  40: -1565.58428432227720   -1.62632e-06   1.18308e-03 DIIS
   @DF-RKS iter  41: -1565.58428601946093   -1.69718e-06   1.18307e-03 DIIS
   @DF-RKS iter  42: -1565.58428358918013    2.43028e-06   1.18308e-03 DIIS
   @DF-RKS iter  43: -1565.58428516316621   -1.57399e-06   1.18308e-03 DIIS
   @DF-RKS iter  44: -1565.58428477465213    3.88514e-07   1.18308e-03 DIIS
   @DF-RKS iter  45: -1565.58428440906641    3.65586e-07   1.18308e-03 DIIS
   @DF-RKS iter  46: -1565.58428495098451   -5.41918e-07   1.18308e-03 DIIS
   @DF-RKS iter  47: -1565.58428459385368    3.57131e-07   1.18308e-03 DIIS
   @DF-RKS iter  48: -1565.58428468662214   -9.27685e-08   1.18308e-03 DIIS
   @DF-RKS iter  49: -1565.58428476502786   -7.84057e-08   1.18308e-03 DIIS
   @DF-RKS iter  50: -1565.58428464436361    1.20664e-07   1.18308e-03 DIIS
   @DF-RKS iter  51: -1565.58428472528931   -8.09257e-08   1.18308e-03 DIIS
   @DF-RKS iter  52: -1565.58428470326794    2.20214e-08   1.18308e-03 DIIS
   @DF-RKS iter  53: -1565.58428468647639    1.67915e-08   1.18308e-03 DIIS
   @DF-RKS iter  54: -1565.58428471334310   -2.68667e-08   1.18308e-03 DIIS
   @DF-RKS iter  55: -1565.58428469500745    1.83356e-08   1.18308e-03 DIIS
   @DF-RKS iter  56: -1565.58428470022113   -5.21368e-09   1.18308e-03 DIIS
   @DF-RKS iter  57: -1565.58428470380500   -3.58386e-09   1.18308e-03 DIIS
   @DF-RKS iter  58: -1565.58428469783030    5.97470e-09   1.18308e-03 DIIS
   @DF-RKS iter  59: -1565.58428470197896   -4.14866e-09   1.18308e-03 DIIS
   @DF-RKS iter  60: -1565.58428470074773    1.23123e-09   1.18308e-03 DIIS
   @DF-RKS iter  61: -1565.58428469998626    7.61474e-10   1.18308e-03 DIIS
   @DF-RKS iter  62: -1565.58428470131480   -1.32854e-09   1.18308e-03 DIIS
   @DF-RKS iter  63: -1565.58428470037711    9.37689e-10   1.18308e-03 DIIS
   @DF-RKS iter  64: -1565.58428470066701   -2.89901e-10   1.18308e-03 DIIS
   @DF-RKS iter  65: -1565.58428470082890   -1.61890e-10   1.18308e-03 DIIS
   @DF-RKS iter  66: -1565.58428470053059    2.98314e-10   1.18308e-03 DIIS
   @DF-RKS iter  67: -1565.58428470074432   -2.13731e-10   1.18308e-03 DIIS
   @DF-RKS iter  68: -1565.58428470067770    6.66205e-11   1.18308e-03 DIIS
   @DF-RKS iter  69: -1565.58428470064223    3.54703e-11   1.18308e-03 DIIS
   @DF-RKS iter  70: -1565.58428470070658   -6.43468e-11   1.18308e-03 DIIS
   @DF-RKS iter  71: -1565.58428470065996    4.66116e-11   1.18308e-03 DIIS
   @DF-RKS iter  72: -1565.58428470067474   -1.47793e-11   1.18308e-03 DIIS
   @DF-RKS iter  73: -1565.58428470068384   -9.09495e-12   1.18308e-03 DIIS
   @DF-RKS iter  74: -1565.58428470067020    1.36424e-11   1.18308e-03 DIIS
   @DF-RKS iter  75: -1565.58428470067929   -9.09495e-12   1.18308e-03 DIIS
   @DF-RKS iter  76: -1565.58428470067497    4.32010e-12   1.18308e-03 DIIS
   @DF-RKS iter  77: -1565.58428470067543   -4.54747e-13   1.18308e-03 DIIS
   @DF-RKS iter  78: -1565.58428470067815   -2.72848e-12   1.18308e-03 DIIS
   @DF-RKS iter  79: -1565.58428470067679    1.36424e-12   1.18308e-03 DIIS
   @DF-RKS iter  80: -1565.58428470067565    1.13687e-12   1.18308e-03 DIIS
   @DF-RKS iter  81: -1565.58428470067543    2.27374e-13   1.18308e-03 DIIS
   @DF-RKS iter  82: -1565.58428470067543    0.00000e+00   1.18308e-03 DIIS
   @DF-RKS iter  83: -1565.58428470067588   -4.54747e-13   1.18308e-03 DIIS
   @DF-RKS iter  84: -1565.58428470067724   -1.36424e-12   1.18308e-03 DIIS
   @DF-RKS iter  85: -1565.58428470067588    1.36424e-12   1.18308e-03 DIIS
   @DF-RKS iter  86: -1565.58428470067497    9.09495e-13   1.18308e-03 DIIS
   @DF-RKS iter  87: -1565.58428470067747   -2.50111e-12   1.18308e-03 DIIS
   @DF-RKS iter  88: -1565.58428470067452    2.95586e-12   1.18308e-03 DIIS
   @DF-RKS iter  89: -1565.58428470067406    4.54747e-13   1.18308e-03 DIIS
   @DF-RKS iter  90: -1565.58428470067588   -1.81899e-12   1.18308e-03 DIIS
   @DF-RKS iter  91: -1565.58428470067543    4.54747e-13   1.18308e-03 DIIS
   @DF-RKS iter  92: -1565.58428470067429    1.13687e-12   1.18308e-03 DIIS
   @DF-RKS iter  93: -1565.58428470067679   -2.50111e-12   1.18308e-03 DIIS
   @DF-RKS iter  94: -1565.58428470067543    1.36424e-12   1.18308e-03 DIIS

Other things I’ve tried →

damping_percentage → 20-40%
diis_max_vecs → 50-100 (info I found here)
changing initial guess → huckel
making the grid finer →

set dft_spherical_points 434
set dft_radial_points 85
set dft_pruning_scheme robust

Rollin_King · May 8, 2024, 5:51pm

So the density-fitted SCF procedure is not converging. Perhaps (I’m not sure this is the parameter that is the problem) try changing from the default:

CHOLESKY_TOLERANCE = 1e-4

susilehtola · May 9, 2024, 2:29pm

Psi4 1.5 doesn’t have ADIIS which is the solution to most of this type of problems. Updating to a newer release should fix the issue. Plain DIIS simply isn’t good enough an algorithm to land on the minimum.

sbembenek · May 17, 2024, 1:50am

Hi @Rollin_King

I did try changing the value:

CHOLESKY_TOLERANCE = 1e-3

but it did not help. Thanks.

sbembenek · May 21, 2024, 2:49am

Thanks for the info @susilehtola

I have run the same calculation under 1.6, which does have the ADIIS, but it is still failing to converge:

   @DF-RKS iter  82: -1566.28338767164814   -2.87132e+00   3.60315e-04 DIIS/ADIIS
   @DF-RKS iter  83: -1566.28188741430313    1.50026e-03   3.60649e-04 DIIS/ADIIS
   @DF-RKS iter  84: -1566.29166893139882   -9.78152e-03   3.56495e-04 DIIS/ADIIS
   @DF-RKS iter  85: -1567.22658715057037   -9.34918e-01   5.34598e-05 DIIS
   @DF-RKS iter  86: -1567.22951375307025   -2.92660e-03   3.98180e-05 DIIS
   @DF-RKS iter  87: -1567.22971632027429   -2.02567e-04   4.10351e-05 DIIS
   @DF-RKS iter  88: -1566.31295608924188    9.16760e-01   3.47884e-04 DIIS/ADIIS
   @DF-RKS iter  89: -1567.23100097849124   -9.18045e-01   3.45396e-05 DIIS
   @DF-RKS iter  90: -1566.28146330815116    9.49538e-01   3.60723e-04 DIIS/ADIIS
   @DF-RKS iter  91: -1566.27857066020761    2.89265e-03   3.61701e-04 DIIS/ADIIS
   @DF-RKS iter  92: -1566.27152156198849    7.04910e-03   3.64526e-04 DIIS/ADIIS
   @DF-RKS iter  93: -1566.27398687179038   -2.46531e-03   3.63504e-04 DIIS/ADIIS
   @DF-RKS iter  94: -1566.27030109443081    3.68578e-03   3.64873e-04 DIIS/ADIIS
   @DF-RKS iter  95: -1566.27428672143651   -3.98563e-03   3.63757e-04 DIIS/ADIIS
   @DF-RKS iter  96: -1566.31729719175746   -4.30105e-02   3.47321e-04 DIIS/ADIIS
   @DF-RKS iter  97: -1567.23142540517097   -9.14128e-01   3.21628e-05 DIIS
   @DF-RKS iter  98: -1566.31314783377638    9.18278e-01   3.48605e-04 DIIS/ADIIS
   @DF-RKS iter  99: -1566.32974119762684   -1.65934e-02   3.48657e-04 DIIS/ADIIS
   @DF-RKS iter 100: -1566.34244197116936   -1.27008e-02   3.50019e-04 DIIS/ADIIS

PsiException: Could not converge SCF iterations in 100 iterations.

  Failed to converge.

I could run it for longer than a 100 iterations but it does not seem like this is going to help.

susilehtola · May 26, 2024, 9:42am

The monomer not converging is already very suspicious, and it’s no wonder that the non-bonded dimer also doesn’t converge. Are you sure you have the right charge state and multiplicity? It is hard to comment on the issue without knowing more about the system. You could try to run a CASSCF calculation and see whether the system has polyradical character. Depending on the problem, changing the functional to one that has 100% long-range exchange like wB97M-V might also help.

sbembenek · May 26, 2024, 7:34pm

thanks for the info. Actually, the monomer does converge. And the other tautomer of the dimer (where the acid is protonated and the amine is deprotonated) also converges. It’s this tautomer of the dimer, where the acid has been deprotonated giving it a -1 charge and the amine has been protonated to give it a +1 that is not converging. The system would have an overall charge of 0 and the multiplicity would 1. I’ll try the new basis set. thanks again!

jmisiewicz · May 29, 2024, 2:53am

I have a very bad feeling. Could you please converge with HF and report the molecular orbital energies for LUMO-5 to HOMO+5?

sbembenek · May 29, 2024, 3:49am

Something like this (which would be comparable to the inputs I’m using):

set scf_type df
set basis 6-31G(d)
set reference rhf
set s_tolerance 1e-9

set optking { opt_coordinates = cartesian }
set ensure_bt_convergence true

mol.update_geometry()
mol.symmetrize(1e-3)

e, wfn = optimize('HF', return_wfn=True)

mol.print_out()
print_variables()

jmisiewicz · May 29, 2024, 11:30am

I just need an HF orbital energies at a geometry where B3LYP cannot converge SCF. This is not about geometry convergence.

sbembenek · May 29, 2024, 7:07pm

Got it.

Here’s the HF orbital energies:

119A     -0.291266   
120A     -0.256456  
121A     -0.245427   
122A     -0.238691   
123A     -0.225714  
124A     -0.219924  (HOMO)

    Virtual:                                                              

125A     -0.007588  (LUMO)
126A      0.016528   
127A      0.023819  
128A      0.068447   
129A      0.070236   
130A      0.083415

jmisiewicz · May 29, 2024, 11:20pm

That’s bad. I’ll consult with a colleague tomorrow, but I fear that the answer here is “diagonalization based SCF solvers will not be able to converge this; Psi4 needs more SCF technology.”

sbembenek · May 29, 2024, 11:41pm

Help me understand what you’re seeing as ‘bad’.
thanks.

jmisiewicz · May 31, 2024, 11:45am

Negative LUMO is always questionable, and it’s known that unusual protonation states and DFT convergence issues can go together.

Could you post the HOMO LUMO for a single monomer (in the proper protonation state) as well? Once with HF and once with B3LYP.

jmisiewicz · June 4, 2024, 3:18am

And now that the preprint is up, this scenario is what I’m afraid of.

susilehtola · June 16, 2024, 11:00am

… and this is exactly what OpenOrbitalOptimizer is meant for!

sbembenek · June 20, 2024, 2:10am

Thanks for the preprint. Indeed, this seems like the issue.

My level of theory was b3lyp-d3bj/aug-cc-pVTZ. And even though ADIIS was able to force a convergence on the density for the geometry, it looked ‘too shaky’ to me. However, I still went ahead and calculated some properties of interest, like the polarizability.

I did see some of issues the preprint notes. In particular, the unphysically small HOMO-LUMO gap. I also got a ridiculously high polarizability.

In general, I am always concerned about optimizing systems with strong Coulomb interactions in vacuum. I see the preprint noted that using implicit solvent might help. Honestly, this would have been my preferred conditions, but I find the solvent method in psi4 too computationally expensive – perhaps I am doing something wrong here(?)

After reviewing the preprint and some other references, and noting the suggestion by @susilehtola to use a range-separated functional, I went with CAM-B3LYP-D3BJ/aug-cc-pVTZ. This converged the density very smoothly over the geometry – even with DIIS alone.

Look forward to your thoughts!