Sudden jump in binding energy with separation distance

[kayahans]

Describe the bug
Binding energy of 3x3x1 bilayer graphene jumps by 0.6 eV/atom, suddenly at interlayer separation of 5.0 A.

This behavior is shown below for VMC (no jastrow, PBE orbitals):

By also performing calculations at 4.99 and 5.01 A, the source of the jump can be isolated to the e-e interaction. Units below are eV/atom.

Component	E(5.0)	(E(4.99)+E(5.01))/2	Difference
Total energy	-149.10(2)	-149.79(3)	0.69(3)
Kinetic	117.80(6)	117.77(9)	0.03(9)
Ion-Ion	1169.7759	1169.7769	0.001
Local ECP	-2701.1(2)	-2701.1(3)	0.3(3)
Non-local ECP	11.09(4)	11.08(6)	0.01(8)
Elec-Elec	1253.4(1)	1252.7(2)	0.7(2)

Error shows insensitivity to the variations we have tried:

• Use of hybrid rep
• Lr_dim cutoff (25-100)

The error also appears consistently at every statistical block. Plot below shows averages over 10 blocks for clarity. Therefore, the error could be more likely to reside in the construction of the potential, rather than a feature of the configuration sampled during the random walk.

DFT curve is smooth at these interlayer separations and occupations are either 0 or 1 down to machine precision.

To Reproduce
Steps to reproduce the behavior:

1. release version or git commit hash being built :
   Git branch: develop
   Last git commit: 89af1d7-dirty
   Last git commit date: Fri Nov 4 10:37:03 2022 -0400
   Last git commit subject: Merge pull request Orbital rotation test with hcp Be #4304 from markdewing/orb_rot_hcpBe
2. cmake command: Andes build script

Files needed to reproduce are located at: /gpfs/alpine/mat151/proj-shared/qmcpack_issue_4549
All the results are from a complex build

Expected behavior
Binding curve should be smooth.

System:

• system name: Andes

[prckent]

Yikes! This is a fairly old build but I don't recall us fixing anything that would cause such behavior. Indeed it really looks like the e-e potential is simply computed consistently higher. It will be interesting to see if the problem is reproducible on another computer system. Andes should not be special in any way.

[prckent]

I was not able to reproduce this problem by rerunning your inputs with the current develop version of the code. However I got the -197.3 Ha energy for 4.99, 5.00, 5.01. I also spotted an issue with qmca. Please can you rerun with either latest release 3.16.0 from 31 January or (preferred) the current development version? I reran on nitrogen2 using 32 omp threads and increased the input to 4 steps/block to get the same weight per block as the files you put in proj-shared.

[kayahans]

Thanks Paul for rerunning my inputs. -197.3 Ha makes about 149.13 eV/atom in this supercell, which is looks like within the statistical uncertainty of the value I have found at d=5.0 A. The values I calculated at d=4.99 and 5.01 are lower than 149.13 eV/atom. I will repeat the test with the newer versions of the code, and get back.

[prckent]

Could be that a bug has been fixed or that there is somehow an issue on Andes...

[kayahans]

I am posting the updated binding curve with runs performed on Cades and Andes using QMCPACK versions of 3.15.9 and 3.16.9. On the left is the full binding curve and on the right is the part zoomed in 4.99-5.01 interval.

The files needed to reproduce these energies are located at /ccs/home/kayahan/shared/qmcpack_issue_4549_2/ in Andes/Summit filesystem.

[jtkrogel]

Correction on the file location:

/gpfs/alpine/mat151/proj-shared/qmcpack_issue_4549_2

[prckent]

Full set of data from nitrogen2 run with develop (amdclang cpu, openblas). The statistics are different from than the previous runs. The results are clearly different from Andes, with different trends and one or two unphysical discontinuities. I will rerun with more statistics and different seeds to check for reproducibility.

[ye-luo]

I dug into the *bandinfo.dat which the spline builder outputs. There are 4 states having very close eigenvalues in DFT.
band 0-71 are occupied by the builder by default.

at 5.00A

spin 0
#  Band    State   TwistIndex BandIndex Energy      Kx      Ky      Kz      K1      K2      K3    KmK 
      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
      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
spin 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
      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

TwistIndex BandIndex 4 7 and 8 7 are picked

at 4.99A

spin 0
      70       70        4        7    -0.063064  0.7791  0.4498  0.0000  0.3333  0.3333  0.0000      1
      71       71        4        8    -0.063063  0.7791  0.4498  0.0000  0.3333  0.3333  0.0000      1
      72       72        8        7    -0.063064  1.5582  0.8996  0.0000  0.6667  0.6667  0.0000      1
      73       73        8        8    -0.063063  1.5582  0.8996  0.0000  0.6667  0.6667  0.0000      1
spin 1
      70       70        4        7    -0.063063  0.7791  0.4498  0.0000  0.3333  0.3333  0.0000      1
      71       71        4        8    -0.063063  0.7791  0.4498  0.0000  0.3333  0.3333  0.0000      1
      72       72        8        7    -0.063063  1.5582  0.8996  0.0000  0.6667  0.6667  0.0000      1
      73       73        8        8    -0.063063  1.5582  0.8996  0.0000  0.6667  0.6667  0.0000      1

TI BI 4 7 and 4 8 are picked

distance 5.01A

spin 0
      70       70        4        7    -0.063060  0.7791  0.4498  0.0000  0.3333  0.3333  0.0000      1
      71       71        4        8    -0.063060  0.7791  0.4498  0.0000  0.3333  0.3333  0.0000      1
      72       72        8        7    -0.063060  1.5582  0.8996  0.0000  0.6667  0.6667  0.0000      1
      73       73        8        8    -0.063060  1.5582  0.8996  0.0000  0.6667  0.6667  0.0000      1
spin 1
      70       70        4        7    -0.063060  0.7791  0.4498  0.0000  0.3333  0.3333  0.0000      1
      71       71        4        8    -0.063059  0.7791  0.4498  0.0000  0.3333  0.3333  0.0000      1
      72       72        8        7    -0.063060  1.5582  0.8996  0.0000  0.6667  0.6667  0.0000      1
      73       73        8        8    -0.063059  1.5582  0.8996  0.0000  0.6667  0.6667  0.0000      1

TI BI 4 7 and 4 8 are picked.

I suggested @kayahans to promote electron from 71-72 on both spin channels at distance 5.00A.
In such a way, consistent orbitals are picked to form the single det.

[kayahans]

With @ye-luo's suggestion I used swapped the occupations on 5.00A distance calculation:

            <slaterdeterminant>
                <determinant id="updet" group="u" sposet="spo_u" size="72">
                <occupation mode="excited" spindataset="0" format="band" pairs="1" >
                    8 7 4 8
               </occupation>
               </determinant>
               <determinant id="downdet" group="d" sposet="spo_d" size="72">
                <occupation mode="excited" spindataset="1" format="band" pairs="1" >
                    8 7 4 8
                </occupation>
               </determinant>
           </slaterdeterminant>

The initial figure is updated as below:

Energies at 4.99, 5.00 and 5.01 become statistically indistinguishable.

[jtkrogel]

Something is strange here. In each case (4.99, 5.00, 5.01), the 4 7 and 8 7 orbitals are lowest in KS energy. Only at 5.00 are they picked. Why isn't QMCPACK occupying by lowest energy?

[ye-luo]

Something is strange here. In each case (4.99, 5.00, 5.01), the 4 7 and 8 7 orbitals are lowest in KS energy. Only at 5.00 are they picked. Why isn't QMCPACK occupying by lowest energy?

Good question.