VASP Forum

Sorry for another question here, but I'm not able to get convergence of my HSE calculation. I've tried different mixing parameters and some other parameters but doesn't work. The input and output file are all attached below. Here are a few points:

• I'm, using 5.4.4, because version 6 is not available on the HPC machine I'm using, where it has more number of cores. HSE is too slow so I'm running it there. I do PBE with same KPOINTS first, and then HSE reads the wafefunction.

• My system contains f-electron. Physically speaking, HSE+U is not quite physical because they have same effect. But in meta-GGA case, I tested that if I don't include U, f-electron will come to fermi level, which is wrong. I didn't test it with HSE here, but I guess it would be same.

• I'm studying 3 cases, spin along 3 directions. Spin directions are defined with SAXIS, and MAGMOM of magnetic atom is in the form of 0 0 m, where m should be the starting magnitude of spin.

• 001 direction convergence is fine. However, for the other 2, convergence is not achieved after thousands of steps. When it hits NELM, I continue the calculation by reading WAVECAR. For 110 direction, the energy oscillates around 0.5E-05, and -0.15E-03 for 111

• Following the page Troubleshooting electronic convergence, I tried some mixing parameters (AMIX, AMIX_MAG, BMIX, BMIX_MAG) in the previous PBE step to obtain smaller number of electronic steps first, then run HSE. Also, LHFCALC this pages says we could reduce TIME. So I tried several values, 0.6 0.5 0.4 0.3 for different calculations but doesn't work.

• In 110 case, because it oscillates around 0.5E-05, I tried to reduce EDIFF to 1E-05 at the end. But it still cannot converge.

• There's a symmetry warning, but that doesn't show up in vasp 6 case (also attached, just when calculation starts). I indeed created gamma-centered kpoints. Don't know why the warning is there.

• ALGO=All seems to be worse. I tested it in 110 case. So I think ALGO=Damped should still be used here.

 -----------------------------------------------------------------------------
|                                                                             |
|           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   !!!           |
|                                                                             |
|      Your generating k-point grid is not commensurate to the symmetry       |
|      of the lattice.  This can cause   slow convergence with respect        |
|      to k-points for HF type calculations                                   |
|      suggested SOLUTIONS:                                                   |
|       ) if not already the case, use automatic k-point generation           |
|       ) shift your grid to Gamma (G) (e.g. required for hex or fcc lattice) |
|                                                                             |
 -----------------------------------------------------------------------------

Dear Yuija,

Don't worry, this forum is for questions, so feel free to ask.

I strongly urge you to compile VASP 6.5.1 on the HPC facility with more cores yourself. Please follow the installation guide, and search the wiki and the respective section of the forum if you have issues. I would also urge you to contact your HPC center and look at their documentation to find infos about the modules and the best toolchain to use, etc.

The symmetry routines have been reworked since VASP 5, so I hope you understand that it does not make much sense to troubleshoot your problem on 8 year old software that is outdated, especially since you have access to the current release.

General points to mention are:

• Read the section on f-electrons in DFT on the wiki. If at all possible, use the Ho_3 potential to get rid of the wrong representation of the f-electrons. Then the Hubbard U is not needed.

• Use a KPOINTS_OPT file!

• You can try to start from a WAVECAR file that has been preconverged with PBE.

• You should probably set GGA_COMPAT=F.

• Note that SAXIS is in cartesian coordinates. Your cell is not cubic, so make sure that the direction you want to set is correct.

Please feel free to post in the installation part of the forum if you have problems there that you cannot find elsewhere. And if you still have issues with this topic in VASP 6, reply here again. I will keep the topic open until you report positive or negative results.

Cheers, Michael

Dear Michael,
Currently I'm not able to install 6.5.1 on the HPC cluster with more cores due to some reasons. There's 6.1.2 there, and I'm not sure if it will make any difference. Many new changes are in later versions, like KPOINTS_OPT. So I still ran 6.5.1 on the usual cluster for now. And I understand that it doesn't make sense to check the symmetry problem in 5.4.4.

I just got one result with only 60 electronic steps and it's far from convergence, same as the case in 5.4.4. The input/output are attached. All the HSE calculation I did starts from preconverged PBE WAVECAR. Except kpoints, all other setting are same as 5.4.4 that I attached before. I don't know how to get the convergence here. I've tried those usual methods that I said in my original post but that doesn't work.

I tried Ho_3. It indeed converges soon, but the result is not what I want. I tried GGA_COMPAT=F in vasp 5.4.4 but there's no change.

Besides the convergence problem, I noticed that there may be a bug, the output is also attached. I see that this bug was fixed in line 2 in 6.5.0. Maybe that is just for usual calculation?

-----------------------------------------------------------------------------
|                                                                             |
|     EEEEEEE  RRRRRR   RRRRRR   OOOOOOO  RRRRRR      ###     ###     ###     |
|     E        R     R  R     R  O     O  R     R     ###     ###     ###     |
|     E        R     R  R     R  O     O  R     R     ###     ###     ###     |
|     EEEEE    RRRRRR   RRRRRR   O     O  RRRRRR       #       #       #      |
|     E        R   R    R   R    O     O  R   R                               |
|     E        R    R   R    R   O     O  R    R      ###     ###     ###     |
|     EEEEEEE  R     R  R     R  OOOOOOO  R     R     ###     ###     ###     |
|                                                                             |
|     A calculation of the band-structure with hybrid-functionals on the      |
|     k-points read from KPOINTS_OPT file was requested with NCORE/=1.        |
|     This is currently not supported.                                        |
|                                                                             |
|       ---->  I REFUSE TO CONTINUE WITH THIS SICK JOB ... BYE!!! <----       |
|                                                                             |
 -----------------------------------------------------------------------------

Best,
Yujia

Dear Yujia Teng,

Good that you have tried 6.5.1.

I have played a bit with your calculation, although I have reduced the KPOINTS mesh from 6x6x6 to 3x3x3 to make the tests more tractable and disabled the PBE+U part. I think it is valuable to see if the hybrid functional is enough to depin the Ho f-states from the Fermi level. I also switched to ALGO=All with ISEARCH=1, since this is usually very stable and converges well. I also focused on the main convergence, not the bandstructure calculation, so I did not provide a KPOINTS_OPT file.

What I noticed as well is that in your HSE calculation, the magnetic moment disappears. This is probably a bad sign.

I just got one result with only 60 electronic steps and it's far from convergence, same as the case in 5.4.4

I have looked at your OSZICAR, and it does not look too bad. The total energy is decreasing constantly. I think if you just allow the calculation to do more steps, you could converge it.

I tried Ho_3. It indeed converges soon, but the result is not what I want. I tried GGA_COMPAT=F in vasp 5.4.4 but there's no change.

Indeed, if you want to observe the magnetism in the 4f electrons, this potential will not work.

Indeed, I got your system to converge after 157 steps using the setup mentioned above. As you can see from my attached OUTCAR file, the magnetization is not stable along the axis you prescribed. It moves from (0 0 1) to (1 1 0) gradually.

I am not sure why it is shown as zero at the end of the OUTCAR, I will discuss this tomorrow with a collegue who is more knowledgable concerning non-collinear magnetism.

Ho, and all f-electron systems are very tricky, which you complicate by doing a non-collinear calculation and HSE. It is quite common that such setups do not converge in the default 60 NELM steps, so please try again with NELM=200. And make sure your magnetic setup is correct.

Cheers, Michael

Dear Michael,
Thanks very much for your test. Following your test on no U & HSE scf only, I did the usual HSE band calculation with also smaller parameters to reduce time (6x6x6 to 3x3x3 and ENCUT from 600 to 300). It converged at 111 steps and the resulted band looks reasonable (in/output and result attached). So looks like hybrid functional is able to push Ho-f away, similar effect as Hubbard U, and HSE+U shouldn't be used. I will do more tests to make sure everything is correct, which would take a while.

My magnetization setting of the attached file here is along 1 1 0 Cartesian direction (SAXIS = 1 1 0), not along 0 0 1. Maybe I errorly described it at first. In OSZICAR it is indeed still along 1 1 0 direction. However, I also observe that in OUTCAR it's all 0. By looking at 'grep magnetization OUTCAR' , everything is fine. But when I look at the block below, I found it suddenly drops to 0 right after 'aborting loop because EDIFF is reached', both in your OUTCAR and mine here. I also checked the spin along 001 case, but there's no such problem there (file hse 001 attached).

 magnetization (x)

# of ion       s       p       d       f       tot
--------------------------------------------------
    1        0.000   0.000   0.000   0.000   0.000
    2        0.000   0.000   0.000   0.000   0.000
    3        0.000   0.000   0.000   0.000   0.000
--------------------------------------------------
tot          0.000   0.000   0.000   0.000   0.000

And just want to confirm about magnetization output. In OSZICAR (component meaning should be same in OUTCAR) last line below, does the 3 components of mag mean Cartesian coordinate or components relative to SAXIS? PBE and HSE gives different result.

PBE: 1 F= -.23633979E+02 E0= -.23618885E+02  d E =-.301877E-01  mag=     0.0000     0.0000     3.6828

HSE: 1 F= -.35803301E+02 E0= -.35802468E+02  d E =-.166490E-02  mag=    -2.8448    -2.8450    -0.0008

Best,
Yujia

Just attaching one more result here. I did same thing in 6.1.2 (should almost same as 5.4.4), and magnetization is same in OUTCAR and OSZICAR there. The magnetization is still along z component, while in 6.5.1 it's along x and y. Maybe there's some change in new version that cause mag be 0 in OUTCAR and mag component definition different?

And if it doesn't converge within NELM, for hybrid functional or meta-GGA, is using ISTART=1 to read WAVECAR to continue the job reasonable? I found that convergence becomes harder again when I go back to 6*6*6 kpoints.

Dear Yuija Teng,

This problem is more complex than I thought, and I can now almost certainly say that things are buggy. Or at least they appear so at the moment.

Restarting from an existing WAVECAR is fine, but you should also read the magnetization density from CHGCAR.
So copy both files over and set ISTART=1 and ICHARG=1

For now, I advise you not to use SAXIS= 1 1 0, but keep the default SAXIS=0 0 1, and set MAGMOM= 5 5 0 0 0 0 0 0 0 instead of 0 0 10 0 0 0 0 0 0, both for the PBE+U and the HSE calculation. (Note that this should be physically, but not numerically equivalent, especially for a lower value of ENCUT, due to different FFT grids.)

For me, this results in much better convergence (30 steps for the HSE calculation at the lower k-density and lower ENCUT), and the magnetic moments are correctly reported in the OUTCAR file as well.

I will investigate further what goes wrong here, but I suspect that SAXIS might not work correctly for non-local functionals. Nothing definitive yet, however.

I hope I can get back to you soon with a better explanation, but this solution should give you quicker convergence for now.
Cheers, Michael

Dear Michael,
I did a test with your recommended setting and it also works well. For k-666, the first 2 iterations gives magnetization along 110 and this is correct (OUTCAR OSZICAR consistent). So at this point, this setting should give the correct result, and currently I just need to wait for the result to verify it, which may take 1-2 weeks.

I also doubted the problem may related to SAXIS, because spin along 001 case works fine in both new and old versions. So I tested the changed MAGMOM setting in 6.1.2 few days ago (it's faster there), but it doesn't work there. Even constraining magnetic moment with I_CONSTRAINED_M=1 doesn't work. One strange thing is that if I still use previous SAXIS = 1 1 0 in 6.1.2, k-333 works, but for k-666 the same problem shows up. Seems that there're more problems in previous versions. Just tested it because it's on HPC with more cores, so that I can get results to check in a few hours rather than days.

Best,
Yujia

Dear Yujia,

We double checked, and there is indeed a bug in the code for ISYM=3 (set by LHFCALC=T), and SAXIS≠0 0 1. The symmetry routines cannot handle different magnetization axes.
I added a warning to the SAXIS page on the wiki to reflect that, and also made a note in the known issues list, referencing this thread.

In the next version of VASP, the code will exit with an error if SAXIS is not the default, LHFCALC=T, and ISYM>=0.

Note that, as far as I know, this behaviour was not different in older VASP versions. SAXIS≠0 0 1 and LHFCALC=T never worked correctly together, I think.

Turning off symmetry operations with ISYM=-1 will result in correct behaviour, but will increase computation time by a lot, potentially.

Thanks for bringing this problem to our attention! Please let me know if you have any more questions on the topic, or if I can lock the thread.
Cheers, Michael

Dear Michael,
Thanks for confirming the bug. My problem is solved here with the correct method, and I also want to further clarify the bug here. I did more tests, in old version 6.1.2 (should have same result as 5.4.4) and 6.5.1 with old method. The results is attached below.

From my tests, except the bug with SAXIS, I think there's also a problem with the old method, i.e. k-points from IBZKPT + 0-weighted kpoints along band path. In previous version where only old method is avaliable, even setting MAGMOM = 1 1 0 and using default SAXIS doesn't give correct result (different magnetization direction in OUTCAR and OSZICAR). Constraining magnetization also doesn't work. In newest version 6.5.1 with old method, the result is also wrong with above setting.

So I think we can make a conclusion here about SAXIS ≠ 001:

• With symmetry on, we must use default SAXIS, setting MAGMOM, and using new KPOINTS_OPT method for hybrid functional band structure calculation. The old method would give wrong result.

• In previous versions where only old method is available, turning off symmetry with ISYM=-1 is the only way to obtain correct result.

Besides the SAXIS problem, I mentioned another problem before. I copy it below. Is this different from the BUGFIX here?

Besides the convergence problem, I noticed that there may be a bug, the output is also attached. I see that this bug was fixed in line 2 in 6.5.0. Maybe that is just for usual calculation?

Code: Select all

-----------------------------------------------------------------------------
| |
| EEEEEEE RRRRRR RRRRRR OOOOOOO RRRRRR ### ### ### |
| E R R R R O O R R ### ### ### |
| E R R R R O O R R ### ### ### |
| EEEEE RRRRRR RRRRRR O O RRRRRR # # # |
| E R R R R O O R R |
| E R R R R O O R R ### ### ### |
| EEEEEEE R R R R OOOOOOO R R ### ### ### |
| |
| A calculation of the band-structure with hybrid-functionals on the |
| k-points read from KPOINTS_OPT file was requested with NCORE/=1. |
| This is currently not supported. |
| |
| ----> I REFUSE TO CONTINUE WITH THIS SICK JOB ... BYE!!! <---- |
| |
-----------------------------------------------------------------------------

Best,
Yujia

I did some more tests with ISYM=-1 and looks like it's more problematic than we thought. Although some of them are not converged, I don't think reaching convergence would change the result here.

I found that in previous version with old method, even turning off symmetry doesn't give correct result. The magnetization from OSZICAR is always tilted a little bit.

For newest version 6.5.1, using SAXIS = 1 1 0 with symmetry off seems still isn't correct. There's about 0.01 \(\mu_B\) in other 2 directions.

For 6.5.1 with rotated MAGMOM, it seems correct, but there's still little difference between magnetization x and y component in OSZICAR. The difference is less with less k points (k-333) but increased with more k points (k-666), while it's exactly same in new method.

Dear Yujia,

I apologize for the late reply. I somehow did not get a notification or an email about your last posts, when, according to my forum settings, I should have gotten both.

I will look into this ASAP and get back to you as fast as I can. Thanks for further digging into this, and sorry for missing this for a while.

Cheers, Michael

Dear Yujia,

OK, I have some trouble bringing your comments in line with the files you uploaded, since you did so many different tests. So just to recap and to see if I got all the remaining issues sorted:

• You are seeing different magnetization reported in the OSZICAR and OUTCAR files.

• You do get correct results only when using KPOINTS_OPT, not when using the "old" method for HF bandstructure calculations.

• Even when turning off symmetry via ISYM=-1, you see variations between the x and y components in the magnetization when using the old method.

• You get an error message when using NCORE≠1

Different magnetizations reported in OSZICAR and OUTCAR files:
This is normal. The magnetization in OSZICAR is calculated by subtracting the total charge densities (up-down; more complicated in the NCL case), so the whole cell contributes. OUTCAR reports the local moments that are calculated based on a projection scheme set by the LORBIT tag. See the wiki documentation for more details. So in OUTCAR, you are missing the interstitial regions; thus, a difference in the values is to be expected.

You do get correct results only when using KPOINTS_OPT, not when using the "old" method for HF bandstructure calculations:
From looking at your files, I think your KPOINTS files are wrong. I am not talking about the added points with zero weight, which I did not check, but the weighted ones seem to have to much symmetry.
e.g. in the magmom-110_k-333_0.001_mag_difference folder (which should correspond to a 3x3x3 mesh I would think, you have only 4 kpoints with weights:

0.080   3   3   3    4  0.040   43    2   20   23           # Parameters to Generate KPOINTS (Do NOT Edit This Line)
    47
Reciprocal lattice
    0.00000000000000    0.00000000000000    0.00000000000000     1
    0.33333333333333    0.00000000000000    0.00000000000000     8
    0.33333333333333    0.33333333333333    0.00000000000000     6
   -0.33333333333333    0.33333333333333    0.00000000000000    12
    0.50000000000000    0.50000000000000    0.50000000000000     0
    0.47368421052632    0.47368421052632    0.47368421052632     0
    :
    :
    :
    

With ISYM=-1, the 3x3x3 mesh translates to 27 kpoints, all with the same weight of 1:

Automatically generated mesh
      27
Reciprocal lattice
    0.00000000000000    0.00000000000000    0.00000000000000             1
    0.33333333333334    0.00000000000000    0.00000000000000             1
   -0.33333333333333    0.00000000000000    0.00000000000000             1
    0.00000000000000    0.33333333333334    0.00000000000000             1
    0.33333333333334    0.33333333333334    0.00000000000000             1
   -0.33333333333333    0.33333333333334    0.00000000000000             1
    0.00000000000000   -0.33333333333333    0.00000000000000             1
    0.33333333333334   -0.33333333333333    0.00000000000000             1
   -0.33333333333333   -0.33333333333333    0.00000000000000             1
    0.00000000000000    0.00000000000000    0.33333333333334             1
    0.33333333333334    0.00000000000000    0.33333333333334             1
   -0.33333333333333    0.00000000000000    0.33333333333334             1
    0.00000000000000    0.33333333333334    0.33333333333334             1
    0.33333333333334    0.33333333333334    0.33333333333334             1
   -0.33333333333333    0.33333333333334    0.33333333333334             1
    0.00000000000000   -0.33333333333333    0.33333333333334             1
    0.33333333333334   -0.33333333333333    0.33333333333334             1
   -0.33333333333333   -0.33333333333333    0.33333333333334             1
    0.00000000000000    0.00000000000000   -0.33333333333333             1
    0.33333333333334    0.00000000000000   -0.33333333333333             1
   -0.33333333333333    0.00000000000000   -0.33333333333333             1
    0.00000000000000    0.33333333333334   -0.33333333333333             1
    0.33333333333334    0.33333333333334   -0.33333333333333             1
   -0.33333333333333    0.33333333333334   -0.33333333333333             1
    0.00000000000000   -0.33333333333333   -0.33333333333333             1
    0.33333333333334   -0.33333333333333   -0.33333333333333             1
   -0.33333333333333   -0.33333333333333   -0.33333333333333             1

If you don't sample your Brillouin zone evenly, the results will be wrong. Note that you can use the --dry-run feature to generate an IBZKPT file quickly for your current INCAR settings.

Even when turning off symmetry via ISYM=-1, you see variations between the x and y components in the magnetisation when using the old method:
So I think this is a combination of using the wrong KPOINTS data for uniform sampling and too high expectations. If you turn off symmetry operations, the magnetic moments can easily differ by a couple of per cent in x and y, even if they should be equivalent by symmetry. The magnetic subsystem has little effect on the energetics, and fluctuations are normal. Everything below 0.1muB can be pretty safely ignored in most cases.

You get an error message when using NCORE≠1:
This is expected behaviour. The bug fix was to add the error message. In the past, the results were just wrong because KPOINTS_OPT was never properly supported for hybrid functionals with NCORE≠1. Note that you have to use HFRCUT=-1 for HF calculations with KPOINTS_OPT, as written in the documentation.

Please let me know if these explanations have covered all your issues, or if there is still something unclear. If there is, please point me to the exact input and output files corresponding to the issue.

Cheers, Michael

Dear Michael,
Thanks very much for the detailed reply. Now I see what I did wrong. Setting ISYM=-1 means sampling all the Brillouin zone, not just adding this tag to INCAR. I usually make this mistake because I use vaspkit to generate the KPOINTS file with old method. The default setting of vaspkit is just sampling the irreducible BZ. I need to change it's setting to make it correct. I checked the k-333 case and got same 27 k-points as you provided.

Also, if we sample all the BZ (like 27-kpts in k-333 case), seems that ISYM=-1 or default doesn't change the result, because looks like the main point of turning off symmetry is sampling all the BZ here. I tested it in mBJ case and confirmed this if there're enough k points, but not sure about HSE.

Different magnetizations reported in OSZICAR and OUTCAR files:
This is normal. The magnetization in OSZICAR is calculated by subtracting the total charge densities (up-down; more complicated in the NCL case), so the whole cell contributes. OUTCAR reports the local moments that are calculated based on a projection scheme set by the LORBIT tag. See the wiki documentation for more details. So in OUTCAR, you are missing the interstitial regions; thus, a difference in the values is to be expected.

Yeah I'm aware that the difference is due this reason. What I mean here is that in one file it's 001 direction while the other is like 0.5 0 1 direction, not sth like 0 0 1 and 0 0 0.98.

You do get correct results only when using KPOINTS_OPT, not when using the "old" method for HF bandstructure calculations.

Even when turning off symmetry via ISYM=-1, you see variations between the x and y components in the magnetization when using the old method.

These two should due to the wrong KPOINTS file. I will recheck those.

You get an error message when using NCORE≠1:
This is expected behaviour. The bug fix was to add the error message. In the past, the results were just wrong because KPOINTS_OPT was never properly supported for hybrid functionals with NCORE≠1. Note that you have to use HFRCUT=-1 for HF calculations with KPOINTS_OPT, as written in the documentation.

Thanks for the info. Looks like I didn't use HFRCUT=-1 in my calculation, but there's no error. I'm not quite sure if the result is still correct. I would use it for future calculations.

I would do some further tests and reply back. There should be no problem this time.

Best,
Yujia

Dear Michael,
I did a few more test with correct KPOINTS file and there's no issue now.

• With new method, we can either rotating MAGMOM and with symmetry on or turning off symmetry with ISYM=-1 with SAXIS.

• With old method, we must turn off symmetry to obtain correct result (sampling all BZ).

Thanks again for the clear explanation. So the solution is indeed the one provided by you earlier or written on the 'known issue' page.

I just have some other thoughts below that may save some computation time in some cases, which may not strictly follow your solution. I'm not sure if they are correct.

For old method, we need the k points sampling all BZ here. I didn't use ISYM=-1, because I directly generate k points with vaspkit, which is same as using ISYM=-1 for an auto generated k mesh like 3*3*3 in KPOINTS and do a dry run to obtain IBZKPT.

• Here, the reason that we need to sample k points in all BZ I think is because magnetism lowers symmetry? So irreducible BZ changed and sampling irreducible BZ then map to whole BZ with symmetry of the crystal lattice is wrong (e.g. cubic lowered to monoclinic, but still use cubic symmetry). If the Bravais lattice is same with some magnetic direction, just sampling the irreducible BZ in old method should work? or it doesn't work due to some other reason so we always need to sample all BZ in old method?

• For ISYM=-1, in my several tests, I feel like the main role is let k points sample all the BZ. If I already did this in KPOINTS, then using this tag or not doesn't change result. I'm aware that it turns off symmetry operations, but in my tests, looks like it just affects k points. Because the result is same, using default ISYM could save a little bit time. But I'm not sure if this is always a 100% safe option.

Best,
Yujia