Numerical problems with relativistic zora calculations

[moritzgubler]

Hi

I recently noticed, that I am not able to converge zora calculations with tight precisions. The problem is the following: consider the input file

world_prec = 1.0e-6
world_unit = angstrom
WaveFunction {
  method = lda
  relativity = zora
}
Molecule {
$coords
H 0.0 0.0 0.0
H 0.9 0.0 0.0
$end
}

This does not converge. It is close though, it almost reaches the threshold after 10 SCF iterations. After that the MO residual stays slightly above the convergence threshold for the next 20 SCF iteration. Then the simulation goes really wrong and the residual starts to diverge. Also, the memory required by mrchem continuously grows and after 47 iterations i stopped the program because it used up the entire memory of my workstation (32 GB).

I used the newest version of mrchem and tested the functionals lda, pbe and pbe0 where the problem is the same. It is also present in all the other molecules I tested: H_2O, NH_4, it is even more severe in the NH_4 molecule. When the world_prec is increased to 1e-4 or the zora method is not used, the calculation converges. Below you find the output of the SCF cycle from the input file from above.

I tried to play around with some SCF keywords to get these calculations to converge, but was not able to do so. Do you have an idea what could go wrong here? Thanks a lot for looking into this.

===========================================================================
 Iter           MO residual             Total energy                Update
---------------------------------------------------------------------------
    0          1.000000e+00          -1.126549655620         -1.126550e+00
    1          5.588320e-02          -1.137476166839         -1.092651e-02
    2          5.683084e-03          -1.137719316527         -2.431497e-04
    3          1.426500e-03          -1.137732864726         -1.354820e-05
    4          8.197909e-04          -1.137724341320          8.523405e-06
    5          7.077225e-04          -1.137733969255         -9.627935e-06
    6          5.908577e-05          -1.137733876118          9.313768e-08
    7          6.918668e-05          -1.137733974038         -9.792065e-08
    8          3.802536e-05          -1.137733939618          3.442023e-08
    9          4.411585e-05          -1.137733973516         -3.389800e-08
   10          2.451983e-05          -1.137733968648          4.867959e-09
   11          4.849671e-05          -1.137733930095          3.855352e-08
   12          5.043205e-05          -1.137733966139         -3.604488e-08
   13          2.160990e-05          -1.137733974271         -8.131049e-09
   14          5.883780e-05          -1.137733878856          9.541477e-08
   15          6.316945e-05          -1.137733974749         -9.589299e-08
   16          1.297690e-05          -1.137733971889          2.859604e-09
   17          3.504807e-05          -1.137733951260          2.062898e-08
   18          3.231851e-05          -1.137733973203         -2.194333e-08
   19          1.141251e-05          -1.137733973642         -4.380603e-10
   20          1.059654e-05          -1.137733974027         -3.849538e-10
   21          1.189108e-05          -1.137733971784          2.242024e-09
   22          2.156554e-05          -1.137733970533          1.251466e-09
   23          2.344845e-05          -1.137733967163          3.370509e-09
   24          3.120666e-05          -1.137733957546          9.616949e-09
   25          3.934829e-05          -1.137733952101          5.444069e-09
   26          3.975666e-05          -1.137733953968         -1.866136e-09
   27          5.046707e-05          -1.137733923887          3.008081e-08
   28          6.774338e-05          -1.137733892736          3.115061e-08
   29          9.342978e-05          -1.137733801676          9.106051e-08
   30          1.294986e-04          -1.137733711446          9.023012e-08
   31          1.567931e-04          -1.137733572297          1.391487e-07
   32          1.830440e-04          -1.137733321878          2.504186e-07
   33          2.267721e-04          -1.137732863967          4.579118e-07
   34          3.672284e-04          -1.137730986239          1.877727e-06
   35          5.095956e-04          -1.137730116688          8.695510e-07
   36          5.585423e-04          -1.137728828327          1.288361e-06
   37          6.677414e-04          -1.137724757660          4.070667e-06
   38          9.842135e-04          -1.137711255339          1.350232e-05
   39          1.440041e-03          -1.137693119222          1.813612e-05
   40          1.688429e-03          -1.137696282315         -3.163093e-06
   41          1.908598e-03          -1.137656400602          3.988171e-05
   42          3.059050e-03          -1.137439252897          2.171477e-04
   43          4.264359e-03          -1.137195631600          2.436213e-04
   44          4.850040e-03          -1.136858487554          3.371440e-04
   45          5.472838e-03          -1.135977026135          8.814614e-04
   46          8.564126e-03          -1.128026930496          7.950096e-03
   47          1.256601e-02          -1.125935204092          2.091726e-03

[moritzgubler]

Here is the NH_4 example where the problem is much more severe:

world_prec = 1.0e-5
world_unit = bohr

WaveFunction {
  method = pbe
  relativity = zora
}

Molecule {
$coords
  N   0.068688       -0.001101        0.003535
  H   0.016955       -1.773313       -0.722080
  H   1.587247        0.909119       -0.752035
  H  -1.479628        0.865089       -0.692427
$end
}

Once again, I killed the program after it used more than 30 GB of memory. Here this happened already after 7 scf cycles.

===========================================================================
 Iter           MO residual             Total energy                Update
---------------------------------------------------------------------------
    0          2.236068e+00         -56.099856264610         -5.609986e+01
    1          2.608680e-01         -56.528455582620         -4.285993e-01
    2          1.092690e-01         -56.560817280211         -3.236170e-02
    3          5.336892e-02         -56.564123235751         -3.305956e-03
    4          1.669575e-02         -56.564797846569         -6.746108e-04
    5          8.402572e-03         -56.564231974860          5.658717e-04
    6          8.388056e-03         -56.564867508999         -6.355341e-04
    7          1.053038e-01         -56.323330957508          2.415366e-01

[ilfreddy]

Just to avoid misunderstandings. You say you run NH_4 but your input clearly shows you are running NH_3. That is OK because NH_4 would be an open shell system, which is not implemented for ZORA + DFT.

As for your convergence issues, we are aware that ZORA is not yet very stable as it is, but I'm indeed surprised to see that you get issues already for such small systems. As for the memory issues, if we exclude a possible memory leak, the issue is related to the additional cusp introduced by the $$\kappa$$ function, coming from the ZORA kinetic energy.
In conclusion we need to look into this.

For the time being what I can suggest to make the SCF converge in your case is the following:

1. Try to experiment with the KAIN history
2. Change the polynomial order
3. Use the low-prec calculations as a starting guess for higher precision calculations.
4. Change functional

If none of these works we definitely need to look into this issue.

[moritzgubler]

Hi Luca

Thanks for your quick reply:) Sorry, I meant to say NH3.

I will try to change the settings and see if it gets better.

About the memory problem: It did not really feel like a memory leak, when I monitored the memory usage, it was constant for quite some time followed by a quick increase when the residual started to rise quickly.

Best,
Moritz

[moritzgubler]

Hi Luca

We tried everything you mentioned except suggestion 3 since I do not know how to do this. Your suggestions help a bit but we still cannot run all of our systems with the desired precision (which is quite high, but this is the reason we would like to use MRCHEM).

[Dheasra]

Hi Moritz, sorry I'm late to the party, but I hope I might be of help here.
I experienced similar issues you have (almost converged but then divergence into oblivion).
To get reliable-ish convergence using ZORA, you need to start with low precision first and use it as a starting guess.
It's quite easy, here are code snippets that gives you the inputs necessary to get an initial guess and then use it for a higher precision job.
First, get the low precision initial guess by adding the following in your input file.

SCF {
  path_orbitals = orbitals              # Path to orbital files
  write_orbitals = true                # Save converged orbitals to file
}

The "path_orbitals" keyword can be omitted and it will just give you the orbitals in a folder called "orbitals" within the directory of the input file, but I added it here for completeness.

Then to use your initial guess, you have to use the following keywords in your actual run's input file

SCF {
  guess_type = mw                  # mw means it starts from a pre-existing function tree file 
}
Files {
  guess_phi_p = PATH/TO/PAIRED/ORBITALS/phi_p     # Path to where you stored your initial guess, don't forget the /phi_p
}

As one can expect, you can repeat this procedure ad nauseam, starting with an extremely low precision initial guess and slowly increasing the precision to whatever your heart desires.

Should that not suffice (and sometimes it doesn't), I noticed that sometimes the default world_size parameter is too large which causes enough numerical instability to cause similar problems to what you experienced. In that case, try lowering fro the default 7 (-ish) to 6 or even 5. It helped a couple of times.

Another source of issues might be the localisation of orbitals. If all else fails, try this:

Response {
  localize = false                      # Use canonical or localized  orbitals
}

Hope that helps with your issue :)

[ilfreddy]

One more comment from my side: are you using point charges as nuclei or smeared nuclei? For 4d elements you should definitely go for smeared nuclei.

[moritzgubler]

We used the default setting which I think is point charges. I did not even see the smearing option, that sounds promising.