Questions about mBJ band and fermi surface

[yujia_teng]

Dear admin,
I have two questions about mBJ band and fermi surface.

For mBJ band:
When computing mBJ band, I found that using 'ALGO = Damped' and 'All' gives slightly different band. Both calculation are converged, and 'All' takes much more steps to converge. I'm wondering why this happens. By comparing to paper, in my case, 'Damped' gives more matched band. The input and OUTCAR of 2nd step below are attached.

And just to make sure the step for mBJ is correct:
1. do PBE scf calculation with modified KPOINTS file (IBZKPT + 0 weight k-points along high symmetry line, can be directly obtained from vaspkit function 251 https://vaspkit.com/tutorials.html#pre- ... functional.), write WAVECAR.
2. Add relevant MBJ tags to INCAR and read WAVECAR, with same KPOINTS file.

Code: Select all

ISTART = 1
#mBJ
METAGGA = MBJ
LASPH = .TRUE.
ALGO = Damped/All
#TIME = 0.05
NELM = 300

For mBJ fermi surface:
Can we still obtain the fermi surface with mBJ correction? Looks like the difference with usual fermi surface calculation is at KPOINTS and reading WAVECAR not CHGCAR. And step is still the two above?

For vaspkit, https://vaspkit.com/tutorials.html#fermi-surface, its method already contains generating a list of points, which should be IBZKPT of scf calculation. So compared to mBJ band KPOINTS, it just doesn't have 0-weighted points along high symmetry path. Compared to usual fermi surface calculation, this seems reasonable.

For pyprocar, https://pyprocar.readthedocs.io/en/stab ... /vasp.html, the step is same as dos calculation. So similarly, I should use the IBZKPT file generated from regular k-mesh like 9*9*9 to compute mBJ dos?

[michael_wolloch]

Dear yujiua_teng,

Thanks for the detailed explanation of your problem and for sending the relevant files and results.

I think the major problem with your setup is the usage of mBJ to analyze the band structure of a metallic system. This functional was developed to correct the underestimation of the bandgap that is systematic in LDA/GGA DFT functionals. It was not developed to give better band curvatures or a better Fermi surface.
I would thus strongly recommend switching to another functional.

I see a few additional issues with your bandstructure calculations that are relevant to some extent for all calculations and/or metaGGAs:

• INCARs: I see that your INCAR files are very long with a lot of commented-out tags. This is not recommended since errors can happen quickly, and it is hard to find if a tag is defined twice or commented out (or in) by mistake.

• POTCARs: Please be sure that the Ho_3 POTCAR you are using is appropriate for the valence that Ho has in your system. Also check if the Bi potential is appropriate. We recommend the Bi_d POTCAR in most cases.

• Magnetism: Are you sure that your magnetic settings are correct? You first select a non-spin polarized calculation by setting ISPIN=1, but then you are specifying LSORBIT=T, which internally sets LNONCOLLINEAR=T and makes VASP ignore ISPIN. I guess this is what you want, but then you also set all magnetic moments to zero using MAGMOM. Are you sure that all moments are zero? It is also recommended to check results with respect to switching symmetry operations off, i.e., ISYM=-1. Maybe your system would not even be metallic with another magnetic structure?

• Electronic optimization settings: Your settings are not optimized for what you are going to do. For direct optimization it is recommended to set PREC=Accurate. Furthermore, ALGO=A can optimize the timestep automatically, which you do not use. If you do not set TIME, the calculation will converge in fewer steps, but still to a different result than ALGO=Damped.

I would also urge you to upgrade to a more modern VASP version if possible. The line search for the conjugate gradient (ALGO=All) algorithm was recently much improved! And there are, of course, a large number of other new features and improvements as well.
Especially useful for bandstructure calculations is the KPOINTS_OPT file.

As for the Fermi surface calculation. You are correct that you will have to sample the Brillouin zone uniformly, and do not need the 0 weigth k-points from the band structure calculations. However, you will have to converge your KPOINTS mesh with respect to the results. Probably you will have to use a pretty dense mesh. And again: mBJ is probably not the correct functional to use.

Please understand that we cannot give any support for tools like vaspkit and pyprocar.

We will look a bit more closely into mBJ and ALGO=Damped, since it has not been tested much. I will get back to you if there is something worth reporting. If any of my comments are unclear, please feel free to respond to this post as well.

Cheers, Michael

[yujia_teng]

Dear Michael,
Thanks for the detailed explanation. I was not aware that mBJ was not developed for metal. But I used that because there're paper that uses mBJ to correct the band of semi-metal from LDA (Actually for this kind of half-Heusler material, all those paper use mBJ to correct band althought it's not insulator). https://journals.aps.org/prb/abstract/1 ... .82.235121 and https://journals.aps.org/prb/abstract/1 ... .82.125208. Also, this paper https://journals.aps.org/prb/abstract/1 ... 108.165154 uses VASP on one of the compound I'm studying. I'm following these paper and am able to reproduce the result. Also the experiment result matches the mBJ band. So I'm not sure if it's still a good idea to switch to another functional.

And thanks for your suggestions on input files. I forgot why I didn't use Bi_d, I will correct that and check it.

• As for the magnetism, I'm just testing non-magnetic state here, so I use Ho-3 and magnetism should be 0. For usual case with Ho-f on valence state, the effect on band near E_F is just spin-split band and is still metal. In soc case, ISPIN is ignored, so I used to don't care about it and just adding LSORBIT tag to INCAR. So is it a better idea to comment out ISPIN in soc case? I set MAGMOM = 9*0 because if I don't set it, in the mBJ steps, the magnetic moment won't converge to 0. For this point, I'm attaching another material YPtBi result here.

  An update that I just found after running VASP 6.5.1 below. So according to the following warning, I think even for non-magnetic materials, it's better to specify MAGMOM=0 not comment it out.

  Code: Select all

   -----------------------------------------------------------------------------
  |                                                                             |
  |           W    W    AA    RRRRR   N    N  II  N    N   GGGG   !!!           |
  |           W    W   A  A   R    R  NN   N  II  NN   N  G    G  !!!           |
  |           W    W  A    A  R    R  N N  N  II  N N  N  G       !!!           |
  |           W WW W  AAAAAA  RRRRR   N  N N  II  N  N N  G  GGG   !            |
  |           WW  WW  A    A  R   R   N   NN  II  N   NN  G    G                |
  |           W    W  A    A  R    R  N    N  II  N    N   GGGG   !!!           |
  |                                                                             |
  |     You requested a magnetic or noncollinear calculation, but did not       |
  |     specify the initial magnetic moment with the MAGMOM tag. Note that      |
  |     a default of 1 will be used for all atoms. This ferromagnetic setup     |
  |     may break the symmetry of the crystal, in particular it may rule        |
  |     out finding an antiferromagnetic solution. Thence, we recommend         |
  |     setting the initial magnetic moment manually or verifying carefully     |
  |     that this magnetic setup is desired.                                    |
  |                                                                             |
   -----------------------------------------------------------------------------
  

• For optimization settings, I will try newer version of VASP. Also would try KPOINTS_OPT method

So at this point, looks like using vasp6 and ALGO=A(or All, IALGO=58 should be same, and only specify 1 of these tags is enough?) is best for mBJ band calculation. I will try it to see if there's any difference.

[yujia_teng]

Dear Michael,
I did some tests with ALGO in newest version 6.5.1 and have some questions.

1. So 'Damped' method still gives matched band with paper, as can be seen in the figure I put in attachment below. But 'All' gives an artifact near the touching point near Gamma point.

Is this because the following reason from vaspwiki? https://vasp.at/wiki/IALGO#Direct_optimization It says that 'All' is best for insulator while 'Damped' for metal.

The preconditioned conjugate gradient (IALGO= 58, ALGO= A) algorithm is recommended for insulators. The best stability is usually obtained if the number of bands equals half the number of electrons (non-spin-polarized case). In this case, the algorithm is fairly robust and foolproof and might even outperform the mixing algorithm.

For small gap systems and for metals, it is however usually required (metals) or desirable (semiconductors) to use a larger value for NBANDS. In this case, we recommend using the damped MD algorithm (IALGO=53, ALGO=Damped) instead of the conjugate gradient.

2. And I'm still not sure if optimization settings is correct. From your suggestion, ALGO=All shouldn't have TIME set. In my test (file attached), indeed it's 10min faster than setting TIME, although it has more steps. I also tried the HFRCUT=-1 tag for ALGO=All, suggestedin step 3 here https://vasp.at/wiki/Band-structure_cal ... unctionals, and there isn't much different in resulted band.

3. The mBJ band from 6.5.1 is different from 5.4.4. See the figures attached, the main difference is at the red box, where 6.5.1 gives a crossing there. Because the algorithm is improved a lot in 6.5.1 compared to 5.4.4 (much much less steps in 6.5.1), and the band in ref paper is 15 years old, the band from 6.5.1 should be more accurate?

Best,
Yujia

[fabien_tran1]

Hi Yujia,

A few points:

mBJ:
As Michael mentioned, mBJ was originally developed for band gap calculations. However, mBJ has also been applied to semi-metals, as you saw. It is difficult to make a general statement about the adequacy of mBJ for metals. This will depend on the studied system and property. But why not if the papers that you mentioned report good results with mBJ and your calculations are of similar type (system and property).

ISPIN:
If LSORBIT=.TRUE., then ISPIN is irrelevant. However, for a clean INCAR I would comment out ISPIN.

MAGMOM:
By default MAGMOM=3NIONS1.0. Therefore, to increase the chance to get a non-magnetic state, one should of course set MAGMOM=3NIONS0.0.

HFRCUT:
HFRCUT concerns the treatment of the singularity in the exact Fock exchange. Thus, this concerns only the Hartree-Fock and hybrid functionals that you are not using. Note that these functionals are very expensive and not recommended for metals.

VASP version:
It is strongly recommended to use a recent version rather than the old 5.4.4.

Commented lines in INCAR:
For modern versions of VASP, the lines should be commented out with # and not with tt.

Difference between ALGO=All and Damped:
There is indeed a (small) artifact in the ALGO=All plot. This would be interesting to see what you obtain with the KPOINTS_OPT file method mentioned by Michael. Maybe using EDIFF=1.0E-7 or 1.0E-8 would remove this artifact?

Difference between 5.4.4 and 6.5.1:
There is indeed a difference. I would, in principle, trust the results obtained with the more recent version of VASP. As mentioned above, this would be interesting to see what KPOINTS_OPT gives.

SCF convergence:
Have you tried ALGO=Normal (the default)? This is a quite general observation that convergence is difficult when LSORBIT=.TRUE.

[yujia_teng]

Dear Fabien,
Thanks very much for your suggestions. This make things quite clear.

And I will switch to newest version 6.5.1. I did the band with the new KPOINTS_OPT method. But I'm not able to extract the data to plot it. The external packages I used before was developed for vasp5, so it's not compatible with the new method. The band extracted from KPOINTS_OPT method with those packages doesn't make sense. And currently the version on my cluster is not compiled with HDF5, so py4vasp cannot be used to plot band. I've asked the admin to compile it. Once it's done I will report results here.

Best,
Yujia

[yujia_teng]

Dear Fabien,
I did some more tests and previous questions should be clear now.

SCF convergence:
Have you tried ALGO=Normal (the default)? This is a quite general observation that convergence is difficult when LSORBIT=.TRUE.

In 5.4.4, the default algorithm diverges in mBJ calculation, so I didn't use it there. But in 6.5.1, it converges well with both old and new KPOINTS_OPT method, and total job time is much faster than Damped and All. The resulted band is also correct, same as Damped result. From wiki, looks like ALGO=All converges faster, but I'm not sure in my case why the default ALGO is faster while All is the slowest.

Difference between ALGO=All and Damped:
There is indeed a (small) artifact in the ALGO=All plot. This would be interesting to see what you obtain with the KPOINTS_OPT file method mentioned by Michael. Maybe using EDIFF=1.0E-7 or 1.0E-8 would remove this artifact?

The artifact in ALGO=All plot is removed by using the new KPOINTS_OPT method, with same EDIFF=1.0E-6. Increasing precision to 1.0E-8 doesn't change band, so the artifact is still present in old method even with EDIFF=1.0E-8. For ALGO=Damped or Normal, old and new method give same band (figure attached below).

But for ALGO=All, the way how 2 bands touch around Gamma (top band minimum slighted shifted away) is still slightly different compared to ALGO=Damped/Normal (parabola min/max head-to-head). Probably this is just due to algorithm difference.

Difference between 5.4.4 and 6.5.1:
There is indeed a difference. I would, in principle, trust the results obtained with the more recent version of VASP. As mentioned above, this would be interesting to see what KPOINTS_OPT gives.

In 6.5.1, the old method gives a crossing between lowest 2 bands shown below at Gamma. With new KPOINTS_OPT method, the gap is reopened a little bit, but actually doesn't change much. The key feature, i.e. how band touches at Gamma point, is same for both versions.

So in conclusion, there's some difference of band among different algorithm and version, but the key feature is same. And we should use the newest version.

Thanks again for your detailed explanation. It's very helpful.

Best,
Yujia

[fabien_tran1]

Hi Yujia,

Thank you for these new tests. Now it seems clear which version (6.5.1) and method for bland plotting (the new KPOINTS_OPT) should be used. With this choice, the differences between the different algorithms for ALGO seem rather unimportant. Actually, were the three band structures shown in your last message obtained with 6.5.1 and KPOINTS_OPT?

Have you compared the PBE band structures obtained with "All" and "Damped"? Are they exactly the same or is there also some difference?

As a side remark, by looking at the "total charge" at the end of the OUTCAR files, one can also notice differences between the "All" and "Damped" algorithms, consistent with the observed differences in the band structures.

[yujia_teng]

Dear Fabien,
Yeah the bands in last message is from 6.5.1 with KPOINTS_OPT method.

Have you compared the PBE band structures obtained with "All" and "Damped"? Are they exactly the same or is there also some difference?

I just tried it. 'All' gives exactly same result as 'Normal'. The "total charge" is same for both. However, Damped is very difficult to converge. I tried tuning mixing parameter and reading WAVECAR but still hard to converge. So I don't want to spend more time on this.

I also looked at the "total charge" at the end of OUTCAR, it's indeed different between All and Damped. And same for Damped and Normal, this should be why Damped band is same as Normal one. Below is the result from 6.5.1 mBJ new KPOINTS_OPT method:

All:

Code: Select all

 total charge

# of ion       s       p       d       f       tot
--------------------------------------------------
    1        2.207   6.354   1.053   0.050   9.664
    2        1.422   1.581   9.946   0.000  12.948
    3        0.704   0.567   8.173   0.000   9.443
--------------------------------------------------
tot          4.332   8.502  19.172   0.050  32.056

Damped:

Code: Select all

 total charge

# of ion       s       p       d       f       tot
--------------------------------------------------
    1        2.195   6.320   1.057   0.049   9.620
    2        1.323   1.629   9.947   0.000  12.900
    3        0.582   0.532   8.158   0.000   9.271
--------------------------------------------------
tot          4.100   8.480  19.162   0.049  31.791

Normal:

Code: Select all

 total charge

# of ion       s       p       d       f       tot
--------------------------------------------------
    1        2.195   6.320   1.057   0.049   9.620
    2        1.323   1.629   9.947   0.000  12.900
    3        0.582   0.532   8.158   0.000   9.271
--------------------------------------------------
tot          4.100   8.480  19.162   0.049  31.791

So for usual band calculation the usual setting is fine. But for mBJ, it's better to use the newest version and new method.

[fabien_tran1]

Thanks for these additional data.

The most probable reason to explain the slightly different mBJ results obtained with the "ALL" method is the fact that "ALL" relies much more than the others (e.g. ALGO=Normal) on the consistency between the energy and potential. With METAGGA=MBJ this consistency is not satisfied, since mBJ is only a potential (with no existing corresponding energy expression), while LDA is used for the energy as an arbitrary choice. Thus, the use of ALGO=All with mBJ should be done with great care.

[yujia_teng]

Sorry that I'm still confused and have some questions when I turn on the magnetization of my Ho compound system. The results and input/output files of 2nd case (exp lattice) is attached.

So I'm comparing the band from 5 methods for 3 different lattice constants of same system: 5.4.4 Damped, 6.5.1 old Normal, 6.5.1 new Normal,6.5.1 old Damped,6.5.1 new Damped. 'New' means the new KPOINTS_OPT method, 'old' means the previous IBZKPT + 0-weighted kpoints along high-sym path. Because from previous results, and some in the current attached results, All doesn't give correct band, so I didn't test it here. 3 lattice constants are: DFT-relaxed > exp lattice > 2% strained exp lattice.

In the pdf file, I made a table with checks and crosses. Crosses means I don't think the band make sense. I have some observations from these tests:

• Damped works well in all cases. Old method gives quite similar result between 5.4.4 and 6.5.1, but is different from new method.

• Normal sometimes works but sometimes not. I don't understand why. In page 1, DFT-relaxed lattice constant, 'old Normal' doesn't make sense because the top bands is weird. But 'new Normal' band is almost same as 'new Damped'. Their total charge differs only 0.001.

  However, for page2 exp lattice, 'old Normal' gives quite similar band as '5.4.4 Damped' and '6.5.1 old Damped' while 'new Normal' band doesn't make sense, because of top bands become 3.

  In page3, both Normal doesn't work.

• In my calculation, sometimes 'ALGO=Normal' can't converge if it reads WAVECAR, but can converge if start from scratch. I also don't understand this.

So from above result, and a paragraph in Direct optimization:

IALGO=53-58: Treat total free energy as variational quantity and minimize the functional completely selfconsistently[3][4][5].
These algorithms have been carefully optimized and should be selected for Hartree-Fock type as well as meta-GGA functionals. The present version is rather stable and robust even for metallic systems.

looks like 'ALGO=Damped' should be used. Also, as mentioned by you before, it's recommended to use 6.5.1 with new KPOINTS_OPT method. But I have some questions below:

• I want to ask if is there a reliable way to know which method should be used. From above tests and my next question, I'm still not 100% sure which method should I use.

• I need to plot the fermi surface of the mBJ band. The way to I used to compute fermi surface is same as the old method (KPOINTS is IBZKPT from scf, without 0-weighted kpoints along high-sym path), and the new KPOINTS_OPT method is a little bit different. So looks like for this purpose, I should use the old method for band?

• For Fermi level, I always take the fermi level from dos calculation to plot the usual PBE band, because the fermi calculated from denser grid in dos calculation is more accurate. But in mBJ band, or in general meta-GGA/hybrid band, is the fermi level computed from itself is accurate enough? If I still take fermi from dos, then the band would shift a lot and doesn't make sense.

[fabien_tran1]

Hi,

I have repeated the two calculations of your last zip file. I used EDIFF=1.0E-8 for a better SCF convergence, but I reduced the number of k-points and ENCUT to make the calculations faster. ALGO=Normal and Damped give exactly the same total energies and band structures (with KPOINTS_OPT and a recent VASP version). I did first the calculation with ALGO=Normal, and then used the WAVECAR to start the calculation with ALGO=Damped. I suggest that you repeat your two calculations (and maybe some others) with EDIFF=1.0E-8. In principle, you should also get agreement between Normal and Damped.

[yujia_teng]

Hi Fabien,
Thanks very much for your tests. I found the weird Normal band is because it reads WAVECAR from step 1 below. If I don't read it in step 2, then it gives matched result with Damped. Increasing EDIFF to 1.0E-08 can improve the band a little bit and I think should be used here, but is not the reason why they are so different. I still couldn't get exactly same band and energy, maybe I'm wrong in some steps. But they are quite close and the difference is in acceptable range.

I used the following two steps:
1. Run PBE scf calculation with KPOINTS_OPT to generate WAVECAR.
2. Add meta-GGA related tags and Read WAVECAR from step 1.

So for Normal, we should only add the tags but keep ISTART=0. Or in other words, ALGO=Normal can directly use step2, without a previous step.

[fabien_tran1]

Yes you are right, EDIFF is not really important. I could reproduce what you described: the results with mBJ+U-Normal depend on how the calculation is started: from PBE+U WAVECAR or with ISTART=0. What matters is the way the calculation is started. The space of electronic configurations of magnetic systems, in particular those with f-electrons, can have multiple local minima and the initial conditions of the calculation can play an important role for determining in which minimum the calculation finishes. This is sometimes a question of luck.

Not related: It is very common to combine PBE with U, but I have rarely seen mBJ combined with U. Is it the case in the literature that you have read for this study?

[yujia_teng]

Hi Fabien,
Thanks for the explanation of different behavior of ALGO=Normal, it's very clear.

And it's true that usually people don't use mBJ+U. Same for HSE+U or other hybrid/meta-GGA+U I think. The paper I'm referring to are mostly about non-magnetic materials. So I tested mBJ without U, on my magnetic Ho-compound system. I tried two ways:
1. no U starting step 1 above
2. step 1 still has U but no U in step 2. (still PBE+U in 1st step and 2nd step is mBJ)

However, the resulted band below is clearly wrong from any method above, because the flat f-electron band comes to the fermi level. So for my system here, according to the test, U should be included in mBJ band. But maybe not be included in other cases or in general.

(The proj band is obtained with old method. So the two bands looks slightly different. But it doesn't change the fact that f-electrons comes to E_F.)