# Curse of the nano-Hartree

5 years 1 week ago #538 by Nike
Replied by Nike on topic Curse of the nano-Hartree
Dear Mihaly,
Thanks for replying so fast!

1) Perhaps MOLPRO is also using the AO Fock Matrix, because it agrees with CFOUR to ~1 nHartree. The ROHF in MOLPRO was converged with an energy tolerance of 1e-16 ( here ), so the SCF cannot really be converged any better, and the AO Fock Matrix is probably yielding the correct MO Fock matrix and correct correlation energies to sub-nHartree precision.

With CFOUR the SCF was converged using SCF_CONV=10, which can be achieved with 14,372 iterations, with "largest density difference" being only: 0.9669940548D-10. The energy was converged to 0.1 pHartree at iteration 14, and the next 14,000 iterations were just jiggling around the orbitals until the largest density difference went below 10^-10:
14     -14.864012799596406 0.3904492329D-04
14371  -14.864012799596448 0.2061507832D-07
14372  -14.864012799596420 0.9669940548D-10

This makes me think that converging the SCF in CFOUR won't be able to change very much.

I tried to converge with the largest density difference being < 10^-11, but after 49060 iterations I gave up.

The reason why I didn't post the ROHF results with CFOUR was because it *looks* like the energy doesn't agree with MOLPRO and MRCC, but actually this is because of the Bohr -> Angstrom conversion in CFOUR:
2014 CODATA: 1 Angstrom  = 1/0.529,177,210,67(12)  Bohr (4.17 Ang = 7.880,157,943,159,143  Bohr)
MOLPRO:      1 Angstrom  = 1/0.529,177,209         Bohr (4.17 Ang = 7.880,157,968,027,6825 Bohr)
MRCC:        1 Angstrom  = 1/0.529,177,20859       Bohr (4.17 Ang = 7.880,157,974,133,132  Bohr)
CFOUR:       1 Angstrom  = 1/0.529,177,249         Bohr (4.17 Ang = 7.880,157,372,374,11   Bohr)   

I thought that a difference of 0.0000006 Bohr would not make a difference, but actually if I use MOLPRO's geometry of 7.880,157,968,027,6825 Bohr, we get perfect agreement between CFOUR and MOLPRO at ROHF level:
CFOUR with 7.880,157,372,374 Bohr: -14.864,012,799,596
CFOUR with 7.880,157,968,027 Bohr: -14.864,012,797,989
MOLPR with 7.880,157,968,027 Bohr: -14.864,012,797,989

So actually the 25 pHartree difference between MOLPRO and MRCC at ROHF level might not have anything to do with the Boys function (why would it?), it might be because of the 0.00000001 Bohr difference in geometry due to using difference conversion factors. This was something that I never would have expected, especially for such a shallow state like Li2,1^3\Sigma_u^+, which has only a -333cm-1 depth and an r_e = 4.17 Angstroms.

Anyway, based on the convergence I've shown for the CFOUR ROHF with SCF_CONV=10, do you still think the difference in MO Fock Matrix can yield a 9 nano-Hartree difference in CCSD? This is 1000 times bigger than the difference between CFOUR's ROHF at iteration 14 and iteration 14,000.

49060 iterations was not enough to get convergence with SCF_CONV=11, so maybe we just have to compare the MO Fock matrices somehow. A Root Mean Square difference between matrix elements with SCF_CONV=10, 9, 8..., along with differences in CCSD for those cases might show how much a small change in matrix elements changes the CCSD energy. This would tell us if SCF_CONV=11 will help or not. Is it possible to do this somehow? Otherwise I might guess that there might be some truncation happening when writing the fort.55 just as there was in ovirt.f, but perhaps you have already checked that.

3) Thank you for looking into the AO-MO transformation performance!

With best wishes!!
Nike

Please Log in or Create an account to join the conversation.

5 years 1 day ago - 5 years 1 day ago #548 by Nike
Replied by Nike on topic Curse of the nano-Hartree
Greetings,
I have not been able to get SCF to converge in CFOUR with "largest density difference" smaller than 1e-11 with 49060 iterations, and I don't think it will be possible since it took 14,372 iterations to get it to become smaller than 1e-10. But I can live with the 9 nH difference this seems to cause from the MO integrals giving the wrong MO Fock matrix.

However at ROHF level for the Li atom, I'm getting a rather large error in MRCC vs MOLPRO for 6Z, and it's relevant compared to the size of our discrepancy with experiment. For some reason 5Z and 7Z are perfectly fine.
aug-cc-pCV5Z (< 1 pico-Hartree Boys function error):
-------------------------------------------------------------
ROHF/MRCC (canonical) -7.432 722 692 551 (default)
ROHF/MRCC (semi-can)  -7.432 722 692 551 (default)
ROHF/MOLPRO:          -7.432 722 692 551 (itol=18,scftol=16)
-------------------------------------------------------------
aug-cc-pCV6Z (62,000 pico-Hartree Boys function error ??):
-------------------------------------------------------------
ROHF-MRCC (semi-can)  -7.432 726 059 513 (itol=18, scftol=13)
ROHF-MRCC (canonical) -7.432 726 059 513 (itol=18, scftol=13)
ROHF-MOLRPO           -7.432 726 121 753 (itol=18, scftol=16)
-------------------------------------------------------------
aug-cc-pCV7Z (3 pico-Hartree Boys function error):
-------------------------------------------------------------
ROHF/MRCC (semi-can)  -7.432 726 687 162
ROHF/MRCC (canonical) -7.432 726 687 162
ROHF/MOLRPO           -7.432 726 687 159 

For the Li2 molecule, this huge discrepancy doubles in size:
aug-cc-pCV5Z (25 pico-Hartree Boys function error):
-------------------------------------------------------------
ROHF/MRCC             -14.864 012 797 964 (itol=18,scftol=14)
ROHF/MOLPRO           -14.864 012 797 989 (itol=18,scftol=16)
-------------------------------------------------------------
aug-cc-pCV6Z (125,000 pico-Hartree Boys function error ??):
-------------------------------------------------------------
ROHF/MRCC (semi-can)  -14.864 019 536 409 (itol=18, scftol=13)
ROHF/MRCC (canonical) -14.864 019 536 410 (itol=18, scftol=13)
ROHF/MOLPRO           -14.864 019 661 611 (itol=18, scftol=16)
-------------------------------------------------------------
aug-cc-pCV7Z (28 pico-Hartree Boys function error):
-------------------------------------------------------------
ROHF/MRCC (semi-can)  -14.864 020 697 720 (default settings)
ROHF/MRCC (canonical) -14.864 020 697 720 (default settings)
ROHF/MOLPRO           -14.864 020 697 748 (itol=18,scftol=16)

Clearly there's more than just Boys function error in the 6Z case, so I compared the basis set exponents and contraction coefficients but I really so no difference in them (MOLPRO was copied and pasted from the output file and MRCC should have been taken from the MOLDEN file but it was in scientific notation so to make it easier to compare I just took the numbers from the GENBAS file):
MOLPRO:                           MRCC:
=======                           =======
s 103207.000000 0.000004 0.000001 103207.00000 0.000004 0.000001
s  15664.800000 0.000029 0.000005  15664.80000 0.000029 0.000005
s   3600.840000 0.000149 0.000023   3600.84000 0.000149 0.000023
s   1026.350000 0.000627 0.000098   1026.35000 0.000627 0.000098
s    335.925000 0.002274 0.000355    335.92500 0.002274 0.000355
s    121.409000 0.007358 0.001159    121.40900 0.007358 0.001159
s     47.337200 0.021450 0.003373     47.33720 0.021450 0.003373
s     19.630900 0.055599 0.008978     19.63090 0.055599 0.008978
s      8.547060 0.124328 0.020423      8.54706 0.124328 0.020423
s      3.853030 0.228994 0.040705      3.85303 0.228994 0.040705
s      1.780350 0.323811 0.064625      1.78035 0.323811 0.064625
s      0.837931 1.000000              0.837931 0.000000 0.000000
s      0.390527 1.000000              0.390527 0.000000 0.000000
s      0.097776 1.000000              0.097776 0.000000 0.000000
s      0.044933 1.000000              0.044933 0.000000 0.000000
s      0.020602 1.000000              0.020602 0.000000 0.000000
s      0.006100 1.000000              0.006100 0.000000 0.000000
s     49.390700 1.000000             49.390700 0.000000 0.000000
s     26.340100 1.000000             26.340100 0.000000 0.000000
s     14.047200 1.000000             14.047200 0.000000 0.000000
s      7.491400 1.000000              7.491400 0.000000 0.000000
s      3.995100 1.000000              3.995100 0.000000 0.000000

p     32.644300 0.000229             32.644300 0.000229
p      7.758170 0.001768              7.758170 0.001768
p      2.461690 0.007560              2.461690 0.007560
p      0.869615 0.023322              0.869615 0.023322
p      0.327205 1.000000              0.327205 0.000000
p      0.134815 1.000000              0.134815 0.000000
p      0.059995 1.000000              0.059995 0.000000
p      0.028101 1.000000              0.028101 0.000000
p      0.013320 1.000000              0.013320 0.000000
p      0.003900 1.000000              0.003900 0.000000
p     51.457700 1.000000             51.457700 0.000000
p     23.052200 1.000000             23.052200 0.000000
p     10.327000 1.000000             10.327000 0.000000
p      4.626300 1.000000              4.626300 0.000000
p      2.072500 1.000000              2.072500 0.000000

d      0.825300 1.000000              0.825300
d      0.437300 1.000000              0.437300
d      0.231700 1.000000              0.231700
d      0.122800 1.000000              0.122800
d      0.065000 1.000000              0.065000
d      0.034500 1.000000              0.034500
d     34.914700 1.000000             34.914700
d     15.004700 1.000000             15.004700
d      6.448300 1.000000              6.448300
d      2.771200 1.000000              2.771200

f      0.679100 1.000000              0.679100
f      0.358800 1.000000              0.358800
f      0.189500 1.000000              0.189500
f      0.100100 1.000000              0.100100
f      0.052900 1.000000              0.052900
f     23.047600 1.000000             23.047600
f      8.994200 1.000000              8.994200
f      3.509900 1.000000              3.509900

g      0.419900 1.000000              0.419900
g      0.230900 1.000000              0.230900
g      0.127000 1.000000              0.127000
g      0.069900 1.000000              0.069900
g     15.913400 1.000000             15.913400
g      5.499300 1.000000              5.499300

h      0.337100 1.000000              0.337100
h      0.183500 1.000000              0.183500
h      0.099900 1.000000              0.099900
h     10.203500 1.000000             10.203500

i      0.292700 1.000000              0.292700
i      0.146200 1.000000              0.146200

Actually now I have directly compared the exponents and contraction coefficients from the MOLPRO output with the MOLDEN file generated in MRCC and after putting the two into arrays, Matlab says the two arrays are the same: isequal(MRCC,molpro) gives a value of 1. So unfortunately the basis sets seem to be exactly the same. The 1's and 0's for the contraction coefficients in the GENBAS are correct, and the likewise the numbers at the top of the GENBAS describing how many exponents and contraction coefficients.

I cannot imagine any other reason why the 6Z results are so much worse than the 5Z and 7Z ones, apart from maybe one code being more prone to errors due to linear dependencies in the basis set (but I think this would hurt the 7Z results more than the 6Z results and this isn't the case here).

The discrepancy is relevant on the scale of our discrepancy with experiment.
After adding an X2C+SFDC+Breit correction of 0.255cm-1 to each De value, we get, before any CBS extrapolations:
1) -333.58 cm-1 (FCIest./aCV7Z) -> Discrepancy with expt = 0.07 cm-1
2) -333.59 cm-1 (FCIest./aCV6Z) -> Discrepancy with expt = 0.06 cm-1
3) -333.65 cm-1 (empirical with finite-mass effects removed)
4) -333.74 cm-1 (FCIest./aCV5Z) -> Discrepancy with expt = 0.09cm-1

But the 125 nHartree difference between ROHF/MRCC and ROHF/MOLPRO is 0.025cm-1 which is more than 40% the size of the 0.06cm-1 discrepancy with experiment.

The 6Z MRCC output is here (we're just comparing ROHF, because the norm goes bad when MRCC tries to do correlation on a basis set this big): github.com/ndattani/AI_ENERGIES/blob/mas..._betterIntegrals.txt

The MOLPRO output is here:
github.com/ndattani/AI_ENERGIES/blob/mas...etterConvergence.txt

Is there anything else that might be causing this big discrepancy for the 6Z case?

With best wishes!
Nike Dattani
Last edit: 5 years 1 day ago by Nike.

Please Log in or Create an account to join the conversation.

5 years 1 day ago #551 by kallay
Replied by kallay on topic Curse of the nano-Hartree
Nike,
Can you please post the GENBAS file with the aug-cc-pCV6Z basis you use?

Best regards,
Mihaly Kallay

Please Log in or Create an account to join the conversation.

5 years 1 day ago #552 by Nike
Replied by Nike on topic Curse of the nano-Hartree
Dear Mihaly,
Here is the GENBAS file I used, which contains the basis named "LI:aCV6Z-KOPUT":
github.com/ndattani/AI_ENERGIES/blob/master/GENBAS

With best wishes!
Nike

Please Log in or Create an account to join the conversation.

5 years 1 day ago #555 by kallay
Replied by kallay on topic Curse of the nano-Hartree
Dear Nike,
I have tested the 6Z calculation with all the integral algorithms implemented in mrcc. I have also run the calculation with cfour. I am getting -7.432 726 059 511 in every case. So something might be wrong with your molpro calculation.

Best regards,
Mihaly Kallay

Please Log in or Create an account to join the conversation.

5 years 1 day ago - 5 years 1 day ago #556 by Nike
Replied by Nike on topic Curse of the nano-Hartree
Dear Mihaly,
Thank you so much for taking the time to test it.
The agreement with CFOUR seems good indication that MOLPRO might be the problem.

MOLPRO had a lower energy at iteration 4, and then never went up:
From github.com/ndattani/AI_ENERGIES/blob/mas...-cc-pCV6Z/MOLPRO.txt:
IT  DDIFF   GRAD         ENERGY     2-EL.EN. DIIS ORB.
1  0.0D+0 0.0D+00     -7.43042075 4.450025   0  start
2  0.0D+0 0.8D-03     -7.43214903 4.508094   1  diag,B
3  0.1D-2 0.4D-03     -7.43272528 4.560007   2  diag,B
4  0.1D-2 0.2D-04     -7.43272612 4.561845   3  diag,B
5  0.7D-4 0.2D-05     -7.43272612 4.561833   4  diag,B
6  0.4D-5 0.2D-06     -7.43272612 4.561833   5  diag,B
7  0.6D-6 0.4D-07     -7.43272612 4.561833   6  fixocc
8  0.1D-6 0.3D-08     -7.43272612 4.561833   7  diag,B
9  0.2D-7 0.3D-09     -7.43272612 4.561833   8  diag,B
10  0.2D-8 0.2D-10     -7.43272612 4.561833   0  orth

All I told it to do was converge integrals to 1e-18 and energy to 1e-16:

geometry={Li}
gthresh,oneint=1e-18,twoint=1e-18
THRESH,ENERGY=1d-16
wf,charge=0,spin=1
rhf

I will check with the MOLPRO developers.
With best wishes!
Nike
Last edit: 5 years 1 day ago by Nike.

Please Log in or Create an account to join the conversation.

Time to create page: 0.030 seconds