My Community

Dear VASP community,

I hope you could help me out with my problem with van der Waals calculations (relaxation of Bi2Se3) that has been bugging me for long. But first a quick introduction to Bi2Se3 structure: Bi2Se3 is composed of quintuple layers that are separated by van der Waals gaps. Here is an example structure (Se is yellow and Bi is purple):



Bi2Se3 belongs to a peculiar R3̅m space group, which has both rhombohedral (5 atoms in the unit cell) and hexagonal (15 atoms in the unit cell) representations. The hexagonal representation is more illustrative, and the above picture is a 3x3x1 times repeated hexagonal unit cell.

Computationally, it is expected that with PBE the vdW gap is overestimated, and when vdW functionals are applied, the gaps shrink close to the experimental value. This is also reported in practice in many papers. However, in my vdW simulations the energy is minimized by boundless increase of the vdW gaps (!), and also the a-paramater (see the above figure) undergoes a mystical increase from about 4.15Å to about 4.35Å so that it seems that the vdW interactions are strongly repulsive for me. This is for optPBE-vdW, optB88-vdW, optB86b-vdW and vdW-DF2 functionals. But for Tkatchenko-Scheffler the results have neglible or nearly neglible difference to PBE, so there are two extreme cases of things going wrong here.

In my relaxation method I have made a grid with different values for cell parameters a and c, and then relaxed the ions for each structure with ISIF=2 and finally made a 2D surface fit to find the correct a and c. This is because due to the vdW gaps the total energy dependence on the c-parameter was very weak for PBE (the PBE results were compatible with other literature for me). I thought that due to this ISIF=3 or ISIF=4 would be less accurate (the forces are very small even quite far from the energy minimum).

INCAR. Spin-orbit coupling (SOC) is on since Bi has strong SOC and it is what makes Bi2Se3 a topological insulator (the reason why I'm trying to study this material). In this INCAR I have used the optPBE-vdW functional. KPOINTS for rhombohedral calculations: POSCAR with R3̅m symmetry for rhombohedral calculations, with very large vdW gaps (a=4.36Å, c=35.9Å). The hexagonal version of this POSCAR is visualized in the picture above. After two ionic steps I obtained energy -17.750816eV. By the way, why does VASP like to break the R3̅m symmetry in the ionic optimization? My CONTCAR ends up having a lesser symmetry (with very small deviations from the symmetric situation).
POSCAR with R3̅m symmetry for rhombohedral calculations, and with cell parameters that according to this paper are the relaxed parameters. I obtained energy -17.108817eV (again after two ionic steps so that the latter one breaks the symmetry), which is more than 0.6eV higher than the surely incorrect structure above.
A picture of hexagonal 3x3x1 times repeated version of this POSCAR (so that you can see that the vdW gaps should in reality be rather small:



Things I have tried without success:

• More accurate calculations (particularly larger ENCUT)
• Both rhombohedral and hexagonal calculations
• Different POTCARS (Bi and Bi_d) for bismuth

Thank you in advance!

Okay, there was a problem with the figures, I used straight links to my Dropbox, and even though there should have been a permission for anyone to view them, they apparently don't work when embedded like this. For me it was possible to see them until I logged out from Dropbox. And apparently it's not possible to edit posts, so here are the figures, this time uploaded to internet:





Sorry about this.