Inconsistent twist directions between electron gas and einspline wavefunctions

[Paul-St-Young]

Sympton: electron gas momentum distribution n(k) is wrong at fractional twist (i.e. not 0 or 0.5)
Direct cause: MomentumEstimator::putSpecial subtracts the twist vector from kpoint grid, whereas HEGGrid::create_kpoints adds the twist vector to the kpoint grid.
Root cause: einspline wavefunction puts negative sign in front of kgrid?

I found this bug when trying to calculate the momentum distribution n(k) of the homogeneous electron gas (HEG) at a fractional twist e.g. (0.3, 0.0, 0.0). I realized that the MomentumEstimator subtracts the twist vector from the Gamma-point kgrid
ikpt[0]=i-twist[0]; ikpt[1]=j-twist[1]; ikpt[2]=k-twist[2];
whereas the electron gas orbital builder adds the twist to the Gamma-point kgrid
PosType k(i0+tw[0],i1+tw[1],i2+tw[2]);

If I add the twist to kgrid in MomentumEstimator, then n(k) becomes correct. (see figure below).
I have also attached a tiny example (2 electrons) short.zip, which could be added as a short test for the bug in question.

Of course, the current MomentumEstimator is correct when used with einspline wavefunction, so I think the bigger question is: why does einspline wavefunction negate the twist direction?

Unfortunately, I am not sure how exactly the einspline wavefunction sets up the kgrid and the twist vectors, because the wavefunction construction is quite complicated. File potentially involved:

• EinsplineSetBuilderCommon.cpp
• EinsplineSetBuilder_createSPOs.cpp
• EinsplineSetBuilderESHDF.fft.cpp
• EinsplineSetBuilder.h
• EinsplineSetBuilderOld.cpp
• EinsplineSetBuilderReadBands_ESHDF.cpp
• EinsplineSet.cpp
• EinsplineSet.h
• Bspline*?

In my opinion, adding twist to a positive kgrid is the intuitive option, so I am tempted to change MomentumEstimator and the einspline wavefunction. On the other hand, changing HEGGrid is a lot easier. How should I proceed?

[prckent]

Thanks for catching this.

I agree with your intuition. We'll have to investigate how convoluted the spline wavefunction builder is though. There is also work going on with periodic Gaussians where this should be consistent.

Is anyone familiar with the twist vector handling in the splines?

[jtkrogel]

Wouldn't be surprised if it was a funny sign definition issue.

[jtkrogel]

Relationship between twist angle and phase in einspline:

EinsplineSetBuilderOld.cpp:608:          double phi = -2.0*M_PI*dot (ru, TwistAngles[ti]);

[jtkrogel]

It looks like the sign is fixed once it is read in from h5, so probably the origin of the sign convention is propagated back into pw2qmcpack.

[ye-luo]

QMCPACK xml input twist for the splines is opposite to the QE input.

[jtkrogel]

@ye-luo Do you know where the inconsistency originates?

[ye-luo]

It was wfconvert and then carried over by pw2qmcpack if I remember correctly. @lshulen could you confirm?

[jtkrogel]

So is there just a line somewhere in pw2qmcpack where we can apply a sign change to fix this? Any anticipated consequences to just changing the sign in the converter so QE and QMCPACK will match?

[Paul-St-Young]

The reduced_k fields in the wavefunction hdf5 file is typically positive, which is consistent with QE convention. The sign change actually first happened in EinsplineSetBuilderESHDF.fft.cpp:303. Here is the surrounding excerpt:

    // Early versions from wfconv had wrong sign convention for
    // k-points.  EinsplineSet uses the opposite sign convention
    // from most DFT codes.
    if (Version[0] >= 2)
      for (int dim=0; dim<OHMMS_DIM; dim++)
        TwistAngles[ti][dim] *= -1.0;

Is it enough to reverse EinsplineSetBuilderESHDF.fft.cpp:303 together with EinsplineSetBuilderOld.cpp:608?

[mcbennet]

I'd like to re-awaken this discussion.

I imagine this inconsistency confuses a lot of newcomers. Is it possible to make QE and QMCPACK match via EinsplineSetBuilderESHDF.fft.cpp and EinsplineSetBuilderOld.cpp as suggested by Paul?

If so, I will have a bit of free time next week and would enjoy this as a toy project, if someone wants to assign me

[jtkrogel]

I'm aware that both MomentumEstimator and MomentumDistribution rely on the current inverted sign, and so would also need to be updated.

From the discussion above, it seems pw2qmcpack would also need an update. Correct @ye-luo?

@anbenali does the LCAO PBC + twists code follow the same inverted convention?

[ye-luo]

@jtkrogel I don't think pw2qmcpack is wrong. If I add [0.0 0.5 0.25] to nscf, after pw2qmcpack I got

$ h5ls -d pwscf_output/pwscf.pwscf.h5/electrons/kpoint_1/reduced_k
reduced_k                Dataset {3}
    Data:
        (0) 4.16333634234434e-17, 0.5, 0.25

I think the issue on the qmcpack side.
In the printout.

Super twist #0:  [   0.00000   0.00000   0.00000 ]
Super twist #1:  [  -0.00000   0.50000  -0.25000 ]

and I need to input [ -0.00000 0.50000 -0.25000 ] or [ -0.00000 -0.50000 -0.25000 ] to hit the Super twist 1. This is bad.

I also noticed that I need to input

twistnum="-1" twist="0.0 -0.5 -0.25"

to specific twist. I'm in favor of removing twistnum and default twist="0.0 0.0 0.0" if there is no input.

[jtkrogel]

Thanks @ye-luo. So the sign inversion is local to QMCPACK. This problem really should be fixed.

An interesting feature is that the Bloch state must be constructed inside QMCPACK with an inverted wavevector, so this is more than just a simple labeling issue. If it was just a labeling issue and a Bloch state w/ the correct wavevector was actually created, it would show up in the momentum distribution since n(k) is a projection of the state onto planewaves. The QMCPACK B-spline states are absent at the correct wavevector, as @Paul-St-Young noticed.

I'm curious to know what the situation is with LCAO's. @anbenali, can you comment?