Energy inconsistency in Quantum Package to QMCPACK converter

[djstaros]

Describe the bug When a single determinant of natural orbitals (SD-NO's) obtained with the software Quantum Package 2 is converted into QMCPACK input, the VMC energy of this set of orbitals without a Jastrow factor is ~2 Ha higher than the variational PT2 energy of the input SD-NO's using QMCPACK 3.16.9 with batched drivers.

To Reproduce NOTE: The SBATCH files may need to be changed, depending where this code is rerun. Step 0. Download the attached ZIP file (cri3_issue.zip) and uncompress. Step 1. From the top directory, execute 1_save_for_qmcp.sbatch.in or working analogue. Step 2. From the top directory, execute 2_conv4qmc.sbatch.in or working analogue. Step 3. From the top directory, execute 3_vmc.sbatch.in or working analogue. Step 4. From the top directory, execute the command: qmca -q ev vmc.s000.scalar.dat Step 5 From the top directory, execute the command: grep -e "* Energy of state 1" fci.out -m 1

Expected behavior Upon executing step 1, three output files will be generated, the contents of which should match those in /cri3_issue/x_output/1_save_for_qmcpack. Upon executing step 2, several output and XML files will be generated, the contents of which should match those in /cri3_issue/x_output/2_conv4qmc. Lastly, execution of step 3 will give output files that should match those in /cri3_issue/x_output/3_vmc.

Note that the contents of the QP2QMCPACK. files generated in step 2 (with the exception of QP2QMCPACK.h5) should NOT match the vmc. files provided for running step 3; the QMC input files generated in step 2 were renamed, and their parameters manually changed to perform the desired VMC calculation.

Step 4 should yield something like the following... LocalEnergy Variance ratio QP2QMCPACK series 0 -237.342646 +/- 0.004025 47.095968 +/- 0.197978 0.1984 while Step 5 will yield: Energy of state 1 -239.1695046989506. Both numbers are in units of Ha.

System:

• Conversion performed on Cori, VMC calculation performed on Perlmutter
• Other systems where this is reproducible: Polaris

Additional context This followed a workflow in which PySCF orbitals were expanded with CIPSI and diagonalized back to NO's, on which CIPSI was run again. The NO's provided here were the 3rd iteration of such a procedure. The CIPSI calculations weighted ground state- and excited state-contributing determinants equally.

Additionally, it was noticed that the QP2QMCPACK converter generated an input vmc.structure.xml which erroneously has number of up-spin electrons = 41 and number of down-spin electrons = 29, in contrast to the original 38 up-spin and 32 down-spin electrons. This was kept for reproducibility, and it was checked that correcting these two numbers still yields an erroneously high VMC energy.

Several crucial files for reproducing this were NOT uploaded because of their size. I have added text placeholders, so please ask me for files as you may need them and we can transfer them some other way. cri3_issue.zip

[kgasperich]

If the beginning of yout workflow is a PySCF DFT calculation, then I think a good place to start would be to compare single-det energies in the following cases:

1. NSCF energy of single det of KS orbitals in PySCF
2. VMC (no jastrow) energy of the same single determinant (converted directly from PySCF -> QMCPACK)
3. single-det energy in QP2 (e.g. from the beginning of a CIPSI calculation in those orbitals)
4. VMC (no jastrow) energy of the same single determinant (this time converted from QP2 -> QMCPACK)

I'm not sure whether there's a simple function that already exists to obtain the NSCF energy in PySCF, so I would start with the other three cases. Someone might have this already as part of some existing workflow, but if not then I can help with that step.

[prckent]

Is the kinetic energy available? This is a single particle quantity, fairly easy compute, and if it does not agree it is already game over.

[djstaros]

Dear Kevin and Paul,

Thank you for the input. It was noticed that the VMC structure and wavefunction files obtained from the converter had the wrong number of up/down spin electrons, and fixing this still resulted in a 2 Ha disagreement.

An open shell restricted Hartree Fock calculation in PySCF gave a ground state energy of -239 Ha, following which VMC using the convert4qmc wavefunction still yielded a number around -237 Ha. The kinetic energy from PySCF was 104.616 Ha whereas that from VMC was 105.834 Ha.

[prckent]

Thanks for the update. We should discuss how to stop this happening again. Does the converter need an update?

[prckent]

Wondering about the status. Have correct energies been obtained?

[djstaros]

Not yet, but correct energies have been obtained for a planar CrI3 molecule in open boundary conditions. It looks like the issue lies somewhere in the use of periodic boundary conditions; this is still being teased out.

[kgasperich]

I've attached an example workflow that shows agreement between PySCF and QP2, and disagreement between ElecElec energies between QMCPACK (VMC with no Jastrow) and the other two codes. cri3-pyscf-krohf.tar.gz I only included input and text-based output. I can also share the hdf5 output containing orbitals and integrals (~250MB) upon request.

It would probably be good to have a smaller test system that is more well-behaved. This took something like 170 iterations for the ROHF to converge (the attached PySCF output (scf.out) is from a restarted calculation, so it only shows 70 iterations, but the attached input (scf.py) has a larger maximum number of SCF cycles, so that should converge without needing to restart). We've also already seen it converge to a couple of different ROHF minima. We should be able to identify some system with fewer electrons (so that we can get better-converged VMC error bars more quickly) that consistently converges to the same ROHF state in a reasonable number of iterations.

Here are the energy components from this set of calculations: The stdout from run_test.sh should look something like this

QP2 single-det energy components
 Vnn  =    479.051850583662     
 Ven  =   -1473.88579263905     
 Vee  =    649.402408799605     
 Vecp =   0.000000000000000E+000
 T    =    106.342171565587     

PySCF KROHF energy components
E_nuc         =  477.66663853604416
E_madelung    =  1.3852120476176977
E_kin   = 106.34217156558691
E_ne    = -1506.0985507680837
E_ecp   = 32.21275812903489
E_veff  = 650.7876208472238

and here is the output from qmca: (I replaced the IonIon line to show more digits)

vmc/cri3-pbe  series 0 
  LocalEnergy           =         -239.6108 +/-           0.0061 
  Variance              =             25.82 +/-             0.24 
  Kinetic               =           106.339 +/-            0.026 
  LocalPotential        =          -345.950 +/-            0.027 
  ElecElec              =           650.312 +/-            0.030 
  LocalECP              =         -1517.419 +/-            0.067 
  NonLocalECP           =            43.490 +/-            0.033 
  IonIon                =            477.66662609 +/-             0.00 
  LocalEnergy_sq        =          57439.17 +/-             2.87 
  BlockWeight           =           8000.00 +/-             0.00 
  BlockCPU              =           3.72112 +/-          0.00080 
  AcceptRatio           =          0.672895 +/-         0.000029 
  Efficiency            =             32.64 +/-             0.00 
  TotalTime             =            372.11 +/-             0.00 
  TotalSamples          =            800000 +/-                0 

Key takeaways: PySCF and QP2 agree on everything except for where to put the E_madelung term (0.5 * nelec * madelung_const). In PySCF, they add this as a shift in the exchange part of V_eff. In QP2, we just treat it as a constant shift and lump it in with the IonIon energy (probably not the correct way to report it, but this has been the behavior at least as long as I've been working on the project, and nobody has complained about it, so we haven't really had a reason to change it yet). In the PBC version of QP2, our ECP term (including local and nonlocal parts) is combined with our Ion-Elec Coulomb potential, so Ven from QP2 should be compared to E_ecp + E_ne from PySCF.

The VMC kinetic energy matches the PySCF (and QP2) value. The sum of LocalECP and NonLocalECP matches E_ne + E_ecp from PySCF (and Ven from QP2). (I assume that the Ion-Elec Coulomb term is included in LocalECP, but I haven't actually verified this yet)

The IonIon and E_nuc terms don't agree exactly, but the difference is on the order of tens of microHartrees, so I'm not sure it's worth worrying about. We might be losing some precision in a converter somewhere. (There's also a comment in the PySCF function that does the Ewald sum mentioning that a formula in Martin is incorrect; I haven't done any more digging into this but I'd be surprised if we were using the wrong version and getting this close to agreement.)

So the only significant disagreement is in the ElecElec term. The VMC value falls between the PySCF one (which includes the Madelung term) and the QP2 one (without the Madelung term). The difference between these in either direction is an order of magnitude larger than the error bars, so it's not close to agreeing with either of them. This is where it would be useful to have a smaller test system where we could get tighter error bars on this term. It might be off by some amount that would be recognizable if we could get a more precise estimate.