Energy difference between total energies for wavefunction using hybrid C2R and C2C representation of spline

[fyzibxm]

Describe the bug
Difference in total energies calculated using hybrid C2R and C2C representation of spline of wavefunction. First observed for rather large bilayer heterostructure (GrMoS2) calculating expected energy difference in DMC between different states, GS using real hybrid and EX using complex hybrid representation comparing to non hybrid differences, subsequently reduced just to small graphene supercell (12 unit cell) and comparing just GS for real and complex with and without hybrid rep for VMC energies.

To Reproduce
All runs are copied to OLCF /lustre/orion/mat201/world-shared/brndiar/hybrid-bug
qmca energies (eV):

difference for hybrid rep
complex/hybrid/qmc series 1 -3699.309368 +/- 0.018067 1866.704581 +/- 1.300452 0.5046
complex/hybrid-patched/batch-4.2/qmc series 1 -3699.476471 +/- 0.012745 2128.588869 +/- 8.075983 0.5754
complex/hybrid-patched/qmc series 1 -3699.489952 +/- 0.012274 2118.320382 +/- 2.598431 0.5726
real/hybrid/qmc series 1 -3699.855098 +/- 0.016589 2335.605611 +/- 5.976172 0.6313

no difference without hybrid rep
complex/qmc series 1 -3699.733354 +/- 0.011550 1827.619036 +/- 1.790855 0.4940
complex/patched/qmc series 1 -3699.956620 +/- 0.010278 2091.714849 +/- 12.463772 0.5653
real/batch-4.2/qmc series 1 -3699.934053 +/- 0.011538 2081.696419 +/- 5.166362 0.5626
real/qmc series 1 -3699.951925 +/- 0.011387 2077.847896 +/- 3.287861 0.5616

energy for real with and without hybrid rep are not expected to be equal because hybrid parameters was just set and not properly chosen to minimize energy and variance

parameter for hybrid for carbon was

            <parameter name="cutoff_radius" > 1 </parameter>
            <parameter name="lmax" >          4 </parameter>

meshfactor 0.5
totalwalker count was set to 112000 to directly see difference in 1sigma of variance and not just statistically.

Expected behavior
No energy difference within error bars between real and complex representation for system where state respect k -k symmetry

System:
Leonardo CPU partition
Tested for 3.17.1 legacy and 4.2 batch versions.
I have put all calculation on OLCF storage /lustre/orion/mat201/world-shared/brndiar/hybrid-bug

Additional context
in qmca results there are calculation called patched. During performing complex calculation for graphene I have found that threshold for state ordering operator in src/QMCWaveFunctions/BandInfo.h is to large to occupied K K' valence band properly (see spline bandinfo file where ordering is K VB K CB K'VB K'CB expected is K VB K'VB K CB K'CB, bands 46 to 49) that's why all complex calculation with large threshold do not calculated GS as real calculation where K VB K' VB is replaced with twice K VB due to k -k symmetry. And for such small system, gap on graphene is not properly closed due to finite size, that's why there is difference. For patched calculations I have changed the comparison threshold to 1e-12 which produced proper GS. I'm not sure if this is also a significant new issue.

[ye-luo]

-3699.309368 Hartree in energy with 1866.704581 variance indicates really poor quality of wavefunction.
1866.704581/-3699.309368 is around 0.5 that is really bad. I would control the ratio at least below 0.05.
Most of the calculation I worked with have ratios around 0.02 using DFT orbitals and J1 and J2.
Is your calculation using both one and two-body Jastrow? If the answer is yes, you'll need to double check the planewave cutoff in DFT calculation. I'm afraid that you had it too low.
When WF is in poor quality, hybrid rep can potentially cause issues making fixed node and fixed phase approximation showing large discrepancy. So please try to see if orbitals can be improved helping both accuracy and efficiency.

[ye-luo]

Regarding, a tighter threashold results different state selection and hence the energy value. I'll take a closer look.

[fyzibxm]

-3699.309368 Hartree in energy with 1866.704581 variance indicates really poor quality of wavefunction. 1866.704581/-3699.309368 is around 0.5 that is really bad. I would control the ratio at least below 0.05. Most of the calculation I worked with have ratios around 0.02 using DFT orbitals and J1 and J2. Is your calculation using both one and two-body Jastrow? If the answer is yes, you'll need to double check the planewave cutoff in DFT calculation. I'm afraid that you had it too low. When WF is in poor quality, hybrid rep can potentially cause issues making fixed node and fixed phase approximation showing large discrepancy. So please try to see if orbitals can be improved helping both accuracy and efficiency.

all is in eV energy and also variance, transferring to ha gives 0.02 ...

[fyzibxm]

So as I wrote in last comment my qmca results were in eV to better see difference. Transferring to ha gives normal ratio 0.02. So orbitals from DFT are proper (cutoff 300 Ry for carbon). As I understand from code tracing there is no large difference in hybrid C2R and C2C. I suppose that arithmetic for Y is correct so I have tried to look to radial parts (I have extracted splineData_r and splineData_i) and found two strange things:
One is that for C2R orbitals in bandgroup in initialize_spline_pio_gather is not energy ordered but there is this parts with BandIndexMap and also Dets is invariant wrt to order of columns/rows unless wrong occupation happened so I suppose this is not a problem for any evaluation of energy parts.
Second is I think the smoking gun of difference between C2R and C2C but I do not find the reason. All radial parts for bands (exactly matching twists due to fact that in C2C bandgroup is energy ordered) from C2R in C2C are same but -k parts completely differ. So k -k pairs is not equivalent in radial parts after transforming for C2C. On plotting they look smooth and reasonable but I was unable to find, if it is correct, to which bands or Ylm-s inside it could belongs with respect to this miss order of C2R case of bandgroup.
So I will try to dig it deeper but it is complicated without proper knowledge of code structure and flow ...

[fyzibxm]

Update to k -k problem. I have changed the direction of for cycles which choose distinct twists in EinsplineSetBuilderCommon.cpp to be able to pick k or -k groups in C2R to check them. Their radial parts are as in C2C and different between k and -k series. But what was also different was the total energy of these two C2R runs. Of course without hybrid representation there is no energy difference between those two different C2R runs as obvious ... I have increased 16 times the walker count to be able to see also energy parts and found that there is significant correlation in kinetic energy difference wrt total energy difference. Also variance have changed. Here is results for cutoff_radius 1.2, lmax 4 and meshfactor 0.5 where difference in much large as was for 1, 4 and 0.5 combination (all in eV ! again):

C2R -k
Using 1 copies of twist angle [ 0.000,  0.000,  0.000]
Using 1 copies of twist angle [ 0.500,  0.500,  0.000]
Using 2 copies of twist angle [ 0.667,  0.667,  0.000]
Using 2 copies of twist angle [ 0.833,  0.833,  0.000]
Using 1 copies of twist angle [ 0.000,  0.500,  0.000]
Using 2 copies of twist angle [ 0.167,  0.667,  0.000]
Using 1 copies of twist angle [ 0.500,  0.000,  0.000]
Using 2 copies of twist angle [ 0.667,  0.167,  0.000]
LocalEnergy           =        -3699.6262 +/-           0.0032   1.715e-03    1.0
Variance              =           3334.53 +/-             3.14   1.671e+03    1.0
Kinetic               =          2866.933 +/-            0.045   1.178e-01    3.0
LocalPotential        =         -6566.559 +/-            0.045   1.167e-01    3.0

C2R k
Using 1 copies of twist angle [ 0.500,  0.000,  0.000]
Using 2 copies of twist angle [ 0.333,  0.833,  0.000]
Using 1 copies of twist angle [ 0.000,  0.500,  0.000]
Using 2 copies of twist angle [ 0.833,  0.333,  0.000]
Using 1 copies of twist angle [ 0.500,  0.500,  0.000]
Using 2 copies of twist angle [ 0.333,  0.333,  0.000]
Using 2 copies of twist angle [ 0.167,  0.167,  0.000]
Using 1 copies of twist angle [ 0.000,  0.000,  0.000]
LocalEnergy           =        -3700.0379 +/-           0.0032   1.507e-03    1.2
Variance              =           2463.03 +/-             1.65   3.746e+02    1.2
Kinetic               =          2866.462 +/-            0.054   1.009e-01    4.9
LocalPotential        =         -6566.500 +/-            0.052   9.762e-02    4.7

C2C
qmc series 1
LocalEnergy = -3698.3454 +/- 0.0039 1.461e-03 1.7
Variance = 2380.26 +/- 7.35 8.865e+03 1.0
Kinetic = 2868.264 +/- 0.032 8.117e-02 2.0
LocalPotential = -6566.609 +/- 0.031 8.078e-02 2.0

[prckent]

Are these determinant only no-drift VMC runs? If so, analysing the components of the total energy might reveal where the difference lies. e.g. If we are lucky and the e-e terms match, we would know that the underlying probability distribution is correct and presumably the k.e. is wrong, implying the derivatives in the hybrid rep are different. Most likely everything is bad/different implying the wavefunction values are different. To save cost, potentially these differences can be studied with greatly reduced electron count (same wavefunctions otherwise).

[fyzibxm]

I think I have figured out why this is happening. It looks like this is not an error but a "feature". And the reason is implementation in QMCPACK, but the source of the problem is not an error in implementation (in transformation of orbitals, or in energy evaluation) but lies in what is transformed. pw.x do not strictly take care of "physical gauge fixing" for single particle orbitals like complex conjugation for k -k SPO sets or any other space symmetries. So quality of transformation could vary because one lmax or cutoff could be sufficient for one series of SPO (e.g k) but insufficient for other SPO series (e.g -k) and only for large lmax or small cutoff limit hybrid transformation is robust enough (see my previous results where after increasing cutoff and fixing lmax difference was larger and always disappeared for very large lmax, not shown here). So any calculations using hybrid representation which somehow use different SPO series or deal with R2C and C2C, like in twist averages or energy differences, could be done with care or we need to develop a code to find unitary matrix to produce SPO in "physical gauge" which will satisfy all symmetries which we want to use. So in my case all calculations have to be done in C2C even for GS where k or -k set of SPO could save time and memory, comparing to excitations (where only limited could be done in C2R)