Dear Colleagues,
I am a beginner in the field of computational chemistry, and my interests are primarily focused on the calculation of magnetic states and exchange interactions for 3d- and 4f-block metal coordination compounds. In particular, I am interested in the broken-symmetry DFT calculation for a binuclear Gd(III) complex using ORCA 5 software. The overall procedure seems clear to me, yet the following issue arises.
First of all, I generate a "fragmented" initial guess orbitals via separate DFT calculations for gadolinium centers and organic ligands, and then combine them together through the orga_mergefrag routine. Nevertheless, when I use this guess for the further DFT calculation of the overall molecule (no matter, HS or BS state), the SCF convergence exhibits a sudden increase in energy (ca. 0.0006 Eh) after resetting KDIIS, and eventually, the SCF converges to a higher energy than one could expect.
I use DKH approximation with the appropriate basis sets (DKH-DEF2-SVP and SARC2-DKH-QZVP for organic ligands and Gd, respectively) and PBE0 functional. Additionally, the RIJ-COSX approximation was implied.
Which set of ORCA parameters would be appropriate in my case to achieve a stable and reliable SCF convergence for binuclear gadolinium complexes? Perhaps I should use !SlowSCF or !defgrid3 instructions?
I'd suggest to optimize the geometry in high spin state (UKS) first and follow the following steps.
1. Check if the high spin is correctly assigned on the metal centers (or ligands) by analyzing atomic Mullikan spin density values.
2. Use the high spin wave-function guess (Keywords: MORead, %moinp, FlipSpin) to get the broken symmetry solution (single point calculation.) It should converge to the similar high spin density profile, except for the different signs at two (or more) centers. The energy however is not the same as UKS solution, due to exchange and correlation effects of the paired vs. unpaired electrons. Unless the energy difference is significant (> 5 kcal/mol), we should not hastily arrive at the conclusion regarding most stable spin state.
3. Finally, use the wave-function from the single-point broken-symmetry solution to begin geometry optimization. The optimized geometry is the broken symmetry solution geometry.
4. SCF stability analysis: this is a diagnostics and it is better to check. Follow section section 9.10 on the orca 5 manual. If the analysis deems unstable wave-function, go to step-3.
I also do not understand why you chose basis sets (DKH-DEF2-SVP and SARC2-DKH-QZVP for organic ligands and Gd, respectively). Because you are using RIJCosX approximation, you'll need to need to define def2/JK and SARC2-DKH-QZVP/JK auxiliary basis on your input file. Otherwise by default, the code chooses AutoAux.
I hope my suggestions help.
Dear Bimal Pudasaini
Thank you for your response! Following your suggestions, I have recalculated the complex, yet the problem still retains.
Regarding the question why I use DKH-DEF2-SVP and SARC2-DKH-QZVP basis sets, it is necessary to account for relativistic effects in the case of heavy (Gd) metal centers.
I have followed the same calculation procedure you have suggested with a preliminary preparation of fragmented guess orbitals, yet the SCF convergence behaves strangely in both HS and BS cases.
Namely, the convergence looks like this:
ITER Energy Delta-E Max-DP RMS-DP [F,P] Gradient Norm
15 -24319.4524280248 -0.000000465316 0.00037270 0.00000147 0.0000980 0.000030780
16 -24319.4524282064 -0.000000181626 0.00018519 0.00000076 0.0000676 0.000023497
17 -24319.4524281782 0.000000028245 0.00009817 0.00000054 0.0000419 0.000014140
18 -24319.4524281501 0.000000028027 0.00008242 0.00000039 0.0000225 0.000009913
19 -24319.4524281358 0.000000014355 0.00032630 0.00000073 0.0000167 0.000007415
*** Restarting incremental Fock matrix formation ***
****Resetting KDIIS****
20 -24319.4517824195 0.000645716289 0.00002772 0.00000038 0.0010639 0.000425113
21 -24319.4517836071 -0.000001187560 0.00008226 0.00000068 0.0009575 0.000382645
22 -24319.4517846322 -0.000001025117 0.00015731 0.00000125 0.0007660 0.000306284
At a certain step (highlighted in bold), the energy drastically increases and does not decrease back. What seems to be the problem and how could I tune SCF settings to achieve a stable convergence to the energy minimum?
Although I am not an expert of Orca, there should some flag to set after which step the KDIIS is reset. You might try to increase it to reach convergence before the reset. Or you might try to reduce it to see if the system explores better the energy surface.
On a general discussion however, can I point out that you are possibly over-converging your calculation? The relative error in the calculation of the total energy is around 10^-9 at the reset step, and the step before was along 10^-14. Your presumed numerical accuracy is well below the 50 meV that is the standard for quantum chemistry. So it is completely possible that your SCF calculation was driving you toward a state, but the reset got caught by another local minimum. Their energy difference — 7x10^-4 (eV or Ry?)— is again below chemical accuracy. For large complexes, this should be more than enough and admittedly I would be suspicious on more converged results as they would be likely non reproducible.
Dear Roberto D'Agosta
Thank you so much, your suggestion about adjusting the KDIIS reset seems to be a highly plausible solution.
Actually, in lanthanide chemistry, "sub-chemical" accuracy is mandatory, since the exchange parameters (ca. 0.01 - 1 cm-1) are far below the 50 meV cutoff, and thus I am concerned about the correctness of my calculation.
Thank you once again, I will try to adjust the KDIIS resetting!
How you achieve that depends in too many factors: computing time available, hardware involved, OS involved, even the point at the universe where you are located. Try by example, to run your calculation using Jaguar and an octocore Mac Pro ("trash can") model as regerence or standard. If you have an AppleSeed cluster, that time will be decresed following the Moore laws. If you have a Beowulf cluster, that time will be increased or viceversa. Hope this long reasoning could help you!!!