Calculating Dielectric Tensor of Bi2Te3 (HSE06 + SOC): Unphysical results and SCF oscillations

[michal_hartstein]

I am calculating the macroscopic static dielectric tensor of a Bi2Te3 rhombohedral primitive cell using the HSE06 hybrid functional alongside Spin-Orbit Coupling (SOC). The calculation employs the PEAD method (LCALCEPS = .TRUE.).

I encounter two main issues depending on the k-point mesh density:

1. Unphysical Results (2x2x2 Gamma-centered k-mesh): The calculation finishes successfully, but the resulting macroscopic static dielectric tensor is highly unphysical (diagonal values exceeding 1150).
2. SCF Convergence Oscillations (3x3x3 and denser k-meshes): I increased the k-mesh density. However, the SCF calculation now fails to converge, exhibiting severe energy oscillations during the conjugate gradient steps.

I have already tried switching from ALGO = N to ALGO = All and ALGO = Damped along with TIME = 0.2, but the severe oscillations persist in all cases.
I am looking for advice on how to stabilize the SCF convergence or recommended best practices for narrow-gap/topological insulators within this framework.

Technical details:

1. VASP version: 6.4.2
2. Pseudopotentials: Bi_d (15 electrons) and Te (6 electrons).

I have attached a zip file includes the INCAR, POSCAR, KPOINTS, POTCAR, OUTCAR, and OSZICAR files for the failing 3x3x3 calculation.

[alexey.tal]

Dear michal_hartstein,

I was able to reproduce the problem. Here is what you can try:

1. With sparse k-point sampling, you should use IPEAD = 1 because 3x3x3 mesh is not sufficiently dense to compute higher-order finite differences in PEAD.
   With it the calculation converges. For dense meshes you can keep the default.

2. The large dielectric constant is most likely caused by the too sparse k-point sampling. Note that in your POSCAR the c-axis points almost along x, so the large xx
   component is the out-of-plane value. You should try to converge the dielectric constant with respect to the k-mesh (e.g. 4x4x4, 6x6x6, 8x8x8 with PBE+SOC) before running HSE.

3. The HSE calculation is quite slow in this system, so to make the HSE calculations faster:

   • Set MAGMOM = 15*0 and remove SAXIS. Otherwise the symmetry is lowered and you get
     27 instead of 18 irreducible k-points in the field calculations.
   • EDIFF = 1.e-6 is probably sufficient.
   • PRECFOCK = Fast usually gives about a 2x speedup. Please check on one calculation that it does not change the gap noticeably.
   • Restart the LCALCEPS run from a converged zero-field WAVECAR on the same k-mesh.

Best regards,
Alexey