Odd permutation applied to orbitals for complex to real B-splines

[jtkrogel]

Uncovered while investigating #4549.

In #4549, a sudden jump in the energy is observed for bilayer graphene, which is largest in the elec-elec energy. @rcclay observed that the charge density in the problem case is asymmetric between the layers.

Here, I printed orbital values directly from QMCPACK to disk (no QMC) and integrated the per orbital charge densities to get per orbital per layer charge contributions.

Similar to @rcclay, I found one layer has more charge (74 elec on one layer, 70 on the other). The charge imbalance is observed with the complex build of QMCPACK. The real build (same inputs!) shows no such charge imbalance at the orbital level.

This suggests a possible bug in how occupations are handled in the complex code while running with real valued orbitals (gamma point).

Complex code results:

orbital_index  orb_charge_layer1  orb_charge_layer2  cum_charge_layer1  cum_charge_layer2
68  1.0 1.0    69.0  69.0
69  1.0 1.0    70.0  70.0
70  2.0 0.0    72.0  70.0
71  2.0 0.0    74.0  70.0   (charge imbalance)
--------------------------------  <- Fermi level
72  0.0 2.0    74.0  72.0
73  0.0 2.0    74.0  74.0
74  1.0 1.0    75.0  75.0
75  1.0 1.0    76.0  76.0

Real code results:

orbital_index  orb_charge_layer1  orb_charge_layer2  cum_charge_layer1  cum_charge_layer2
68  1.0 1.0    69.0  69.0
69  1.0 1.0    70.0  70.0
70  1.0 1.0    71.0  71.0
71  1.0 1.0    72.0  72.0  (charge balance)
-------------------------------- 
72  1.0 1.0    73.0  73.0
73  1.0 1.0    74.0  74.0
74  1.0 1.0    75.0  75.0
75  1.0 1.0    76.0  76.0

[jtkrogel]

Investigating further, I computed overlaps between the outputted complex orbitals and the real ones.

Below is a listing of the complex orbitals (index in first column) along with which real orbitals (indices in third column) they have finite overlap with. Blank lines separate different primitive cell k-points.

0 1.0 [2 3]
1 1.0 [2 3]
2 1.0 [4 6]
3 1.0 [4 6]
4 1.0 [5 7]
5 1.0 [5 7]
6 1.0 [8 9]
7 1.0 [8 9]
8 1.0 [10 12]
9 1.0 [10 12]
10 1.0 [11 13]

11 1.0 [11 13]
12 1.0 [14 15]
13 1.0 [14 15]
14 1.0 [16 18]
15 1.0 [16 18]
16 1.0 [17 19]
17 1.0 [17 19]
18 1.0 [20 21]
19 1.0 [20 21]
20 1.0 [22 24]
21 1.0 [22 24]

22 1.0 [23 25]
23 1.0 [23 25]
24 1.0 [26 27]
25 1.0 [26 27]
26 1.0 [28 30]
27 1.0 [28 30]
28 1.0 [29 31]
29 1.0 [29 31]
30 1.0 [32 33]
31 1.0 [32 33]
32 1.0 [34 36]

33 1.0 [34 36]
34 1.0 [35 37]
35 1.0 [35 37]
36 1.0 [40 41]
37 1.0 [40 41]
38 1.0 [42 44]
39 1.0 [42 44]
40 1.0 [43 45]
41 1.0 [43 45]
42 1.0 [46 47]
43 1.0 [46 47]

44 1.0 [48 50]
45 1.0 [48 50]
46 1.0 [49 51]
47 1.0 [49 51]
48 1.0 [52 53]
49 1.0 [52 53]
50 1.0 [54 56]
51 1.0 [54 56]
52 1.0 [55 57]
53 1.0 [55 57]
54 1.0 [58 59]

55 1.0 [58 59]
56 1.0 [60 62]
57 1.0 [60 62]
58 1.0 [61 63]
59 1.0 [61 63]
60 1.0 [68 69]
61 1.0 [68 69]
62 1.0 [70 72]
63 1.0 [70 72]
64 1.0 [71 73]
65 1.0 [71 73]

66 1.0 [74 75]
67 1.0 [74 75]
68 1.0 [79 80]
69 1.0 [79 80]
70 1.0 [81 83]
71 1.0 [81 83]
72 1.0 [82 84]
73 1.0 [82 84]
74 1.0 [85 86]
75 1.0 [85 86]
76 1.0 [87 89]

77 1.0 [87 89]
78 1.0 [88 90]
79 1.0 [88 90]
80 1.0 [91 92]
81 1.0 [91 92]
82 1.0 [93 95]
83 1.0 [93 95]
84 1.0 [94 96]
85 1.0 [94 96]
86 1.0 [97 98]
87 1.0 [97 98]

88 1.0 [0]
89 1.0 [1]
90 1.0 [38]
91 1.0 [39]
92 1.0 [64]
93 1.0 [65]
94 1.0 [66]
95 1.0 [67]
96 1.0 [76]
97 1.0 [77]
98 1.0 [78]

As can be seen, there is a large permutation applied to the orbital list (in addition to sqrt(2)/2 orbital pair mixing for complex to real). In the complex code, the gamma point orbitals appear at the very end of the list.

Possibly this permutation is the root cause of the occupation (and charge imbalance) difference between the real and complex codes.

Since the permutation will in general negatively impact usage of the orbital sets as a basis (as in the 1RDM), I'm labeling this as a bug, even in absence of full proof that the permutation causes the charge imbalance in #4549.

[ye-luo]

Looked at /gpfs/alpine/mat151/proj-shared/qmcpack_issue_4549_2/orbital_images/6-real-stat/vmc_perf_pbe_u_None_3x3x1_1x1x1_4990_ccecparep_legacy_500_500_1.5_50_True_hf_april/einspline.tile_300030001.spin_0.tw_0.g0.bandinfo.dat

at 5.00A real code.

#  Band    State   TwistIndex BandIndex Energy      Kx      Ky      Kz      K1      K2      K3    KmK 
      38       68        8        6    -0.063721  1.5582  0.8996  0.0000  0.6667  0.6667  0.0000      2
      39       70        8        7    -0.063056  1.5582  0.8996  0.0000  0.6667  0.6667  0.0000      2
-------------------------------------------------- above occupied.
      40       72        8        8    -0.063053  1.5582  0.8996  0.0000  0.6667  0.6667  0.0000      2
      41       74        8        9    -0.062382  1.5582  0.8996  0.0000  0.6667  0.6667  0.0000      2

complex code
/gpfs/alpine/mat151/proj-shared/qmcpack_issue_4549_2/4-3.16.9-andes/vmc_perf_pbe_u_None_3x3x1_1x1x1_5000_ccecparep_legacy_500_500_1.5_50_True_hf_april_orbitals/einspline.tile_300030001.spin_0.tw_0.g0.bandinfo.dat

#  Band    State   TwistIndex BandIndex Energy      Kx      Ky      Kz      K1      K2      K3    KmK 
      68       68        4        6    -0.063721  0.7791  0.4498  0.0000  0.3333  0.3333  0.0000      1
      69       69        8        6    -0.063721  1.5582  0.8996  0.0000  0.6667  0.6667  0.0000      1
      70       70        4        7    -0.063056  0.7791  0.4498  0.0000  0.3333  0.3333  0.0000      1
      71       71        8        7    -0.063056  1.5582  0.8996  0.0000  0.6667  0.6667  0.0000      1
-------------------------------------------------- above occupied.
      72       72        4        8    -0.063053  0.7791  0.4498  0.0000  0.3333  0.3333  0.0000      1
      73       73        8        8    -0.063053  1.5582  0.8996  0.0000  0.6667  0.6667  0.0000      1
      74       74        4        9    -0.062382  0.7791  0.4498  0.0000  0.3333  0.3333  0.0000      1
      75       75        8        9    -0.062382  1.5582  0.8996  0.0000  0.6667  0.6667  0.0000      1

So far from the builder output, I think the code picked consistent orbitals between real and complex builds.
@jtkrogel I'm a bit lost in the files. Could you make a _3 folder containing only the runs that you used for investigation.