QMCPACK
Sudden jump in binding energy with separation distance
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:
- release version or git commit hash being built : Git branch: develop Last git commit: 89af1d7f49ee536b8046cce24da561858da10bc9-dirty Last git commit date: Fri Nov 4 10:37:03 2022 -0400 Last git commit subject: Merge pull request #4304 from markdewing/orb_rot_hcpBe
- 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
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.
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.
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.
Could be that a bug has been fixed or that there is somehow an issue on Andes...
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.
Correction on the file location:
/gpfs/alpine/mat151/proj-shared/qmcpack_issue_4549_2
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.
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.
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.
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?
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.