Discrepancy in Z₂ Invariant Results for Bi₂Se₃ Using WannierTools

[alok_shukla1]

Dear Experts,
I have performed Z₂ invariant calculations for Bi₂Se₃ using WannierTools (with VASP). For reference, I followed the example provided in the official documentation: https://wannier-tools.readthedocs.io/en ... i2Se3.html

In that example, the Z₂ invariants for the six time-reversal invariant momentum (TRIM) planes are reported as:

k₁ = 0.0, k₂-k₃ plane: z₂ = 1

k₁ = 0.5, k₂-k₃ plane: z₂ = 0

k₂ = 0.0, k₁-k₃ plane: z₂ = 1

k₂ = 0.5, k₁-k₃ plane: z₂ = 0

k₃ = 0.0, k₁-k₂ plane: z₂ = 1

k₃ = 0.5, k₁-k₂ plane: z₂ = 0

However, in my case, following the same procedure, I obtained:

k₁ = 0.0, k₂-k₃ plane: z₂ = 1

k₁ = 0.5, k₂-k₃ plane: z₂ = 0

k₂ = 0.0, k₁-k₃ plane: z₂ = 1

k₂ = 0.5, k₁-k₃ plane: z₂ = 0

k₃ = 0.0, k₁-k₂ plane: z₂ = 0

k₃ = 0.5, k₁-k₂ plane: z₂ = 1

So, for the k₃ direction, my results are exactly the opposite of those reported in the example.
Here is my POSCAR file:
Bi2Se3
1.0
-2.069 -3.583614 0.000000
2.069 -3.583614 0.000000
0.000 2.389075 9.546667
Bi Se
2 3
Direct
0.3990 0.3990 0.6970
0.6010 0.6010 0.3030
0 0 0.5
0.2060 0.2060 0.1180
0.7940 0.7940 0.8820

I'm unable to figure out why this discrepancy is occurring. I am new to these calculations, and my projection part is coming completely fine (perfect overlapping). I am not getting what could be the possible region of this. If additional files are required, please let me know.

[pedro_melo]

Dear alok_shukla1,

This seems more likely to be a question for the WannierTools community. I will try to help in what VASP is concerned, but I cannot promise that a solution will be found.

Where did you obtain the input files for this example? Did you take them from some reference paper/website and are running a calculation from the very beginning (i.e. DFT ground state), or are you running the example that come with the WannierTools distribution? If it's the former, please paste here the INCAR, POSCAR, POTCAR, and KPOINTS file. I will then have a look to see if there are any issues with generating the file wannier90_hr.dat.

If it's the later and you are just running the examples already provided by WannierTools, then I am afraid I cannot help you, as we do not develop WannierTools.

Kind regards,
Pedro

[alok_shukla1]

Dear Expert,
Thank you very much for your response. Please see and finds the attached files below:

(a) Where did you obtain the input files for this example?
I obtained the POSCAR file from the official WannierTools GitHub repository:
https://github.com/quanshengwu/wannier_tools.

(b) Did you take them from some reference paper/website and are running a calculation from the very beginning (i.e., DFT ground state), or are you running the example that comes with the WannierTools distribution?
Yes, I only took the POSCAR file from the GitHub repository mentioned above. I performed all the subsequent steps from scratch — including structural relaxation, self-consistent field (SCF) calculation, Wannier90 calculations, and Z₂ topological invariant calculations - using WannierTools in combination with VASP.
So, I am not directly running the pre-prepared example from WannierTools; all steps were carried out starting from the POSCAR file.

(c) If it's the former, please paste here the INCAR, POSCAR, POTCAR, and KPOINTS file. I will then have a look to see if there are any issues with generating the file wannier90_hr.dat.
Yes, I am attaching the requested files: INCAR, POSCAR, POTCAR, and KPOINTS. I tried with SOC these calcultions and two different cases I have considered: (a) ISYM = 0, and (b) ISYM = -1.
https://drive.google.com/drive/folders/ ... sp=sharing

Please let me know if any further information is needed.

[pedro_melo]

Dear alok_shukla1,

Sorry, I forgot to ask: which version of VASP are you using?

Best,
Pedro

[alok_shukla1]

Dear Expert,
We attempted these calculations using two different versions of VASP
1. In vasp.5.4.4, after the follwoing line in the OUTCAR, "------------------ aborting loop because EDIFF is reached --------", the job terminates prematurely on its own.

2. In vasp.6.4.3 (with both ISYM = 0 and ISYM =-1), calculations are complted sucessfully but with the reported discrepancy in Z2 values.

Also, I am sharing some addtional settings that were used in our calculations which might be helpful. Please see below:

begin projections
Bi : px; py; pz
Se : px; py; pz
end projections

num_bands = 88
num_wann = 30

num_iter=2000
dis_num_iter=30000
dis_win_min = -4.7375
dis_win_max = 12.0
dis_froz_min = -4.7375
dis_froz_max = 5.26
guiding_centres = .true.
dis_conv_tol = 1.0E-8
conv_tol = 1.0E-5
#hr_plot =.true.
write_hr = true

# This part band related
bands_plot = true
bands_num_points = 50
begin kpoint_path
G 0.0000 0.0000 0.00000 Z 0.00000 0.00000 0.50000
Z 0.00000 0.00000 0.500 F 0.50000 0.00000 0.5000
F 0.50000 0.0000 0.5000 G 0.0000 0.000 0.0000
G 0.0000 0.0000 0.0000 L 0.50000 0.000 0.0000
end kpoint_path

Please let us know if any further informations are needed.

[pedro_melo]

Dear alok_shukla1,

I tried repeating your ground-state calculation. I did not get to the point of computing the Z2 invariant because I think that the issue already shows up when computing the band structure. I attach what I obtained with the wannier90 interface. You can see that there is no signal of an inverted band gap at \(\Gamma\). I tried with different values of num_iter, and there were no significant changes in the bands. Is this also the case in your calculation?

Note: the energies are all shifted downwards so that the Fermi energy is 0 eV.

Best,
Pedro

[alok_shukla1]

Dear Expert,
Thank you so much for performing the calculations.
Yes, indeed I am getting the exact same band structure from the Wannier90 interface in my calculations as well.

However, since Bi₂Se₃ is a well-known topological insulator, I was wondering — shouldn't we see a clear band inversion at Γ? Could it be that the inversion is not captured properly in the Wannierization?

I’d really appreciate your thoughts on this.
Thank you so much once again.

[pedro_melo]

Dear alok_shukla1,

Yes, that was precisely the point I was trying to make. Bi2Se3 is a strong TI, so the signature of band-gap inversion should be present in the bands at \(\Gamma\). The picture I posted here was using num_iter = 0 in the INCAR, so wannier90 should have reproduced the bands exactly has they are from VASP, with no minimisation steps. I will investigate this further, but I would also advise you to check that the bands you get from VASP show the inverted band gap.

Best,
Pedro

[alok_shukla1]

Dear Expert,

Thank you very much for your response.
I have uploaded all my plotted results at the following link:
https://drive.google.com/drive/folders/ ... sp=sharing

Upon reviewing my DFT band structure, I noticed that the inverted band gap at the Γ point is indeed missing.
I would greatly appreciate any further insights or suggestions.